Density-functional model for van der Waals interactions: Unifying many-body atomic approaches with nonlocal functionals

Noncovalent van der Waals (vdW) interactions are responsible for a wide range of phenomena in matter. Popular density-functional methods that treat vdW interactions use disparate physical models for these intricate forces, and as a result the applicability of these methods is often restricted to a subset of relevant molecules and materials. Aiming towards a general-purpose density functional model of vdW interactions, here we unify two complementary approaches: nonlocal vdW functionals for polarization and interatomic methods for many-body interactions. The developed nonlocal many-body dispersion method (MBD-NL) increases the accuracy and efficiency of existing vdW functionals and is shown to be broadly applicable to molecules, soft and hard materials including ionic and metallic compounds, as well as organic/inorganic interfaces.

Van der Waals (vdW) interactions originate from nonlocal correlations in the quantum motion of electrons and give rise to a wide spectrum of physical phenomena from attraction between two atoms (London, 1930) to the macroscopic Casimir effect (Jaffe, 2005). As a result, vdW interactions are one of the prime targets in material modeling, which has led to a plethora of approaches that either treat vdW forces in the same way as the rest of electron correlation, or model them with effective potentials (Klimeš and Michaelides, 2012;Grimme et al., 2016;Hermann et al., 2017). They include quantum Monte-Carlo (QMC) (Ambrosetti et al., 2014a), coupled cluster methods (Yang et al., 2014), random-phase approximation (Lu et al., 2009), nonlocal density functionals (Dion et al., 2004;Vydrov and Van Voorhis, 2009), and coarse-grained approaches ranging from pairwise (Grimme et al., 2010;Becke and Johnson, 2007;Tkatchenko and Scheffler, 2009) to many-body models Silvestrelli, 2013;Caldeweyher et al., 2019).
From a theoretical perspective, this status quo is undesirable, because different models offer often disparate pictures of the nature of vdW forces, leading to incoherent understanding of vdW interactions in molecules and materials. From a practical perspective, the three main characteristics of a method are its generality, accuracy, and computational efficiency, and so far, no single method has satisfied all three requirements while being applicable to all relevant types of matter. For instance, QMC and coupled cluster are limited by computational efficiency, pairwise approaches and two-point vdW functionals lack in accuracy for nanostructured and supramolecular compounds, and atomic models have qualitative problems with ionic and hybrid metal-organic systems.
In this work, we present a unified density-functional model of vdW interactions that couples polarizability density functionals and atomic models, inheriting broad applicability of the former and excellent accuracy of the latter. We integrate the polarizability functional of Vydrov and Van Voorhis (2010b) (VV), normalization to exact free-atom vdW parameters of the Tkatchenko-Scheffler (TS) model (Tkatchenko and Scheffler, 2009), normal- * Emails: science@jan.hermann.name, alexandre.tkatchenko@uni.lu ization to jellium via a zero-gradient limit from the VV10 functional (Vydrov and Van Voorhis, 2010a), and the Hamiltonian for the dispersion energy of the many-body dispersion (MBD) model . Compared to the range-separated selfconsistently screened (rsSCS) variant of MBD (Ambrosetti et al., 2014b), the VV polarizability functional enables consistent treatment of ionic compounds, normalization to the free-atom reference balances the accuracy of the VV polarizability across the periodic table, and normalization to jellium enables effective modeling of metals and their surfaces (Ruiz et al., 2012). The new model involves a similar level of empiricism as MBD@rsSCSwe remove the tabulated vdW radii and short-range screening, while introducing a mechanism to avoid double counting of electron correlation in near-uniform density regions. We demonstrate on a series of benchmark calculations that the new model enables for the first time a consistent treatment of vdW interactions in molecular, covalent, ionic, metallic, and hybrid metal-organic systems.
Some of the problems of MBD@rsSCS have been previously treated by Gould et al. (2016). Their fractionally ionic variant of MBD@rsSCS uses iterative Hirshfeld partitioning in combination with a piecewise linear dependence of atomic polarizability on charge, together with a rescaling scheme for the diverging MBD Hamiltonian in highly polarizable systems or with dipole smearing (Kim et al., 2020). Our approach is instead based on a general polarizability functional.
Atomic models, such as MBD, require an atomic response model in the form of static polarizabilities 0, ≡ (0) and 6, coefficients. In the new model, dubbed MBD-NL, we parametrize the response of atoms by coarse-graining the VV polarizability density to atomic fragments (Hirshfeld, 1977;Nakai, 2009, 2010). The VV polarizability functional is a semilocal functional of the electron density ( ), which models the local dynamic polarizability density (Vydrov and Van Voorhis, 2010b), where i is imaginary frequency and is an empirical param-  (Tkatchenko and Scheffler, 2009). In contrast, the present model is exact by construction (Eq. 4).
eter. The atomic dynamic polarizabilities are obtained by partitioning the polarizability density with Hirshfeld weights H ( ) = free ( )∕ ∑ free ( ), The 6 coefficients are then calculated directly from (i ) via the Casimir-Polder formula (McLachlan, 1963), Unlike approaches that use Hirshfeld fragments to define atomic volumes, MBD-NL is independent of the choice of a particular atomic partitioning, because this influences only local redistribution of the polarizability between atoms, conserving the total polarizability. Our approach is also different from that of Silvestrelli (2008), where the electron density is coarse-grained first and a polarizability functional is evaluated over the fragment densities. Already this bare combination of MBD and the VV polarizability functional substantially improves the description of ionic systems compared to MBD@rsSCS, because the VV functional gives a good estimate of the ionic polarizabilities, unlike the Hirshfeld volume scaling used in MBD@rsSCS. However, this bare combination suffers from two fundamental shortcomings. First, the polarizability functional is not evenly accurate across the periodic table. Second, when combined with semilocal densityfunctional theory (DFT), it suffers from double counting of electron correlation in regions of slowly-varying electron density. To solve these two challenges, we normalize the atomic VV polarizabilities and 6 coefficients to exact values for free atoms (Tkatchenko and Scheffler, 2009), and then normalize MBD-NL to give zero vdW energy for jellium by subtracting the portion of the polarizability from slowly-varying electron-density regions.
The VV polarizability functional is approximate, which is manifest already for free-atom polarizabilities and 6 coefficients, where accurate reference values are known ( Figure 1). It especially underestimates the vdW parameters of metallic elements. To mitigate this, we normalize the VV atomic quantities with the ratio of the respective free-atom values obtained from accurate reference calculations and from the VV functional, Many exchange-correlation (XC) functionals are exact for jellium by construction, even though the portion of electron correla-tion from the nonlocal plasmons is long-ranged and should not be included in semilocal XC functionals. As a result, most XC functionals describe accurately the electron correlation within slowlyvarying density regions, such as found in metals, and those cases require no addition of vdW forces. This is different in most general systems, in which semilocal functionals neglect long-range vdW interactions. At the same time, these metallic-density regions contribute dominantly to the VV polarizability, and hence to the vdW energy in any vdW model that would use the VV functional directly. When combined with semilocal DFT, this would result in overpolarization and overbinding of bulk metals, and of adsorbates on metallic surfaces. To avoid this double counting, the VV10 expression for the vdW energy subtracts the limit of VV10 as the density gradient approaches zero (Vydrov and Van Voorhis, 2010a) Such an approach cannot be used directly in a many-body model such as MBD, because unlike in a pairwise model the many-body vdW energy is not linear in the polarizability.
To ensure the correct zero limit of MBD-NL for uniform densities, we smoothly cut off the contribution of jellium-like regions to the polarizability. These regions can be distinguished with the combination of two local electron-density descriptors: the local ionization potential (Gutle et al., 1999) and the iso-orbital indicator (Becke and Edgecombe, 1990;Kümmel and Perdew, 2003;Sun et al., 2013), is the positive kinetic energy density of occupied orbitals , which for single-orbital densities reduces to the von Weizsäcker kinetic energy density, W [ ] = | | 2 ∕8 , and for jellium to unif [ ] = 3(3 2 ) 2∕3 5∕3 ∕10. The local ionization potential is a form of a reduced gradient with the units of energy, which attempts to model the electronic gap locally. The density gradient alone is insufficient to characterize metallic densities. In particular, both ∼ 0 and ∼ 1 must be true for density to be metallic, whereas ∼ 0 and ∼ 0 corresponds to centers of covalent bonds, and ∼ 0 and ≫ 1 signifies overlap of electron-density tails between noncovalently bound systems. Since the normalization of VV10 to jellium uses only the density gradient, it partially omits contributions from covalent bonds. By using also the iso-orbital indicator, we make MBD-NL more precise in this regard. In practical calculations, the evaluation of the kinetic energy density is the computationally most demanding part of MBD-NL, but this means that its cost is only a fraction of a single self-consistent cycle of a regular meta-GGA KS-DFT calculation. Figure 2a presents polarizability density distributions of and in three benzene compounds and in a set of simple solids (Zhang et al., 2018). In an organic molecule such as benzene (Figure 2a), the vast majority of the polarizability comes from electron density with > 5 eV, with a small part from low-gradient regions with < 1. The intermolecular interactions in the benzene dimer and crystal add a significant amount of polarizability in regions with ≫ 1, despite the electron density being low there. A richer spectrum of patterns is found in simple solids (Figure 2b). Most similar to the benzene compounds are semiconductors. In contrast, the polarizability in main-group metals is dominated by jellium-like regions near ( , ) = (0, 1). In transition metals, the  (Zhang et al., 2018). Each distribution is normalized to 62 (a. u.), the VV polarizability of a benzene monomer, for a single color scale with a.
polarizability is distributed in a wider range of the local gap along the 1 < < 2 strip, with a larger part still in the low-gradient regions. In simple ionic solids, most of the polarizability comes from single-orbital regions ( < 1).
To avoid the double counting of vdW interactions of lowgradient densities, we smoothly cut off their contribution to the polarizability functional, We impose two simple requirements on this cutoff. First, the density regions with a local gap lower than the work function of conductors should not contribute to the calculated vdW energy, because those are assumed to be covered by a semilocal XC functional. We chose the cutoff value of 5 eV, which is around the peak of the work function of elemental metals. Second, the VV polarizability of simple covalent compounds (exemplified by a benzene molecule) should not be influenced by the cutoff. The following function satisfies these two requirements:  Figure 3 | Distributions of relative changes in atomic static polarizabilities and 6 coefficients from monomers to dimers. The distributions are calculated over all atoms from all complexes in the S66 data set (Řezáč et al., 2011).
Function consists of a logistic function centered at 5 eV and of function taken from the SCAN functional (Sun et al., 2015), where it also interpolates between = 0 and = 1 (see Appendix for a plot of ( , )). We find that = 1∕10 ensures that the effect of the cutoff on the VV polarizability of a benzene molecule is negligible (< 2%). The performance of the resulting model depends only weakly on the precise values of the chosen parameters, as long as the local gap cutoff sufficiently covers the work function of a given conductor. Nevertheless, a more rigorous formulation of the model in this direction would be desirable. Apart from avoiding the double counting of long-range electron correlation, the cutoff function removes another deficiency of the VV polarizability functional. When molecules form vdW-bound compounds, the introduction of density-tail overlaps significantly increases the VV polarizability compared to the monomers (Figure 2a). This effect is an artifact of the VV functional that causes increasingly large vdW-bound systems to be overbound, and cutting off the polarizability of low-gradient regions with > 1 eliminates this issue without affecting the polarizabilities of isolated monomers (Figure 3).
Finally, the static polarizabilities and 6 coefficients calculated as described above are directly used in the MBD Hamiltonian to obtain the vdW energy. This Hamiltonian describes a system of charges in harmonic potentials-Drude oscillatorscharacterized by their static polarizabilities 0, and resonance frequencies = 4 6, ∕3 2 0, , and interacting via a long-range dipole potential lr ( ) ≡ lr ( ) ( ) (Lucas, 1967;, where ≡ √ are displacements of the charges weighted with masses (which have no effect on the energy). The interaction energy of this model system-the vdW energy-is obtained by direct diagonalization yielding a set of coupled oscillation fre-quencies̃ , In MBD-NL, we use the same long-range coupling lr as in the MBD@rsSCS variant (Ambrosetti et al., 2014b), lr ( ) = 1 1 + e −6( ∕ ( vdW + vdW )−1) (10) but with a simplified definition of the vdW radii. Rather than tabulated vdW radii, we use the quantum-mechanical formula for vdW radii of free atoms from Fedorov et al. (2018), and scale them similarly to the vdW parameters as in (4) We optimize the damping parameter of lr ( ) on the S66 data set (Řezáč et al., 2011), as was done for MBD@rsSCS, and find the optimal values of 0.81 and 0.83 for the XC functionals PBE (Perdew et al., 1996) and PBE0 (Adamo and Barone, 1999), respectively, only slightly smaller than the values of 0.83 and 0.85 for MBD@rsSCS. Next, we briefly describe several benchmark tests of MBD-NL (see Appendix for a more detailed description of the calculations and for additional results). On a set of small organic dimers (S66, Řezáč et al., 2011), MBD-NL performs nearly identically to MBD@rsSCS (Figure 4a), which is already excellent for a DFT+vdW approach. In contrast, the errors in lattice energies of 64 hard solids (Zhang et al., 2018) are re-duced drastically when MBD@rsSCS is replaced with MBD-NL (Figure 4b). This improvement comes mainly from the errors for metals and ionic solids, which MBD@rsSCS overbinds substantially, whereas plain PBE performs reasonably well, and MBD-NL retains this good performance. MBD-NL still somewhat overbinds the metals compared to PBE, as could be expected, because bare PBE does not underbind the metals despite the missing (non-jellium) long-range vdW interactions. Ionic solids are underbound by 4% with PBE, which is reduced nearly to zero when the missing nonlocal correlation is added by MBD-NL, whereas MBD@rsSCS overbinds some of them substantially. The performance of MBD-NL on semiconductors is similar to MBD@rsSCS. On a set of organic molecular crystals (X23, Reilly and Tkatchenko, 2013), MBD-NL performs nearly identically to MBD@rsSCS, with a similar tendency to underbind (2%) as MBD@rsSCS has to overbind (3%). On a set of supramolecular complexes (S12L, Grimme, 2012), the accuracy of MBD-NL is reduced compared to MBD@rsSCS, from 5% to 9% in terms of the mean absolute relative error (MARE), but the accuracy of the two methods is equal with the PBE0 functional, with MBD-NL having a smaller mean relative error compared to MBD@rsSCS.
Compounds with small or zero electronic gap pose the hardest problem for DFT+vdW approaches, because such systems require in principle long-range coupling of delocalized electronic fluctuations. Despite that, MBD-NL reaches the accuracy of established effective models for hybrid interfaces of metallic surfaces and organic molecules, such as the MBD@rsSCS[surf] method (Ruiz et al., 2016), with a difference in the binding energy between the two methods below 10% for a benzene molecule on a silver (111) surface. This is only possible because the long-wavelength electronic fluctuations in the metal have no correlation counterpart in the adsorbed molecule, so a fully delocalized treatment of the fluctuations is not necessary in this case.
In contrast, the delocalized fluctuations cannot be effectively neglected in layered vdW materials with small band gaps, such as the transition-metal dichalcogenides (TMDCs), which comprise 23 of the benchmark set of 26 layered materials (here dubbed "26", Björkman, 2012). MBD@rsSCS and VV10 overbind the "26" set by 10% and 52%, respectively, indicating that both models overpolarize these small-gap layered compounds. In contrast, the nonlocal part of the polarizability from low-gradient density regions is removed in MBD-NL, resulting in its underbinding of the "26" set by 21% (the accuracy of the reference calculations is 10%-20%). Of the three methods, the three non-TMDC layered materials in the "26" set (graphite, BN, PbO) are described most accurately by MBD-NL (MARE of 7%, compared to 27% for MBD@rsSCS and 53% for VV10).
Before concluding, we discuss some open questions regarding MBD-NL. First, the VV polarizability functional is semiempirical and it can be improved by including nonlocal density information, for example by developing a meta-GGA polarizability functional. Such an extension could improve the overall accuracy significantly, but requires nontrivial advances in the general theory of polarizability functionals. Another possibility would be to normalize the vdW parameters not only to free atoms, but also to ions (Gould et al., 2016), which is comparably more straightforward. Second, although MBD-NL can effectively treat hybrid interfaces between organic and metallic compounds, it does not capture the truly nonlocal electronic fluctuations that can be found in conductors (Dobson, 2014). Incorporating such a mechanism would not only enable MBD-NL to treat long-range interactions between fully metallic bodies, but also increase its accuracy for interacting systems of small-gap compounds, such as TMDCs. How to do this in practice is at the moment unclear, and we see this as the largest remaining theoretical gap in the general understanding of vdW interactions in materials. Third, MBD-NL uses an empirical range-separating function parameterized for a given XC functional. XC functionals differ substantially in their mid-range behavior (unlike the PBE and PBE0 functionals used here), and developing seamless range-separation approaches that couple semilocal XC functionals and vdW methods in a universal and transferable way remains an open challenge (Hermann and Tkatchenko, 2018).
In conclusion, we have developed a vdW model that unifies atomic many-body methods and nonlocal vdW functionals. By normalizing to free atoms and jellium, we have retained the accuracy of best DFT+vdW approaches while extending applicability to ionic and metallic compounds, and hybrid metal-organic interfaces. Our approach enables efficient, accurate, and consistent modeling of many-body vdW interactions in a substantially broader range of systems than previously possible. were used for all DFT and MBD calculations. For hard solids, we have used the -point density from (Zhang et al., 2018). All molecular and crystal geometries were taken directly from the respective benchmark sets without any relaxation. Table 1 reports the performance of MBD-NL, MBD@rsSCS, and VV10, in combination with the PBE and PBE0 functionals, on the set of organic molecular crystals (X23, Reilly and Tkatchenko, 2013), a set of supramolecular complexes (S12L, Grimme, 2012), and a set of 26 layered materials (dubbed "26", Björkman, 2012). Of the standard vdW data sets, S12L is the only one where MBD-NL achieves a different performance with the PBE and PBE0 functionals. This results mostly from PBE binding the π-π complexes somewhat more than PBE0. The proper balance between semilocal DFT and long-range vdW models in the case of large π-π complexes is unclear (Hermann and Tkatchenko, 2018). On the "26" set, the MBD@rsSCS Hamiltonian has negative eigenvalues for 20 of the 26 compounds. To obtain finite energies nevertheless, we use the eigenvalue rescaling as proposed by Gould et al. (2016). Figure 5 compares the binding energy curve of a hybrid organic/inorganic interface as calculated by PBE-NL and surface variants of the MBD@rsSCS and TS methods by Ruiz et al. (2016). Table 2 presents timings of illustrative DFT+MBD calculations and their individual components for a simple inorganic material and a large organic complex. In both cases, the evaluation of the MBD energy is only a small fraction of the cost of the DFT calculation, even when the evaluation of the kinetic energy density needed for parametrization of MBD-NL is included in the cost of the MBD calculation.
The function ( , ) from eq. (6) of the main text is visualized in Figure 6.