Accurate Ground State Energies of Solids and Molecules from Time Dependent Density Functional Theory
We demonstrate that ground state energies approaching chemical accuracy can be obtained by combining the adiabatic connection fluctuation-dissipation theorem (ACFDT) with time-dependent density functional theory (TDDFT). The key ingredient is a renormalization scheme, which eliminates the divergence of the correlation hole characteristic of any local kernel. This new class of renormalized kernels gives a significantly better description of the short-range correlations in covalent bonds compared to the random phase approximation (RPA) and yields a four fold improvement of RPA binding energies in both molecules and solids. We also consider examples of barrier heights in chemical reactions, molecular adsorption and graphene interacting with metal surfaces, which are three examples where RPA has been successful. In these cases, the renormalized kernel provides results that are of equal quality or even slightly better than RPA, with a similar computational cost.
pacs:31.15.E-, 31.15.ve, 31.15.vn, 71.15.Mb
The calculation of ground state energies with a precision on the order of kT at room temperature (around 1 kcal/mol or 43 meV) is a long standing challenge in computational chemistry and materials science. This degree of precision, sometimes referred to as ”chemical accuracy”, can in principle be achieved using quantum chemistry wave function methods Bartlett and Muslat (2007). However, the extreme scaling of such methods limits their usage to relatively small systems, and their application to periodic systems in general and metals in particular, seems highly non-trivial Booth GH (2013); Shepherd and Grüneis (2013).
Recently, the adiabatic-connection fluctuation-dissipation theorem (ACFDT) has become a popular method for obtaining electronic ground state energies form first principles. With this approach the electronic correlation energy is obtained from an approximation to the dynamic density response function , which can be obtained from Time-Dependent Density Functional Theory (TDDFT) through the Dyson equation
Here is the non-interacting Kohn-Sham response function and with the Coulomb interaction and the xc-kernel. In static Density Functional Theory (DFT), one needs a rather involved approximation for the xc-energy in order to provide a decent description of the ground state properties of a particular system and no approximation has so far been able to provide accurate results across different binding regimes. In contrast, the ACFDT and Eq. (1) gives a simple framework for calculating correlation energies in terms of the excited states of a non-interacting auxiliary system and is expected to provide a high degree of accuracy with simple approximations for the xc-kernel. In particular, if we use we obtain the Random Phase Approximation (RPA) Bohm and Pines (1951), which has been shown to give an accurate account of dispersive interactions, static correlation, weak covalent bonds, and is presently considered state of the art in ab initio electronic structure theory involving solid state systems Furche (2001); Eshuis et al. (2012); Ren et al. (2012); Olsen et al. (2011); Olsen and Thygesen (2013,); Yan et al. (2013); Harl and Kresse (2009); Harl et al. (2010); Schimka et al. (2010). Nevertheless, RPA suffers from large self-correlation errors and predicts too weak binding of solids and molecules, which have severely limited the universal applicability of the method. Furthermore, extending RPA with semi-local approximations for the xc-kernel has been shown to fail dramatically, and progress in ”beyond RPA” methods within the framework of TDDFT has so far been very limited Lein et al. (2000); Furche and Voorhis (2005); Dobson and Wang (2000); Jung et al. (2004); Gould (2012); Heßelmann and Görling (2010, 2011). A somewhat orthogonal approach for improving RPA ground state energies, is based on eliminating the RPA self-correlation energy within many-body perturbation theory and is referred to as second-order screened exchange (SOSEX). The SOSEX correlation energy vanishes by construction for any one-electron system and has been shown to improve the accuracy of covalent bonds slightly, while it completely deteriorates the good description of static correlation and barrier heights in RPA Grüneis et al. (2009); Ren et al. (2012). In addition, SOSEX scales as with system size and therefore quickly becomes much more computationally demanding than RPA, which scales as .
In this Letter, we apply a renormalization scheme that solves the divergence problem associated with semi-local adiabatic TDDFT in the ACDFT. We have implemented the xc-kernels for the two most commonly used xc-approximations (LDA and PBE) Perdew et al. (1996) within this renormalization scheme and show that the results obtained with adiabatic renormalized versions are superior to RPA in all cases. In particular, with the renormalized adiabatic PBE we obtain four-fold improvement of binding energies of molecules and solids, while it maintains the scaling of RPA.
The construction of the renormalized kernel is motivated by the observation that the Fourier transformed correlation hole in the homogeneous electron gas (HEG) is closely reproduced by the ALDA exchange kernel for wave vectors . For larger the exact hole is essentially zero. In contrast, while the ALDA hole vanishes exactly for , it attains a finite value for larger and decays to zero much too slowly. In fact, this long tail produces a divergence of the correlation hole at . For the ALDA correlation hole vanishes, since , and it is tempting to define a new renormalized kernel as (note that in this work we only consider the exchange part of the LDA and PBE kernels). Alternatively, one can generalize the LDA energy functional to include non-local exchange-correlation effects, by replacing the local density by an averaged quantity . This substitution will not alter the ground state energy or potential of the homogeneous electron gas, but it will introduce non-locality into the kernel. In fact, it seems physically reasonable that the exchange-correlation energy density should include contributions from the density in the vicinity of the xc hole, which will be the case if has a width similar to the xc hole. This is accomplished by choosing as the Fourier transform of the step function and we obtain the renormalized kernel as the adiabatic kernel derived from the generalized LDA energy functional. In the supplementary material we provide more details on this derivation.
In Refs. Olsen and Thygesen (2012, 2013) we investigated this rALDA kernel in detail and we refer to those papers for more details. The idea is readily generalized to any adiabatic semi-local kernel, AX, and we obtain the following expression for the renormalized kernel in real space:
Here is the cut-off wave vector where the Fourier transform of the correlation hole of the HEG obtained with becomes zero (for X=LDA the cut-off becomes ). To generalize the construction to inhomogeneous systems we apply the substitutions
Thus both and become functions of and . Importantly, the delta function of the original adiabatic kernel acquires a finite (density-dependent) width which ensures a finite value of the correlation hole for and entails a better description of the short-range correlation effects. We will only discuss non-self-consistent applications of the kernel. This implies that calculated energies will depend on the input density, orbitals and eigenvalues. This situation is well known from RPA calculations. We note, however, that in contrast to RPA where no obvious starting point exists, it is natural to base an rAX calculation on a DFT-X calculation. We regard this is as a fundamental merit of the method, since it eliminates the arbitrary choice of input orbitals in non-self-consistent RPA calculations.
We have previously shown that the class of kernels (Accurate Ground State Energies of Solids and Molecules from Time Dependent Density
Functional Theory) significantly improves the correlation energies of the Homogeneous electron gas. Furthermore, using the exchange part of rALDA, reduces the correlation energy of a Hydrogen atom to -0.1 eV. This is already a major improvement compared to RPA which gives eV, but the value still comprises a significant self-correlation error. However, inclusion of gradient corrections in the form of PBE exchange, reduces the correlation energy to less than one meV (less than the resolution of the implementation). In the following we will refer to this kernel as rAPBE and demonstrate the superiority of the resulting kernel over the RPA and other ”beyond RPA” methods
Fig. 1 shows the deviation from experimental values of molecular atomization energies. RPA has a well-known tendency to underbind and performs somewhat worse than PBE. We also observe a significant difference between RPA@LDA and RPA@PBE with RPA@PBE being the more accurate. The renormalized kernels rALDA and rAPBE shows a striking improvement compared to both RPA and PBE with most errors being on the order of 1-3 kcal/mol. In Fig. 2 we compare the mean absolute relative errors of different methods. It is seen that already the rALDA is better than the SOSEX and approaches the accuracy of r2PT (SOSEX corrected for single excitations Ren et al. (2011)) and the hybrid PBE0. However, the gradient-corrected rAPBE reduces the error of rALDA by a factor of two and outperforms both PBE0 and rPT2. Since RPA already performs very well for geometry optimization, we have not performed a detailed comparison of bond lengths with the renormalized kernels. We have only considered the F and N molecules where we find deviations from the experimental bond lengths of and with rAPBE (RPA), respectively.
In contrast to the case of atomization energies, RPA is highly accurate for molecular barrier heights where it performs much better than both the SOSEX methods and hybrid functionals Ren et al. (2012). Here we have just tested a single reaction barrier with the rAPBE functional namely: H+NOOH+N. For this case we obtain errors of 0.9(0.4) kcal/mol and 0.2(2.4) kcal/mol for the forward and backward reactions, respectively, using the rAPBE(RPA) functional. For comparison we obtain errors of 3.2(5.9) and 12.2(8.7) using PBE0(B3LYP). Another important property of the renormalized kernels, is the fact that they describe the dissociation of the H molecule correctly in the strict atomic limit Olsen and Thygesen (2013). This is in contrast to most MBPT extensions of RPA, where one has to sacrifice the good description of static correlation in molecules Ren et al. (2012).
While the hybrid functionals may be a good choice for a decent accuracy in molecular atomization energies, these methods fail completely for the cohesive energies of solids. The RPA also performs poorly, and the PBE functional seems to be the best choice for this problem. Nevertheless, from Figs. 3 and 4 it is clear that the rALDA functional outperforms RPA and produces an accuracy similar to PBE. Inclusion of gradient corrections through the rAPBE functional reduces the error by a factor of four and is thus much more accurate than any of the other functionals considered. We note in passing that SOSEX has been shown to produce results of similar accuracy as the rAPBE for a small test set of five semiconductors Grüneis et al. (2009), but scales as with system size and therefore becomes significantly more computationally demanding for solid state systems than RPA. We have confirmed that the rAPBE functional inherits the good description of lattice constants within RPA Harl et al. (2010). In particular, we have calculated the lattice constants for bulk C, Si, Na and Pd. The results are provided in the supplementary. The deviations from experimental values corrected for zero point anharmonic effects are , , , and for rAPBE(RPA) respectively.
One might speculate that the rAPBE accuracy of atomization energies and cohesive energies of solids are simply due to a better description of the isolated atoms compared to RPA. Of course, the already much improved correlation energy of the HEG compared to RPA indicates that this is not the case. As another demonstration of this fact we have calculated the formation energy of MgO defined as the energy of the reaction Mg(s) + 1/2O(g) MgO(s). The result is shown in Tab. 1. As can be seen, rAPBE reduces the RPA error by eV corresponding to a four-fold decrease in relative error. Thus, the renormalized kernel does not only improve the description of isolated atoms relative to RPA, but gives rise to a universal improvement of correlation energies that manifests itself in any calculated quantity.
One of the great success stories of RPA, is that it solves the famous ”CO puzzle”. Most GGA functionals predict the wrong adsorption site for CO on metal surfaces, whereas RPA predict the correct adsorption order Schimka et al. (2010). Furthermore, it is possible to choose a semi-local functional that gives the correct adsorption energy, but this will be at the cost of a highly inaccurate metal surface energy. In contrast, RPA produces an accurate adsorption energy and still yields a reasonable metal surface energy. As shown in Fig. 5, this trend is inherited by the rAPBE functional which produce results very similar to those of RPA. For comparison we also show the results for various GGAs Schimka et al. (2010) and van der Waals functionals Wellendorff et al. (2012).
Another successful application of the RPA method is the adsorption of graphene on metal surfaces Olsen et al. (2011); Mittendorfer et al. (2011); Olsen and Thygesen (2013,). Graphene interacts with metals through long range dispersive- and short range weak covalent interactions. The adsorption geometry and binding energy is determined by a detailed balance between these contributions which are of equal magnitude. Semi-local functionals cannot account for the dispersive interactions and in many cases do not provide the required accuracy for weak covalent interactions either. On the other hand, most van der Waals functionals account well for the dispersive interactions at large distances, but fail at shorter distances due to incorrect description of the exchange interactionWellendorff et al. (2012); Mittendorfer et al. (2011). The RPA predict metal-graphene binding distances in overall good agreement with experiments. In particular, for the case of Ni(111) and Co(0001) RPA yields two distinct minima around Å and Å corresponding to (weak) chemisorption and physisorption, respectively. The chemisorption minimum is slightly deeper in good agreement with experimental findings of Å. However, due to the inaccurate description of covalent bonds in molecules and solids one could question the accuracy of the RPA description of the chemisorption minimum. So far, RPA has been the best possible method to analyze these systems, since quantum chemistry methods are out of the question for this class of systems. We have calculated the binding energy curve for graphene on Ni(111) using the rAPBE kernel and compared it to various other methods, see Fig. 6. It is seen that rAPBE and RPA produce very similar results. As one could perhaps expect from the general tendency of RPA to underbind, the rAPBE kernel lowers the chemisorption minimum relative to the physisorption minimum. While this change is in the right direction compared to experiments, it does not change the qualitative picture of the binding obtained from RPA.
In conclusion, we have shown how to construct xc-kernels within TDDFT that extends the RPA in a natural way and improves its description of short-range correlations while retaining the good description of dispersive interactions and static correlation with similar computational cost. The proposed renormalization procedure can be applied to any known semi-local xc-functional and thus defines an entirely new class of adiabatic non-local xc-kernels which could pave the way for achieving chemical accuracy in solid state calculations.
The authors acknowledge support from the Danish Council for Independent Research’s Sapere Aude Program, Grant No. 11-1051390. The Center for Nanostructured Graphene is sponsored by the Danish National Research Foundation, Project DNRF58.
Appendix A The renormalized kernel from a weighted density approximation
Here we will show how the renormalized kernels arise naturally as a generalization of (semi-)local functionals, if one replaces the local density by a modified density obtained as an average of the original density over the extend of the xc-hole. In the case of the homogeneous electron gas, the LDA exchange-correlation functional gives the exact ground state energy (by construction). However, the adiabatic kernel obtained from LDA gives rise to poor correlation energies and has the unphysical property that the pair-correlation function diverges at the origin. Nevertheless, it has previously been shown that the main problem with the exchange-correlation kernel is not the lack of frequency dependence, but rather the local nature of the exchange-correlation kernel. Thus, the adiabatic approximation in itself does not seem to pose a fundamental problem for ground state energies within the ACDFT framework. It is therefore natural to ask whether it is possible to generalize LDA such that it gives a more consistent description of exchange-correlation in the sense that both the energy functional itself and the energy obtained from the ACDFT using the derived adiabatic kernel yield accurate ground state energies of the HEG.
To that end, we consider replacing the density in by a locally averaged density of the form
where is a suitably chosen weighting function. On physical grounds, should integrate to 1 and have a width of around which is the characteristic size of the xc-hole in the HEG. Apart from this, could be considered as a free parameter to be optimized such as to yield the best ACDFT correlation energy of the HEG. It is clear that the substitution (1) does nothing for a constant density and thus the LDA remains exact for homogeneous systems. It does, however, affect the xc-kernel which becomes non-local. Specifically we get
For inhomogeneous systems the substitution (7) will of course modify the energy and potential functionals, but the HEG represents a fixed point for which the energy and potential remains the same. Therefore, one could argue that an LDA based on (7) may be equally justified and that the usual LDA is just a special case corresponding to the choice . In fact, for an inhomogeneous system it seems physically reasonable to use the energy density of a HEG with the locally averaged density rather than the strict local density because the xc-hole is not local anyway, but extends over . Moreover, in the context of the ACFDT formalism, the use of a non-local will produce a non-local kernel, which is known to be a crucial property.
The double convolution in Eq. (8) becomes a simple product in Fourier space. A particular simple form for the kernel in Eq. (8) is obtained by choosing to be the Fourier transform of a step function of the form . To ensure that the renormalized kernel is continuous in , we must fix by the condition . We note that in the case of pure LDA exchange, this condition becomes meaning that will have a characteristic width of corresponding to the size of the xc-hole in the HEG.
The above discussion, in particular the derivation of the renormalized kernel , was mainly focused on the HEG. To obtain a kernel appropriate for inhomogeneous systems it is natural to make a local density approximation for the renormalized kernel, i.e.
This is the kernel we have investigated in the present work. It should be noted that we implicitly assumed that was independent of the density when deriving (8). Nevertheless, the physical contents and the idea of deriving non-local kernels from LDA type functionals should be clear from the above discussion.
Appendix B Implementation
The response function and exchange-correlation kernel has been implemented in a plane wave basis set. The kernel Eq. (1) is invariant under simultaneous lattice translations in and and its plane wave representation takes the form
Here and are reciprocal lattice vectors, is the number of sampled unit cells, belongs to the first Brillouin zone, and
where we have introduced the lattice point difference and used that each of the sampled unit cell integrals in Eq. (10) can be transferred into a single unit cell by letting . We also used that . The function is thus periodic in both and . In principle the sum over lattice point differences should be converged independently, but in all the calculations in the present work we have set the number of sampled unit cells in the kernel equal to the k-point sampling.
b.1 Computational details
The correlation energy is given by
where the interacting response function is calculated from the Dyson equation Eq. (1). We have only considered exchange kernels, which have the property . In principle, the method could easily be applied to kernels including correlation, since the coupling constant scaling is well-known, but this would then require an evaluation of the kernel at all sampled -points. The coupling constant integration in Eq. (12) is evaluated numerically using 8 Gauss-Legendre points. The frequency integration is also evaluated numerically using 16 Gauss-Legendre points distributed from 0-800 eV. We have checked that these samplings produce converged results for a wide range of different electronic systems.
In all calculations the total number of bands used to evaluate was set equal to the number of plane waves in the calculation. This ensures a convergence of correlation energies with respect to the plane wave cutoff . For all systems we have evaluated the correlation energies by calculating the expression Eq. (12) at a sequence of correlation energies up to 300-400 eV and then extrapolated energy differences to the infinite cutoff limit. This procedure has previously been shown to be highly accurate for RPA and the convergence properties of the rAX are very similar to that of RPA. The initial input orbitals was calculated with a cutoff of 600 eV and the Hartree-Fock energies were all calculated with 800 eV.
For most molecular and atomic systems the correlation energies were computed in a periodic unit cell with at least 6 Å between neighboring cells. For all metal atoms as well as the Cl and P atoms a separation of 8 Å was required and for Na we used a separation of 10 Å to ensure decoupling of periodic images. The Hartree-Fock and hybrid calculations were converged separately, since much larger unit cells were required for those calculations.
For the bulk solid state systems, we have used k-point samplings of 12x12x12 for Ag, Cu, Rh, Pd, Na, Al, MgO, GaN, AlN, LiF, BN, C, Si, Ge, SiC and 8x8x8 for the remaining. The Hartree-Fock calculations were converged separately with much higher k-point samplings. For the CO@Pt(111) system we used 6x6x1 k-points for the correlation energies and 16x16x1 k-points for the Hartree-Fock energy. For graphene@Ni(111) we used 16x16x1 k-points for both Hartree-Fock and correlation. All calculations were performed with a Fermi-smearing of 0.01 eV.
The calculations were performed within the projector-augmented wave framework. For the cohesive energy of solids we found that it is important to include semi-core states of Ag, Pd, Rh and Mg. More details on this point can be found in Ref. .
Appendix C Binding energies
Appendix D Lattice constants
- For simplicity, we have neglected the variation in the potential with respect to the density gradient and thus use to construct the rAPBE kernel.
- R. J. Bartlett and M. Muslat, Rev. Mod. Phys. 79, 555 (2007).
- K. G. A. A. Booth GH, GrÃneis A, Nature 493, 365 (2013).
- J. J. Shepherd and A. Grüneis, Phys. Rev. Lett. 110, 226401 (2013).
- D. Bohm and D. Pines, Phys. Rev. 82, 625 (1951).
- F. Furche, Phys. Rev. B 64, 195120 (2001).
- H. Eshuis, J. E. Bates, and F. Furche, Theor. Chem. Acc. 131, 1084 (2012).
- X. Ren, P. Rinke, C. Joas, and M. Scheffler, J. Mater. Sci 47, 7447 (2012).
- T. Olsen, J. Yan, J. J. Mortensen, and K. S. Thygesen, Phys. Rev. Lett. 107, 156401 (2011).
- T. Olsen and K. S. Thygesen, Phys. Rev. B 87, 075111 (2013,).
- J. Yan, J. S. Hummelshøj, and J. K. Nørskov, Phys. Rev. B 87, 075207 (2013).
- J. Harl and G. Kresse, Phys. Rev. Lett. 103, 056401 (2009).
- J. Harl, L. Schimka, and G. Kresse, Phys. Rev. B 81, 115126 (2010).
- L. Schimka, J. Harl, A. Stroppa, A. Grüneis, M. Marsman, F. Mittendorfer, and G. Kresse, Nature Materials 9, 741 (2010).
- M. Lein, E. K. U. Gross, and J. P. Perdew, Phys. Rev. B 61, 13431 (2000).
- F. Furche and T. V. Voorhis, J. Chem. Phys. 122, 164106 (2005).
- J. F. Dobson and J. Wang, Phys. Rev. B 62, 10038 (2000).
- J. Jung, P. Garcia-González, J. F. Dobson, and R. W. Godby, Phys. Rev. B 70, 205107 (2004).
- T. Gould, J. Chem. Phys. 137, 111101 (2012).
- A. Heßelmann and A. Görling, Mol. Phys. 108, 359 (2010).
- A. Heßelmann and A. Görling, Phys. Rev. Lett. 106, 093001 (2011).
- A. Grüneis, M. Marsman, J. Harl, L. Schimka, and G. Kresse, J. Chem. Phys. 131, 154115 (2009).
- A. Karton, E. Rabinovich, J. M. L. Martin, and B. Ruscic, J. Chem. Phys. 125, 144108 (2006).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- T. Olsen and K. S. Thygesen, Phys. Rev. B 86, 081103(R) (2012).
- T. Olsen and K. S. Thygesen, Phys. Rev. B 88, 115131 (2013).
- J. A. Alonso and L. A. Girifalco, Phys. Rev. B 17, 3735 (1978).
- J. Enkovaara et al., J. Phys.: Condens. Matter 22, 253202 (2010).
- J. Yan, J. J. Mortensen, K. W. Jacobsen, and K. S. Thygesen, Phys. Rev. B 83, 245122 (2011).
- P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
- X. Ren, A. Tkatchenko, P. Rinke, and M. Scheffler, Phys. Rev. Lett. 106, 153003 (2011).
- J. Paier, M. Marsman, K. Hummer, G. Kresse, I. C. Gerber, and J. G. Angyan, J. Chem. Phys. 124, 154709 (2006).
- J. Wellendorff, K. T. Lundgaard, A. Møgelhøj, V. Petzold, D. D. Landis, J. K. Nørskov, T. Bligaard, and K. W. Jacobsen, Phys. Rev. B 85, 235149 (2012).
- F. Mittendorfer, A. Garhofer, J. Redinger, J. Klimes, J. Harl, and G. Kresse, Phys. Rev. B 84, 201401 (2011).