Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data

We present a parameter-free method for an accurate determination of long-range van der Waals interactions from mean-field electronic structure calculations. Our method relies on the summation of interatomic C6 coefficients, derived from the electron density of a molecule or solid and accurate reference data for the free atoms. The mean absolute error in the C6 coefficients is 5.5% when compared to accurate experimental values for 1225 intermolecular pairs, irrespective of the employed exchangecorrelation functional. We show that the effective atomic C6 coefficients depend strongly on the bonding environment of an atom in a molecule. Finally, we analyze the van der Waals radii and the damping function in the C6R 6 correction method for density-functional theory calculations.

Noncovalent forces, such as hydrogen bonding and van der Waals (vdW) interactions, are crucial for the formation, stability, and function of molecules and materials.At present, ubiquitous vdW interactions [1] can only be accounted for properly by high-level quantum-chemical wave function or by the Quantum Monte Carlo (QMC) method.In contrast, the correct long-range interaction tail, e.g., for separated molecules, is absent from all popular local-density or gradient corrected exchange-correlation (xc) functionals of density-functional theory (called DFT-XCA in what follows), as well as from the Hartree-Fock (HF) approximation [2,3].
Many encouraging concepts and methods have been proposed to include vdW interactions in DFT calculations [4][5][6][7][8][9].Ultimately, a functional which is able to account for vdW interactions in a ''seamless'' manner is desirable.The Chalmers-Rutgers approach [4,10] is a step in that direction, but at this point, the performance is not certain and, for example, errors are as large as 70% for the binding energy of rare-gas dimers-prototypical vdW systems [10].
A popular remedy for the missing vdW interaction in present-day DFT consists of adding a pairwise interatomic C 6 R À6 term (E vdW ) to the DFT energy [5,6,[11][12][13][14], where R AB is the distance between atoms A and B, C 6AB is the corresponding C 6 coefficient, R 0 A and R 0 B are the vdW radii.The R À6 AB singularity at small distances is eliminated by the short-ranged damping function f damp ðR AB ; R 0 A ; R 0 B Þ. Obviously, at the medium and short range, this approach can only work together with xc functionals that underestimate the binding energy.In particular, Grimme has proven such schemes to be accurate for a range of molecular systems [6,12].A serious shortcoming of the C 6 R À6 schemes is their empirical nature, since the parameters do not depend on the electronic structure, but are rather ob-tained by fitting to experimental C 6 coefficients and/or post-Hartree-Fock binding energy data.Furthermore, we note that the damping function will also correct (or affect) other properties of the employed xc functional at short distances.Though a correction of present-day xc functionals is necessary, it is not satisfactorily handled by such an approach.Several methods exist to determine the C 6 coefficients either from ground-state orbitals [8,9,15] or timedependent DFT (TDDFT) [16,17].Unfortunately, the errors are quite large (15%-20% on average and maximum deviation of 40%-60%) [8,15,16].The origin of such errors is related to the grossly overestimated polarizability in DFT-XCA [18,19].
In this Letter, we develop and assess a scheme to determine the C 6 coefficients and vdW radii from the mean-field ground-state electron density (DFT-XCA or HF) for molecules and solids.It is largely independent of the employed DFT-XCA approximation [tested here for local-density approximation (LDA), Perdew-Burke-Ernzerhof (PBE), and Becke-Lee-Yang-Parr (BLYP) functionals], and it shows a mean absolute error of 5.5% for intermolecular C 6 coefficients on a database of experimental dipole oscillator strength distribution (DOSD) data of Meath and coworkers for 1225 complexes (see, e.g., Refs.[5,8,20,21]).The DOSD for a given molecule is the (differential) dipole oscillator strength df=dE, as a function of excitation energy E, from the electronic absorption threshold E 0 to very high energies.Many important molecular dipole properties can be evaluated as integrals involving the DOSD [22].The critical idea of our method is to use the electron density to compute the relative and not the absolute polarizability of an atom in a molecule.The method includes charge polarization effects in a transparent way, as shown for different atomic hybridization states and hydrogen bonding.
In order to develop our method, we start with the exact expression (Casimir-Polder integral) for the leading isotropic C 6 term describing the vdW interaction between two atoms or molecules A and B (Hartree atomic units used The American Physical Society throughout) [23], where A=B ði!Þ is the frequency-dependent polarizability of A and B evaluated at imaginary frequencies.We retain only the leading term in the Pade ´series [24] for where 0 A is the static polarizability of A and A is an effective frequency.Substituting 1 ði!Þ of Eq. ( 3) for ði!Þ in Eq. ( 2) yields the London formula [25], For A ¼ B, we obtain (5) with a corresponding expression for B .Combining Eqs. ( 4) and ( 5), we arrive at a formula for C 6AB which depends only on homonuclear parameters C 6AA , C 6BB , 0 A , and 0 B , For the free-atom reference values of 0 A and C 6AA , we rely on the database of Chu and Dalgarno [16], which reports self-interaction corrected TDDFT values scaled to reproduce accurate all-order many-body calculations for rare gases, alkalis, and alkaline earth atoms.The accuracy of the free-atom results is presumed to be better than 3% for 0 A and C 6AA for nonmetallic elements (1% for rare gases, alkalis, and alkaline earth atoms).Using these homonuclear values along with Eq. ( 6), we obtain a mean absolute relative error (MARE) of just 2.7% on a database of 70 heteronuclear C 6 coefficients between light elements, rare gases, alkalis, and alkaline earth atoms from accurate many-body calculations [26][27][28].The almost perfect correlation is shown in Fig. 1.
Let us now define the C 6 coefficient for an atom inside a molecule or a solid.This requires the definition of the effective volume, referenced to the free atom in vacuo.We take advantage of the direct relation between polarizability and volume [29], and employ the Hirshfeld partitioning of the electron density for the latter [8,30], where i A is the proportionality constant between volume and polarizability for the free-atom and atom-in-amolecule, w A ðrÞ is the Hirshfeld atomic partitioning weight for the atom A, r 3 is the cube of the distance from the nucleus of an atom A, nðrÞ is the total electron density, n free A ðrÞ is the electron density of the free atom A, and the sum goes over all atoms in the system.Both nðrÞ and n free A ðrÞ are calculated from DFT-XCA.The effective coefficient C eff 6AA for an atom in a molecule is determined in the following way from Eqs. ( 4) and ( 5): We assume the proportionality constant Since the C 6 coefficients are additive [22,31], the intermolecular C 6 coefficient, C mol 6 , is given by the sum of all interatomic contributions where M 1 and M 2 refers to the first and the second molecule, respectively.
To assess the accuracy of our scheme, it was benchmarked on a database of 1225 intermolecular C 6 coefficients, derived from pseudo DOSD data of Meath and coworkers (see, e.g., Refs.[5,8,20,21]).The database con- tains the C 6 coefficients for the interaction between 8 atoms and 42 molecules (organic and inorganic, from small dimers to C 8 H 18 ).The geometry of every molecule was fully relaxed using the FHI-AIMS [32] code together with LDA, PBE, and BLYP functionals.The atomic C eff 6 coefficients were calculated from Eq. ( 9) using the DFT density and then the molecular C mol coefficients were computed using Eq.(10).The correlation between the calculated C 6 coefficients and reference results is shown in Fig. 1 and compared with other methods for obtaining the C 6 coefficients from the ground-state density.With a MARE of 5.5%, we note that our scheme is a factor of 2-3 more accurate than existing methods [8,9,15].Furthermore, if the H 2 molecule is excluded from the database, the MARE drops down to 4.5%.The C 6 coefficients vary only slightly for different exchange-correlation (xc) functionals.The MARE in the C 6 coefficients between LDA and PBE is 1.4% (maximum deviation of 3.8%), while for BLYP and PBE it is 0.68% (maximum deviation of 2.1%).Clearly, the difference in the electron density between xc functionals is compensated when computing the ratio in Eq. (7).
In Table I, we show the values of the C 6 coefficients for various atoms in the DOSD database.It is encouraging that the C 6 coefficients correspond very closely to those empirically fitted by Wu and Yang [5] for different hybridization states of C and O atoms.This further indicates that our scheme correctly accounts for different atomic environments without the need of empirical parameters.The largest variations occur for carbon and silicon, due to various possible hybridization states for these elements.
In order to further illustrate the change in the atomic C 6 coefficients as a function of molecular bonding and geometry, we show the case of the hydrogen-bonded water dimer in Fig. 2. As the H 2 O molecules approach to form a hydrogen bond, the polarizabilities and the C 6 coefficients of all atoms are modified.The hydrogen involved in the hydrogen bond becomes significantly more polarizable along with the donor oxygen.On the other hand, the acceptor oxygen becomes less polarizable along with the attached hydrogens.The plot in Fig. 2 was done for a fixed water dimer geometry to illustrate that the change in the C 6 coefficients is a purely electronic effect.Relaxing the geometries for every O-O distance yields a similar plot for distances larger than the equilibrium one.
We carried out preliminary tests of our scheme for solids, calculating the C 6 coefficient of the carbon atom in a graphene sheet and in a diamond crystal at experimental geometries.For carbon in a graphene sheet, we get a value of 33.0, close to the expected value of 30.3 for the sp 2 hybridization in benzene.For diamond, we get a value of 38.6, significantly larger than 24.1 for the sp 3 hybridization in methane.
Our method can in principle treat all elements, including ions, on the same footing.Most problematic cases are those where the concept of atoms-in-molecules cannot be applied.Clearly, a pairwise summation of C 6 R À6 interactions can fail qualitatively for metallic low-dimensional systems due to nonadditive higher-order effects [33].Our scheme can be further generalized to higher-order van der Waals coefficients (C 8 , etc.) since approximations similar to the London formula are known [22].The Axilrod-Teller-Muto three-body term [22] can also be calculated since it involves an integral similar to the Casimir-Polder one.
Let us now briefly discuss the coupling of the above long-range scheme for correcting DFT-XCA calculations for shorter distances.For this, the damping function 1) must be defined.The vdW radii, R 0 A=B , are not experimental observables, unlike the C 6 coefficients.However, a rigorous theoretical definition does exist: the vdW radius corresponds to half of the distance between two atoms where the Pauli repulsion balances the London dispersion attraction [34].Using the definition of the effective atomic volume in Eq. ( 7), the vdW radius of an atom in a molecule becomes According to the above definition, the free-atom vdW radii, R 0 free , for rare-gas atoms correspond to the equilibrium distance of rare-gas dimers.Unfortunately, the vdW radii of Bondi [34] for other elements cannot be used since they correspond to atoms-in-molecules case.The simplest an-TABLE I. Free-atom C 6 coefficients (hartree Á bohr 6 ) from Chu and Dalgarno [16] along with atom-in-a-molecule (minimum and maximum) C 6 coefficients for various atoms computed for molecules from the DOSD database.For CNOH atoms, these values are compared to empirical results of Wu and Yang [5]   satz for defining consistent free-atom vdW radii comes from the electron density for the (spherical) free atoms.The electron density contour value corresponding to the vdW radius can be determined for the rare-gas atoms and then used to define R 0 free for other elements in the same row of the periodic table.We used the above ansatz for the light elements (CNOH), to illustrate the performance of our nonempirical C 6 scheme with DFT.We use a Fermi-type damping function [5], where s R are free parameters.The d parameter adjusts the damping function steepness.Our analysis indicates that the results change negligibly for 12 < d < 45.A small value of d affects covalently bonded systems, whereas a large value yields kinked binding energy curves.The choice of d ¼ 20 turned out to satisfy both constraints as tested for binding energy curves of raregases and vdW-bonded organic molecule dimers (CH 4 and benzene), as also shown by Grimme in previous work [12].Therefore, the scaling coefficient s R remains as a single empirical parameter which determines the onset of the vdW correction for a particular xc functional in terms of the distance [14].We use the S22 database of Jurecka et al. [35] to obtain the s R parameter [36].The database reports converged CCSD(T) binding energies for 22 different dimers with varying interaction strength, from a weakly vdW-bonded CH 4 dimer (23 meV) to a hydrogen-bonded uracil dimer (0.9 eV).The performance of our scheme, when coupled with the PBE functional, is significantly better than for highly empirical C 6 R À6 approaches [6,12,14].The mean absolute error (MAE) of our approach on the S22 database is 13 meV (20, 13, and 6 meV on hydrogen-bonded, vdW-bonded and mixed systems, respectively).This can be compared to 20 meV (29, 20, and 9 meV) when using empirical C 6 coefficients and vdW radii from Ref. [14].It is especially encouraging that the overestimation of the hydrogen-bonded systems is significantly reduced due to larger effective vdW radii for atoms involved in hydrogen bonds.Since the atomic C 6 coefficients are functionals of the electron density [Eq.( 9)], the potential due to the energy expression in Eq. ( 1) should enter the Kohn-Sham equations in DFT calculations.However, we do not expect a self-consistent treatment to give major changes as also noticed in Ref. [10].This aspect will be addressed in future work.
In summary, we have presented an accurate nonempirical method to obtain molecular C 6 coefficients from ground-state electron density and reference values for the free atoms.Our scheme can also be used to improve the description of weakly bonded systems in DFT for a range of xc functionals.
A. T. acknowledges the Alexander von Humboldt (AvH) foundation for funding.
this choice to be remarkably good for a large variety of molecules.Indeed, this approximation breaks down only for the smallest H 2 molecule, with deviation of 44% for the C 6 coefficient.Already for molecules such as N 2 and CO 2 , the error is less than 10%.Since the static polarizability of a molecule cannot be expressed as a linear combination of atomic polarizabilities in general, and thus free A for different hybridization states.
FIG. 2 (color online).Dependence of the atomic C 6 coefficients on the O-O distance in the water dimer.Note the two different scales on both sides for H and O atoms.