Correcting density-driven errors in projection-based embedding

Density-functional theory (DFT) is an extremely powerful tool for exploring a range of chemical problems, including reactivity,^1–3 molecular spectroscopy,^4,5 the study of inorganic crystals,⁶ enzyme catalysis⁷ and simulations of materials,⁸ amongst many others. The attraction of the method is that it can produce accurate results with a relatively low scaling—cubic or lower—of computational cost with system size.

However, results can often be rather sensitive to the choice of the approximate exchange-correlation functional, leading to uncertain conclusions,^2,3,9 particularly when there are no experimental data to help inform this choice. The correlated post-Hartree-Fock (HF) method of coupled-cluster theory with singles, doubles and pertubative triples (CCSD(T)) is known to provide much more consistently reliable results, at least for systems whose electronic structure is dominated by a single reference determinant.¹⁰ 

The $O(N7)$ computational scaling of this method prevents application to large chemical systems without further approximations, and two basic strategies have evolved to address this: first, one can construct approximations that reduce computational scaling for the calculation on the whole system, as done in various local-correlation approaches.^11–19 A second strategy is to develop multiscale techniques in which a high-level calculation on a chemically important subsystem is embedded into a more approximate model describing the environment.^7,20–24

In common with many other groups, we find that DFT provides an excellent framework for constructing embedding theories,^25–30 based on a simple additive partition of the density into subsystem contributions,

The energy of the total system can then be calculated as

where the first two terms are the energies of the isolated subsystems and the third includes all nonadditive contributions. Most contributions to $ΔE[ρA,ρB]$ are straightforward, but the nonadditive-kinetic-energy term $ΔTs[ρA,ρB]$ is difficult. This term is often approximated^{27,28,31–36} but can also be calculated exactly through potential inversion methods,^31,37–40 although this remains numerically challenging.

Projection-based embedding avoids this issue by defining subsystem densities in terms of disjoint subsets of orthogonal orbitals such that $ΔTs[ρA,ρB]=0$⁠.³⁰ The fact that there is no non-additive kinetic energy to approximate, coupled with the simplicity of other terms, makes it straightforward to achieve the exactness property that a DFT-in-DFT calculation is equivalent to performing DFT on the whole system. An accurate wavefunction (WF) method (e.g., CCSD(T)) can then be applied to subsystem A, with interactions between subsystems and the energy contribution of subsystem B handled at the DFT level.

The coulomb and exchange-correlation potential of the subsystem B is added to the core Hamiltonian $H^A-in-B$⁠, which describes all interactions between subsystems A and B. It enforces orthogonality by using a level shift projector that shifts all orbitals of subsystem B towards an infinite energy (in practice, a large finite value is used). The total energy is then calculated by

where $ΨA$ is the wavefunction for subsystem A and E^B is the DFT energy of subsystem B. Further details can be found in our previous work.^30,41,42

Projection-based CCSD(T)-in-DFT embedding calculations for open-shell systems are performed with the following procedure:

1. perform a restricted open-shell DFT calculation on the entire system (HF or in general any self-consistent field (SCF) method can be used);

2. localize the occupied molecular orbitals;

3. select atoms and the number of electrons in subsystem A and determine the partitioning from the local molecular orbitals. All open-shell orbitals are taken to be in subsystem A, resulting in a closed-shell subsystem B;

4. form the embedded Hamiltonian $H^A-in-B$⁠;

5. perform restricted open-shell HF followed by a restricted open-shell CCSD(T) calculations on subsystem A.

This approach has been used to model adsorption of cobalt onto graphene analogues;⁴² in conjunction with combined quantum mechanics/molecular mechanics (QM/MM) to accurately model an enzyme active site;⁴³ in many-body expansions of neutral water clusters;^41,44 to investigate new cobalt-based catalysts for hydrogen evolution;⁴⁵ and to explore electrochemical stability and solvation effects in battery electrolytes.⁴⁶ 

Combinations of different methods, as in projection-based CCSD(T)-in-DFT, are often accurate but can begin to breakdown in some regimes. For example, it is possible that the chemical process in subsystem A induces a response in subsystem B that is incorrectly handled by the low level of theory in that region. Main sources of error in projection-based error have been analysed by Goodpaster et al., who found that in cases where embedding performs less well, the error is often dominated by the non-additive exchange-correlation energy.⁴⁷ Other errors such as those that occur from the use of DFT on subsystem B, or in the embedding potential, were observed to be relatively minor for the systems studied by Goodpaster et al.⁴⁷ In an earlier incarnation of our embedding work, we observed some cases where errors primarily arose from inaccuracies in the density computed using the local-density approximation, in calculations of spin-state splittings of the hexaaquairon(II) cation.⁴⁸ 

Here we investigate cases where large errors occur for CCSD(T)-in-DFT embedding but not for CCSD(T)-in-HF. We show that these errors arise from an inaccurate description of the density by DFT, which originate from the incomplete cancellation of the self-interaction component of the Coulomb and exchange-correlation energies. This is a well-known problem for DFT, leading to issues such as the overstabilization of hemibonded structures^49–57 or in modelling the interaction of charge-transfer complexes⁵⁸ and in systems with a non-integer number of electrons.⁵⁹ The problems can be reduced through self-interaction corrections^60–62 or by including a portion of exact exchange in a hybrid functional.⁶³ However, for many chemical problems, density errors driven by self-interaction are small and are less important than errors in the evaluation of the energy.⁵⁷ 

Burke et al. show that in some chemical systems such as odd-electron radical species, the self-interaction error from DFT is large enough that it dominates over energetic errors arising from using an approximate exchange-correlation functional.^57,64 They demonstrate that in such cases accuracy can be greatly improved by evaluating the approximate DFT energy using the qualitatively accurate HF density. These “density-corrected” DFT calculations produce potential-energy surfaces that are much more accurate than those produced by approximate DFT alone.^57,64

Here we show that a similar self-interaction effect occurs for some embedding calculations, which rely on an accurate density to perform an effective partitioning of the system into subsystems A and B. This could be of particular importance for calculations on transition metal clusters, which often have a large number of unpaired electrons.

It is well-known that pure generalized-gradient approximation (GGA) functionals incorrectly predict the water dimer cation to have the hemibonded [$H2O⋯OH2$]⁺ structure;^51,65,66 more sophisticated theoretical methods^{49,56,66–70} and vibrational spectroscopy⁷¹ show that the correct geometry is the contact-pair $H3O+⋯∙OH$⁠. The overstabilization of the hemibonded structure by GGAs is due to the self-interaction error, which leads to the spurious delocalization of the electron hole.^66,72,73 Where the self-interaction error is reduced, by inclusion of exact exchange^50,66,69,73 or self-interaction corrections,^67,74 DFT approximations typically predict the correct geometry.

This problem is known to extend into larger cationic water clusters, where many possible isomers can form,^70,73,75,76 representing various solvation structures for the dimer cation in the aqueous solution. For $(H2O)n+$ clusters with $n≤4$⁠, the lowest-energy structures have an adjacent hydronium ion and a hydroxyl radical, which is termed as a contact-pair isomer. For $n>4$⁠, a separated-pair isomer is preferred, with at least one water monomer separating the ion and radical.⁷⁶ For the purposes of using these systems to assess errors in embedding, it is useful to choose the contact-pair isomers, as the hydroxyl-hydronium ion pair is contained within a contiguous region. Therefore, the two tetramer (n = 4) contact-pair isomers were studied, with the structures shown in Figure 1.

Herr et al. found structure (1) to be lower in energy than (2) by 4.85 kcal mol⁻¹ using the EOM-IP-CCSD(dT)/CBS^DT method.⁷⁶ We found restricted open-shell CCSD(T)/aug-cc-pVTZ to be in excellent agreement, producing a relative energy of 4.80 kcal mol⁻¹ to use as a benchmark. The following DFT and coupled-cluster calculations are all restricted open-shell calculations to coincide with this benchmark number for the case where subsystem A encompasses the whole system.

Here we perform a range of embedding calculations on this system (computational details are given in the  Appendix). Figure 2 shows the error in relative energy with respect to CCSD(T) for various methods. The relative energies from the non-embedded SCF methods (HF and all six DFT functionals) vary in accuracy, with the energy of −0.6 kcal mol⁻¹ predicted by the local-density approximation (LDA) being in error by 4.2 kcal mol⁻¹ from CCSD(T). The hybrid functionals B3LYP (20% exact exchange), PBE0 (25%), and BH&HLYP (50%) perform better than the pure functionals LDA, BLYP, and PBE, showing the importance of exact exchange in the treatment of this system. As with predictions of the ground-state structure of the dimer cation,^50,66,69,73 the improvement seen in relative energy here with the inclusion of exact exchange is likely due to a reduction of self-interaction error.

CCSD(T)-in-HF shows excellent agreement with CCSD(T), with an error of just 0.5 kcal mol⁻¹. However, for these open-shell examples, CCSD(T)-in-DFT errors are strongly dependent on the exchange-correlation functional. This contrasts with our previous studies which showed a great reduction of sensitivity to the choice of functional in CCSD(T)-in-DFT calculations. When a pure GGA is used in subsystem B, the relative energies computed using embedding showed large errors, with a maximum of −10.8 kcal mol⁻¹ for LDA—significantly larger than the errors from DFT alone. The hybrid functionals produced much better embedding results, with CCSD(T)-in-DFT using PBE0, B3LYP, and BH&HLYP showing errors of less than 0.5 kcal mol⁻¹ from CCSD(T).

These results clearly show that whilst HF is reliable for embedding, DFT requires exact exchange to achieve a similar level of accuracy. This suggests that the self-interaction error of LDA and GGAs may play a role. To investigate this, we performed embedding calculations using the density-corrected DFT scheme suggested by Burke and co-workers.^57,64 These calculations are performed precisely as described in the Introduction, except that the molecular orbital for the whole system is computed using the HF theory, and the DFT calculation is performed without iterative update of these orbitals. HF-density-corrected DFT (which we denote HF-DFT) calculations produced mildly improved relative energies, though errors of up to 2.5 kcal mol⁻¹ from CCSD(T) remain (not shown in Figure 2). However, CCSD(T)-in-(HF-DFT) calculations are strikingly more accurate than CCSD(T)-in-DFT, with all relative energies within 0.4 kcal mol⁻¹ of CCSD(T). This confirms that embedding calculations are particularly sensitive to errors in the mean-field electronic density.

As further evidence that density errors underpin the poor performance of CCSD(T)-in-DFT for these systems, we show spin-density differences (relative to HF) in Figure 3. This illustrates the spurious delocalization of spin density from the hydroxy radical M2 to the oxygen atom of water-molecule M1. It can also be seen that the effect is significantly smaller for BH&HLYP than for BLYP, a difference that is reflected in the relative performance of these two functionals in embedding calculations, as shown in Figure 2.

It appears that the delocalization of spin density from monomer M2 to M1 decreases with increasing fraction of exact exchange. To confirm this, we investigated the degree of this delocalization in calculations using hybrid analogues of BLYP with

as a function of the fraction of exact exchange. Figure 4 shows the population of spin density on monomer M1 as a function of a. Populations were computed using the active-orbital charges of the intrinsic-bond-orbital (IBO) method—this method determines the fraction of the orbital that resides on each atom of the species and has been shown to be much less sensitive to the basis set than the Mulliken population method.⁷⁹ As expected, increasing the proportion of exact exchange in the functional decreases the degree of delocalisation onto the M1 oxygen systematically. The effect is much larger for structure (2) than structure (1), and the delocalization effect in the former dominates in the energy error.

Figure 5 clearly highlights the linear relationship between the degree of delocalisation of the singly occupied orbital (measured as the fractional population of the unpaired electron on monomer M1) and the error in the relative energy. Meta-GGAs such as the Minnesota functionals from the Truhlar group^9,80–82 are a popular choice of functional for many DFT studies. CCSD(T)-in-M06-L performs significantly better than other non-hybrid GGAs, with an absolute error below 2 kcal mol⁻¹. This error is again further reduced with the inclusion of exact exchange in M06, M06-2X, and M06-HF.

For this ionised water system, the self-interaction error causes DFT to produce a poor density. This causes some error in the DFT energies, but these remain at least in reasonable agreement with more accurate methods. However, embedded Hartree-Fock (and thus coupled-cluster) theory appears to be much more sensitive to deficiencies in the density, presumably due to the role that the density plays in determining the subsystem partition, and this can cause large errors in relative energy for CCSD(T)-in-DFT. The reduction of self-interaction (and delocalisation) by using hybrid functionals improves the density and the relative energies. However, an alternative solution is to use a density correction (HF-DFT), which dramatically reduces sensitivity on the choice of functional, as shown in Figure 2.

For the cases considered so far, the density errors in DFT arise from spurious delocalisation of unpaired electrons. Therefore, we now investigate the role of this error in open- and closed-shell systems with otherwise similar chemistry. The open-shell system here is the coupling of methyl and pentyl radicals to form hexane, and the closed shell analogue is the formation of a dative bond between borane and butylamine (Figure 6).

Subsystem A was initially chosen to include the two reacting atoms (radical carbons in hexane, boron, and nitrogen in the amine borane) and connected hydrogen atoms, thus including 18 electrons. The nature of the system allows the size of subsystem A to be systematically increased to include successive –CH₂– fragments. With this increase in the size of subsystem A, the reaction energy should converge to the CCSD(T) solution, which is reached exactly when subsystem A encompasses the whole system.

The results are shown in Figure 7. The left panel shows embedding results for the formation of hexane from open-shell reactants, and CCSD(T)-in-HF can be seen to converge to the CCSD(T) solution faster than CCSD(T)-in-BLYP, with respect to the size of subsystem A. This is most pronounced with the smallest subsystem A for this system (18 electrons, two chain atoms), where CCSD(T)-in-HF gives an error in relative energy of −0.2 kcal mol⁻¹, while CCSD(T)-in-BLYP gives an error of −7.3 kcal mol⁻¹.

As in Section II, this error for CCSD(T)-in-BLYP is observed to be a result of spurious delocalisation of the singly occupied orbital in the pentyl radical. HF gives an orbital charge on the first carbon atom outside subsystem A to be 0.96, compared to 0.92 for BLYP densities. Again, the density correction removes this delocalisation error, with the CCSD(T)-in-(HF-BLYP) error reduced to −2.0 kcal mol⁻¹. However, unlike the corresponding calculations in Section II, these density-corrected embedding calculations are still not as accurate as CCSD(T)-in-HF. This suggests that other non-density errors such as those described by Goodpaster et al.⁴⁷ come into play for this system, and that those errors are larger with DFT in subsystem B than with HF.

It should be noted that where three chain atoms are included in subsystem A, CCSD(T)-in-BLYP is considerably more accurate than non-embedded BLYP, regardless of the density error. Although density-correction should be used where relevant, we usually recommend that subsystem A is defined to include atoms that neighbour the reacting atoms (i.e., three chain atoms in this system) as a minimum to ensure reliable results are achieved.

The right-hand panel in Figure 7 shows the same embedding results for the formation of the amine-borane system, which involves only closed-shell species. CCSD(T)-in-HF is again seen to converge to CCSD(T) faster than CCSD(T)-in-BLYP. However, in this case, the error for CCSD(T)-in-BLYP using the smallest subsystem A is much smaller: 1.7 kcal mol⁻¹ compared to 0.6 kcal mol⁻¹ for CCSD(T)-in-HF. Furthermore, CCSD(T)-in-(HF-BLYP) shows no significant improvement on CCSD(T)-in-BLYP. This suggests that the errors observed for this system are not due to inaccuracy of the BLYP density, which is expected as there are no unpaired electrons. This is confirmed by the reverse calculation, where a BLYP density is used in an HF embedding calculation (CCSD(T)-in-(BLYP-HF)). An error of 0.7 kcal mol⁻¹ is produced, which is in very close agreement with the CCSD(T)-in-HF calculation. Therefore, the error presumably comes from another source and is again seen to be somewhat larger with BLYP than with HF.

Projection-based embedding³⁰ is an effective method for performing wavefunction-in-density-functional-theory (WF-in-DFT) calculations, and even when subsystem A is small, coupled-cluster-type accuracy is typically approached. Previous work has shown that the main sources of embedding error are either small or can be improved with correction terms,^30,47 and WF-in-DFT can be used to greatly reduce dependence on the choice of approximate density functional.⁴³ 

For some open-shell systems, the self-interaction error in approximate density functionals leads to qualitative errors in the electronic density. In subsequent projection-based embedding calculations nothing that is done in subsystem A can compensate for the global inaccuracy of the density, so large energetic errors arise. The inclusion of exact exchange reduces the delocalization error of approximate DFT calculations, but a more systematic and satisfactory solution is obtained by using the density-corrected DFT.^57,64 The use of density-corrected DFT for embedding eliminates the major pathology introduced by the delocalization error, producing more consistently reliable accuracy across closed- and open-shell systems. This approach can be expected to improve the reliability of embedding calcualtions in cases where delocalization error arises, for example, in dissociative processes involving an odd number of electrons;^59,72 dissociation of metal halides, even when the electron number is even;⁸³ and in the computation of transition-state energies.

We are grateful to EPSRC for funding: RCRP is supported through the Doctoral Training Grant, S.J.B. through EP/K018965/1, EP/M022129/1, and the BBSRC through BB/L018756/1. T.F.M. acknowledges support from the National Science Foundation (NSF) CAREER Award under Grant No. CHE-1057112 and from the U.S. Army Research Laboratory under Grant No. W911NF-12-2-0023.

All data produced in this work are available in an open-access repository.¹⁰² 

Geometries for the ionised water clusters (Section II) were taken from the CCSD(T)/aug-cc-pVTZ optimised structures of Herr et al.⁷⁶ Geometries for the species involved in the reactions of hexane and butylamine borane (Section III) were optimised using the Gaussian09 Revision D.01 software package⁸⁴ with B3LYP-D3^85–87 in the 6-31G(d)^88,89 basis set. Single-point energies were obtained using HF, DFT (LDA,⁹⁰ BLYP,^86,91 PBE,⁹² B3LYP, BH&HLYP,^86,93 PBE0,^94,95 M06-L,⁸⁰ M06,⁸² M06-2X,⁸² M06-HF⁸¹), and CCSD(T),^96,97 as well as CCSD(T)-in-SCF embedded single-point energy calculations using the SCF methods listed above. All of these single-point calculations were performed using the MOLPRO 2015.1 software package.^98,99 Whilst the basis-set truncation is available to reduce the computational cost of projection-based embedding,^41,44 the calculations here were small enough to be performed in the full basis. All single point energy calculations were performed using the aug-cc-pVTZ^100,101 basis set. The IBO method with the zero bond-dipole scheme⁷⁹ was used to produce local orbitals for defining subsystems A and B.