Density functional theory embedding for correlated wavefunctions: Improved methods for open-shell systems and transition metal complexes

# Density functional theory embedding for correlated wavefunctions: Improved methods for open-shell systems and transition metal complexes

## Abstract

Density functional theory (DFT) embedding provides a formally exact framework for interfacing correlated wave-function theory (WFT) methods with lower-level descriptions of electronic structure. Here, we report techniques to improve the accuracy and stability of WFT-in-DFT embedding calculations. In particular, we develop spin-dependent embedding potentials in both restricted and unrestricted orbital formulations to enable WFT-in-DFT embedding for open-shell systems, and we develop an orbital-occupation-freezing technique to improve the convergence of optimized effective potential (OEP) calculations that arise in the evaluation of the embedding potential. The new techniques are demonstrated in applications to the van-der-Waals-bound ethylene-propylene dimer and to the hexaaquairon(II) transition-metal cation. Calculation of the dissociation curve for the ethylene-propylene dimer reveals that WFT-in-DFT embedding reproduces full CCSD(T) energies to within 0.1 kcal/mol at all distances, eliminating errors in the dispersion interactions due to conventional exchange-correlation (XC) functionals while simultaneously avoiding errors due to subsystem partitioning across covalent bonds. Application of WFT-in-DFT embedding to the calculation of the low-spin/high-spin splitting energy in the hexaaquairon(II) cation reveals that the majority of the dependence on the DFT XC functional can be eliminated by treating only the single transition-metal atom at the WFT level; furthermore, these calculations demonstrate the substantial effects of open-shell contributions to the embedding potential, and they suggest that restricted open-shell WFT-in-DFT embedding provides better accuracy than unrestricted open-shell WFT-in-DFT embedding due to the removal of spin contamination.

## I Introduction

The demand for accurate and efficient descriptions of complex molecular systems requires development of quantum embedding methods for electronic structure in which a small subsystem is treated with a high level of theory while the remainder of the system is treated at a more affordable level. Widely used examples of quantum embedding include QM/MM,(1); (2); (3); (4); (5); (6) ONIOM,(7); (8) and fragment molecular orbital (FMO) approaches,(9); (10); (11) which have led to significant advances in the simulation of condensed-phase and biomolecular systems. However, such methods generally rely on empirical models for the subsystem interactions, including link-atom approximations for embedding across covalent bonds(13); (12); (14); (15) and point-charge electrostatic descriptions of the environment,(4); (9) that are difficult to systematically improve and that can fail in practical applications.(3); (6); (16); (17)

Density functional theory (DFT) offers an appealing framework for addressing this challenge.(18); (19); (20); (21); (22); (23); (24); (25); (26); (27); (28); (29); (32); (33); (30); (31); (36); (37); (38); (39); (41); (34); (35); (40); (42) DFT embedding provides a formulation of electronic structure theory in which subsystem interactions depend only on their electronic densities, including non-additive contributions due to the electrostatic, exchange-correlation (XC), and kinetic energy terms. In the WFT-in-DFT embedding approach, the DFT embedding potential is included as an external potential for WFT calculations, providing a WFT-level description for one (or more) subsystem while the remaining subsystems and their interactions are seamlessly treated at the DFT level of theory.

Several groups, including this one, have recently demonstrated that non-additive kinetic energy contributions to the embedding potential can be exactly computed(27); (28); (33); (26); (35) with the use of optimized effective potential (OEP) methods.(43); (44); (45); (46); (47); (48); (49) In this paper, we introduce a simple technique to improve the robustness of OEP calculations in systems that exhibit small HOMO-LUMO gaps, such as transition metal complexes. In addition, we derive spin-dependent embedding potentials to enable the accurate description of open-shell systems in the WFT-in-DFT embedding framework. Numerical applications to the van-der-Waals-bound ethylene-propylene dimer and to the hexaaquairon(II) transition-metal cation illustrate the applicability of these new techniques and demonstrate the accuracy of the WFT-in-DFT approach in systems for which conventional density functional theory methods exhibit substantial errors.

## Ii Theory

Like Kohn-Sham (KS)-DFT, DFT embedding provides a formally exact framework for the ground-state electronic structure problem. Here, we review DFT-in-DFT embedding and its basis for WFT-in-DFT calculations.

### ii.1 DFT-in-DFT embedding

We begin by considering a closed-shell system in which the total electronic density consists of two subsystems, . The corresponding one-electron orbitals for and obey the Kohn-Sham Equations with Constrained Electron Density (KSCED),(22)

 [−12∇2+veff[ρA,ρAB;r]]ϕAi(r)=ϵAiϕAi(r), (1) [−12∇2+veff[ρB,ρAB;r]]ϕBj(r)=ϵBjϕBj(r), (2)

where , , and and are the number of electrons in the respective subsystems. is the effective potential for the coupled one-electron equations,

 veff[ρA,ρAB;r] = vKSeff[ρA;r]+vemb(r), (3)

where is the standard KS potential for subsystem A, and

Here, is the nuclear-electron Coulomb potential from the nuclei contained in subsystem B, is the Hartree potential, and is the XC potential. In addition to these familiar terms from conventional KS-DFT calculations, DFT embedding introduces the non-additive kinetic potential (NAKP) which properly enforces Pauli exclusion between the subsystem densities. It is obtained from

where , and is the non-interacting kinetic energy functional. The total energy functional for the full system is then

 E[ρAB] = Ts[ρA]+Ts[ρB% ]+Tnads[ρA,ρB]+ (6) Vne[ρAB]+J[ρAB]+Exc[ρAB],

where the last three terms on the right-hand side (RHS) are the nuclear-electron Coulomb energy, the Hartree energy, and the XC energy computed over the total density. Enforcing to be identical for all subsystems (see Sec. III B) leads to a unique partitioning in the DFT embedding formulation, such that the specification of the nuclei and the integer number of electrons in subsystem A and B fully specify the density partitioning.(26); (35) Eqs. 1-6 are easily generalized to the description of multiple embedded subsystems.

We have previously demonstrated that by using OEP methods to calculate the NAKP, DFT-in-DFT embedding can accurately describe both weakly and strongly interacting subsystems, including subsystems connected by covalent bonds;(27); (28) and we have shown that this method is computationally feasible for large systems by way of localized approximations to the NAKP.(28) More recently, we have introduced a projection approach that completely avoids the NAKP calculation in exact DFT embedding,(29) which appears worthy of further investigation. The OEP-based approach employed here is appealing because it provides a local embedding potential that is a functional of only the subsystem electronic densities.

In practice, the KSCED equations (Eq. 1 and Eq. 2) are solved by simply modifying the core Hamiltonian in the self-consistent field (SCF) calculation to include the additional embedding terms. The embedding potential (Eq. 4) can be written in the atomic orbital (AO) basis as

 vemb = vBne+J[γ% B]+vxc[γAB]− (7) vxc[γA]+v%nad[γA,γAB],

where the various terms on the RHS of this expression correspond to those in Eq. 4. The subsystem and total AO density matrices in Eq. 7 satisfy . It follows that the Fock matrix for subsystem A can be expressed as

 fA in B = hA in B+J[γA]+vxc[γA], (8)

where

 hA in B=hA+vemb, (9)

and is the core Hamiltonian for subsystem A (the kinetic energy plus external potential due to the nuclei in subsystem A). The Fock matrix for subsystem B, , is analogously defined.

### ii.2 WFT-in-DFT embedding

The embedding potential in Eq. 4 describes the subsystem interactions in terms of their corresponding electronic densities. However, the subsystem densities can be computed with any level of theory, thus allowing for the description of one subsystem at the (single- or multi-reference) WFT level, while the remaining environment is treated at the DFT level.(23); (29); (34); (35); (50); (51); (52); (53); (54) Closed-shell WFT-in-DFT embedding simply involves performing a WFT calculation on a given subsystem using the modified core Hamiltonian, in Eq. 9, that contains the embedding terms due to the environment of the other subsystem. The WFT-in-DFT energy is then obtained by modifying the DFT energy with respect to subsystem contributions at the WFT level,(35)

 Etot[ρWFTA,ρDFTB] = EDFTAB[ρDFTAB] (10) −(EDFTA[ρDFT%A]+∫vemb(r)ρDFTA(r)dr) +(EWFTA[ρWFT%A]+∫vemb(r)ρWFTA(r)dr).

This expression is easy to evaluate since the terms in the parentheses are just the DFT and WFT energies of subsystem A performed using the modified core Hamiltonian, . Just as DFT-in-DFT embedding is an exact theory for the case of an exact DFT XC functional, WFT-in-DFT embedding is an exact theory for the case of an exact DFT XC functional and a full configuration interaction (FCI) WFT description.(52)

## Iii Methods of Implementation

Here, we describe techniques to improve the accuracy and convergence of both DFT-in-DFT and WFT-in-DFT calculations. First, a description of open-shell DFT embedding is developed to incorporate the effects of spin-dependence in the embedding potential. Then, implementation of the OEP calculation is discussed, and an orbital occupation constraint is introduced to enable robust DFT-in-DFT and WFT-in-DFT embedding calculations for systems with low-lying virtual orbitals, such as transition metal complexes.

### iii.1 Embedding for open-shell systems

For an open-shell embedded subsystem, the and electrons generally experience different embedding potentials due to differing non-additive XC and NAKP contributions. Previous WFT-in-DFT implementations for open-shell systems have in practice neglected this difference, effectively assuming that the spin polarization is localized within the WFT subsystem.(34); (53) In this study, we show that effects due to spin-dependent potentials are substantial and easily included via separate and embedding potentials. We develop approaches to utilize both restricted and unrestricted open-shell orbital formulations in WFT calculations.

#### Open-shell DFT-in-DFT embedding

We begin by considering an open-shell system for which the total electronic density is comprised of the and density of the two subsystems, . The effective potential for the unrestricted spin orbitals(27) is

 vαeff[ραA,ρβA,ραB,ρβB;r] = vKS,αeff[ραA,ρβA;r]+vαemb(r), (11)

where is the standard KS effective potential for the unrestricted (U)KS orbitals, and is a spin-dependent embedding potential applied only to the -spin electrons,

The corresponding quantities for the -spin electrons are analogously defined. The total energy functional for the full open-shell system is then

Separate OEP calculations are performed over the and spin-densities for the exact calculation of the NAKP, which allows for numerically exact unrestricted open-shell DFT embedding (U-DFT-in-DFT).(27)

In practice, we solve for the unrestricted spin orbitals by adding the spin-dependent embedding potentials to the and Fock matrices. The and embedding potential can be written in the AO basis as

 vξemb = vBne+J[γ% B]+vξxc[γαAB,γβAB]− (14) vξxc[γαA,γβA]+vξnad[γξA% ,γξAB],

where , and the corresponding Fock matrices are

 fξ,A in B = hA+J[γA]+vξxc[γαA,γβA]+vξemb. (15)

The unrestricted spin orbitals for subsystem A are then obtained by separately diagonalizing and in the usual way.

Practical implementations for performing OEP calculations using restricted open-shell orbitals have yet to be developed. We thus introduce a simple, approximate scheme for restricted open-shell DFT embedding (RO-DFT-in-DFT). In this approach, a U-DFT-in-DFT calculation is first performed, and the embedding potentials and are constructed using Eq. 14. Then, and are constructed using Eq. 15, and the usual RO approach is employed to obtain subsystem orbitals that are spatially identical for the and electrons. Specifically, is diagonalized to obtain a set of occupied spin orbitals, {}, and is then projected into the space of the occupied spin orbitals using

 ~fβ,A in B=cTfβ,A in Bc, (16)

where is the matrix with columns comprised of the AO coefficients for {}. Finally, the projected Fock matrix, , is diagonalized to obtain the set of RO orbitals, {}, with the first orbitals doubly occupied and with orbitals singly occupied, where and indicate the number of and electrons in subsystem A. Although the second and third terms on the RHS of Eq. 15 are updated at each iteration of the RO-DFT-in-DFT calculation, we leave the embedding potentials unchanged to avoid performing OEP calculations using restricted open-shell orbitals. The RO-DFT-in-DFT energy for the total density is calculated using Eq. 13.

Several different schemes have been proposed to calculate the embedding potential for open-shell subsystems while neglecting its spin-dependence.(34); (35); (53) These approaches generally assume that interactions between the subsystems can be described by a single embedding potential. For example, in systems with an even number of electrons, the embedding potential, in Eq. 7, can be obtained assuming that each embedded subsystem is closed-shell, and then the open-shell subsystem is calculated using . In this approach, Eq. 15 is solved self-consistently while and are held fixed, and the final DFT-in-DFT energy is calculated using Eq. 13. This spin-independent description of the embedding potential can be used with either an unrestricted treatment of the open-shell subsystem (U-DFT-in-DFT-CS) or a restricted treatment of the open-shell subsystem (RO-DFT-in-DFT-CS); we later employ the approach to compare with the previously described methods (U-DFT-in-DFT and RO-DFT-in-DFT) that include spin-dependence in the embedding potential.

#### Open-shell WFT-in-DFT embedding

Unrestricted open-shell WFT-in-DFT (U-WFT-in-DFT) calculations are performed by first computing unrestricted Hartree-Fock (UHF) orbitals in the spin-dependent embedding potential, and then using these orbitals for a post-HF WFT calculation. The and Fock matrices for the calculation of the UHF orbitals are

 fξ,A in B = hA+J[γA]+Kξ[γξA]+vξemb, (17)

where , is the HF exchange matrix, and the embedding potentials (Eq. 14) are obtained from a U-DFT-in-DFT embedding calculation. The total energy is evaluated using

 Etot[ρWFTA,ρDFTB] = EDFTAB[ρDFTAB] (18) −⎛⎝EDFTA[ρDFT%A]+∑ξ∈{α,β}∫vξemb(r)ρξ,DFTA(r)dr⎞⎠ +⎛⎝EWFTA[ρWFT%A]+∑ξ∈{α,β}∫vξemb(r)ρξ,WFTA(r)dr⎞⎠.

Restricted open-shell WFT-in-DFT (RO-WFT-in-DFT) calculations are performed by solving for restricted open-shell HF (ROHF) orbitals in the spin-dependent embedding potential. Just as in RO-DFT-in-DFT embedding, a U-DFT-in-DFT calculation is first performed, and the embedding potentials and are constructed using Eq. 14. Then, the Fock matrices in Eq. 17 are constructed and the usual approach is employed to obtain RO orbitals; the second and third terms on the RHS of Eq. 17 are updated at each iteration while the embedding potential is left unchanged. The ROHF orbitals are used in the post-HF WFT calculation, and the total energy is then evaluated using Eq. 18. We note that for the RO-WFT-in-DFT energy calculation, the term , is evaluated using the ROKS-DFT energy.

### iii.2 Optimized effective potential

As seen in Eqs. 5 and 6, DFT embedding requires computation of both the kinetic energy, , and its functional derivative. However, since the explicit functional form for the kinetic energy is unknown, OEP methods are needed to obtain these terms exactly.

The OEP is the local potential for which solution of the one-electron equations

 [−12∇2+vOEP(r)]ϕi = ϵiϕi (19)

yields orbitals that correspond to a given target density while minimizing the non-interacting kinetic energy. A variety of methods for determining such potentials from an input target density have been developed.(43); (44); (45); (46); (47); (48); (49) Calculations reported here employ the direct optimization procedure developed by Wu and Yang,(43); (44) in which the kinetic energy is obtained via the unconstrained maximization

 Ts[ρin]=maxvOEP(r){Ws[Ψdet,vOEP(r)]}, (20)

where

 Ws[Ψdet,vOEP(r)] = 2N2∑i⟨ϕi|−12∇2|ϕi⟩ (21) +∫(ρOEP(r)−ρin(r))vOEP(r)dr − ζ||∇vλ(r)||2,

and

 vOEP(r) = vKSeff[ρin;r]+vλ(r). (22)

In these equations, , {} comprise an auxiliary basis set for the potential, are the corresponding expansion coefficients, and is a regularization parameter.(44) Maximization of utilizes the Newton method for optimization with a back-tracking line search in the expansion coefficients,(55)

 b(i+1)=b(i)+τH−1g, (23)

where is the iteration number, and are the Hessian and gradient of , respectively, and [0,1] is the step-size in the line search.

In practice, to obtain the embedding potential, we do not explicitly calculate the NAKP for each subsystem. Instead, for closed shell subsystems, we directly update the embedding potential (Eq. 4) at each iteration of the KSCED equations using(26)

 v(i+1)emb(r)=v(i)emb(r)−θvλ(r), (24)

where [0,1] is a damping coefficient. By construction, the embedding potential for each subsystem is identical at every iteration. The KSCED equations are initialized using , such that the initial guess for the NAKP for subsystem A exactly cancels the remaining terms in Eq. 4 (and the initial guess for the NAKP for subsystem B likewise cancels the corresponding terms). Upon convergence of the KSCED equations, , , and . Enforcing the embedding potential to be identical for all subsystems leads to a unique partitioning of the subsystem densities.(26); (35)

For open-shell calculations, we similarly update the spin-dependent embedding potential (Eq. 12); the OEP obtained for a given spin density is

 vξOEP(r)=vKS,ξeff[(ραA+ραB),(ρβA+ρβB);r]+vξλ(r), (25)

and as in Eq. 24, is used to update at each iteration.

Finally, we note that XC functionals that include a fraction of the exact exchange can be employed in DFT embedding via the OEP calculation. The HF exchange matrix, , is evaluated using , the OEP density matrix in the AO basis. For DFT-in-DFT embedding, the exchange energy is calculated using

 EX[γOEP,γin] = −tr(γOEPK[γ%OEP]) (26) +tr((γOEP−γin)K[γOEP]),

where the second term on the RHS corrects the exchange energy for small numerical differences between and . For calculations on the low-spin state of the hexaaquairon(II) cation, this correction is found to be as large as 20 kcal/mol; however, the correction is not required for the evaluation of the WFT-in-DFT energy (Eq. 10), since is obtained directly from a KS-DFT calculation on the full system.

### iii.3 Orbital-occupation freezing

For to be a concave function of , it is necessary(43) that the orbitals used to construct in Eq. 21 correspond to the lowest eigenvalues of Eq. 19. However, this can be problematic for systems with small energy differences between the occupied and virtual orbitals, where small changes in can alter the relative ordering of the orbitals.

To illustrate this issue, Fig. 1 shows the line search for an illustrative Newton step in an OEP calculation for the low-spin hexaaquairon(II) cation. is plotted as a function of , where is the step-size in Eq. 23. For any step-size larger than in this case, the orbital occupancy changes from one in which only t-like orbitals are occupied to one in which e-like orbitals are occupied. In traditional back-tracking line searches, any step which increases would be accepted, including the step indicated with the red arrow. However, this step is problematic since the Hessian and gradient of for the next Newton step would be evaluated using a density that corresponds to the wrong orbitals. The net results are poor convergence and incorrect solutions for the OEP.

We introduce a simple method to alleviate this problem by modifying the back-tracking line search. Reference (=0) orbitals are computed from Eq. 19 using , and for any proposed step-size , the corresponding orbitals are computed using . The proposed step is rejected if the overlap between these two sets of orbitals is less than 0.5, regardless of the change in ; otherwise, it is subjected to the usual criteria of the back-tracking line search algorithm. Upon rejection, the step-size is reduced by a factor of 2. This technique ensures that the correct orbitals remain occupied throughout the maximization of . In Fig. 1, the proposed step indicated by the red arrow is rejected, whereas the proposed shorter step indicated by the black arrow is accepted; not only is the value of increased, but the correct orbitals remain occupied. By utilizing this technique, we found that the maximization of typically requires less than 20 Newton steps for the low spin state of the hexaaquairon(II) cation, whereas the optimization failed to converge without the use of orbital-occupation freezing.

### iii.4 Computational Details

The DFT embedding methods employed here are all implemented in the development version of the Molpro software package.(56) All calculations employ the supermolecular basis set convention, in which the molecular orbitals for each subsystem are described in the AO basis for the full system.(57) Calculations on the ethylene-propylene dimer use the aug-cc-pVTZ orbital basis set for the carbon atoms and the aug-cc-pVDZ orbital basis set for the hydrogen atoms. Calculations on the hexaaquairon(II) cation use the aug-cc-pVTZ orbital basis set for the iron atom and the aug-cc-pVDZ orbital basis set for the hydrogen and oxygen atoms. For the auxiliary basis set used in the OEP calculations, we employ atom-centered Gaussian basis functions (, where is the normalization constant) for which the coefficient assumes values of , where . Calculations on the ethylene-propylene dimer employ the basis set for which the -type functions for the carbon and hydrogen atoms span {} = {-4, 4}, and the -type functions for the carbon and hydrogen atoms span {-2, 2}. Calculations for the hexaaquairon(II) cation employ the basis set for which the -type functions for the iron atom span {-4, 6}, the -type functions for the iron atom span {-4, 6}, the -type functions for the iron atom span {-2, 2}, the -type functions for the oxygen atoms span {-4, 6}, the -type functions for the oxygen atoms span {-2, 4}, the -type functions for the hydrogen atoms span {-4, 4}, and the -type functions for the hydrogen atoms span {-2, 2}. For all systems, the finite auxiliary basis set for the OEP calculations was confirmed to introduce a difference of less than 1 kcal/mol between the total energy computed using KS-DFT and either closed-shell or unrestricted open-shell DFT-in-DFT embedding. The regularization parameter used in the OEP calculations is set to ; smaller values were tested on the ethylene-propylene dimer and the hexaaquairon(II) cation and were found to have only a small ( effect on the total DFT-in-DFT energy.

The KSCED equations are initialized with subsystem densities comprised of the superposition of HF atomic densities and with ; different initial guesses for the embedding potential were tested on the hexaaquairon(II) cation and were found to yield similar final embedding potentials with only small () changes in the total DFT-in-DFT energy.

## Iv Results

### iv.1 The Ethylene-Propylene Dimer: WFT-in-DFT Embedding

The ethylene-propylene dimer is a prototypical system for which quantum embedding methods, such as QM/MM or ONIOM, may be employed. It exhibits a weak interaction that is difficult to address with conventional KS-DFT methods, while also exhibiting a spectator -CH moiety that contributes little to the interaction energy while substantially increasing the cost of the high-level calculation. However, unlike the QM/MM treatment of subsystems, the interactions between the system and the -CH moiety can be treated seamlessly using WFT-in-DFT embedding, as is now demonstrated.

Fig. 2(a) presents the ethylene-propylene dimer dissociation curve plotted as a function of the distance between the ethylene and propylene bonds, with the equilibrium dimer geometry obtained via minimization at the MP2 level of theory. Other geometries along the curve are obtained by displacing the two molecules along the vector formed between the midpoints of the two C=C bonds, while fixing all other internal coordinates. The relative energies are plotted by aligning each curve at infinite separation. The full CCSD(T) calculation (black) shows a binding energy of 2.0 kcal/mol. KS-DFT calculations using the PBE(58); (59) (green), B3LYP(60) (orange), B-LYP(61); (62) (blue) and B88-P86(61); (63) (cyan) XC functionals illustrate the difficulty in describing dispersion interactions using KS-DFT. The PBE functional underestimates the binding energy by 1.3 kcal/mol, while the rest of the XC functionals fail to capture any of the attractive interactions.

Finally, the red curve in Fig. 2(a) presents the results of WFT-in-DFT embedding, using a subsystem partitioning in which the 32 electrons associated with the system (CH-CH-, red in Fig. 2(b)) are treated at the WFT level of theory and the remaining 8 electrons in the -CH moiety are treated at the DFT level of theory. We employ CCSD(T) for the WFT and the B3LYP XC functional for the DFT (i.e., CCSD(T)-in-B3LYP). Fig. 2(a) shows excellent agreement between the CCSD(T) (black) and CCSD(T)-in-B3LYP (red) calculations; these curves, which are graphically indistinguishable, differ by less than 0.10 kcal/mol through the entire range of distances. We have confirmed that this level of accuracy is maintained with different XC functionals used for the DFT; specifically, CCSD(T)-in-(B-LYP) energies differ from the CCSD(T) results by less than 0.20 kcal/mol throughout the entire curve. These results illustrate that WFT-in-DFT embedding can be used to systematically improve DFT results and to avoid embedding errors while partitioning across covalent bonds.

### iv.2 The hexaaquairon(II) cation

We now present DFT-in-DFT and WFT-in-DFT calculations for the high-spin [] and low-spin [] states of the hexaaquairon(II) cation, a system that presents challenges due to the presence of low-lying unoccupied orbitals, the important role of unpaired electrons, and the relatively large number of electrons (84 e) in the full system. First, we test the accuracy of DFT-in-DFT embedding for the various treatments of the open-shell embedding potential described earlier. We then employ WFT-in-DFT calculations to investigate the low-spin/high-spin energy splitting and the ligation energy for this transition metal complex.

#### DFT-in-DFT embedding

Fig. 3(a) presents the potential energy curve for the simultaneous dissociation of all six HO ligands of the hexaaquairon(II) cation, plotted as a function of the average iron-oxygen distance. The equilibrium geometries for the low-spin [] and high-spin [] states are obtained using KS-DFT energy minimization with the B3LYP XC functional; all other geometries are obtained by uniformly stretching the iron-oxygen distances in the complex, keeping all other internal coordinates unchanged. All KS-DFT and DFT-in-DFT embedding results reported in this section are obtained using the B3LYP XC functional. The curves in the main panel of Fig. 3(a) are vertically shifted to share a common minimum value; they are not horizontally shifted. The high-spin state is lower in energy and exhibits a longer average iron-oxygen distance than the low-spin state.

We perform DFT-in-DFT embedding using a subsystem partitioning in which the 24 electrons associated with the iron center comprise one subsystem (red in Fig. 3(b)) and the remaining 60 electrons associated with the six water ligands comprise a second subsystem (blue in Fig. 3(b)). For the low-spin state, Fig. 3(a) demonstrates good numerical agreement between DFT-in-DFT (red) and KS-DFT (black); the relative energies differ by less than 0.6 kcal/mol throughout the range of reported internuclear distances.

For the high-spin state of the hexaaquairon(II) cation, Fig. 3(a) shows that the UKS-DFT and ROKS-DFT methods are in good agreement with each other, as well as with the corresponding U-DFT-in-DFT and RO-DFT-in-DFT embedding approaches described in Sec. III A 1. The U-DFT-in-DFT calculation accurately reproduces the relative energies obtained from UKS-DFT to within 0.4 kcal/mol throughout the attractive branch of the curve and to within 0.8 kcal/mol at shorter distances. The RO-DFT-in-DFT calculation reproduces the relative energy obtained from ROKS-DFT to within 1.0 kcal/mol throughout the attractive branch of the curve and to within 2.2 kcal/mol at shorter distances.

The inset of Fig. 3(a) shows the various potential energy curves computed for the high-spin state of the hexaaquairon(II) cation, with each curve vertically shifted by only the UKS-DFT minimum energy. This inset demonstrates relatively small differences in the total energies computed with the various embedding and open-shell treatments.

Finally, the dashed black curve in Fig. 3(a) demonstrates the importance of including spin-dependence in the embedding potential. This curve corresponds to the RO-DFT-in-DFT-CS treatment of the high-spin state of the hexaaquairon(II) cation described in Sec. III A 1. It exhibits large relative errors (over 70 kcal/mol) compared to the other treatments of the high-spin state of the hexaaquairon(II) cation, as well as qualitatively incorrectly shortening of the equilibrium internuclear distance. Although this approximation is expected to be more reliable for systems in which the spin-density is strongly localized with a single subsystem, the result demonstrates that substantial errors can emerge due to the neglect of spin-dependence in the embedding potential.

#### WFT-in-DFT Embedding

We now consider WFT-in-DFT embedding for the hexaaquairon(II) cation, employing the same subsystem partitioning as in the DFT-in-DFT embedding calculations (Fig. 3(b)). The hexaaquairon(II) cation is a benchmark system for spin splittings in transition metal complexes.(66) We initially discuss results for MP2 embedding to compare the U-WFT-in-DFT and RO-WFT-in-DFT approaches, and we then present results obtained using CCSD(T) embedding.

Fig. 4 presents results for the low-spin/high-spin energy difference () obtained using MP2, KS-DFT, and MP2-in-DFT embedding; detailed values are reported in Table 1. For KS-DFT calculations of , the energy for the high-spin state of the hexaaquairon(II) cation was obtained at the UKS-DFT level of theory. The WFT-in-DFT embedding energy for the low-spin state of the hexaaquairon(II) cation is obtained using closed-shell WFT-in-DFT (Sec. II B), while the high-spin state is treated using either U-WFT-in-DFT or RO-WFT-in-DFT (Sec. III A 2). The KS-DFT results (red in Fig. 4) exhibit strong dependence on the XC functional, with hybrid functionals underestimating to a somewhat lesser degree than the semi-local functionals.

Fig. 4 clearly illustrates that the RO-MP2-in-DFT results (blue) are in better agreement with the full MP2 calculation than the corresponding U-MP2-in-DFT results (green), particularly for semi-local XC functionals. Removal of spin-contamination in the WFT calculation reduces the energy of the high-spin state RO-WFT-in-DFT calculation with respect to that obtained using U-WFT-in-DFT.

Another important observation from Fig. 4 is that the dependence of on the DFT XC functional is greatly reduced in the embedding calculation, even though only the single transition metal atom is treated at the WFT level. The spread of values obtained at the KS-DFT level of theory is over 6000 cm, which is reduced by a factor of 3 in the RO-MP2-in-DFT embedding calculations.

Fig. 5(a) presents calculations of the low-spin/high-spin splitting obtained using WFT-in-DFT calculations at the RO-CCSD(T)-in-DFT level of theory; detailed values are reported in Table 2. For the reference calculation obtained at the full RO-CCSD(T) level of theory,(67) no T2 amplitudes were found to exceed 0.05, indicating that a single-reference description of the wavefunction is adequate. The general trend for the RO-CCSD(T)-in-DFT calculations is consistent with the results obtained from RO-MP2-in-DFT. It is again seen that the dependence of on the XC functional is substantially reduced using RO-CCSD(T)-in-DFT embedding, and the accuracy of the KS-DFT results are generally improved by treating the transition metal atom at the WFT level. For this system, the embedded RO-CCSD(T) calculation involves correlating significantly fewer electrons than the full RO-CCSD(T) calculation, and we found that the WFT step in the RO-CCSD(T)-in-DFT calculation required approximately 50 times less wall-clock time than the full RO-CCSD(T) calculation.

Fig. 5(b) shows that the LDA functional(68); (69) presents an interesting outlier compared to the other results in Fig. 5(a). Unlike the semi-local and hybrid functionals, RO-CCSD(T)-in-LDA calculations do not exhibit a significant improvement with respect to the corresponding KS-DFT result. We now show that this anomalous result arises from a density-based error in the LDA functional.

Fig. 5(c) and Fig. 5(d) present the charge on the Fe atom from a Mulliken population analysis for the low-spin state of the hexaaquairon(II) cation. Fig. 5(c) shows that the semi-local and hybrid functionals all yield a similar charge for the Fe atom, which is very close to that of the full (relaxed) CCSD density. In contrast, Fig. 5(d) reveals the LDA functional significantly underestimates the Fe atomic charge, which indicates a significant error in the calculation of the ground state density. Although the use of embedded WFT can be expected to overcome the error in the contribution to the spin-splitting energy due to the LDA functional, it can not overcome this error in the actual ground state density due to LDA.

To confirm this interpretation, we show that removing the error in the LDA density leads to improved WFT-in-DFT estimates for the spin-splitting energy, even if the LDA functional is still employed for the DFT contributions to the energy. In Fig. 5(b), the B-LYP+LDA result for WFT-in-DFT embedding (blue, triangle) is obtained by (i) calculating the embedding potential and the subsystem densities using the B-LYP XC functional, (ii) performing the embedded WFT calculation at the CCSD(T) level, and (iii) using the LDA functional and CCSD(T) to evaluate the respective DFT and WFT contributions to the total energy in Eq. 18. The corresponding B-LYP+LDA result for KS-DFT (red, circle) is obtained by calculating the total density using KS-DFT with the B-LYP XC functional and then using the LDA functional to evaluate the KS-DFT energy. As is seen in Fig. 5(d), the B-LYP treatment of the subsystem densities leads to the expected partial charge for the Fe atom; it avoids the error in the electronic density that is introduced using LDA. However, the spin-splitting energy obtained using the B-LYP+LDA result for KS-DFT is essentially no better than that obtained using KS-DFT with the LDA functional (Fig.5(b)), indicating that simply correcting the LDA error in the density is not enough to avoid the LDA error in the energies. Finally, Fig. 5(b) shows that the B-LYP+LDA result for WFT-in-DFT does exhibit a substantial improvement over the corresponding KS-DFT result; this confirms that WFT embedding is able to overcome energy-based errors due to the DFT XC functional, although it is less effective at overcoming density-based errors due to the DFT XC functional.

Although we have shown that WFT-in-DFT embedding with the subsystem partitioning shown in Fig. 3(b) generally leads to improved estimates for the low-spin/high-spin splitting energy over KS-DFT, the same does not hold true for calculated ligation energies of the hexaaquairon(II) cation. Ligation energies calculated using RO-CCSD(T)-in-DFT embedding are essentially unchanged from those obtained using KS-DFT with the corresponding XC functional; indeed, the mean absolute difference between the computed WFT-in-DFT and KS-DFT ligation energy is only 0.6 kcal/mol per ligand across the set of functionals that includes LDA, B-LYP, PBE, PW91, B3LYP, and PBE0. Unlike the spin-splitting energy, which is highly sensitive to the electronic structure of the Fe atom and is thus impacted by the WFT subsystem description, the ligation energy is dominated by interactions between the Fe atom and the water ligands; these inter-subsystem interactions are still treated essentially at the DFT level in WFT-in-DFT embedding. An improved description for the ligation energy could be obtained by simply expanding the number of electrons that are treated at the WFT level of theory, or by including two-body correlation corrections through an embedded many-body expansion description of the system.(29)

## Conclusions

In this work, we have introduced and demonstrated improved methods for the implementation of WFT-in-DFT calculations for open-shell systems and systems with low-lying virtual orbitals. A simple orbital-occupation-freezing technique is introduced to enable robust OEP calculations on systems with small HOMO-LUMO gaps, leading to accurate DFT-in-DFT and WFT-in-DFT embedding calculations on transition-metal complexes. Furthemore, the use of spin-dependent embedding potentials is shown to preserve the accuracy of open-shell DFT-in-DFT calculations in both the restricted and unrestricted orbital formulations, whereas neglect of the spin polarization leads to significant errors in both computed energies and geometries. WFT-in-DFT calculations on the hexaaquairon(II) cation reveal that the treatment of only the single transition metal atom leads to significant improvements in the accuracy of calculated spin-splittings, as well as marked reduction in the dependence of results on the DFT XC functional. Taken together, the exact embedding techniques reported and demonstrated here offer a promising approach to the robust treatment of systems for which the accuracy of WFT is required but for which the cost of the full WFT calculation is not feasible.

## Acknowledgements

This work is supported in part by the Air Force Office of Scientific Research (FA9550-11-1-0288) and the U. S. Army Research Laboratory and the U. S. Army Research Office (W911NF-10-1-0202). TFM and FRM also gratefully acknowledge network funding from the NSF (CHE-1057112) and EPSRC (EP/J012742/1), respectively. Computational resources were provided by the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231.

### Footnotes

1. RO-MP2 yields 16439 cm.
2. U-MP2 yields 17396 cm.
3. RO-CCSD(T) yields 14149 cm.

### References

1. A. Warshel and M. Levitt,  J. Mol. Biol. 103, 227 (1976).
2. P. Sherwood, A. H. de Vries, S. J. Collins, S. P. Greatbanks, N. A. Burton, M. A. Vincent, and I. H. Hillier,  Farad. Discuss. 106, 79 (1997).
3. J. L. Gao, P. Amara, C. Alhambra, and M. J. Field,  J. Phys. Chem. A. 102, 4714 (1998).
4. H. Lin and D. G. Truhlar,  Theor. Chem. Acc. 117, 185 (2007).
5. H. M. Senn and W. Thiel,  Angew. Chem. Int. Ed. 48, 1198 (2009).
6. L. Hu, P. Söderhjelm, and U. Ryde,  J. Chem. Theory. Comput. 7, 761 (2011).
7. S. Dapprich, I. Komáromi, K. S. Byun, K. Morokuma, and M. J. Frisch,  J. Mol. Struct. THEOCHEM 461-462, 1 (1999).
8. F. Maseras and K. Morokuma,  J. Comp. Chem. 16, 1170 (1995).
9. K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M. Uebayasi,  Chem. Phys. Lett. 313, 701 (1999).
10. D. G. Fedorov and K. Kitaura,  J. Chem. Phys. 120, 6832 (2004).
11. D. G. Fedorov and K. Kitaura,  J. Phys. Chem. A. 111, 6904 (2007).
12. M. J. Field, P. A. Bash, and M. Karplus,  J. Comput. Chem. 11, 700 (1990).
13. U. C. Singh and P. A. Kollman,  J. Comput. Chem. 7, 718 (1986).
14. Y. Zhang, T. S. Lee, and W. Yang,  J. Chem. Phys. 110, 46 (1999).
15. Y. Zhang,  J. Chem. Phys. 122, 024114 (2005).
16. G. Chałasiński, M. M. Szczȩśniak, P. Cieplak, and S. Scheiner,  J. Chem. Phys. 94, 2873 (1991).
17. E. B. Kadossov, K. J. Gaskell, and M. A. Langell,  J. Comput. Chem. 28, 1240 (2007).
18. G. Senatore and K. R. Subbaswamy,  Phys. Rev. B. 34, 5754 (1986).
19. M. D. Johnson, K. R. Subbaswamy and G. Senatore,  Phys. Rev. B. 36, 9202 (1987).
20. P. Cortona,  Phys. Rev. B 44, 8454 (1991).
21. T. A. Wesolowski and A. Warshel,  J. Phys. Chem. 97, 8050 (1993).
22. T. A. Wesolowski, in Computational Chemistry: Reviews of Current Trends - Vol. 10, pp. 1-82 (World Scientific, Singapore, 2006).
23. N. Govind, Y. A. Wang, A. J. R. da Silva, and E. A. Carter,  Chem. Phys. Lett. 295, 129 (1998).
24. P. Huang and E. A. Carter,  J. Chem. Phys. 135, 194104 (2011).
25. P. Elliott, K. Burke, M. H. Cohen, and A. Wasserman,  Phys. Rev. A 82, 024501 (2010).
26. J. Nafziger, Q. Wu, and A. Wasserman,  J. Chem. Phys. 135, 234101 (2011).
27. J. D. Goodpaster, N. Ananth, F. R. Manby, and T. F. Miller, III,   J. Chem. Phys. 133, 084103 (2010).
28. J. D. Goodpaster, T. A. Barnes, and T. F. Miller, III,   J. Chem. Phys. 134, 164108 (2011).
29. F. R. Manby, M. Stella, J. D. Goodpaster, and T. F. Miller, III,   J. Chem. Theory Comput. 8, 2564 (2012).
30. O. Roncero, M. P. de Lara-Castells, P. Villarreal, F. Flores, J. Ortega, M. Paniagua, and A. Aguado,  J. Chem. Phys. 129, 184104 (2008).
31. O. Roncero, A. Zanchet, P. Villarreal, and A. Aguado,   J. Chem. Phys. 131, 234110 (2009).
32. S. Fux, K. Kiewish, C. R. Jacob, J. Neugebauer, and M. Reiher,  Chem. Phys. Lett. 461, 353 (2008).
33. S. Fux, C. R. Jacob, J. Neugebauer, L. Visscher, and M. Reiher,  J. Chem. Phys. 132, 164101 (2010).
34. A. S. P. Gomes, C. R. Jacob, and L. Visscher,  Phys. Chem. Chem. Phys. 10, 5353 (2008).
35. C. Huang, M. Pavone, and E. A. Carter,  J. Chem. Phys. 134, 154110 (2011).
36. A. Solovyeva, M. Pavanello, and J Neugebauer,  J. Chem. Phys. 136, 194104 (2012).
37. C. R. Jacob, J. Neugebauer, and L. Visscher,  J. Comput. Chem. 29, 1011 (2008).
38. M. Iannuzzi, B. Kirchner, and J. Hutter,   Chem. Phys. Lett. 421, 16 (2006).
39. J. W. Kaminski, S. Gusarov, T. A. Wesolowski, and A. Kovalenko,   J. Phys. Chem. A. 114, 6082 (2010).
40. Ł. Rajchel, P. S. Żuchowski, M. M. Szczȩśniak, and G. Chałasiński,  Chem. Phys. Lett. 486, 160 (2010).
41. A. W. Götz, S. M. Beyhan, and L. Visscher,  J. Chem. Theory Comput. 5, 3161 (2009).
42. A. S. P. Gomes and C. R. Jacob,  Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 108, 222 (2012).
43. Q. Wu and W. Yang,  J. Phys. Chem. 118, 2498 (2003).
44. F. A. Bulat, T. Heaton-Burgess, A. J. Cohen, and W. Yang,  J. Chem. Phys. 127, 174101 (2007).
45. Q. Wu, P. W. Ayers, and Y. Zhang,  J. Chem. Phys. 131, 164112 (2009).
46. C. R. Jacob,  J. Chem. Phys. 135, 244102 (2011).
47. Q. S. Zhao and R. G Parr,  Phys. Rev. A 46, 2337 (1992).
48. Q. S. Zhao and R. G Parr,  J. Chem. Phys. 98, 543 (1993).
49. Q. S. Zhao, R. C. Morrison, and R. G Parr,  Phys. Rev. A 50, 2138 (1994).
50. S. Sherifzadeh, P. Huang, and E. A. Carter,  Chem. Phys. Lett. 470, 347 (2009).
51. T. Klüner, N. Govind, Y. A. Yang, and E. A. Carter,  J. Chem. Phys. 116, 42 (2002).
52. T. A. Wesolowski,  Phys. Rev. A 77, 012504 (2008).
53. P. Huang and E. A. Carter,  J. Chem. Phys. 125, 084102 (2006).
54. S. Sherifzadeh, P. Huang, and E. A. Carter,  J. Phys.: Condens. Matter 21, 355501 (2009).
55. A. Mordecai, in Nonlinear Programming: Analysis and Methods,  (Dover Publishing, 2003).
56. H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby, M. Schütz et al. Molpro, version 2012.1; Cardiff University: Cardiff, U. K.; Universitt Stuttgart: Stuttgart, Germany, 2012. See www.molpro.net.
57. T. A. Wesolowski,  J. Chem. Phys. 106, 8516 (1997).
58. J. P. Perdew, K. Burke, and M. Ernzerhof,  Phys. Rev. Lett. 77, 3865 (1996).
59. J. P. Perdew, K. Burke, and M. Ernzerhof,  Phys. Rev. Lett. 78, 1396 (1997).
60. A. D. Becke,  J. Chem. Phys. 98, 5648 (1993).
61. A. D. Becke,  Phys. Rev. A 38, 3098 (1988).
62. C. Lee, W. Yang, and R. G. Parr,  Phys. Rev. B 37, 785 (1988).
63. J. P. Perdew,  Phys. Rev. B 33, 8822 (1986).
64. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais,  Phys. Rev. B. 46, 6671 (1992).
65. C. Adamo and V. Barone,  J. Chem. Phys. 110, 6158 (1999).
66. A. Fouqueau, S. Mer, M. E. Casida, L. M. L. Daku, A. Hauser, T. Mineva, and F. Neese,  J. Chem. Phys. 120, 9473 (2004).
67. P. J. Knowles, C. Hampel, and H.-J. Werner,  J. Chem. Phys. 99, 5219 (1993).
68. J. C. Slater,  Phys. Rev. 81, 385-390 (1951).
69. S. H. Vosko, L. Wilk, and M. Nusair,  Can. J. Phys. 58, 1200-1211 (1980).
101706