Communication: Density functional theory embedding with the orthogonality constrained basis set expansion procedure

Embedding theories are of interest in computational chemistry because they alleviate computational expense through the partitioning of a chemical system into subsystems treated at different levels of theory.^1–14 Recent embedding methods have focused on ensuring the orthogonality of subsystem orbitals through a variety of projector methods,^8,11–13 avoiding the need to calculate the non-additive kinetic potential^1,2,15–19 that appears in DFT-in-DFT embedding. While some of these projector methods require a level shifting parameter,⁸ others are parameter free.^11–13 One recent chemical application of a parameter-free approach¹² uses a self-consistent field (SCF) version of the Huzinaga equation²⁰ to enforce subsystem orbital orthogonality. In this implementation, however, the density of the environment needed to be frozen during the SCF procedure to ensure numerical stability and reasonable convergence behavior.¹² Herein we present the orthogonality constrained basis set expansion^21,22 (OCBSE) procedure as a simple alternative to the Huzinaga equation for enforcing orbital orthogonality in the context of density embedding. Previously we used the OCBSE procedure in embedding theory for multicomponent DFT,²³ where electrons and specified protons are treated quantum mechanically on the same level, but herein this approach is applied to conventional electronic DFT. The main advantage of the OCBSE procedure is that excellent convergence behavior is attained for the systems studied without freezing any of the subsystem densities, while retaining a similar or higher level of accuracy as other approaches.

For the case of DFT-in-DFT embedding, the one-particle density matrix for the total system may be expressed in the atomic orbital (AO) basis as the sum of subsystem one-particle density matrices

where $γ AB$ is the total density matrix, $γ A$ is the subsystem A density matrix, and $γ B$ is the subsystem B density matrix. The one-particle density matrices $γ A$ and $γ B$ are constructed from the Kohn-Sham orbitals $ψ i A$ and $ψ i B$⁠, respectively, associated with subsystems A and B, which will be referred to as the “active” subsystem and the “environment,” respectively. The energy for both subsystems can be expressed as

where k = 1, 2 and refers to the level of theory used for each subsystem, $γ$ is a generic one-particle density matrix, and h is the core Hamiltonian. The two-electron terms are contained in G_k according to the relation

where J is the Coulomb functional and E_xc,k is the exchange-correlation functional at the level of theory indicated by the index k. Within the embedding theory framework, the total energy for the system is

For the purposes of this paper, the higher and lower levels of DFT will be denoted by subscripts 1 and 2, respectively, and only the active subsystem A is treated at the higher level of theory. If the $ψ i A$ and $ψ i B$ subsystem molecular orbitals (MOs) are constrained to be orthogonal, the Fock matrices for subsystems A and B are given by

The SCF procedure with the OCBSE protocol is performed as follows. Consider a system with $n tot$ MOs (i.e., n_tot atomic basis functions) partitioned into subsystems A and B with $n occ A$ occupied orbitals in subsystem A and $n occ B$ occupied orbitals in subsystem B. Alternating single SCF iterations are performed on subsystems A and B, and the densities are converged simultaneously. Iterations for subsystem A will be denoted i_A, and iterations for subsystem B will be denoted i_B. For a given iteration i_A or i_B, the Fock matrix F^A or F^B is constructed. To enforce orthogonality between subsystem A and B MOs, the OCBSE procedure transforms the matrix F^A or F^B to a new truncated basis via the transformation matrix W_A or W_B. The transformation matrix W_A has dimension n_tot × n_A where n_A = $n tot − n occ B$⁠, and the transformation matrix W_B has dimension n_tot × n_B where n_B = $n tot − n occ A$⁠. The first $n occ A$ columns of W_A are the occupied orbitals from subsystem A from iteration i_A − 1, and the remaining $n tot − n occ A − n occ B$ columns are the virtual orbitals of subsystem B from iteration i_B − 1. The matrix W_B is constructed analogously, where the first $n occ B$ columns of W_B are the occupied orbitals of subsystem B from iteration i_B − 1, and the remaining $n tot − n occ A − n occ B$ columns are the virtual orbitals of subsystem A from iteration i_A − 1. Iteration i_A or i_B requires performing the transformation

or

respectively. The SCF procedure is performed in the corresponding W basis for both the $F W A$ and $F W B$ matrices, and the new MOs obtained during iterations i_A and i_B in the W basis are back transformed to the AO basis. The $γ A$ and $γ B$ matrices are then constructed and used to compute the energy according to Eq. (4). New F^A and F^B matrices for iterations i_A + 1 and i_B + 1, respectively, are constructed according to Eq. (5), and new transformation matrices W_A and W_B are constructed as outlined above. This process is continued until both subsystems achieve SCF convergence. Note that the matrix transformations in Eqs. (6) and (7) and the subsequent orbital back transformations are of negligible computational expense compared to the calculation of the two-particle integrals or the contributions from DFT exchange-correlation potentials on a grid.

To demonstrate the effectiveness of the OCBSE method, the same test systems were studied as in Refs. 10 and 12, with geometries obtained from Ref. 10, except for the bond length of OH⁻, which is 0.978 Å as given in Ref. 12. These test systems are depicted in Fig. 1. An SCF calculation was performed on the entire system at the low level of theory, and the occupied MOs were localized using the Pipek-Mezey²⁴ method of orbital localization. Occupied MOs were assigned to the active subsystem by summing the Mulliken populations of atoms in subsystem A and including all orbitals for which this sum is greater than a threshold. To enable direct comparison, the threshold for orbital inclusion was chosen to ensure that the number of orbitals included in the active subsystem was the same as that reported in Ref. 12. The remaining occupied orbitals were assigned to subsystem B. These procedures for orbital localization and selection were chosen for simplicity and to allow comparison to previous results. However, the OCBSE method is general and can be used in conjunction with a variety of orbital localization and selection procedures.

The numbers of electrons in the active subsystem and the environment are fixed according to the number of occupied MOs assigned to each subsystem. The subsystem density matrices $γ A$ and $γ B$ were constructed from the orbitals $ψ i A$ and $ψ i B$⁠, respectively, and the Fock matrices for each subsystem were constructed according to Eq. (5). The subsystem densities were converged simultaneously employing the transformation to the W bases as outlined above. For all calculations, the 6-31G* basis set²⁵ was used. The low level of theory was LDA (local density approximation),^26,27 while the high level of theory was either B3LYP^28–30 or PBE (Perdew-Burke-Ernzerhof)³¹ denoted “B3LYP-in-LDA” or “PBE-in-LDA,” respectively. All calculations were performed with an in-house code for the OCBSE protocol, obtaining the electronic exchange-correlation contributions to the energies and Fock matrices via an interface to the GAMESS program.³² 

As in Refs. 10 and 12, the reference reaction energies obtained with the entire system treated at the high level of theory were subtracted from the reaction energies obtained with the embedding method as a function of the number of carbon atoms in the active subsystem to compute the reaction energy errors. If no carbon atoms are in the active subsystem, the reaction energy error is the difference between the reaction energies computed at the high and low levels of theory.^10,12 The Huzinaga data were obtained from Ref. 12. Note that the reference energies in this work are slightly different from those of Ref. 12 due to the use of Cartesian rather than spherical basis functions.

The first test case is a simple substitution reaction in which chlorodecane reacts with a hydroxide anion to form decanol and a chloride anion [Fig. 1(a)]. As shown in Fig. 2, both the Huzinaga and OCBSE methods perform well, with the OCBSE method performing better when only one carbon atom is included in the active subsystem.

For the second case, decanoic acid is deprotonated to form decanoate, and therefore the charge of the reactant is different from the charge of the product [Fig. 1(b)]. This case is more challenging than the previous reaction because of the effects of charge redistribution. These challenges are reflected in the larger errors as a function of the number of carbon atoms in the active subsystem, as shown in Fig. 3. For B3LYP-in-LDA embedding, the OCBSE method is clearly more accurate than the Huzinaga method, with errors significantly below 1 kcal/mol for every size of the active subsystem with the exception of the two carbon atom point. As in previous calculations on this system,^10,12 the error does not decrease monotonically as the active subsystem is enlarged, suggesting that more than two carbon atoms should be included in the active subsystem. For PBE-in-LDA embedding, the Huzinaga method gives smaller errors for one and two carbon active subsystems, but the OCBSE method becomes more accurate from three carbon atoms in the active subsystem onward.

The final test system is the Diels-Alder reaction of 1,3-butadiene and octadecanonaene [Fig. 1(c)]. This reaction has a much higher degree of electron delocalization and conjugation than the previous two systems. For this system, 1,3-butadiene is always part of the active subsystem, and the x-axis of Fig. 4 is the number of carbon atoms of octadecanonaene that are included in the active subsystem, starting with the central C=C bond and enlarging symmetrically. The reaction energy error for the OCBSE method with both the B3LYP-in-LDA and PBE-in-LDA schemes falls below 1 kcal/mol when four or more carbon atoms are included in the active subsystem, clearly demonstrating the effectiveness of the OCBSE method for this system. The reaction energy errors for the OCBSE method are comparable to those for the Huzinaga method.

The goal of this work is to communicate the effectiveness of the OCBSE procedure as an alternative to the Huzinaga equation in density embedding methods. Both of these methods retain orbital orthogonality between the subsystems to avoid the calculation of the non-additive kinetic potential, and both methods circumvent the necessity of a parameter, as required in level shifting, with negligible computational cost. In the OCBSE procedure, the subsystem densities do not need to be frozen during the SCF cycles and are converged simultaneously. Such an approach is less constrained and therefore is likely to be more generally applicable, although it typically requires more SCF iterations. In terms of accuracy, the OCBSE method performs as well as or better than the Huzinaga method for the three systems studied. Thus, the OCBSE procedure has been shown to be a viable alternative to the Huzinaga method for DFT-in-DFT embedding. Moreover, the ability to obtain stable convergence behavior when both subsystem densities are allowed to respond to each other could prove advantageous for the development of other embedding methods in the future, including extensions to embedding approaches combining wavefunction and DFT methods.

This material is based upon work supported by the National Science Foundation Grant No. CHE-1361293 and by the Center for Chemical Innovation of the National Science Foundation (Solar Fuels, Grant No. CHE-1305124). We would like to thank Professor Mihály Kállay for helpful correspondence.