Extension of frozen-density embedding theory for non-variational embedded wavefunctions

Multi-scale simulation methods using a density-dependent embedding potential became popular in recent years.^1–7 Some of them use the formal framework of Frozen-Density Embedding Theory (FDET),^1,3 providing self-consistent expressions for the total-energy, the variationally obtained wavefunction representing embedded N_A electrons (Ψ_A), and the embedding potential. In particular, the key relation of FDET equating the total energy evaluated for the embedded wavefunction at a given $ρB(r→)$ with the total energy expressed by means of the Hohenberg-Kohn density functional [see Eq. (1) below] holds only for embedded wavefunctions derived from variational methods such as Kohn-Sham, generalized-Kohn-Sham, Hartree-Fock, multiconfigurational self-consistent field, and configuration interaction (CI). For non-variational methods, the use of the FDET expression for the total energy (see, for instance, the expressions for the total energy given in Refs. 4–6 and 8), represents an additional approximation. The numerical significance of this approximation might be entangled with other approximations used in a particular method and vary from system to system. The formal link between FDET and the quantities obtained in non-variational calculations remains, however, to be established. In the present work, we propose and derive a formal extension of FDET to non-variational methods to treat embedded N_A electrons.

For a system comprising N_AB electrons in an external potential $vAB(r→)$⁠, the functional $EvABEWFΨA,ρB$ is defined to satisfy

where $EvABHK[ρ]$ is the Hohenberg-Kohn ground-state energy functional and $ρAo(r→)=ΨAo∑iNAδ(r→−r→i)ΨAo$⁠.

By virtue of the second Hohenberg-Kohn theorem, Eq. (1) leads to

where $E0=EvABHK[ρ0]$ and $ρ0(r→)$ is the ground-state energy and density of the total system. Equality is reached for a large class of densities $ρB(r→)$⁠,

Using conventional density functionals representing components of the total energy, and arbitrary partitioning of the external potential $vAB(r→)=vA(r→)+vB(r→)$⁠, leads to the form of $EvABEWFΨA,ρB$ more suitable for further discussions,

where

and $VNANB$ is the interaction energy between the nuclei defining $vA(r→)$ and $vB(r→)$⁠. The non-additive bi-functional $ExcTnadρA,ρB$ is related to the functionals E_xc[ρ] and T_s[ρ] defined in the constrained search formulation of the Kohn-Sham formalism,⁹ 

The functional $ΔFρ$ on the other hand depends on the form of the wavefunction Ψ used in Eq. (1) and also is defined via the constrained search.³ If Ψ_A is a single determinant (Φ), it reads

In the last line, we indicate that $ΔFρ$ is the correlation functional defined in constrained-search formulation of density functional theory.^9,10 For Ψ of the full CI form, $ΔFρ=0$ by definition.

Euler-Lagrange optimisation of Ψ_A leads to a Schrödinger-like equation,

where

with $vxcTnad[ρA,ρB](r→)$⁠, and $vF[ρA](r→)$ being the first functional derivatives of $ExcTnadρ,ρB$ and $ΔFρ$⁠, respectively. The lowest energy solution of Eq. (10) will be denoted as $ΨAEL$⁠. Note that the energy is given not by the Lagrange multiplier λ but in Eq. (4). Any variational method can be used to obtain $ΨAEL$ and the corresponding density $ρAEL(r→)$⁠.

For a given density, $ρ(r→)$⁠, we can define the corresponding wavefunction as

Up to a unitary transformation, it holds that

As a result, the second equality in Eq. (1) holds for variationally obtained wavefunctions $ΨAEL$ and the corresponding density for any embedding potential. We should note that in Eq. (11), $v^emb$ depends on $ρA(r→)$⁠. The first equality in Eq. (1) holds, therefore, only if $v^emb$ and $ΨAEL$ are self-consistent (⁠$v^emb=vemb[ρAEL,ρB,vB](r→)$⁠). For lower-cost variational methods, this can be achieved in practice in an iterative procedure.¹¹ 

Owing to the fact that any N_A electron wavefunction is an implicit (up to unitary transformation) functional of the density, the right-hand-side of Eq. (4) is an implicit functional of $ρA(r→)$⁠,

Let $ΨAnvar$ and $ρAnvar(r→)$ denote the embedded wavefunction and the corresponding density obtained from some non-variational method to solve Eq. (10). Contrary to the variational case,

As a result, the FDET energy expression given in Eq. (14) cannot be applied unless the following constrained search is performed:

Performing such a search is, however, impractical. More importantly, for methods in which the total energy is not evaluated as the expectation value of the Hamiltonian, Eq. (4) cannot be used. Moreover, the density is either not available in a straightforward manner (in coupled-cluster type of approaches, for instance) or it is available in the lower order than energy [Møller-Plesset perturbation theory (MP), for instance].

In the following, we exploit the fact that the right-hand side of Eq. (14) represents an implicit density functional, which can be expanded in a series at $ρAEL(r→)$⁠. For small differences, $Δρ(r→)=ρAnvar(r→)−ρAEL(r→)$ retaining the linear terms leads to

Besides the last two terms, other linear terms in $Δρ(r→)$ do not survive if $ρAEL(r→)$ is obtained variationally and with the self-consistent embedding potential.

Equation (17) is our principal result. It is exact up to first order in $Δρ(r→)$ provided the embedded wavefunction is obtained variationally with the self-consistent embedding potential. Equation (17) can be used within various non-variational treatments of the embedded N_A electrons and approximations to the explicit density functionals needed in practice. In a method, for which a correction to the energy obtained from variational calculations (correlation energy from Møller-Plesset or coupled-cluster calculations, for instance) is available, it can be expected to provide a good approximation to $ΔFρAEL$⁠. If a method yields also $Δρ(r→)$ (at lower order than the energy in MP calculations, for instance), it can be used in the last two terms.

This section illustrates the usefulness of Eq. (17) as a basis for approximate methods. For several weakly bound intermolecular complexes, the total energies are obtained using the following choice of the simplest approximations in Eq. (17). $ΔFρAEL$ is approximated by means of the correlation energy from second-order Møller-Plesset perturbation theory (MP2). For $ExcTnad[ρAEL,ρB]$⁠, local density approximation is used (decomposable approximations using Dirac-Slater functional for exchange energy,^12,13 the parameterisation of the Ceperley-Alder correlation energy¹⁴ by Vosko et al.,¹⁵ and Thomas-Fermi functional for kinetic energy.^16,17 The use of local density approximation for all relevant functionals was motivated by expected cancellation of errors (sharing the same approximation for the averaged exchange-correlation hole) confirmed in numerical tests on a representative set of pairs of densities $ρA(r→)$ and $ρB(r→)$⁠.^18,19 In the absence of well-tested approximations for the second functional derivatives, the $Δρ(r→)$-dependent contributions to energy are neglected. Preliminary numerical results using local density approximation indicate that these terms are indeed negligible. Concerning $ρB(r→)$⁠, it was optimised using the iterative freeze-and-thaw procedure.²⁰ In view of practical applications, we discuss below not the total energies obtained from the defined above approximated version of Eq. (17) but the interaction energies. The FDET interaction energies obtained in this way are compared to conventional counterpoise corrected MP2 interaction energies. Optimisation of $ρB(r→)$ is made in order to attribute the errors in the FDET interaction energies only to approximations made in Eq. (17). To validate numerical robustness of Eq. (17), the reported FDET energies are evaluated using three atomic basis sets (cc-pVDZ, cc-pVTZ, and cc-pVQZ²¹) in two types of multi-center expansion: either localised on individual molecules (monomer expansion, ME) or on all atoms in the system (supermolecular expansion, SE). The FDET interaction energies are calculated as

where $EBemb$ is the total MP2 energy of subsystem B used as $EvBHKρB$ in Eq. (17) and $EBfree$ is the total MP2 energy of the isolated subsystem B. Note that if $ρB(r→)$ is not optimized, the difference between these two energies is zero. $EA,emb(2)$ and $EA(2)$ denote the correlation energies obtained from second-order perturbation theory using embedded and environment-free reference states, respectively.

The results for four representative intermolecular complexes at their equilibrium geometries are reported in Figs. 1 and 2. FDET energies show a remarkable agreement with the reference interaction energies for each basis set. Rather small variation of the FDET interaction energy and its components with the used basis set indicate numerical robustness of the method. Note also that the contributions to the energy, which are approximated using explicit functionals (⁠$Excnad[ρAEL,ρB]$ and $Tsnad[ρAEL,ρB]$⁠), are numerically significant, opposite in sign, and of similar magnitude as the classical electrostatic interactions. This provides an additional proof of the good cancellation of errors in $Excnad(LDA)[ρAEL,ρB]$ and $Tsnad(LDA)[ρAEL,ρB]$⁠. If other approximations for these functionals would be considered, the balance in their quality is more important than improving the accuracy of only one of them. The cause of this remarkable performance is the object of our current research and will be reported in a full publication.

This research was supported by a grant from the Swiss National Science Foundation (Grant No. 200020-172532).