A Simple, Exact Density-Functional-Theory Embedding Scheme.

Density functional theory (DFT) provides a formally exact framework for quantum embedding. The appearance of nonadditive kinetic energy contributions in this context poses significant challenges, but using optimized effective potential (OEP) methods, various groups have devised DFT-in-DFT methods that are equivalent to Kohn-Sham (KS) theory on the whole system. This being the case, we note that a very considerable simplification arises from doing KS theory instead. We then describe embedding schemes that enforce Pauli exclusion via a projection technique, completely avoiding numerically demanding OEP calculations. Illustrative applications are presented using DFT-in-DFT, wave-function-in-DFT, and wave-function-in-Hartree-Fock embedding, and using an embedded many-body expansion.

Introduction

Kohn–Sham (KS) density functional theory (DFT) provides a powerful theoretical framework for performing calculations on a subsystem exactly embedded in its full, quantum-mechanical environment.[1−3] The quantity that mediates the interaction between subsystems is the electronic density, which is partitioned into two terms

The KS energy consists of terms relating to subsystems A and B and an interaction term containing all nonadditive parts of the energy:This last term includes an explicitly known Coulomb contributionand, depending on the details of the subsystem partitioning, simple electrostatic interaction of electrons in A with nuclei in B and vice versa. Exchange-correlation effects between subsystems are included throughwhich of course must be approximated in practical calculations.

If the subsystem densities ρA and ρB are constructed from mutually orthogonal orbitals, the kinetic energy of the whole system is simply given by the sumIn general, however, there will also be a nonadditive termwhich must be approximated[3,4] or, in more recent versions of the theory, computed exactly using optimized effective potential (OEP) methods.[5−9]

The current authors have demonstrated that OEP methods can be used in an iterative scheme to obtain the Kohn–Sham ground state density for molecular systems, including those that exhibit strongly overlapping subsystem densities.[5,6] Reiher and co-workers used OEP methods to analyze and compare accurate embedding potentials in approximate embedding schemes,[7] and OEP methods were used to solve for a unique partitioning of a system by allowing the subsystems to share a common embedding potential.[8,9] Finally, several groups have utilized these DFT embedding methods to describe the interface between high-level and low-level subsystems in wave-function-in-DFT (WF-in-DFT) approaches.[9,11,12]

In what follows, we avoid the complications of OEP through three simple, robust innovations: (1) We replace the iterated DFT-in-DFT with a single conventional KS calculation. (2) We completely avoid the issue of kinetic energy nonadditivity through the use of a level shifting projection operator to keep the orbitals of one subsystem orthogonal to those of another. And, (3) we develop a simple but effective perturbation theory to eliminate practically all dependence on the level shift parameter.

Embedding methods that maintain orthogonality between subsystem orbitals have been in use for decades. What has not been recognized is that these can be used to formulate a formally exact DFT embedding scheme, equivalent to but much simpler than the recently developed OEP-based methods. The Philips–Kleinman pseudopotential approach starts by level shifting the core orbitals[13] to produce valence orbitals implicitly orthogonalized to the core. Even frozen core approximations, which have been in use in quantum chemistry at least since the 1950s,[14] amount to a sort of Hartree–Fock (HF)-based embedding scheme in which the core and valence subsystems are described by mutually orthogonal orbitals. Similar ideas operate in the incremental scheme introduced by Stoll and co-workers,[15] in the region method of Mata and et al.,[16] and in Henderson’s embedding scheme.[17] Further details of the large amount of literature in this area can be found in a recent and comprehensive review.[18]

The Pauli principle is enforced through a level-shift projector in the ab initio model potential (AIMP) method,[19−21] which in turn derives from earlier work by Huzinaga and Cantu.[22] A similar strategy was employed by Rajchel et al. for the description of intermolecular interactions, but there they iterated equations for each subsystem.[10] Projection is sometimes employed in the embedding used in the fragment molecular orbital method (see for example the recent article ref (23) and references therein). The novelty of the current work resides in the important simplification that arises by performing a conventional Kohn–Sham calculation plus localization to determine fragment orbitals, the recognition that this amounts to a formally exact DFT-in-DFT embedding scheme, and the elimination of practically all dependence on the level-shift parameter through perturbation theory (see below).

Simplified Exact Embedding

Self-consistently iterated embedded DFT calculations are exactly equivalent to full KS theory on the whole system. If the only aim is to perform DFT-in-DFT embedding, it is currently more practical to perform KS theory on the whole system. We now demonstrate that this remains true even if the initial mean-field calculation is merely a prelude for WF-in-DFT embedding.

The sum of the atomic-orbital (AO) density matrices from the embedded one-electron equations is equal to the total KS density matrix:[2] γA + γB = γ. The occupied KS orbitals {ϕ} can be rotated among each other, and one rotation produces the subsystem orbitals {ϕA} ∪ {ϕB}. Let us assume that the whole system and both subsystems A and B are closed-shell and that the total number of electrons is partitioned as nA + nB = n; the extension of the following arguments to open-shell subsystems is straightforward.

To produce the correct total electron density, a KS calculation on subsystem A requires an embedding potential that correctly describes the electrostatic and exchange-correlation interactions with the electrons in subsystem B, while also properly enforcing Pauli exclusion between the subsystem densities. This is in fact the role of the nonadditive kinetic potential vnad,[3] but an alternative scheme can be reached by simply level shifting the orbitals of B to very high energy for the calculation on A. Then, a correctly embedded calculation on the nA electrons in subsystem A has the energy expressionwhere h is the core Hamiltonian (kinetic plus external potential) in the AO basis.

We propose a Fock operator for electrons in subsystem A that includes the appropriate derivative of the energy, plus a level-shift term to preserve the orthogonality relation ⟨ϕA|ϕB⟩ = 0:whereThe term μPB elevates the energy of the ith orbital in subsystem B to εB + μ, and provided that μ is sufficiently large and positive, the eigenstates of f A are then mutually orthogonal to the orbitals in subsystem B.

Evaluation of the gradient in eq 2 produces the Fock operator for subsystem A that we employ in the calculations below:In HF theory or in KS theory with hybrid functionals, the additional (possibly scaled) exchange term −K[γA + γB] appears in the usual way.

As in other DFT-in-DFT embedding schemes, this approach provides the flexibility to perform a completely different type of electronic structure calculation for subsystem A, using a core Hamiltonian that contains the additional embedding terms from above:The total energy for a normalized wave function in the A subsystem is then calculated aswhere ĤA in B is the Hamiltonian built with the embedding potential added to the core Hamiltonian, the central terms correct the exchange-correlation energy contribution, and the final term is the energy of the electrons in the B subsystem.

The code is implemented in the development version of the Molpro software package.[24,25] The current version of the code, and all calculations presented here, employ the supermolecular basis set. Future versions will employ subsystem calculations in a spatially local set of atomic orbitals; we see no technical barriers to this extension.

Example Calculations

Accuracy and Dependence on μ

We begin by demonstrating the efficacy of the projection technique described above, using the example of embedding the 10 electrons of the −OH moiety of ethanol in the environment produced by the ethyl subsystem. KS orbitals from a PBE[26] calculation in the 6-31G* basis[27,28] are localized using the Pipek–Mezey localization scheme,[29] and the five orbitals clearly associated with the −OH group are identified. Note that this step is performed automatically by finding those orbitals with a Mulliken population on any “active” atom greater than 0.4. The embedded core Hamiltonian matrix hA in B is then constructed according to eq 3, with μ = 103Eh, and a 10-electron Kohn–Sham calculation is performed.

If the projection technique has functioned perfectly, the resulting orbitals will be orthogonal to those associated with subsystem B, the energy expression in eq 1 will yield precisely the KS energy of the whole system, and μ tr γAPB will exactly vanish.

We find this not quite to be the case (see Table 1). The value of μ tr γAPB is approximately 20 μEh, and the energy is somewhat too low. In the limit of infinite μ, the agreement with Kohn–Sham theory will be exact by construction, so the error arises from the finite value of this parameter. Upon increasing the value of μ, we do indeed find that the error in the energy decreases (Figure 1, uncorrected), although very large values result in numerical noise associated with the double-precision numerics of the program.

Table 1

PBE/6-31G* Energies for Ethanol, with and without Embedding, Where Subsystem A Is Comprised of the Five Pipek–Mezey Local Orbitals Associated with the −OH Moiety

	energy/Eh
EKS	–154.82798488
E[γA;γB]	–154.82800669
μ tr γAPB	0.00002181
E[γA;γB]+μ tr γAPB	–154.82798488

Figure 1

Error in the uncorrected (E[γA;γB]) and corrected (E[γA;γB] + μ tr γAPB) energy expressions relative to full KS on ethanol using PBE/6-31G*, demonstrating that the perturbation theory correction yields essentially exact embedding energies over a wide range of μ.

A more stable approach is to compute the μ → ∞ limit using perturbation theory. We can introduce a perturbed operator f̂ζ that satisfies the requirements f̂0 = f̂ + μP̂ andby adopting the form(Here, for clarity, we have dropped the A and B superscripts.) Expanding the operator f̂ζ as a power series in ζ and making the usual perturbative expansion of energy and wave function, we obtainwhich yields at first-order the contribution E1 = ⟨0|μP̂|0⟩. In the AO basis E1 = μ tr γAPB, and as shown in Table 1, this correction brings the embedded energy into near-perfect agreement with the KS energy. Higher orders in this perturbation expansion could be used to correct the embedded orbitals, or to further correct the energy, although nothing in our calculations so far suggests that this will be necessary.

In fact, it turns out that the energy expression E[γA;γB] + μ tr γAPB agrees with the full KS energy to 7 pEh, which is essentially exact in double precision. This excellent agreement is very insensitive to the choice of μ, as shown in Figure 1. For values μ < 102Eh, the electronic structure is qualitatively incorrect, and for μ > 107Eh numerical noise begins to become an issue. But for all values of μ in between, spanning 5 orders of magnitude, the discrepancy from full KS theory is below 20 nEh.

Similar results are found for other molecules, and the projection approach can also be employed for HF embedding. For example, a HF-in-HF embedded calculation on pyradine with subsystem A composed of the five Pipek–Mezey orbitals associated with the nitrogen atom and with μ = 104Eh reproduces the total HF energy to within 0.1 nEh. The method is also very stable for WF-in-DFT calculations, which we will demonstrate through CCSD(T)-in-PBE calculations on ethanol, using the same subsystem partitioning as before. The calculations are performed (i) without the perturbative correction to the total energy, (ii) with the pertubative correction obtained using the HF subsystem density for the active subsystem, and (iii) with the perturbative correction obtained using the unrelaxed CCSD density for the active subsystem; the differences between the total energies with μ = 106 and μ = 105 are then 470 nEh, 50 nEh, and 0.74 nEh, respectively.

Chemical Reactions

A key motivation for DFT embedding schemes is the possibility of embedding accurate wave function calculations in a DFT environment. Here, we demonstrate the simplicity and robustness of the current scheme with a few, sample calculations.

First, we consider the deprotonation energy of gas-phase ethanol, an example for which the PBE/aug-cc-pVDZ level of theory is in fairly good agreement with CCSD(T)/aug-cc-pVDZ. Two embedded calculations are performed, using either −OH (10 electrons) or −CH2OH (18 electrons) as the active subsystem. The results are shown in Table 2, where it can be seen that embedded CCSD(T)-in-PBE with the 10-electron active subsystem slightly overcorrects the PBE value, by approximately 6 mEh. The calculation with the 18-electron active subsystem produces a result in much better agreement, differing from the total CCSD(T) energy by a little over 1 mEh. These results are not entirely unexpected, since in the first calculation the system is partitioned across a bond which changes polarization during the deprotonation process.

Table 2

Deprotonation Energies for Ethanol Using PBE/aug-cc-pVDZ and CCSD(T)/aug-cc-pVDZ and Combinations of These Using Embedding

	subsystem A	E/mEh
PBE		611.2
CCSD(T)-in-PBE	–OH	627.5
CCSD(T)-in-PBE	–CH2OH	622.8
CCSD(T)		621.3

Second, we examine the symmetric SN2 reaction of Cl– with propyl chloride. The reactant and transition state structure were optimized at the B3LYP/6-311G*++ level[28,30] of theory using Gaussian.[31] Three coupled-cluster, KS, and embedded activation energies are shown in Table 3. Activation energies calculated at the CCSD(T), DFT, and CCSD(T)-in-DFT levels of theory are shown in Table 3; each of these calculations is performed using the cc-pVTZ[32] basis set with aug-cc-pV(T+d)Z on chlorine.[33] DFT calculations with approximate functionals often underestimate barriers, but in the present case, the PBE and BLYP barriers are actually negative. HF theory overestimates the barrier. The embedded results are consistently much closer to the reference coupled-cluster results, with the accuracy improving as more atoms are included in subsystem A.

Table 3

Activation Barrier of the Reaction of Propyl Chloride with Chloride Anion Using HF, KS, CCSD(T), and Combinations Using Embeddinga

		Ea/mEh
	subsystem A	PBE	BLYP	HF
DFT/HF		–1.2	–1.9	14.6
CCSD(T) embedded	–CH2Cl2–	9.8	12.7	13.2
CCSD(T) embedded	–(CH2)2Cl2–	8.8	9.4	8.6
CCSD(T)		7.8	7.8	7.8

The basis set is cc-pVTZ with aug-cc-pV(T+d)Z on chlorine.

Accurate Embedded Many-Body Expansion

Recent research in the group of one of us has focused on the use of embedded many-body expansions to compute accurate structure and energetics of molecular crystals.[34] The idea builds on earlier attempts to capture many-body effects through embedding,[35,36] and it offers an important step toward systematically improving accuracy by improving the embedding scheme, rather than the number of terms in the many-body expansion. We illustrate this approach for the case of a symmetric model water trimer.

For simplicity, the trimer is optimized in a planar geometry, in which each of the three monomers have equivalent geometries and environment (see structure in Figure 2). The geometry is optimized with these contraints using MP2/6-31G*. Denoting the monomer, dimer, and trimer CCSD(T) energies as E1CC, E2CC, and E3CC, we can immediately express the coupled-cluster binding energy (relative to water molecules in their cluster geometry) in the form EbindCC = E3CC – 3E1CC. Calculations are performed using the aug-cc-pVDZ basis throughout.

Figure 2

Binding curve relative to monomers for a planar water trimer as a function of the distance from the center of mass to the oxygen atoms. Shown are the binding curves from HF theory, CCSD(T) many-body expansion truncated at the two-body level (MBE2), embedded two-body expansion (EMBE2), and from full trimer CCSD(T). Note that EMBE2 is indistinguishable from CCSD(T) on the scale shown. In the lower panel, the error for EMBE2 is compared to that of the simple additive scheme HF-Δ12, in which the HF binding curve is corrected using the two-body CCSD(T) correlation energies.

Several important approximations can be defined without embedding, and two of these are shown in Figure 2. Truncating the many-body expansion at the two-body level gives the approximationwhere δE2CC = E2CC – 2E1CC. The HF approximation to the binding energy is defined by EbindHF = E3HF – 3E1HF. Both of these approximations result in errors of around 5 mEh in the total binding energy of around 20 mEh.

Many-body effects can be incorporated through embedding. Since our HF-in-HF embedding calculations require HF on the full trimer, we focus our attention on the correlation contribution to binding, and we use the notation Ecorr ≡ ECC – EHF for correlation energies. The one-body correlation correction without embedding vanishes (by analogy to eq 4), but with embedding, the contribution added to the trimer energy is different from that subtracted as part of the monomer energy. Thus, one obtains the approximationwhere Ẽ1corr is the correlation energy from a CC-in-HF calculation on a single monomer. A further correction can be made by considering the two-body embedded correlation contribution, givingwhere δẼ2corr = Ẽ2corr – 2Ẽ1corr.

This EMBE2 approximation can be seen in Figure 2 to be in essentially perfect agreement with CCSD(T) across all geometries. Also shown in the figure is the error in this scheme (compared to EbindCC), which can be seen to be on the order of 10 μEh in the binding region. The results are strikingly better than those obtained using the simple unembedded HF-Δ12 approach defined as[37]EbindHF-Δ = EbindHF + 3δE2corr, especially at short range.

Discussion

The ideas presented here derive from a long history of ensuring that the Pauli principle is obeyed through level shifting a subset of the orbitals. But the interesting and novel point is that the DFT embedding techniques recently developed by us[5,6] and by others[7−10] can be cast exactly into the very simple and robust framework described here.

It should be noted that this projector embedding method is not without drawbacks. It is limited to applications for which the electronic structure can be reasonably described using KS theory, because the projector can only be formed if KS orbitals are available. This limitation will become problematic in cases for which the DFT description is fundamentally broken, or if it is impractical to perform a full KS calculation on the entire system. Furthermore, the new method describes embedding in the environment by a projection operator, rather than a simple, local potential.

In its favor, however, the current approach relies on existing, stable, well-developed technologies. It requires nothing more than a KS code, an orbital localization scheme, some elementary matrix operations in the atomic orbital basis, and any wave function method that can accept an arbitrary core Hamiltonian. Since these ingredients are available in practically all molecular electronic structure codes, we anticipate no barriers to widespread adoption of the approach.