Accurate and Efficient Method for Many-Body van der Waals Interactions

An efficient method is developed for the microscopic description of the frequency-dependent polarizability of finite-gap molecules and solids. This is achieved by combining the Tkatchenko-Scheffler van der Waals (vdW) method [Phys. Rev. Lett. 102, 073005 (2009)] with the self-consistent screening equation of classical electrodynamics. This leads to a seamless description of polarization and depolarization for the polarizability tensor of molecules and solids. The screened long-range many-body vdW energy is obtained from the solution of the Schrödinger equation for a system of coupled oscillators. We show that the screening and the many-body vdW energy play a significant role even for rather small molecules, becoming crucial for an accurate treatment of conformational energies for biomolecules and binding of molecular crystals. The computational cost of the developed theory is negligible compared to the underlying electronic structure calculation.

van der Waals (vdW) interactions are ubiquitous in nature, playing a major role in defining the structure, stability, and function for a wide variety of molecules and materials.Thus, the accurate description of vdW interactions is essential for improving our understanding of many biological, chemical, and (hard and soft) condensed matter.Many encouraging ideas and methods have been proposed in recent years for approximately including the missing long-range vdW interactions in generalized gradient and hybrid density-functional theory (DFT) methods (see, e.g., Refs.[1][2][3][4] and references therein).However, the efficiency of popular methods for vdW interactions hinges on two major approximations: (1) the neglect or only effective treatment of electrostatic screening and (2) the neglect of nonadditive many-body vdW energy contributions beyond the pairwise approximation.
The recently developed Tkatchenko-Scheffler (TS) vdW method [5] computes the vdW energy for atoms in molecules using the ground-state electron density and accurately includes hybridization effects for the polarizability.However, the TS-vdW scheme does not include long-range electrostatic screening extending beyond the range of the exponentially decaying atomic densities.Neglecting retardation effects, the two-body vdW energy originates from the electrostatic interaction of ''atomic'' dipolar fluctuations.When an atom is embedded in a condensed phase (or in a large molecule), dipolar fluctuations differ from their free atom counterpart not only because of the local environment but also because a fluctuating dipole interacts electrostatically (via a power law decay) with the more distant fluctuating dipoles.Both short-and long-range effects of the environment act to screen the atomic dipolar fluctuation, and it is important to include these effects in first-principles calculations.In this Letter, we propose to extend the original TS-vdW method for molecules and solids by self-consistently including long-range screening effects in the effective atomic polarizabilities.Even with this improvement, vdW energies beyond two-body interactions are not included, and this Letter proposes a scheme based on the coupled fluctuating dipole model (CFDM) [6] to account for many-body effects.It should be mentioned that in practical applications the treatment of effects beyond two-body has been limited mostly to the three-body Axilrod-Teller formula [7] (a procedure that might even give the wrong sign for the many-body vdW energy) [8,9].The adiabatic-connection fluctuation-dissipation theorem (ACFDT) is an exact expression for the total electron correlation energy [10]: Thus, it clearly contains the full many-body vdW energy, and the formulas proposed in this Letter can be derived with approximations from the ACFDT formula [11].
We assume that the system (molecule or condensed matter) has a finite electronic (highest occupied molecular orbital-lowest unoccupied molecular orbital) gap and can be divided into effective atomic fragments, for example, by using the Hirshfeld electron-density partitioning scheme.The effective frequency-dependent dipole polarizability of every atom is defined as [12] (1) where p ½nðrÞ is the static polarizability of an atom p, and !p ½nðrÞ is the corresponding characteristic excitation frequency [5].Since both parameters are referenced to highly accurate free-atom reference data [5], short-range quantum-mechanical exchange-correlation effects are accounted for in p and !p by construction.However, the TS-vdW definition of the effective atomic polarizability lacks the aforementioned long-range electrostatic screening.In fact, atoms located inside molecules and materials experience the dynamic electric field created by the surrounding atoms, which gives rise to polarization and depolarization effects and is largely responsible for the anisotropy in the molecular polarizability [13,14].These effects can be included microscopically by modeling the environment as a dipole field and solving the resulting classical electrodynamics self-consistent screening (SCS) equation [13][14][15] where TS ðr; i!Þ is the sum of the TS-vdW effective atomic polarizabilities, and T ðr À r 0 Þ is the dipole-dipole interaction tensor (Hartree atomic units used throughout).
The SCS equation requires the dipole polarizability TS ðr; i!Þ as a continuous function over space, while the TS-vdW method in principle only yields point polarizabilities.To extend this description, we note that Eq. ( 1) is the frequency-dependent dipole polarizability of the spherical quantum harmonic oscillator (QHO), with the following Gaussian density where R is the width of the Gaussian function.
By representing the N atoms in a given molecular system as a collection of QHOs, integration over r leads to the following discretized form of Eq. ( 2) where T pq ¼ r r p r r q Wðr pq Þ is the dipole interaction tensor, r p and r q are the QHO positions, r pq ¼ jr p À r q j, and RÞ=r pq is the Coulomb potential for the interaction between two spherical Gaussian distributions.The solution of the SCS equation for every frequency of the electric field yields the molecular polarizability tensor SCS ði!Þ and atomic polarizability tensors SCS p ði!Þ, which now contain both short-range (via the TS-vdW scheme) and long-range (via the SCS equation) electrostatic screening.The width R of the Gaussian in Eq. ( 3) is derived from the dipole self-energy (i.e., the zerodistance limit of the classical dipole-dipole interaction) for a given frequency of the electric field [16,17] As an illustration of including the SCS effects, we begin by considering the static polarizability of the N 2 molecule.The N 2 polarizability is anisotropic, with experimental values being k ¼ 16:1 and ?¼ 9:8 bohr 3 . We first extended the TS-vdW scheme to treat atomic anisotropy by integrating the Hirshfeld volume in the x, y, and z directions, obtaining k ¼ 13:7 and ?¼ 11:9 bohr 3 ( iso ¼ 12:5 bohr 3 ) for N 2 .It is clear that the atomic Hirshfeld partition of the molecular electron density is not sufficient to reproduce the anisotropy in the molecular static polarizability.However, the anisotropy is greatly reduced when calculating the C 6 coefficients, because the response becomes increasingly isotropic for larger imaginary frequencies.In fact, this is the reason why the TS-vdW scheme yields significantly better results for the C 6 coefficients (5.5% error) than for the static polarizability ($ 12% error) [5].The SCS equation with the TS-vdW polarizability as input yields k ¼ 16:0 and ?¼ 8:7 bohr 3 ( iso ¼ 11:2 bohr 3 ) for N 2 , significantly increasing the anisotropy.We further calculated the polarizability tensor of 17 small molecules, ranging from H 2 to cyclohexane.The error in the fractional polarizability anisotropy (( ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ) was found to be 23% compared to experimental values [14].However, because the anisotropy is less important for the integrated vdW coefficients, this error is acceptable, as shown below by the performance of TS-vdW þ SCS for intermolecular C 6 coefficients.We note that the error in fractional anisotropy at the TS-vdW level is 80%.Thus, the screening clearly improves the physical description of the polarizability anisotropy.
The TS-vdW method has been shown to yield accurate isotropic C 6 coefficients (mean absolute error of 5.5%) for a large database of 1225 atomic and molecular dimers [5].The TS-vdW þ SCS method slightly increases the error to 6.3% due to linear alkane chains, where the anisotropy along the chain is overestimated.
The importance of long-range screening grows with the system size, becoming the greatest for solids.We illustrate this effect with the example of hydrogen-saturated silicon clusters of increasing size and also with the silicon bulk.The results are shown in Fig. 1 and compared to the timedependent local density approximation (TDLDA) calculations.We measured the accuracy of TDLDA using the experimentally derived dipole-oscillator strength C 6 coefficient for the SiH 4 molecule [18] (13% error) and the C Si-Si 6 coefficient in the silicon bulk (3% error) using the Clausius-Mossotti equation with the experimental dielectric function [19].Therefore, we deem the TDLDA C 6 coefficients as accurate references for the larger silicon clusters.For smaller clusters, the TS-vdW values are in good agreement with TDLDA and experiment.However, the error in the TS-vdW C 6 coefficients increases progressively with the cluster size, with an error of 27% for the largest cluster, Si 172 H 120 .The inclusion of SCS effects coefficient in comparison to the value derived from the experimental dielectric function, whereas the TS-vdW þ SCS method is in excellent agreement with experiment, with an error of just 8%.
We now extend the TS-vdW þ SCS method to include the fully nonadditive many-body vdW energy.For a system with N atoms, the many-body vdW energy will contain terms up to Nth order (i.e., 2-body, 3-body, . .., N-body).First, the fully screened static polarizability SCS p and the characteristic excitation frequency !SCS p are obtained for every atom in the system using the SCS procedure described above [20].The many-body vdW energy is computed using the coupled fluctuating dipole model (CFDM) [6,21] for a collection of coupled isotropic threedimensional QHOs that represent the atoms in the molecular system of interest.The CFDM Hamiltonian can be written as [21] H where p ¼ ffiffiffiffiffiffi ffi m p p p , with p being the displacement of oscillator p from equilibrium, and m p ¼ ð p ! 2 p Þ À1 .In the above equation, the first two terms correspond to the kinetic and potential energy for an individual QHO, while the third term corresponds to the dipole-dipole interaction between QHOs [T 0 pq ¼ r r p r r q W 0 ðr pq Þ, where W 0 ðr pq Þ will be defined below].The SCS superscript is ommited in the notation, but assumed for all and !.The bilinear form of the interaction term in the CFDM Hamiltonian allows for a numerically exact solution of the Schro ¨dinger equation by diagonalizing the 3N Â 3N interaction matrix.The vdW interaction energy between the QHOs is then computed as the difference between the zero-point energies of the coupled and uncoupled QHOs where p are the Hamiltonian eigenvalues.The above expression for the vdW energy can be derived from ACFDT, since in our case the polarizability is an even function of ! and has only simple poles on or below the real !axis.Under these circumstances the ACFDT correlation energy is given as the difference of zero-point energies [22,23].The many-body vdW energy is part of the long-range correlation energy, but in general the correlation energy includes other contributions.In order to have an electronic structure method that treats the full range of exchange and correlation effects, we need to couple the vdW energy to an approximate DFT functional.This coupling typically introduces empiricism into the approach.However, instead of the ad hoc damping functions used in interatomic pairwise dispersion-correction approaches, here the coupling is more naturally introduced through a modified Coulomb potential [24].We take an arbitrary, albeit one of the simplest, forms for the Coulomb potential where is a range-separation parameter and R vdW pq is the sum of the vdW radii for a pair of atoms p and q.The parameter controls how quickly the Coulomb potential reaches the 1=r pq asymptote in the long range.Physically, the dispersion energy is finite in the zero-distance limit [25]; however, approximate DFT functionals should already capture part of this energy.Thus we choose ! 2, which ensures that W 0 ðr pq Þ !0 when r pq !0. In Eq. ( 7) the vdW radius of an atom is defined as in the TS-vdW scheme, R vdW , where the free atom vdW radius is obtained from the coupled-cluster electron density.This accurate reference electron density is required for a consistent definition of the vdW radius with the polarizability [5].
The combination of the TS-vdW scheme with SCS and CFDM gives rise to an efficient method for obtaining the vdW energy, wherein short-and long-range screening effects are seamlessly included and the many-body vdW energy is computed to infinite order in the dipole approximation.In what follows we call this new scheme DFA þ MBD, where DFA is substituted by the name of the functional used to approximate the DFT energy, and MBD means many-body dispersion.We benchmarked the performance of DFA þ MBD on the S22 database of intermolecular interactions [26,27].All DFA calculations employ the FHI-AIMS code [28].The parameter in Eq. ( 7) is determined using the S22 database ( ¼ 2:56 for Perdew-Burke-Ernzerhof (PBE) [29], ¼ 2:53 for PBE0 [30]).The mean absolute relative error (MARE) with respect to the recent coupled-cluster method with singles, doubles, and perturbative triples [CCSD(T)] binding energies computed at the basis-set limit is shown in Fig. 2 for different state-of-the-art methods.The range of binding energies in the S22 database varies from 23 to 895 meV, thus the MARE is the most unambiguous measure of the performance.The MARE is decreased from 9.2% for PBE þ TS-vdW to 5.4% for PBE þ MBD.In particular, the MARE for vdW-bound systems is reduced by 10% compared to PBE þ TS-vdW.We note that both PBE þ MBD and PBE0 þ MBD methods yield remarkable accuracy for all types of intermolecular interactions in the S22 database.
In order to show that the improvement obtained with DFA þ MBD also holds outside of the S22 database, we calculated the relative energies of 27 conformers of alanine tetrapeptide (Ace-Ala 3 -NMe)-a fundamental biomolecular benchmark system [31,32].As a reference, we use recently calculated CCSD(T) relative energies [31,32].The PBE0 þ TS-vdW method yields a mean absolute error (MAE) of 0:52 kcal=mol compared to CCSD(T).However, PBE0 þ TS-vdW predicts the wrong conformation to be the most stable, and there are changes in the conformational hierarchy compared to CCSD(T).PBE0 þ MBD reduces the MAE to just 0:29 kcal=mol, predicting the right conformation to be the most stable one and leading to a conformational hierarchy in excellent agreement with CCSD(T).The differences in conformational energies between PBE0 þ TS-vdW and PBE0 þ MBD reach 1 kcal=mol, illustrating the importance of including SCS and the MBD energy contributions for reaching chemical accuracy in biomolecular simulations.
The screening effects and the many-body vdW energy become even more pronounced for extended systems.Here we use the benzene molecular crystal as an example, since the crystal geometry and the cohesive energy have been accurately measured [33,34], as well as recently computed using the exact-exchange plus RPA correlation energy method [35].We computed the cohesive energy using the PBE þ TS-vdW and PBE þ MBD methods at the experimental geometry.At the PBE þ TS-vdW level, the cohesive energy was computed as 690 meV=molecule, which is a significant overestimation compared to the experimental values of 518-560 meV=molecule (extrapolated to 0 K and with zero-point vibrational energy subtracted) [35].PBE þ MBD significantly improves the prediction, yielding 565 meV=molecule.Notice the large reduction of the binding energy that comes from the collective inclusion of both SCS and MB effects.
In conclusion, we have developed an efficient method for obtaining an accurate description of vdW interactions that includes both screening effects and treatment of the many-body vdW energy to infinite order.This method has a single physically motivated range-separation parameter for a given DFT functional and is applicable to finite-gap molecules and condensed matter systems.The calculation of the MBD energy takes % 3 min for a system of 1000 atoms on a single processor (DFT calculation time not included).Our method significantly improves the binding energies for small molecules compared to high-level quantum chemical methods, with pronounced improvements found for larger, more complex systems.The efficiency of our method paves the way to molecular dynamics simulations with the full many-body treatment of vdW interactions.Work to derive the self-consistent MBD potential and forces is in progress.
overall depolarization for the larger clusters, decreasing the error significantly in comparison to TDLDA.The depolarization is even larger for the Si bulk, in which the TS-vdW scheme yields an overestimation of 68% in the C Si-Si 6