Exchange–correlation functionals via local interpolation along the adiabatic connection
Abstract
The construction of density-functional approximations is explored by modeling the adiabatic connection locally, using energy densities defined in terms of the electrostatic potential of the exchange-correlation hole. These local models are more amenable to the construction of size-consistent approximations than their global counterparts. In this work we use accurate input local ingredients to assess the accuracy of range of local interpolation models against accurate exchange-correlation energy densities. The importance of the strictly-correlated electrons (SCE) functional describing the strong coupling limit is emphasized, enabling the corresponding interpolated functionals to treat strong correlation effects. In addition to exploring the performance of such models numerically for the helium and beryllium isoelectronic series and the dissociation of the hydrogen molecule, an approximate analytic model is presented for the initial slope of the local adiabatic connection. Comparisons are made with approaches based on global models and prospects for future approximations based on the local adiabatic connection are discussed.
I Introduction
Kohn–Sham density-functional theory (KS DFT)Kohn and Sham (1965) is the method most widely used in electronic structure calculations, due to its modest computational cost combined with an accuracy that is often competitive with much more expensive ab initio methods. The accuracy of the method is limited by the quality of approximations required to describe the quantum mechanical exchange and correlation (XC) interactions of electrons. A large number of density functional approximations (DFAs) for the XC–energy have been developed in recent decades.
The simplest DFAs are based on the local density approximation (LDA), as proposed by KS in their 1965 paper,Kohn and Sham (1965) in which the XC–energy is approximated as a functional of the density at a given point in space. The generalised gradient approximations (GGAs)Lee, Yang, and Parr (1988); Becke (1988); Perdew et al. (1992); Perdew, Burke, and Ernzerhof (1996); Zhao and Truhlar (2008) go beyond the LDA by modelling the XC–energy as a functional of the local density and its first derivative. The meta-GGAsTao et al. (2003); Peverati and Truhlar (2014); Sun, Ruzsinszky, and Perdew (2015) are closely related but their functional forms are also dependent on the KS kinetic energy density and/or, less commonly, the Laplacian of the electron density. Further developments led to the introduction of the occupied KS orbitals as ingredients for the XC energy (hybrid functionals,Becke (1993); Perdew, Ernzerhof, and Burke (1996); Zhao, Schultz, and Truhlar (2006) self-interaction corrections,Perdew and Zunger (1981); Vydrov and Scuseria (2004, 2005); Kümmel and Kronik (2008)), and more recently also the virtual KS orbitals (double–hybrid functionals,Grimme (2006); Goerigk and Grimme (2010); Sharkas, Toulouse, and Savin (2011) random-phase approximationsFurche (2001); Hesselmann and Görling (2010); Ángyán et al. (2011)). Local hybrid functionalsJaramillo, Scuseria, and Ernzerhof (2003); Arbuznikov and Kaupp (2007, 2008); Bahmann and Kaupp (2015) are also an interesting alternative approach to construct hybrid methods that are pertinent to the context of this work.
The inclusion of additional dependencies in XC–functionals has often resulted in significant improvements in their accuracy for general calculations. However, these improvements cannot be described as systematic in the same way that the accuracy of an ab initio calculation may be systematically improved by considering a larger number of excited determinants; some DFAs give excellent results for particular systems but perform very poorly otherwise, and vice versa. There also remain many problems that none of the currently available DFAs can accurately model. An important example of this, which is pertinent to this work, are strong correlation effects, commonly found in systems with near–degenerate orbitals such as the – and –block elements, but also in systems where chemical bonds are being broken or formed.
In the present work, the problem of constructing DFAs accurate for systems with and without strong correlation is examined by considering the adiabatic connection (AC)Harris (1984); Langreth and Perdew (1975) at the local level, i.e., in each point of space.Mirtschink, Seidl, and Gori-Giorgi (2012) The AC, discussed in subsection II.1, provides an exact expression for the exchange and correlation energies by considering the changes that occur as the strength of electron interaction is smoothly increased from zero. This formalism has provided the basis for the development of several DFAs,Ernzerhof (1996); Burke, Ernzerhof, and Perdew (1997); Becke (1993); Mori-Sanchez, Cohen, and Yang (2006) which attempt to interpolate the AC between the non–interacting and physical systems in order to estimate the XC–energy. An advantage of the AC formalism crucial to our construction, is that it allows the problem of strong correlation to be addressed in a more direct way, by creating interpolation models that are explicitly dependent on the strongly–interacting limit, in addition to the non–interacting limit, of the AC.
The strongly–interacting limit of the AC has recently become the subject of much interest.Seidl, Gori-Giorgi, and Savin (2007); Gori-Giorgi, Vignale, and Seidl (2009); Malet and Gori-Giorgi (2012); Mirtschink, Seidl, and Gori-Giorgi (2012, 2013); Vuckovic et al. (2015a) The properties of the AC integrand in this limit reveal highly non–local density dependence of correlation effectsSeidl (1999); Seidl, Gori-Giorgi, and Savin (2007); Malet and Gori-Giorgi (2012); Malet et al. (2013); Mendl, Malet, and Gori-Giorgi (2014a) that cannot be obtained from the standard (semi)local or hybrid functionals. Study of the strongly–interacting limit in DFT has focused on strictly-correlated electrons (SCE) functional, where the electrons have an infinite interaction strength. This limit is of particular interest from a theoretical point of view as it can be studied exactly for one–dimensional systemsColombo, De Pascale, and Di Marino (2015) and may be closely approximated in systems with spherical or cylindrical symmetry.Seidl, Gori-Giorgi, and Savin (2007); Vuckovic et al. (2015b) These studies show that in the limit of infinite interaction strength certain integrals of the density appear in the exchange-correlation functionals, revealing a mathematical structure very different from the one of the usual semi-local or orbital-dependent approximations. The nonlocal radius (NLR) functional proposed in ref 43 approximates the SCE functional with a model that retains some of the SCE nonlocality, introducing the integrals of the spherically averaged density around a reference electron. The inclusion of the NLR functional into global and local interpolations along the adiabatic connection has been very recently explored by Zhou, Bahmann and Ernzerhof.Zhou, Bahmann, and Ernzerhof (2015) In another recent work, Kong and Proynov constructed a functional combining the information from the Becke’13 modelBecke (2013) and approximating local quantities along the AC.Kong and Proynov (2015)
The aim of the present work is to start a systematic study of local interpolation models along the adiabatic connection, using at a first stage exact input ingredients, thus disentangling the errors due to the interpolation models from those due to the approximate ingredients. The local AC for several closed–shell atoms has been recently computed Irons and Teale (2015) to high accuracy between the non–interacting and physical systems using the Legendre–Fenchel formulation of DFT due to Lieb,Lieb (1983) and the Lieb maximisation method of refs49; 50; 51; 52. Local information pertaining to the strongly–interacting limit is calculated using the SCE functional, and together these quantities are used to calculate local analogues of some established global AC interpolation functionals. We also discuss how to approximate crucial local ingredients such as the initial slope of the local adiabatic curve.
In section II, relevant theoretical background is given including an overview of the AC formalism and the construction of DFAs from both global and local variants of the AC. Techniques for computing quantities along the AC are discussed, including the determination of the local AC as introduced in ref 47. In section II.3 the construction of a local model for the AC is discussed, considering the non-interacting and strong-interaction limits carefully in this context. The role of the SCE in constructing local AC interpolation models is examined. Finally the forms of some local interpolation models, taken from successful existing global models are introduced. In section IV the performance of these models is assessed for the helium and beryllium isoelectronic series and for dissociation of the H molecule, a system that typifies the failure of present DFAs to properly account for strong–correlation. Directions for future work are outlined in section V.
Ii Theoretical Background
ii.1 The Adiabatic Connection
The AC was proposed in a series of papers,Harris and Jones (1974); Langreth and Perdew (1975); Gunnarsson and Lundqvist (1976); Harris (1984) which suggested that further insight into electronic correlations in DFT may be gained by considering a system at constant electron density as the interaction strength is smoothly scaled between zero, i.e. the KS auxiliary system, and the full physical interaction strength. This scaling of the interaction strength is achieved by the introduction of a simple coupling–constant coefficient , such that the Hamiltonian for any given is written as
(1) |
where is the kinetic energy operator, is the physical electron interaction operator and is the operator representing an external potential that binds the electron density at , such that it is always equal to the density of the physically interacting system (). As the value of is smoothly increased from zero to one, the system evolves adiabatically through a family of –dependent wave functions to the physical system described by .
Given a Hamiltonian , one can define the corresponding –dependant universal density functional as
(2a) | ||||
(2b) |
where eq 2b follows from the application of the Hellmann–Feynman theorem to eq 2a. This allows the well–known AC formula to be derived, yielding the following exact expression for the XC–energy of an electronic system,
(3) |
where is the (global) AC integrand, given by
(4) |
is the ground state wavefunction of in eq 1, and the Hartree (Coulomb) energy.
The AC integrand may be characterized by several features that can be exactly defined: the expansion of in the non–interacting limit is given byGörling and Levy (1993)
(5) |
whilst its expansion in the strongly–interacting limit can be expressed asSeidl (1999); Seidl, Gori-Giorgi, and Savin (2007); Gori-Giorgi, Vignale, and Seidl (2009)
(6) | ||||
Here, the non–interacting terms and are the exchange energy and twice the second–order correlation energy given by Görling-Levy perturbation theory (GL2)Görling and Levy (1993, 1994) (see section III.1.2), respectively. Their analogues at the strongly–interacting limit, and , have been studied in Refs. 38; 33; 34 and will be discussed further in section III.2. In addition to these asymptotic limits, the behaviour of under uniform coordinate scaling is also well–defined, as discussed in ref 57.
ii.2 DFAs based on the global adiabatic connection
To construct practical DFAs one could consider modelling the integrand of eq 4 using a function that interpolates between the limits of eq 5 and eq 6. The SCE limit is of particular importance in the present work, however, one could also consider models that intercept any other known point on the adiabatic connection for . Several attempts to develop DFAs based on these ideas have been put forward in the literature, see e.g. Refs. 30; 31; 58; 59; 60; 33; 34; 61; 62. Each form makes a choice of a simple model function and the parameters on which to base the model. These parameters often include the known exact expressions for the parameters in eq 5: , , since these may be computed from the set of Kohn–Sham orbitals () and orbital energies ().
The calculation of requires the GL2 correlation energy,Görling and Levy (1993, 1994) which leads to a computational cost similar to the second-order Møller–Plesset (MP2) model used in ab initio quantum chemistry.Möller and Plesset (1934) The parameters in the SCE limit entering eq 6 are clearly also of special interest in this context, and they can be computed numerically for atomic systems and molecules with cylindrical symmetry.Seidl, Gori-Giorgi, and Savin (2007); Di Marino et al. ; Vuckovic et al. (2015b) More frequently DFAs are derived for points along the AC with , often by employing scaling relations to derive forms from existing DFAs. A similar strategy can also be used to approximate by a DFA, see for example ref 60.
In tandem with choosing a set of exact or approximate values to parameterize a model for the AC one must also choose an appropriate model function for the AC integrand. A number of these have been suggested and many have been benchmarked in practical applications. One of the simplest (and most often used) is that of a [1/1] Padé, as suggested by Ernzerhof.Ernzerhof (1996) A range of forms were suggested by Cohen et al. and tested using approximate parameterizations,Cohen, Mori-Sánchez, and Yang (2007) leading to the MCY1 functional in which a [1/1] Padé model is employed. Peach et al.Peach, Teale, and Tozer (2007); Peach et al. (2008) attempted to disentangle approximations in the choice of parameters from those in the choice of model AC function by utilizing nearly exact KS orbitals and orbital energies derived from full configuration interaction data to calculate and and the corresponding interacting wave functions to evaluate via eq 4. Our present study follows a similar philosophy, but applied to local, rather than global, interpolations.
Seidl and co workersSeidl, Perdew, and Levy (1999); Seidl, Perdew, and Kurth (2000) were the first to make use of the strong-interaction limit (although approximated at a semilocal level, using the so-called point-charge-plus-continuum, or PC, functional) in constructing a global AC model, known as the interaction strength interpolation (ISI) functional. The revISI modelGori-Giorgi, Vignale, and Seidl (2009) and the models of Liu and BurkeLiu and Burke (2009) were later constructed to take account of the dependence of the second term of eq 6, which was not correctly described by the ISI approach. Teale, Coriani and Helgaker also proposed forms for the AC integrand based on the structure of traditional ab initio methodologiesTeale, Coriani, and Helgaker (2010c) and parameterized these forms to intercept values of the AC at any .
The majority of these models for the global AC suffer from the fact they are not size consistent in practice. This deficiency arises from a non-linear dependence on the parameters , , and a chosen approximation to . When these global parameters enter in a non-linear fashion (often as ratios) then size consistency is difficult to achieve. One route forward is to construct local AC models, which can replace these global parameters with local values defined at each point in space and may be more amenable to the construction of models that recover size-consistency (at least in the usual density-functional senseGori-Giorgi and Savin (2008); Savin (2009)).
ii.3 Constructing a local adiabatic connection
The AC expression for the XC–energy of a system as given in eq 3 describes a global quantity, integrated over the coupling–constant . However, it may equally be written as the spatial integral of a local quantity analogous to the local value of an XC–functional. To this effect, eq 3 may be re–written as
(7) |
where is the energy density at a given . It is well known that the energy density cannot be uniquely defined;Burke, Cruz, and Lam (1998a); Cruz, Lam, and Burke (1998); Tao et al. (2008) an arbitrary number of terms may be added to , yet an identical recovered if their spatial integral is zero. Thus any such energy densities are only defined within a particular gauge, and only energy densities defined in the same gauge may be meaningfully compared.
In the context of the present work, it is both convenient and physically meaningful to define in the gauge of the electrostatic potential of the exchange–correlation hole,
(8) |
where is the exchange–correlation hole,
(9) |
and is the pair–density obtained from the wave function
(10) | ||||
The definition of energy densities in the gauge of the XC–hole is well–established in the literature, and further discussion may be found in refs72; 73; 29. The coupling–constant averaged (–averaged) XC–energy density is defined as
(11) |
As the spatial integral of the product of this quantity and the density yields the XC–energy, the same quantity may be considered as a target to be modelled by XC–functionals,Irons and Teale (2015) although GGAs and metaGGAs often aim at energy densities within different definitions.Burke, Cruz, and Lam (1998b); Perdew, Burke, and Ernzerhof (1996); Tao et al. (2003) Given the invariance of the exchange energy to electron–interaction strength, eq 11 may be trivially resolved into separate exchange and correlation terms as
(12) | ||||
The aim of the local interpolation schemes examined in this work is to approximate and through interpolating the local AC. In principle, this approach is analogous to that of the global AC interpolation schemes previously discussed, but rather than depending on quantities pertaining to the global AC, they are instead constructed from their local equivalents, . Obviously, a local interpolation will only be meaningful if all of the local terms are defined in the same gauge. It is again both convenient and logical to define all local quantities in the gauge of eq 8, as in which highly accurate energy densities in the range have previously been calculated,Irons and Teale (2015) and additionally can be computed for small systems in the limit .Mirtschink, Seidl, and Gori-Giorgi (2012); Vuckovic et al. (2015a)
At , the energy density in the gauge of eq 8 is the exchange energy density , often denoted in the literature (also equal to the non-local Slater potentialSlater (1951)), which is the crucial ingredient of local hybrid functionals. Accurate and efficient computational schemes for this quantity have become available in the recent years.Bahmann and Kaupp (2015); Janesko, Krukau, and Scuseria (2008) In a way, local interpolation models can be viewed as local hybrids that carefully address the gauge problem.
The local equivalent of is not as simple to define, yet is an essential component of AC interpolation schemes as it provides a measure of the departure from exchange-only behaviour, in other words provides the information from which the correlation energy is approximated. Whilst many global models use GL2 theory for this purpose, its dependence on global quantities makes it unclear how it could be applied to a local interpolation scheme. This is discussed in detail in section III.1.2.
ii.4 The Lieb Maximization
In order to assess the quality of our local interpolation functionals, it is necessary to have accurate data of energy–densities, defined in the gauge of the XC–hole. These may be acquired by the method of Lieb maximisation, described in Refs. 50; 51; 52.
The Lieb maximisation is an optimisation algorithm developed using the convex–conjugate functional defined by Lieb in ref 48 as the Legendre–Fenchel transform to the energy,
(13a) | ||||
(13b) |
in which the density and potential are conjugate variables, belonging to the dual vector spaces
(14) |
and is the energy yielded by a given electronic structure calculation at potential . This convex–conjugate formulation follows from the concavity of variationally-determined energy in , from which Lieb showed that must be convex in . Furthermore the conjugate functional to a nonconcave energy, such as that which may result from a non–variational calculation, remains well–defined as it is necessarily convex. A subsequent Legendre–Fenchel transform of yields the concave envelope (least concave upper bound) to , hence unique solutions to can always be obtained.
In the Lieb maximisation, the optimised density is obtained by maximising with respect to variations in the potential , rather than by minimising with respect to as is usually the case. Therefore at convergence, the potential in eq 13b is that which yields . In the present work, Lieb maximisations have been carried–out at a number of points along the AC in the range , hence the density is constrained such that , resulting in a –dependent optimizing potential.
In order to effectively optimize with respect to the potential, we parameterize it by using the method of Wu and Yang (WY) Wu and Yang (2003) as
(15) |
where is the external potential due to nuclei, is a reference potential chosen to ensure that has the correct asymptotic behaviour, and are a set of Gaussian functions with coefficients . In all calculations in this work we choose the potential expansion basis set to be identical to the primary orbital basis set. The reference potential used in this work is the Fermi–Amaldi potentialFermi and Amaldi (1934) and is optimized with respect to the coefficients of the potential basis . Additionally, convergence is accelerated through the use of the Newton method described in Refs. 50; 51; 52. The relaxed–Lagrangian formulation of Helgaker and JørgensenHelgaker and Jørgensen (1989) is used to obtain relaxed densities for non–variational wavefunctions, which serve as input to the Lieb functional and are used in the determination of the derivatives required for its optimization.
In this work, Lieb maximisation calculations are performed using the implementation of Refs. 50; 51; 52 in a development version of the Dalton quantum chemistry software package,Aidas et al. (2014) in which is computed by using coupled–cluster singles and doubles (CCSD)Purvis and Bartlett (1982) and full configuration–interaction (FCI) wavefunctions. At convergence, where the optimising potential is such that , the relaxed –interacting one– and two–particle density matrices are computed, with which the -dependent XC energy densities may be obtained as
(16) |
Iii Modelling the Local AC
iii.1 The local slope in the non-interacting limit
As described in section II.3, the initial slope of the AC is an important part of many global AC models, in which it may be calculated directly by GL2 perturbation theory, however there is no analogous expression that yields the local equivalent and we give such an expression in section III.1.2. Here, the local initial slope of the XC energy density that is given in eq 8 is defined as
(17) |
and is related to the global slope, hence the GL2 correlation energy, by
(18) |
iii.1.1 Numerical calculation of the local slope
In this study is numerically approximated by the method of finite difference, with a series of for .
The resulting local slopes in the molecule with bond length of a.u. and a.u. are plotted along the H–H bond in Figure 1, along with the densities from which they are calculated, at the FCI level of theory and in the uncontracted aug-cc-pCVTZ basis set.Dunning (1989) In both cases, the local slope is greatest in magnitude at the nuclei, as has been seen previously in atoms.Irons and Teale (2015) It can be seen that the magnitude of the local slope is significantly larger in the stretched molecule, mirroring observations previously made of the global AC in the dissociating hydrogen molecule.Teale, Coriani, and Helgaker (2010a)
The local slopes in the He and Be isoelectronic series are plotted in Figures 2 and 3 respectively. It is clear that, with increasing nuclear charge, the charge densities in both series become increaslingly contracted. The –axis in both plots has been scaled by nuclear charge, highlighting a contrast in their behaviour with respect to the uniform scaling condition,
(19) |
which holds for non-degenerate KS systems.Perdew et al. (2008) In Figure 2, it can be seen that the slope of the AC for the He series becomes less negative with increasing , tending to an asymptotic value as , consistent with the scaling relation of eq 19. However, the slope of the AC in the Be isoelectronic series becomes more negative with increasing , indicating that the scaling relation is not satisfied by this series.Colonna and Savin (1999)
iii.1.2 A functional approximation for the local slope
Whilst it is useful to numerically approximate the local slope for the purposes of evaluating local interpolation schemes, such functionals would only be viable for mainstream use in DFT calculations if they can be described by simple functional forms.
In global models, the initial slope can be calculated directly from the occupied and virtual KS orbitals according to GL2 theory,
(20) | ||||
where the indices and pertain to occupied and virtual KS orbitals respectively, is the local KS potential and the non–local Hartree–Fock (HF) exchange potential. The first term in eq 20 is analogous to the correlation energy given by MP2 theory, in which and are canonical HF orbitals and eigenvalues rather than KS ones. The second term accounts for the difference between the KS and HF exchange potentials and has a form similar to a singles term in many–body perturbation theory. Previous studies of GL2 theory have found that the second term, although non–negligible, is small in magnitude relative to the MP2–like term evaluated with the KS orbitals. Toulouse et al. (2011) Therefore, it follows that an approximation to the GL2 correlation energy may be obtained by evaluating the MP2 correlation energyGrimme (2006) with the KS orbitals and eigenvalues, .
Given that an approximation to the global AC slope may be obtained from an MP2–like calculation, it follows that an approximation to the local AC slope may be obtained by deriving a local form of this expression. Whilst MP2 theory treats perturbations of the wavefunction, the analysis may be extended to energy densities in the gauge of the XC–hole by means of eq 10, as the substitution of eq 16 into eq 17 yields the following,
(21) |
where is the derivative of the pair–density at ,
(22) |
Notice that eq 21 ensures that is in the gauge of the electrostatic potential of the xc hole. Given a non–interacting ground–state wavefunction , the perturbed wavefunction for can be appproximated by the series expansion
(23) |
If one assumes that is non–degenerate and has the form of a single Slater determinant, the first–order correction to the wavefunction is given by
(24) |
Restricting the space of to doubly–excited determinants reduces this expression to
(25) |
In MP2 theory, contributions to the correlation energy from singly–excited determinants are necessarily zero due to Brillouin’s theorem. However, this is not strictly true in GL2 theory as the singles term in eq 20 makes a small, but non–zero, contribution to the GL2 correlation energy.Görling and Levy (1993) As such, considering only double–excitations in the model for the local slope can only yield approximations to the local slope; spatial integration of this quantity will not return the exact GL2 correlation energy.
Application of the Slater-Condon rules to eq 25 allows it to be re–written as
(26) |
where the coefficient to may be identified as an MP2 doubles–amplitude ,
(27) |
To obtain , it is necessary to take the derivative of the pair–density corresponding to the perturbed wavefunction, . Substituting this into eq 10 and rearranging the resulting expressions yields the following,
(28) | ||||
where we assume that and are real and is the pair–density operator. Substituting eq 26 into this expression gives,
(29) |
which may then be resolved into the following orbital–explicit expression,
(30) | ||||
where is the Kronecker delta over two spin functions: . Substituting eq 30 into eq 21 finally results in an expression for the local slope,
(31) |
where is the antisymmetrized orbital potential,
(32) | ||||
Multiplying the right–hand side of eq 31 by the density and integrating over all space, we recover twice the MP2 like expression. This is not an exact expression for the local slope, as the second term of eq 20 is not accounted for. However, the omitted term is generally small relative to the MP2–like term, and vanishes entirely for two–electron systems, hence the expression for the local slope in eq 31 should, in principle, be a fair approximation of the exact local slope.
In future work we will implement and test eq 31 against the numerical results in section III.1.1. The doubles amplitudes are readily obtainable from standard quantum chemical packages, the potential can also be readily calculated; however, it would likely be computationally expensive to evaluate on a numerical grid. To reduce this cost a range of techniques, commonly used to accelerate the calculation of integrals in linear–scaling software packages, may be employed.Werner, Manby, and Knowles (2003); Reine et al. (2008); Domínguez-Soria et al. (2009); Bahmann and Kaupp (2015)
We note that the behaviour of the local slopes presented in Figures 1, 2 and 3 may be rationalized in a similar manner to that commonly discussed for global models in terms of eq 20. This is because of the key role of the doubles amplitude in eq 31. The doubles amplitude has a dependence on the orbitals and orbital energies that is similar to that of the GL2 energy in eq 20. We see in Figure 1 that the local slope of the hydrogen molecule displays the minima at the nuclei. Equation 31, which is exact for two–electron systems, can be used to rationalize this observation. For closed shell two–electron systems with only one virtual orbital, eq 31 is simplified as follows:
(33) |
Even if we used a minimal orbital basis for the evaluation of the expression given in eq 33 for the hydrogen molecule, we would see that the local slope is most negative at the two nuclei, for any bond length. Whilst this effect is captured with the minimal basis, the same minimal basis model would incorrectly describe the slope at at the bond midpoint. For example, in the top panel of Figure 1 we see that is less than at the bond midpoint of H at . Within the minimal basis, the local slope would be exactly for any , as the antibonding orbital which enters eq 33 has a node at the bond midpoint.
We also see in Figure 2 that the correlation energy density for the He isoelectronic series scales quickly towards an asymptotic constant as increases. Furthermore, the local slope decays smoothly with distance from the nucleus, reflecting the behaviour of . The behaviour for the Be isoelectronic series in Figure 3 is more complex. The KS HOMO-LUMO gap is known to increaseSavin, Colonna, and Pollet (2003) as increases from to , from which one would expect the correlation energy to become less negative according to the behaviour of . In the core region this behaviour holds, however in the valence region the trend is opposite, with the correlation energy density becoming more negative with increasing . This suggests that the numerator of and the spatial dependence of due to the form of the KS orbitals are dominant in this region, provided that eq 31 is sufficiently accurate for the Be isoelectronic series.
iii.2 The SCE model and the strong interaction limit
In recent years, the exact strong-coupling limit of the AC has been intensively studied.Seidl, Gori-Giorgi, and Savin (2007); Gori-Giorgi, Vignale, and Seidl (2009); Malet and Gori-Giorgi (2012); Mirtschink, Seidl, and Gori-Giorgi (2012, 2013); Vuckovic et al. (2015a) This limit reveals a new structure for the XC functional: instead of the traditional ingredients of DFAs (local density, density gradients, KS kinetic energy density, occupied and unoccupied KS orbitals) it is observed that certain integrals of the density appear in this limit, encoding highly non-local information.Seidl (1999); Seidl, Gori-Giorgi, and Savin (2007); Malet and Gori-Giorgi (2012); Malet et al. (2013); Mendl, Malet, and Gori-Giorgi (2014a)
Tests on model physical and chemical systems (electrons confined in low-dimensional geometries and low-density, ultracold dipolar systems, simple stretched bonds and anions) have shownMalet and Gori-Giorgi (2012); Malet et al. (2013); Mendl, Malet, and Gori-Giorgi (2014a); Malet et al. (2014); Chen, Friesecke, and Mendl (2014); Vuckovic et al. (2015a); Malet et al. (2015) that taking into account this exact behaviour can pave the way for the solution of the strong correlation problem in DFT. However, the exact information encoded in the infinite coupling limit, described by the SCE functional, does not come for free: the SCE problem is ultra non-local, and, although sparse in principle, its non-linearity makes its exact evaluation for general three-dimensional geometry a complex task. A possible route to find suitable algorithms relies on the fact that constructing the exact SCE functional for a given density is equivalent to solving an optimal transport (or mass transportation theory) problem with a cost function given by the Coulomb interaction.Buttazzo, De Pascale, and Gori-Giorgi (2012); Cotar, Friesecke, and Klüppelberg (2013) This equivalence has triggered interest from the applied mathematics community working on optimal transport problems, which has led to the suggestion of several algorithms,Mendl and Lin (2013); Chen, Friesecke, and Mendl (2014); Benamou et al. (2015); Benamou, Carlier, and Nenna (2015) together with very interesting exact results.Colombo, De Pascale, and Di Marino (2015); Di Marino, Gerolin, and Nenna (2015); Colombo and Stra (2015)
So far, the SCE solution is known exactly for one-dimensional systems.Colombo, De Pascale, and Di Marino (2015) For spherically symmetric systems, a conjectured solutionSeidl, Gori-Giorgi, and Savin (2007) that is very close to the exact oneDi Marino et al. (and it is in many cases, but not always,Colombo and Stra (2015) exact) has been proposed and used to address interesting physical problems.Mendl, Malet, and Gori-Giorgi (2014b); Malet et al. (2015) Using algorithms and ideas from optimal transport, the SCE problem for the hydrogen molecule along the dissociation curve has just recently been solved and both the globalChen, Friesecke, and Mendl (2014); Vuckovic et al. (2015a) and localVuckovic et al. (2015a) SCE quantities have been computed. A more practical way to proceed is to build approximations for the SCE functional inspired by its exact form, as it was done in the construction of the already mentioned NLR functional.Wagner and Gori-Giorgi (2014); Zhou, Bahmann, and Ernzerhof (2015)
The SCE system complements the KS system.Seidl (1999); Seidl, Gori-Giorgi, and Savin (2007); Gori-Giorgi, Vignale, and Seidl (2009) It corresponds to the wave function that minimizes the Hamiltonian of eq 1 when . One can argue that the SCE system is a better starting point than the Kohn-Sham system for the description of very strongly correlated systems.Malet et al. (2013); Mendl, Malet, and Gori-Giorgi (2014a); Chen, Friesecke, and Mendl (2014); Vuckovic et al. (2015a)
The SCE functional is defined asSeidl, Gori-Giorgi, and Savin (2007); Malet and Gori-Giorgi (2012); Mirtschink, Seidl, and Gori-Giorgi (2012):
(34) |
The XC part can be easily extracted from , as . The KS SCE approximation, proposed in ref 35, uses the SCE functional to approximate the Hartree and exchange-correlation energy, and it is equivalent to setting for all . It has been shown that KS SCE yields good energies for systems where correlation plays a dominant role, like electrons confined in low-density nanodevices or extremely stretched bonds. Mendl, Malet, and Gori-Giorgi (2014a); Malet et al. (2013); Vuckovic et al. (2015a); Chen, Friesecke, and Mendl (2014); Malet et al. (2013) On the other hand, KS SCE treats moderately and weakly correlated systems very poorly, giving energies that are unacceptably too low.Malet et al. (2014); Vuckovic et al. (2015a) A less drastic approximation is to construct a model, in such a way that its limit is given by the exact or approximate value of , as done in the pioneering work of Seidl et. al.Seidl, Perdew, and Levy (1999); Seidl, Perdew, and Kurth (2000) Analogously, one can also model , imposing that its limit is given by . This latter approach is the main object of the following sections.
In the SCE limit, the electrons are infinitely or perfectly correlated and their positions are given by an infinite superposition of classical configurations. The basic idea is that the electronic positions are all determined by a collective variable , a feature that is captured by the so-called co-motion functions .Seidl, Gori-Giorgi, and Savin (2007) If a reference electron is at , then the position of all the other electrons in the system will be given by .Seidl, Gori-Giorgi, and Savin (2007) Since the electrons are perfectly correlated, the probability of finding the reference electron at has to be the same as the probability of finding the th electron at . Therefore, the co-motion functions have to satisfy the following differential equation:Seidl, Gori-Giorgi, and Savin (2007)
(35) |
For more details on the co-motion functions, including their group properties, see Refs. 29; 91; 99 and 33.
In terms of the co-motion functions, the SCE functional is given byMirtschink, Seidl, and Gori-Giorgi (2012)
(36) |
Despite the high nonlocality of the SCE functional, evident from eq 35, we can easily compute its functional derivative from the following expressionMalet and Gori-Giorgi (2012); Buttazzo, De Pascale, and Gori-Giorgi (2012)
(37) |
Equation 36 suggests the following energy density in the SCE limit:
(38) |
where is the Hartree potential. This expression is indeed in the gauge of the XC hole potential of eq 8, as proven in ref 29. Being derived from a wavefunction, the energy density decays like , similar to the physical () and the exchange () energy densities of eq 16. Its functional derivative, eq (37), has also the correct asymptotic behaviour.
To solve the SCE problem for spherically symmetric systems (the He and Be isoelectronic series considered in this paper) we have used the approach presented in ref 33, which is exact if . For atomic densities with it could be either a very good approximation for the true minimum of eq 34, or again, the exact result.Colombo and Stra (2015); Di Marino et al. For the H molecule we have used the results of ref 37, where the SCE energy density has been computed by obtaining the co-motion function from the dual Kantorovich formulationButtazzo, De Pascale, and Gori-Giorgi (2012); De Pascale (2015) of the SCE problem.
iii.3 Local interpolation models
The local interpolation models tested in this work are largely simple translations of the well–established global interpolation models into a local form. This was done for the model of Seidl, Perdew and Levy (SPL),Seidl, Perdew, and Levy (1999) the “simplified” model of Liu and Burke,Liu and Burke (2009) which will be referred here as the LB model and the Padé model.Ernzerhof (1996); Mori-Sánchez, Cohen, and Yang (2006) Each of the energy densities resulting from the three mentioned models is constructed from three local parameters, , and , which are defined in the gauge of the XC–hole. The functional forms of these three models are summarized in Table 1.
In addition to these, we constructed a local form of the two–legged representationBurke, Ernzerhof, and Perdew (1997) which, given some value of , takes the form
(39a) | ||||
(39b) |
Whenever we used the two–legged representation to model the local AC in this work, we did it by incorporating the interpolated of the LB model: . By doing the local interpolation this way, we use the following three input quantities: , and and circumvent the direct utilization of the full interacting energy density, . In each of these four models, integration of with respect to coupling–constant gives the –averaged energy density which, if spatially integrated according to eq 7, yields the XC–energy .
An important observation in the translation of global to local models is that, whilst the following global inequalities are always satisfied,
(40) |
their local counterparts do not necessarily satisfy these same inequalities. It has previously been observed for the Hooke’s atom series that, in the tail regions of the density, can be less negative than .Mirtschink, Seidl, and Gori-Giorgi (2012) In this work, the crossing of with , and has only been observed in the tail regions of the density and is thought to be an artefact of the numerical instability that occurs where the density is very small.
Iv Results
iv.1 Helium isoelectronic series
Although the helium isoelectronic series is a set of only two–electron systems, it is a useful series to consider in evaluating the local interpolation models as most standard DFAs incorrectly characterize the hydride ion (H), failing to predict its existence as a bound electronic system.Mirtschink et al. (2014); Peach et al. (2008) Here, local interpolation models are constructed from energy densities acquired by the Lieb maximisation at the FCI level, as described in section II.4, in the range and at by evaluating the SCE functional on the density, also at the FCI level of theory.
Z | FCI | local SPL | global SPL | local LB | Padé | local 2–leg |
---|---|---|---|---|---|---|
1 | -0.0409 | -0.0367 | -0.0368 | -0.0398 | -0.0401 | -0.0477 |
2 | -0.0400 | -0.0378 | -0.0380 | -0.0394 | -0.0399 | -0.0435 |
3 | -0.0410 | -0.0393 | -0.0395 | -0.0404 | -0.0409 | -0.0431 |
4 | -0.0416 | -0.0402 | -0.0404 | -0.0411 | -0.0415 | -0.0433 |
5 | -0.0418 | -0.0408 | -0.0409 | -0.0415 | -0.0418 | -0.0433 |
6 | -0.0419 | -0.0410 | -0.0411 | -0.0416 | -0.0418 | -0.0431 |
7 | -0.0414 | -0.0407 | -0.0408 | -0.0412 | -0.0414 | -0.0423 |
8 | -0.0412 | -0.0405 | -0.0406 | -0.0410 | -0.0412 | -0.0420 |
9 | -0.0411 | -0.0405 | -0.0406 | -0.0409 | -0.0411 | -0.0419 |
10 | -0.0411 | -0.0405 | -0.0407 | -0.0408 | -0.0411 | -0.0418 |
In Table 2, the correlation energies given by local forms of the SPL, LB, two-legged representation (the column labelled “2–leg”) and Padé models (the latter parameterized using the accurate values for , in order to compare with models that, instead, use the information) are given, along with that given by the global SPL model and the FCI correlation energy for comparison. This data shows that the local interpolation correlation energies are in close agreement with the FCI reference values; the mean absolute errors (MAE) of the local interpolation models are mH, mH, mH and mH, for the two-legged representation, SPL, LB and Padé models, respectively.
As would be expected, the local Padé is the most accurate of the models, given that it is derived from the full interacting energy density. This data further suggests that the local LB model is marginally superior to the local SPL and the two-legged representation. However, comparing the global and local models shows a slightly lower error for the global model; the local SPL model has an MAE of mH, compared to mH for the global model.
iv.2 Beryllium isoelectronic series
The changes in correlation energy across the beryllium isoelectronic series are somewhat more complicated than those in the helium isoelectronic series, and its explanation involves the interplay of several effects. With increasing nuclear charge, the density becomes increasingly contracted, suggesting that the correlation energy should approach the high–density limit for very large . However, this is accompanied by a changing KS HOMO–LUMO gap, here the energy difference between and orbitals, which increases from before decreasing with higher values.Savin, Colonna, and Pollet (2003)
Z | CCSD | local SPL | global SPL | local LB | Padé | local 2–leg |
---|---|---|---|---|---|---|
4 | -0.0920 | -0.0876 | -0.1049 | -0.0925 | -0.0911 | -0.1046 |
5 | -0.1089 | -0.1041 | -0.1250 | -0.1100 | -0.1076 | -0.1246 |
6 | -0.1244 | -0.1202 | -0.1455 | -0.1271 | -0.1229 | -0.1444 |
7 | -0.1389 | -0.1363 | -0.1668 | -0.1443 | -0.1373 | -0.1645 |
8 | -0.1534 | -0.1532 | -0.1898 | -0.1626 | -0.1517 | -0.1859 |
9 | -0.1683 | -0.1717 | -0.2157 | -0.1826 | -0.1666 | -0.2098 |
10 | -0.1833 | -0.1920 | -0.2447 | -0.2046 | -0.1817 | -0.2361 |
Table 3 shows the reference and interpolated results for the Be isoelectronic series, with in the range . The wave function for values between and has been computed in the same way as for the He isoelectronic series, however at the CCSD level of theory rather than FCI. As for the helium series, the local Padé that uses is the most accurate of the local interpolation models. However, in contrast to the findings for He isoelectronic series, the local interpolation models are much more accurate than the global models. For example, in the case of F the global SPL model has an MAE of mH, whereas the error for the local SPL model is mH. The local two-legged representation interpolation underestimates the correlation energies of the elements of the given series. We discuss in more details this model in the next subsection.
Figure 5 shows the –averaged correlation energy densities for the beryllium atom. The shape of reflects the shell structure of the Be atom.Colonna and Savin (1999); Irons and Teale (2015) The local SPL and LB interpolation models appear to qualitatively capture the shell structure of , however in some regions it overestimates the reference value whilst in other regions the converse is the case. The error cancellation that results from this is the most likely explanation for the superior accuracy of the local models in comparison to the global models.
iv.3 Hydrogen molecule
Despite the development of DFT into the most widely–applied electronic structure method, and the wealth of XC–DFAs that have been developed, there are some systems for which no combination of DFAs provide an accurate description. A well–known example of such a system is the dissociating molecule.Fuchs et al. (2005); Vuckovic et al. (2015b) Standard DFAs become increasingly inaccurate with greater H–H bond length, reflecting a fundamental flaw of DFAs in their inability to properly treat strong correlation.
It has been seen previouslyChen, Friesecke, and Mendl (2014); Vuckovic et al. (2015a) that KS SCE correctly predicts the dissociation of in a spin–restricted formalism, however at equilibrium geometry the energies it predicts are extremely low and the bond lengths predicted are overly short. The overall accuracy of KS-SCE for dissociation can be substantially improved by the addition of nonlocal corrections.Vuckovic et al. (2015a)
Figure 6 shows the dissociation curves for given by the local interpolation models, along with those given by HF, FCI and the PBE functionalPerdew, Burke, and Ernzerhof (1996) for comparison. The computational details are the same as those of the He isoelectronic series, and the PBE, HF and FCI curves have been obtained from the Dalton quantum chemistry packageAidas et al. (2014) all within the uncontracted aug-cc-pCVTZ basis set.Dunning (1989) The SCE energy density has been computed by using the dual Kantorovich method.Vuckovic et al. (2015a)
It can be seen in Figure 6 that all of the interpolation models correctly predict the dissociation of , which follows from their inclusion of . In global AC models, at infinite separation the initial slope diverges as a result of the vanishing HOMO–LUMO gap, and the SPL and LB models reduce to , yielding the exact energies. However, the dissociation curves produced by the local models approach the FCI curve slowly, resulting in an unphysical ‘bump’–like feature. This is a well–known failing of DFT, having been observed with other functionals, such as the random–phase approximation Fuchs et al. (2005) and even the global Padé model with .Peach, Teale, and Tozer (2007) It can be seen in Figure 6 that this is not remedied by the local interpolation approach, as the curve obtained by the local Padé also exhibits this unphysical bump, as does that given by the local SPL model and, to a lesser extent, the local LB model.
To analyse why the intermediate region is less accurately described by the local interpolation methods than the equilibrium and stretched region, we show in Figure 8 the correlation component of the local AC at one of the nuclei of the hydrogen molecule at different bond lengths: a.u. (at equilibrium), a.u. (the intermediate region) and a.u. (stretched bond). The structure of the three local AC curves at one of the nuclei is very similar to the structure of the corresponding global AC curves.Teale, Coriani, and Helgaker (2009) From the given figure we see that at equilibrium the local AC is almost linear, so we can expect that even a single line segment approximation to the local AC: would properly capture the shape of the given local AC curve. The local AC curve at the nuclei of the stretched H exhibits the characteristic ‘L-shape’, which was also observed in the case of the corresponding global AC curve.Teale, Coriani, and Helgaker (2009) We would expect that the two-legged representation would capture the given local AC very well, but even a single line segment approximation: , this time coming from the strong coupling limit, would be highly accurate for the stretched H.Vuckovic et al. (2015a) In contrast to the local AC curves of the stretched and H at equilibrium, the curvature of the local AC curve at the intermediate bond length is highly pronounced. The shapes of the local AC curves at the nuclei mirror the difference in correlation regimes present in the hydrogen molecule at different bond lengths. While in the H at equilibrium and at very stretched bond length, correlation is almost purely dynamical and almost purely static, respectively, in the intermediate dissociation region there is a subtle interplay between the dynamical and static correlation.
In the intermediate region of the dissociation curve, where the unphysical bump is present, the local two-legged representation model is more accurate than the local Padé which we always use here with . This may be understood by comparing the exact local AC data with the interpolated quantities. The top panel of Figure 8 shows the difference between and that of each of the local interpolation models, along the H–H bond at the a.u. geometry, as a function of the distance from the bond midpoint . This difference , is multiplied by the density to represent an energy per volume element. It shows that the local SPL energy density is the one that most overestimates the . The error is smaller for the LB model and even more so for the local Padé model. The error is smallest in the two–legged model, obtained using the of the local LB. Furthermore, there is the error cancellation in the two–legged model, as there are regions where the of this model underestimates .
It can also be seen that the curves shown in the top panel of Figure 8 have a maximum at the nucleus (). Focusing on this region, it appears that the FCI curve meets that of the Padé at , and that the two–legged representation curve meets that of the LB model also at . This follows from the construction of the Padé and two–legged curves from and respectively. All curves, except for that of the two-legged model, lie above the FCI curve. In the case of the two–legged interpolation model, the first line segment is below the FCI curve, as a result of eq 39a and the convexity of the given local AC curve. The second line segment that starts at is given by , and lies above the FCI curve. The resulting error cancellation makes it clear why the two–legged representation appears more accurate than the other models.
V Conclusion and Perspectives
In this work we have studied local interpolations along the adiabatic connection for the He and Be isoelectronic series and the hydrogen molecule, by using accurate input local quantities computed in the gauge of the electrostatic potential of the XC hole, and comparing the results with nearly exact energy densities defined in the same way. In order to obtain approximations to the local AC over the physical regime (), we constructed interpolation models between the weak and strong coupling limits of DFT. The weak coupling energy densities were obtained using the Lieb variation principle, whilst the strong coupling limit energy densities were obtained using the strictly-correlated electrons (SCE) approach. The inclusion of the SCE information in density functional approximations helps to ensure their ability to capture the strong correlation effects.
Unlike previous attempts in this direction that used global (integrated over all space) input quantities to model the AC, the local approach is more amenable to the construction of approximations that do not violate size consistency, at least in the usual DFT sense.Gori-Giorgi and Savin (2008); Savin (2009) Since the aim here is to work in a restricted formalism, avoiding to mimic strong correlation with symmetry breaking, some care must be taken when discussing size consistency. In fact, strictly speaking, in a restricted framework the energy densities of the second-order perturbation theory and exact exchange are not intensive quantities in the presence of near degeneracy,Gori-Giorgi and Savin (2008); Savin (2009) which is the main challenge of capturing strong correlation within DFT.Savin (1996); Cohen, Mori-Sanchez, and Yang (2008); Mori-Sanchez, Cohen, and Yang (2009)
In future work we will test different approximations for the SCE energy densities and the local slope. The development of algorithms for solution of the SCE problem is a very active research field. In spite of the recent improvements, we still lack an algorithm that will solve the SCE problem for general 3D molecular geometries at low computational cost. However, a good candidate to approximate the SCE energy density in the gauge of the XC hole potential is the nonlocal radius functional (NLR),Wagner and Gori-Giorgi (2014) which has been already implemented and used in ref 44.
In addition to numerically exploring the local AC we have also reported the local weak-coupling slope of the adiabatic connection and derived an approximate expression for it in terms of occupied and unoccupied orbitals. This quantity is very important to signal the amount of correlation at each point of space. In our future work we will implement this expression and test it against the results reported here.
Acknowledgements
We are very grateful to Derk Kooi and André Mirtschink for a critical reading of the manuscript and insightful suggestions to improve it. We acknowledge financial support from the European Research Council under H2020/ERC Consolidator Grant corr-DFT (Grant No. 648932) and the Netherlands Organization for Scientific Research (NWO) through an ECHO grant (717.013.004). A. M. T. is grateful for support from the Royal Society University Research Fellowship scheme. A. M. T and T. J. P. I. are grateful for support from the Engineering and Physical Sciences Research Council (EPSRC), (Grant No. EP/M029131/1). We are grateful for access to the University of Nottingham High Performance Computing Facility.
References
- Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A 1133 (1965).
- Lee, Yang, and Parr (1988) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
- Becke (1988) A. D. Becke, Phys. Rev. A 38, 3098 (1988).
- Perdew et al. (1992) 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).
- Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Zhao and Truhlar (2008) Y. Zhao and D. G. Truhlar, J. Chem. Phys. 128, 184109 (2008).
- Tao et al. (2003) J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys. Rev. Lett. 91, 146401 (2003).
- Peverati and Truhlar (2014) R. Peverati and D. G. Truhlar, Phil. Trans. R. Soc. A 372, 20120476 (2014).
- Sun, Ruzsinszky, and Perdew (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
- Becke (1993) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
- Perdew, Ernzerhof, and Burke (1996) J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 105, 9982 (1996).
- Zhao, Schultz, and Truhlar (2006) Y. Zhao, N. E. Schultz, and D. G. Truhlar, J. Chem. Theory Comput. 2, 364 (2006).
- Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
- Vydrov and Scuseria (2004) O. A. Vydrov and G. E. Scuseria, J. Chem. Phys. 121, 8187 (2004).
- Vydrov and Scuseria (2005) O. A. Vydrov and G. E. Scuseria, J. Chem. Phys. 122, 184107 (2005).
- Kümmel and Kronik (2008) S. Kümmel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008).
- Grimme (2006) S. Grimme, J. Chem. Phys. 124, 034108 (2006).
- Goerigk and Grimme (2010) L. Goerigk and S. Grimme, J. Chem. Theory Comput. 7, 291 (2010).
- Sharkas, Toulouse, and Savin (2011) K. Sharkas, J. Toulouse, and A. Savin, J. Chem. Phys. 134, 064113 (2011).
- Furche (2001) F. Furche, Phys. Rev. B 64, 195120 (2001).
- Hesselmann and Görling (2010) A. Hesselmann and A. Görling, Mol. Phys. 108, 359 (2010).
- Ángyán et al. (2011) J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, J. Chem. Theory Comput. 7, 3116 (2011).
- Jaramillo, Scuseria, and Ernzerhof (2003) J. Jaramillo, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 1068 (2003).
- Arbuznikov and Kaupp (2007) A. V. Arbuznikov and M. Kaupp, Chem. Phys. Lett. 440, 160 (2007).
- Arbuznikov and Kaupp (2008) A. V. Arbuznikov and M. Kaupp, J. Chem. Phys. 128, 214107 (2008).
- Bahmann and Kaupp (2015) H. Bahmann and M. Kaupp, J. Chem. Theory Comput. 11, 1540 (2015).
- Harris (1984) J. Harris, Phys. Rev. A 29, 1648 (1984).
- Langreth and Perdew (1975) D. C. Langreth and J. P. Perdew, Solid State Commun. 17, 1425 (1975).
- Mirtschink, Seidl, and Gori-Giorgi (2012) A. Mirtschink, M. Seidl, and P. Gori-Giorgi, J. Chem. Theory Comput. 8, 3097 (2012).
- Ernzerhof (1996) M. Ernzerhof, Chem. Phys. Lett. 263, 499 (1996).
- Burke, Ernzerhof, and Perdew (1997) K. Burke, M. Ernzerhof, and J. P. Perdew, Chem. Phys. Lett. 265, 115 (1997).
- Mori-Sanchez, Cohen, and Yang (2006) P. Mori-Sanchez, A. J. Cohen, and W. T. Yang, J. Chem. Phys. 124, 091102 (2006).
- Seidl, Gori-Giorgi, and Savin (2007) M. Seidl, P. Gori-Giorgi, and A. Savin, Phys. Rev. A 75, 042511 (2007).
- Gori-Giorgi, Vignale, and Seidl (2009) P. Gori-Giorgi, G. Vignale, and M. Seidl, J. Chem. Theory Comput. 5, 743 (2009).
- Malet and Gori-Giorgi (2012) F. Malet and P. Gori-Giorgi, Phys. Rev. Lett. 109, 246402 (2012).
- Mirtschink, Seidl, and Gori-Giorgi (2013) A. Mirtschink, M. Seidl, and P. Gori-Giorgi, Phys. Rev. Lett. 111, 126402 (2013).
- Vuckovic et al. (2015a) S. Vuckovic, L. O. Wagner, A. Mirtschink, and P. Gori-Giorgi, J. Chem. Theory Comput. 11, 3153 (2015a).
- Seidl (1999) M. Seidl, Phys. Rev. A 60, 4387 (1999).
- Malet et al. (2013) F. Malet, A. Mirtschink, J. C. Cremon, S. M. Reimann, and P. Gori-Giorgi, Phys. Rev. B 87, 115146 (2013).
- Mendl, Malet, and Gori-Giorgi (2014a) C. B. Mendl, F. Malet, and P. Gori-Giorgi, Phys. Rev. B 89, 125106 (2014a).
- Colombo, De Pascale, and Di Marino (2015) M. Colombo, L. De Pascale, and S. Di Marino, Can. J. Math. 67, 350 (2015).
- Vuckovic et al. (2015b) S. Vuckovic, L. O. Wagner, A. Mirtschink, and P. Gori-Giorgi, J. Chem. Theory Comput. 11, 3153 (2015b).
- Wagner and Gori-Giorgi (2014) L. O. Wagner and P. Gori-Giorgi, Phys. Rev. A 90, 052512 (2014).
- Zhou, Bahmann, and Ernzerhof (2015) Y. Zhou, H. Bahmann, and M. Ernzerhof, J. Chem. Phys. 143, 124103 (2015).
- Becke (2013) A. D. Becke, J. Chem. Phys. 138, 074109 (2013).
- Kong and Proynov (2015) J. Kong and E. Proynov, J. Chem. Theory Comput. 12, 133 (2015).
- Irons and Teale (2015) T. J. Irons and A. M. Teale, Mol. Phys. , DOI: 10.1080/00268976.2015.1096424 (2015).
- Lieb (1983) E. H. Lieb, Int. J. Quantum. Chem. 24, 24 (1983).
- Wu and Yang (2003) Q. Wu and W. Yang, J. Chem. Phys. 118, 2498 (2003).
- Teale, Coriani, and Helgaker (2009) A. M. Teale, S. Coriani, and T. Helgaker, J. Chem. Phys. 130, 104111 (2009).
- Teale, Coriani, and Helgaker (2010a) A. M. Teale, S. Coriani, and T. Helgaker, J. Chem. Phys. 132, 164115 (2010a).
- Teale, Coriani, and Helgaker (2010b) A. M. Teale, S. Coriani, and T. Helgaker, J. Chem. Phys. 133, 164112 (2010b).
- Harris and Jones (1974) J. Harris and R. Jones, J. Phys. F 4, 1170 (1974).
- Gunnarsson and Lundqvist (1976) O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B 13, 4274 (1976).
- Görling and Levy (1993) A. Görling and M. Levy, Phys. Rev. B 47, 13105 (1993).
- Görling and Levy (1994) A. Görling and M. Levy, Phys. Rev. A 50, 196 (1994).
- Levy (1991) M. Levy, Phys. Rev. A 43, 4637 (1991).
- Seidl, Perdew, and Levy (1999) M. Seidl, J. P. Perdew, and M. Levy, Phys. Rev. A 59, 51 (1999).
- Seidl, Perdew, and Kurth (2000) M. Seidl, J. P. Perdew, and S. Kurth, Phys. Rev. Lett. 84, 5070 (2000).
- Cohen, Mori-Sánchez, and Yang (2007) A. J. Cohen, P. Mori-Sánchez, and W. Yang, J. Chem. Phys. 127, 034101 (2007).
- Liu and Burke (2009) Z. F. Liu and K. Burke, Phys. Rev. A 79, 064503 (2009).
- Teale, Coriani, and Helgaker (2010c) A. M. Teale, S. Coriani, and T. Helgaker, J. Chem. Phys. 132, 164115 (2010c).
- Möller and Plesset (1934) C. Möller and M. S. Plesset, Phys. Rev. 46, 618 (1934).
- (64) S. Di Marino, A. Gerolin, L. Nenna, M. Seidl, and P. Gori-Giorgi, “in preparation,” .
- Peach, Teale, and Tozer (2007) M. J. G. Peach, A. M. Teale, and D. J. Tozer, J. Chem. Phys. 126, 244104 (2007).
- Peach et al. (2008) M. J. G. Peach, A. M. Miller, A. M. Teale, and D. J. Tozer, J. Chem. Phys. 129, 064105 (2008).
- Gori-Giorgi and Savin (2008) P. Gori-Giorgi and A. Savin, J. Phys.: Conf. Ser. 117, 012017 (2008).
- Savin (2009) A. Savin, Chem. Phys. 356, 91 (2009).
- Burke, Cruz, and Lam (1998a) K. Burke, F. G. Cruz, and K.-C. Lam, J. Chem. Phys. 109, 8161 (1998a).
- Cruz, Lam, and Burke (1998) F. G. Cruz, K.-C. Lam, and K. Burke, J. Phys. Chem. A 102, 4911 (1998).
- Tao et al. (2008) J. Tao, V. N. Staroverov, G. E. Scuseria, and J. P. Perdew, Phys. Rev. A 77, 012509 (2008).
- Burke, Cruz, and Lam (1998b) K. Burke, F. G. Cruz, and K.-C. Lam, J. Chem. Phys. 109, 8161 (1998b).
- Armiento and Mattsson (2002) R. Armiento and A. E. Mattsson, Phys. Rev. B 66 (2002).
- Slater (1951) J. C. Slater, Phys. Rev. 81, 385 (1951).
- Janesko, Krukau, and Scuseria (2008) B. G. Janesko, A. V. Krukau, and G. E. Scuseria, J. Chem. Phys. 129, 124110 (2008).
- Fermi and Amaldi (1934) E. Fermi and E. Amaldi, R. Accad. d’Italia. Memorie 6, 119 (1934).
- Helgaker and Jørgensen (1989) T. Helgaker and P. Jørgensen, Theor. Chim. Acta 75, 111 (1989).
- Aidas et al. (2014) K. Aidas, C. Angeli, K. L. Bak, V. Bakken, R. Bast, L. Boman, O. Christiansen, R. Cimiraglia, S. Coriani, P. Dahle, E. K. Dalskov, U. Ekström, T. Enevoldsen, J. J. Eriksen, P. Ettenhuber, B. Fernández, L. Ferrighi, H. Fliegl, L. Frediani, K. Hald, A. Halkier, C. Hättig, H. Heiberg, T. Helgaker, A. C. Hennum, H. Hettema, E. Hjertenaes, S. Høst, I.-M. Høyvik, M. F. Iozzi, B. Jansík, H. J. A. Jensen, D. Jonsson, P. Jørgensen, J. Kauczor, S. Kirpekar, T. Kjaergaard, W. Klopper, S. Knecht, R. Kobayashi, H. Koch, J. Kongsted, A. Krapp, K. Kristensen, A. Ligabue, O. B. Lutnaes, J. I. Melo, K. V. Mikkelsen, R. H. Myhre, C. Neiss, C. B. Nielsen, P. Norman, J. Olsen, J. M. H. Olsen, A. Osted, M. J. Packer, F. Pawlowski, T. B. Pedersen, P. F. Provasi, S. Reine, Z. Rinkevicius, T. A. Ruden, K. Ruud, V. V. Rybkin, P. Sałek, C. C. M. Samson, A. S. de Merás, T. Saue, S. P. A. Sauer, B. Schimmelpfennig, K. Sneskov, A. H. Steindal, K. O. Sylvester-Hvid, P. R. Taylor, A. M. Teale, E. I. Tellgren, D. P. Tew, A. J. Thorvaldsen, L. Thøgersen, O. Vahtras, M. A. Watson, D. J. D. Wilson, M. Ziolkowski, and H. Ågren, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 4, 269 (2014).
- Purvis and Bartlett (1982) G. D. Purvis and R. J. Bartlett, J. Chem. Phys. 76, 1910 (1982).
- Dunning (1989) T. H. Dunning, J. Chem. Phys. 90, 1007 (1989).
- Perdew et al. (2008) J. P. Perdew, V. N. Staroverov, J. Tao, and G. E. Scuseria, Phys. Rev. A 78, 052513 (2008).
- Colonna and Savin (1999) F. Colonna and A. Savin, J. Chem. Phys. 110, 2828 (1999).
- Toulouse et al. (2011) J. Toulouse, K. Sharkas, E. Brémond, and C. Adamo, J. Chem. Phys. 135, 101102 (2011).
- Werner, Manby, and Knowles (2003) H.-J. Werner, F. R. Manby, and P. J. Knowles, J. Chem. Phys. 118, 8149 (2003).
- Reine et al. (2008) S. Reine, E. Tellgren, A. Krapp, T. Kjærgaard, T. Helgaker, B. Jansik, S. Hoøst, and P. Salek, J. Chem. Phys. 129, 104101 (2008).
- Domínguez-Soria et al. (2009) V. D. Domínguez-Soria, G. Geudtner, J. L. Morales, P. Calaminici, and A. M. Köster, J. Chem. Phys. 131, 124102 (2009).
- Savin, Colonna, and Pollet (2003) A. Savin, F. Colonna, and R. Pollet, Int. J. Quantum. Chem. 93, 166 (2003).
- Malet et al. (2014) F. Malet, A. Mirtschink, K. Giesbertz, L. Wagner, and P. Gori-Giorgi, Phys. Chem. Chem. Phys. 16, 14551 (2014).
- Chen, Friesecke, and Mendl (2014) H. Chen, G. Friesecke, and C. B. Mendl, J. Chem. Theory Comput 10, 4360 (2014).
- Malet et al. (2015) F. Malet, A. Mirtschink, C. B. Mendl, J. Bjerlin, E. O. Karabulut, S. M. Reimann, and P. Gori-Giorgi, Phys. Rev. Lett. 115, 033006 (2015).
- Buttazzo, De Pascale, and Gori-Giorgi (2012) G. Buttazzo, L. De Pascale, and P. Gori-Giorgi, Phys. Rev. A 85, 062502 (2012).
- Cotar, Friesecke, and Klüppelberg (2013) C. Cotar, G. Friesecke, and C. Klüppelberg, Comm. Pure Appl. Math. 66, 548 (2013).
- Mendl and Lin (2013) C. B. Mendl and L. Lin, Phys. Rev. B 87, 125106 (2013).
- Benamou et al. (2015) J.-D. Benamou, G. Carlier, M. Cuturi, L. Nenna, and G. Peyré, SIAM J. on Sci. Comput. 37, A1111 (2015).
- Benamou, Carlier, and Nenna (2015) J.-D. Benamou, G. Carlier, and L. Nenna, arXiv preprint arXiv:1505.01136 (2015).
- Colombo, De Pascale, and Di Marino (2015) M. Colombo, L. De Pascale, and S. Di Marino, Canad. J. Math 67, 350 (2015).
- Di Marino, Gerolin, and Nenna (2015) S. Di Marino, A. Gerolin, and L. Nenna, arXiv preprint arXiv:1506.04565 (2015).
- Colombo and Stra (2015) M. Colombo and F. Stra, arXiv preprint arXiv:1507.08522 (2015).
- Mendl, Malet, and Gori-Giorgi (2014b) C. B. Mendl, F. Malet, and P. Gori-Giorgi, Phys. Rev. B 89, 125106 (2014b).
- De Pascale (2015) L. De Pascale, arXiv preprint arXiv:1503.07063 (2015).
- Mori-Sánchez, Cohen, and Yang (2006) P. Mori-Sánchez, A. J. Cohen, and W. Yang, J. Chem. Phys. 125, 201102 (2006).
- Mirtschink et al. (2014) A. Mirtschink, C. J. Umrigar, J. D. Morgan, and P. Gori-Giorgi, J. Chem. Phys. 140, 18A532 (2014).
- Fuchs et al. (2005) M. Fuchs, Y. M. Niquet, X. Gonze, and K. Burke, J. Chem. Phys. 122, 094116 (2005).
- Savin (1996) A. Savin, in Recent Developments of Modern Density Functional Theory, edited by J. M. Seminario (Elsevier, Amsterdam, 1996) pp. 327–357.
- Cohen, Mori-Sanchez, and Yang (2008) A. J. Cohen, P. Mori-Sanchez, and W. Yang, Science 321, 792 (2008).
- Mori-Sanchez, Cohen, and Yang (2009) P. Mori-Sanchez, A. J. Cohen, and W. Yang, Phys. Rev. Lett. 102, 066403 (2009).