Reducing charge delocalization error of density functional theory
11 Reducing charge delocalization error of density functional theory
Emil Proynov and Jing Kong* Department of Chemistry and Center for Computational Sciences, Middle Tennessee State University, 1301 Main St., Murfreesboro, TN 37130, USA * Correspondence: [email protected]
Abstract
Proper treatment of charge delocalization, besides nondynamic correlation, has been a major challenge to density functional theory. Contemporary functionals are known to undershoot the dissociation of symmetric charged dimer A , a simple but stringent test, and predict a spurious barrier and improperly delocalize charges for charged molecular clusters. We extend a functional designed to treat nondynamic correlation to treat also the charge delocalization with a modification of the secondary nondynamic correlation for parallel spins. The modified functional is shown to eliminate those problems and reduces the multi-electron self-interaction error. Furthermore, its results are the closest to that of CCSD(T) in the whole range of the dissociation compared with contemporary functionals. It also correctly localizes the net positive charge in (CH ) n+ clusters and predicts a nearly constant ionization potential as result. Overall, we show the feasibility of treating charge delocalization together with nondynamic correlation. Density Functional Theory (DFT) is the principal method of choice in theoretical chemistry and material science [1, 2]. Contemporary density functional approximations (DFA) perform well when the quantum state is dominated by a single electronic configuration. They are less effective when the nondynamic correlation (NDC) becomes significant and requires a multireference treatment. Multireference quantum effects are associated with degeneracies and near-degeneracies of the ground state and are critical to a large variety of molecular properties, such as transition-metal-ligand interactions, diradicals and molecular magnets. The simple example is the dissociation of H , where Hartree-Fock (HF) method overshoots the dissociation energy by about 190 kcal/mol and contemporary mainstream DFAs overshoot by over 50 kcal/mol [3, 4]. On the other hand, new methods have emerged that make strides in treating the NDC within the single-determinant framework of Kohn-Sham DFT. In particular, Becke has proposed the Becke’05 (B05) functional for nondynamic correlation designed to compensate the delocalization of the exact exchange, a basic characteristic of systems with nondynamic correlation [5, 6]. B05 was extended by Becke in 2013 [7, 8] to include the strong correlation occurring in the dissociation limit of covalent bonds, resulting in the Becke’13 (B13) method. We have modified B13 by introducing a model for the adiabatic connection for nondynamic correlation and reduced the number of empirical parameters significantly as a result [4]. The method is referred as Kong-Proynov’16/Becke’13 (KP16/B13) here. B13 and KP16/B13 were shown to recover most of the correlation in cases where the correlation is extremely strong. B13 and KP16/B13 reduced the errors for dissociation energies of some covalent-bonds to single digits in kcal/mol [4, 7, 8], and KP16/B13 successfully predicted the Mott metal-insulator transition [9] where other mainstream contemporary functionals failed. Both methods were also shown to be competitive in the prediction of various properties in regular systems [10]. Another major challenge to DFAs is the charge delocalization error, also known as fractional charge error. This error has a broad impact on the accuracy of approximate functionals in various applications: incorrect dissociation limits that involve fractional charge, underestimated reaction barriers and ionization potentials, unrealistic charge distributions, to mention just a few [3, 11-18]. The simplest example of this error is the dissociation of one-electron molecule like H , where most common DFAs yield errors larger than 30 kcal/mol, while HF is exact. The dissociation curves of H and other A with ‘A’ being typically a noble gas atom, using various DFAs show two main common failures: a) Contemporary DFAs undershoot the dissociation, i.e. the calculated energy of the asymptotically fractionally charged species A +0.5 is too low whereas HF overshoots it; and b) Most contemporary DFAs predict a spurious barrier along the dissociation path, whereas HF does not. A closely-related phenomenon is that local and semilocal DFAs delocalize charges when the system is not geometrically symmetric. For instance, when an electron is removed from a cluster of well separated molecules such as methanes, the positive charge should reside on one of the molecules, but DFAs would distribute the positive charge across the whole cluster [19]. Similar erratic charge distribution with DFAs is also observed when the geometric symmetry of A is allowed to break [12, 20], while HF distributes the charge properly. The charge delocalization error of DFAs has been associated with the self-interaction error [12, 21, 22], which occurs when the approximate self-exchange component of a DFA does not cancel exactly the self-Coulomb interaction. HF, on the other hand, is free of one-electron self-interaction error, although HF still contains multielectron self-interaction error [12, 22, 23]. Long-range corrected (LRC) DFAs have been a popular remedy of the charge-delocalization error, in which the long-range potential of a DFA is replaced by the HF exchange potential [18, 24, 25]. LRC functionals have been applied successfully to study charge-transfer excitations using time-dependent DFT. They remove the spurious barrier in the dissociation of A but still undershoot the dissociation curve and delocalize charges in asymmetric cations such as methane clusters [19]. That is because they still have the self-interaction error carried by the DFAs used for the short-range interaction. Another obstacle of the LRC approach is the lack of universal range-separation parameter. Other DFAs have also been designed in combination with 100% HF exchange in order to alleviate the charge delocalization error [3, 26]. However, increasing the contribution of HF exchange has the danger of playing a ‘zero-sum’ game. Janesko et al compared many contemporary DFAs for the dissociation limits of H and H among other properties and showed that DFAs that perform well for charge delocalization tend to perform poorly for NDC, and vice versa [3]. In this regard, it is worth noting the localized orbital scaling method developed by Yang’s group [11, 23] that attempts to correct both the NDC and the charge delocalization problems by straightening the exchange-correlation functionals with respect to the effective fractional orbital occupancies. In this work we apply the KP16/B13 method to the problem of charge delocalization, using A (A being a noble gas atom) and a cluster of well-separated methanes as test systems. KP16/B13 is designed for the left-right NDC, without explicit consideration of the charge delocalization. But it may provide a good foundation for performing well for the charge delocalization problem and thus break from the aforementioned zero-sum dilemma since it includes the full exact exchange. Indeed, the method is reduced to HF for one-electron system and thus exact for H . B05 functional [6], a predecessor to B13 and KP16/B13, shows a reduced charge delocalization error for some properties like H and He dissociation [27], dipole moments [28] and charge-transfer excitations [16]. In due course of applying KP16/B13 to systems with pronounced charge delocalization, we find that the part of the model that treats the NDC for parallel spins overshoots in some regions, becoming about the same or even larger in magnitude than the opposite-spin NDC. This is not physically reasonable, because the parallel-spin NDC is a higher order effect. We explore here a simple modification of the functional (named ‘mKP16/B13’) that corrects to a large extent this drawback and improves significantly on the charge delocalization problem. The KP16/B13 correlation functional models the adiabatic connection (AC) for nondynamic correlation, solving locally the integration over the AC coupling constant at each reference point with a -dependent correlation energy density [4]: d ( ) d c cndc E W w r d r (1) The resulting correlation energy expression is bz bz ACndc ndc c ndcbz e eE u d r q u d re bz , (2) where ndc u is the nondynamic correlation energy density and z is a factor of NDC strength [4]: ndcdync uz u , (3) with dync u being the dynamic correlation energy density at . The factor z is small in regions where the dynamic correlation dominates, and large where the nondynamic correlation dominates. This transition is modulated by the empirical parameter b in Eq.(2). The real space function ACc q takes into account the correlation kinetic of NDC. The form of NDC energy density ndc u used in Eq.(3) is ( ) opp parndc ndc ndc u u c u , (4) with the opposite-spin NDC energy density given by [6, 7]: ( )( )( ) ( ) ( )( ) ( ) opp ex exndc opp X X u f u u rrr r rr r (5) where opp f is the B05 local nondynamic correlation factor with the modifications described in ref.[29]: ( ) min( ( ), ( ),1) , opp f f f r r r (6) effeff( ) XX Nf N (7) Here, eff X N is the B05/B13 relaxed normalization of the exchange hole within a region of roughly atomic size. f measures how much nondynamic correlation of spin needs from the opposite spin ( ). The min() function in Eq.(6) ensures that the NDC is provided at the level of mutual need between and electrons. The scaling of the exact exchange by opp f in Eq.(5) adds the correlation effect between electrons of opposite spins at when the exchange hole is delocalized as indicated by eff X N . It counts for the main nondynamic correlation for situations such as stretched covalent bonds. eff X N should not exceed 1 in theory, but practically it had to let be as large as 2 for numerical and other accuracy related reasons, particularly in regard to the accuracy of calculated energy differences [4, 29]. The second term in Eq.(4) counts the nondynamic correlation of electrons with parallel spins. It was found that a secondary compensation for the delocalized exchange is often needed for open-shell systems, with the worst case being triplet oxygen molecule [5, 6]. The factor c in Eq.(4) is an empirical global parameter introduced in KP16/B13 to rescale the original B13 parallel-spin NDC energy density. This energy density has the form [4, 6]: (1) parndc u A M r r r (8) where A are second-order parallel-spin NDC factors, (1) M is the first-order moment of the Becke-Roussel (BR) relaxed exchange hole ( n = 1 in Eq.(9)): ( ) 20 ( ) 4 ( , ) n n effX M s h s ds r r (9) eff eff( )1 2 1 2(2) X opp X
N f N DA A A A AM (10) ii D r r (11) Here, is the KS kinetic energy density by definition of Becke [30], (2) M is the second-order moment of the Becke-Roussel relaxed exchange hole ( n = 2 in Eq.(9) ). The NDC energy density ndc u is further modified in KP16/B13 in order to decrease the fractional spin error:
11 / exp / /2 ndc ndc u u z D (12) where the parameter is determined empirically, and z is given by Eq.(3). The well-known function D is used to eliminate the one-electron self-interaction error [31]. Equation (12) boosts the recovery of nondynamic correlation towards the dissociation limit. The final KP16/B13 energy expression includes the exact exchange x E (calculated in HF like fashion) and the dynamic correlation energy. For the latter we use the -integrated form given by Eqs.(37) and (38) of ref. [7]:
13 30 q AC dynxc x c ndc
E E u dr w d dr (13) The total number of empirical parameters in this functional is three: b in Eq.(2), c in Eq.(4) and in Eq.(12). The KP16/B13 optimal values of these parameters are b c . While KP16/B13 performs well for various chemical systems with correlation strength from moderate to strong, we find that it does not perform well for the charge delocalization problem (results are discussed below.). The main issue is that the nondynamic correlation for the parallel spins, parndc u (Eq.(8)), overcompensates the delocalized exchange in this case. When the unpaired electron is delocalized in A , it does not incur significant nondynamic correlation even though the exchange hole is delocalized due to the delocalization of this electron in the front orbital. We modify here the parallel-spin NDC energy density by imposing the following local lower bound to it: opp parndc ndc u u r r (14) This modification is not expected to impact much systems such as oxygen molecule where the nondynamic correlation of the opposite spins, the main part of the nondynamic correlation, is already substantial. Indeed, the binding energy of oxygen molecule is improved somewhat with this modification from 4.79 eV to 5.41eV compared with the experimental value of 5.21eV. In order to avoid discontinuities of the SCF potential, we implement this lower bound in a smooth fashion, similar to that used in our B05 SCF implementation [29]: (0.5 )0.25( ) ( ) par par opp oppndc ndc ndc ndc opp parndc ndcopp parndc ndc pt u ut u u u u u e u (15) Here, the smoothing parameter p is adjusted to some large value that still would not affect visibly the SCF convergence pattern ( p = 500 in this case). Equation (15) provides a smooth cut-off of the parallel-spin NDC energy density parndc u , keeping its values withing the physically meaningful range. MKP16/B13 are compared with KP16/B13, HF, B3LYP and ωB97-X methods for the dissociation of A molecules (A=He, Ne, Ar) and the ionization potential and charge distribution of (CH ) n chains (n=1,2,3). B3LYP is used to represent mainstream functionals without long-range corrections. ωB97-X represents long-range corrected ones. All calculations, except with ωB97-X and CCSD(T), are done with our home made DFT code xTron using ultra-fine unpruned grid with G3LARGE [32]. ωB97-X and CCSD(T) calculations are done with the Q-Chem program [33] using G3LARGE basis set. HF and B3LYP orbitals are used as the initial guess for KP16/B13 and mKP16/B13 SCF results respectively. The dissociation curves are depicted in Figs. 1 to 3, and vertical ionization potentials and charge distributions of the methane clusters are listed in Table 1. As discussed earlier, dissociation of A molecules is an exemplary case of charge delocalization. Figs 1-3 show that mKP16/B13 improves significantly over KP16/B13, and perform well compared to all other methods. Its results are also the closest to that of CCSD(T) from near the equilibrium to the dissociation limit. Becke’88 (B88) exchange functional in B3LYP, like other local and semilocal exchange functionals, overestimates the exchange in stretched A by assuming the delocalized, parallel-spin electrons fully local, causing the underestimation of the total energy. The potential of those exchange functionals also drops too fast asymptotically, causing the spurious barrier in the dissociation curves. HF assumes the correct number of the parallel-spin electrons locally when A is well stretched Figure 1.
Dissociation of He with several functionals and CCSD(T). Figure 2.
Dissociation of Ne with several functionals and CCSD(T). Figure 3.
Dissociation of Ar with several functionals and CCSD(T). but lacks the correlation energy, and thus overestimates the total energy. Its potential, however, has the correct asymptotic behavior. LRC functionals such as ωB97-X eliminate the spurious barrier by switching to the HF exchange for long range, but still remain overestimating the exchange effect for the short range and thus undershoot the dissociation. Unmodified KP16/B13 is the closest to the CCSD(T) result at the equilibrium, showing its capability of recovering the majority of the correlation. However, it still undershoots the dissociation with a spurious barrier. KP16/B13 employs the full HF exchange and it seems surprising that it behaves more like B3LYP, albeit to a less extent. This means the correlation is overestimated, Eq.(13). Between the nondynamic and dynamic correlations, it is the nondynamic correlation that contributes to the overestimation since the dynamic correlation is of short-range. The nondynamic correlation energy density for opposite spins in this case, oppndc u in Eq.(5), is very small because all the beta electrons are localized (assuming the delocalized electron is ), i.e. eff X N is near one everywhere, and the correlation factors f and thus opp f are near zero per Eqs.(6) and (7) as a consequence. Thus, the NDC energy density for parallel spins, parndc u (Eq.(8)), is the main reason for the observed overestimation of the correlation energy. This conclusion is confirmed with numerical calculations. The reason for this overestimation is that A in Eq.(10) is larger than zero because eff X N is smaller than one due to the delocalization of the electron in the front orbital. A is near zero for stretched H and He because D is zero for regions of one-electron per spin, which explains the closeness of HF, KP16/B13 and mKP16/B13 curves. D is larger than zero in general for Ne and Ar because of multielectrons per spin in regions of the stretched molecules where density is significant, which leads to a positive A and thus a negative parndc u . Furthermore A increases as the bond is stretched because eff X N becomes smaller, causing the spurious barrier. The modification with mKP16/B13 in Eq.(14) provides a simple and significant correction to parndc u . It adheres to the original motivation of the nondynamic correlation for parallel spin as a secondary correction to the nondynamic correlation of the opposite spin, and produces better results. Furthermore, its curves for Ne and Ar are closer to the zero line than HF in the asymptotic region, correcting some of the multielectron self-interaction error. Overall, mKP16/B13 curves are the closest to that of CCSD(T), indicating its capability of recovering correlation at all ranges. The undershooting of the dissociation curve of A also means the over-delocalization of the charge since A +0.5 is predicted to have a lower energy than the average energy of A and A + . This feature is easier to demonstrate through an ensemble of well separated methane molecules with +1 overall charge [19]. The correct picture is that the positive charge is almost exclusively localized on one of the molecules, leading to an almost constant ionization potential of the neutral cluster when the number of molecules of the cluster changes. As shown in Table 1, model functionals that undershoot the dissociation, namely KP16/B13, B3LYP and ωB97-X, predict a distribution of the positive charge between the molecules, and varying ionization potential when the size of the cluster changes. On the other hand, mKP16/B13 and RHF correctly localize the positive charge correctly on one of the molecules. In conclusion, proper treatment of charge delocalization, besides nondynamic correlation, has been a major challenge to DFT. Contemporary functionals are known to undershoot the dissociation of symmetric charged dimer A , a simple but stringent test, and predict a spurious barrier. They over-delocalize charge distributions for charged molecular clusters. The popular long-range correction approach is only a partial remedy. We extend KP16/B13, a functional designed to treat nondynamic correlation, to treat the charge delocalization. A modification is made to the model to cap appropriately the secondary nondynamic correlation for parallel spins while keeping the major nondynamic correlation for opposite spins intact. The modified KP16/B13 is shown to eliminate the spurious barrier and the undershooting of the dissociation curve for A , reduce the multielectron self-interaction error, and be the closest to the result of CCSD(T) for the whole range of dissociation when compared with contemporary functionals. It is also shown to localize the net positive charge to one molecule for a series of (CH ) n+ clusters and predicts a nearly constant ionization potential as result. Overall, we show the feasibility of treating charge delocalization together with nondynamic correlation. Acknowledgement : This work received support from the National Science Foundation (Grant No. 1665344). We are grateful for Dr. Fenglai Liu’s assistance. References
1. Becke, A.D.,
Perspective: Fifty years of density-functional theory in chemical physics.
J. Chem. Phys., 2014. : p. 18A301. 2. Burke, K.,
Perspective on density functional theory.
J. Chem. Phys., 2012. : p. 150901. 3. Janesko, B.G., et al.,
Practical Density Functionals beyond the Overdelocalization−Underbinding Zero-Sum Game.
J. Phys. Chem. Lett., 2017. : p. 4314. 4. Kong, J. and E. Proynov, Density functional model for nondynamic and strong correlation.
J. Chem. Theory Comput., 2016. : p. 133. 5. Becke, A.D., A real-space model of nondynamical correlation.
Journal of Chemical Physics, 2003. : p. 2972.
Method Charges Ionization potential M1 M2 M3 CH (CH ) (CH ) wB97X +0.492 +0.003 +0.505 14.287 14.259 13.933 B3LYP +0.359 +0.283 +0.358 14.169 13.984 12.523 KP16/B13 +0.505 +0.476 +0.024 14.224 13.590 13.619 mKP16/B13 +0.995 +0.004 +0.001 14.380 14.426 14.415 HF +0.998 +0.002 +0.0005 13.321 13.447 13.293 Table1 . The charges are the charge distribution patterns in the (CH ) linear cluster with a distance of 5 Å between the monomers. Mi is the i th molecule on the linear chain. The ionization potentials are vertical first ionization potentials for the denoted neutral clusters.
6. Becke, A.D.,
Real-space post-Hartree-Fock correlation models.
Journal of Chemical Physics, 2005. : p. 64101. 7. Becke, A.D.,
Density functionals for static, dynamical, and strong correlation.
J. Chem. Phys, 2013. : p. 074109. 8. Becke, A.D.,
Communication: Calibration of a strong-correlation density functional on transition-metal atoms.
J. Chem. Phys., 2013. : p. 161101. 9. Kong, J., E. Proynov, and J. Yu,
Describing a Strongly Correlated Model System with Density Functional Theory.
J. Phys. Chem. Lett., 2017. : p. 3142. 10. Wang, M., et al., Performance of new density functionals of nondynamic correlation on chemical properties.
J. Chem. Phys., 2019. : p. 204101. 11. Li, C., et al.,
Local Scaling Correction for Reducing Delocalization Error in Density Functional Approximations.
Phys. Rev. Lett., 2015. : p. 053001. 12. Ruzsinszky, A., et al.,
Density functionals that are one- and two- are not always many-electron self-interaction-free, as shown for H+2, He+2, LiH+, and Ne+2.
J. Chem. Phys., 2007. : p. 104102. 13. Vydrov, O.A., G.E. Scuseria, and J.P. Perdew,
Tests of functionals for systems with fractional electron number.
J. Chem. Phys., 2007. : p. 154109. 14. Otero-de-la-Roza, A., G.A. DiLabio, and E.R. Johnson,
Exchange−Correla(cid:415)on Effects for Noncovalent Interactions in Density Functional Theory.
JCTC, 2016. : p. 3160. 15. Otero-de-la-Roza, A., L.M. LeBlanc, and E.R. Johnson, Dispersion XDM with Hybrid Functionals: Delocalization Error and Halogen Bonding in Molecular Crystals.
JCTC, 2019. : p. 4933. 16. Becke, A.D., S.G. Dale, and E.R. Johnson, Communication: Correct charge transfer in CT complexes from the Becke’05 density functional.
J. Chem. Phys., 2018. : p. 211101. 17. Dreuw, A., J.L. Weisman, and M. Head-Gordon,
Long-range charge-transfer excited states in time-dependent density functional theory require non-local exchange.
J. Chem. Phys., 2003. : p. 2943-2946. 18. Janesko, B.G., et al.,
Long-range-corrected Rung 3.5 density functional approximations.
J. Chem. Phys., 2018. : p. 104112. 19. Whittleton, S.R., et al.,
Density-functional errors in ionization potential with increasing system size.
J. Chem. Phys., 2015. : p. 184106. 20. Lin, Y.-S., et al.,
Long-range corrected hybrid metageneralized-gradient approximations with dispersion corrections.
J. Chem. Phys., 2012. : p. 154109. 21. Perdew, J.P. and A. Zunger,
Self-interaction correction to density-functional approximations for many-electron systems.
Phys. Rev. B, 1981. : p. 5048. 22. Cohen, A.J., P. Mori-Sanchez, and W. Yang, Development of XC functionals with minimal mayelectron self-interaction error.
J. Chem. Phys. , 2007. : p. 191109. 23. Chen, L., et al.,
Localized orbital scaling correction for systematic elimination of delocalization error in density functional approximations.
Nat. Sci. Review, 2017. : p. 203. 24. Chai, J.-D., Martin Head-Gordon, Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections.
Phys. Chem. Chem. Phys., 2008. : p. 6615. 25. Vydrov, O.A., et al., Importance of short-range versus long-range Hartree-Fock exchange for the performance of hybrid density functionals.
J. Chem. Phys., 2006. : p. 074106. 26. Zhao, Y. and D.G. Truhlar,
The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals.
Theor. Chem. Acc., 2008. : p. 215. 27. Liu, F., et al.,
Comparison of the performance of exact-exchange-based density functional methods.
Journal of Chemical Physics, 2012. : p. 114104.
28. Dale, S.G., E.R. Johnson, and A.D. Becke,
Interrogating the Becke’05 density functional for non-locality information.
J. Chem. Phys., 2017. : p. 154103. 29. Proynov, E., et al.,
Improved self-consistent and resolution-of-identity approximated Becke'05 density functional model of nondynamic electron correlation.
Journal of Chemical Physics, 2012. : p. 034102. 30. Becke, A.D. and M.R. Roussel,
Exchange holes in inhomogeneous systems: A coordinate-space model.
Phys. Rev. A, 1989. : p. 3761. 31. Becke, A.D., Correlation energy of an inhomogeneous electron gas: A coordinate-space model.
Journal of Chemical Physics, 1988. : p. 1053. 32. Curtiss, L.A., K. Raghvachari, P. C. Redfern, V. Rassolov, J. A. Pople, Gaussian-3 (G3) theory for molecules containing first and second-row atoms.
J. Chem. Phys., 1998. : p. 7764. 33. Shao, Y., et al.,
Advances in molecular quantum chemistry contained in the Q-Chem 4 program package.
Molecular Physics, 2015. (2): p. 184-215.(2): p. 184-215.