A Variational Approach to London Dispersion Interactions without Density Distortion.

We introduce a class of variational wave functions that captures the long-range interaction between neutral systems (atoms and molecules) without changing the diagonal of the density matrix of each monomer. The corresponding energy optimization yields explicit expressions for the dispersion coefficients in terms of the ground-state pair densities of the isolated systems, providing a clean theoretical framework to build new approximations in several contexts. As the individual monomer densities are kept fixed, we can also unambiguously assess the effect of the density distortion on London dispersion interactions; for example, we obtain virtually exact dispersion coefficients between two hydrogen atoms up to C10 and relative errors below 0.2% in other simple cases.

Within the Born–Oppenheimer approximation, nonrelativistic quantum mechanics of Coulomb systems (atoms, molecules) predicts that everything is attracted to everything, because it is always possible to find a relative orientation that lowers the total energy of two neutral molecules separated by a large distance R by at least a term −C6R–6, with C6 > 0.[1,2] This universal attractive interaction between quantum Coulomb systems, often called London dispersion interaction, is very small compared to the energy scale of a typical chemical (covalent) bond. Yet, it is omnipresent in chemistry, biology, and physics, determining protein folding, the structure of DNA, and the physics of layered materials, just to name a few examples. An accurate, efficient, and fully nonempirical treatment of these forces remains an open challenge in many respects, and it is the focus of several ongoing research efforts.[3−9]

The physical origin of London dispersion interactions is often discussed in terms of “coupling between instantaneous dipoles”,[2,10] which are actually correlations in the interfragment part of the pair density (the diagonal of the two-body reduced density matrix), as illustrated, for example, in ref (11). The textbook treatment is based on perturbation theory,[2,10,12] in which excitations on both monomers are coupled, which allows rewriting the dispersion coefficients in terms of dynamical polarizabilities of the individual systems. In this Letter, inspired by the work of Lieb and Thirring (LT),[1] we attack the dispersion problem from a variational perspective, using only ground-state properties of the monomers. We construct a class of variational wave functions for the long-range interaction between neutral Coulomb systems that are explicitly forbidden to deform the diagonal of the full many-body spatial density matrix (and thus the density) of each individual monomer. This constraint has several advantages: First, it is conceptually appealing, as dispersion becomes a monomer–monomer interaction accompanied by a change in kinetic energy only, because all the remaining intra monomer interactions are kept unchanged. Second, it highly simplifies the optimization procedure, which can be implemented in a very efficient way and leads to expressions for the dispersion coefficients in terms of the ground-state pair densities of the individual monomers, giving a neat framework to build new approximations. For example, in a density functional theory (DFT) setting, one can use approximations for the exchange–correlation holes of the fragments, providing justification and a route for improving models such as the exchange-hole dipole moment (XDM).[13] It also provides a prescription to obtain atomic dispersion coefficients inside molecules and solids, which are a crucial ingredient in the construction of atomistic force fields and dispersion corrections to DFT calculations (DFT+D).

One could easily object that this construction cannot give the exact answer. In a widely quoted paper, Feynman[14] wrote that van der Waals’ forces arise from charge distributions with higher concentration between the nuclei, with a permanent dipole moment of order R–7 being induced on each atom, such that “each nucleus is attracted by the distorted charge distribution of its own electrons”.[14] Feynman’s conclusion was proven by using the electrostatic Helmann–Feynman theorem for both atoms and molecules[15] and has been confirmed with various accurate calculations.[16,17] Therefore, we know that a small density distortion must be there. Nonetheless, it has been observed several times that very accurate dispersion coefficients can be obtained with methods that describe very poorly this density distortion, including cases in which the electrostratic Helmann–Feynman theorem gives zero force at order R–7.[18] The issue has been often debated in the literature, and there is the general feeling that the density distortion, which actually requires expensive calculations to be accurately captured, should have very little influence on dispersion energetics, which is our goal here. Moreover, our construction provides (if combined with accurate pair densities for the monomers) a rigorous upper bound to the best possible variational interaction energy when the density distortion is explicitly forbidden, allowing us to assess its importance in an unambiguous manner, which is per se also conceptually interesting.

For illustrative purposes, we discuss first the case of two systems A and B, each with N = N = 1 electron, generalizing it right after to the many-electron case. The two atoms are separated by a large distance R, and our wave function readswhere r1 and r2 are electronic coordinates centered in nucleus A and nucleus B, respectively; Ψ0(r) is the ground-state wave function of system A (or B) alone, and J, with |J(r1, r2) | ≪ 1, preserves the unperturbed densities at all orders in R–1 via the constraintswhich are imposed with an expansion in terms of one-electron functions b(r) and b(r) and variational parameters cwhere(and similarly for b), with ρ0(r) = |Ψ0(r)|2. The f(r)’s can be chosen in different ways, with the choice influencing the convergence. Our wave function does not include antisymmetrization between the electron on A and the one on B, which affects the asymptotic (large R) interaction energy only with terms that vanish exponentially with R. This wave function correlates the positions of the two electrons in the two different atoms without changing their one-electron densities.

We now consider the full Hamiltonian of the problem, Ĥ = Ĥ(r1) + Ĥ(r2) + Ĥint(r1, r2), where Ĥ and Ĥ are the Hamiltonians of each isolated atom and Ĥint contains all the interactions between the particles in A (electron and nucleus) with those in B, and evaluate its expectation value on our wave function, focusing on the interaction energy Eint = E0 – E0 – E0. Thanks to the density constraint of eqs and 3, there is no change in the expectation value of the nuclear–electron potential, implying that the interaction energy is entirely determined by the change in kinetic energy T̂ with respect to the one of the isolated atomsplus the expectation value of Ĥint, which splits automatically into the sum of an electrostatic part Wint,0(R) and the correlation part Wint,c({c}, R),As usual,[2] we perform a multipolar expansion of Ĥintwhere d stands for dipole, q for quadrupole, o for octopole, etc. The variational problem then takes a simplified form, as detailed in the Supporting Information,[19] with each order R–2 in the energy having its optimal c = cR–. The physics described by this wave function can be understood by looking at the leading order: our wave function gives a ΔT({c}) which is positive (as it costs kinetic energy to correlate the electrons) and quadratic in the variational parameters c, while Wint, ({c}, R) is linear in the c, negative (as correlation occurs exactly to lower the interaction term), and has a prefactor R–3 from eq . This guarantees that there are always optimal variational parameters c, all proportional to R–3, which minimize Eint, whose minimum becomes the familiar −C6/R6. Notice that ΔT is quadratic in the c only by virtue of the density constraint. This structure repeats at each order in 1/R. In LT[1] and subsequent works in a similar spirit,[20,21] the energy cost to correlate the electrons, quadratic in the variational parameters, comes from the full Hamiltonian of the monomers, while with our construction we force it to be of kinetic energy origin only.

Solving the variational equations leads to explicit expressions for all the even dispersion coefficients. For example, C6 is given bywhere the vector w has contracted indices I = ijwith e = x, y, z, and h = (1, 1, −2) when the bond is along the z-axis. The matrix L has two contracted indices I = ij and K = klwithand similarly for B. Analogous expressions hold for all the even dispersion coefficients (see ref (19)).

We have started with a very simple choice for the f(r)’s (just powers of the distance r from the nucleus times spherical harmonics[19]), finding that the convergence is already fast, as shown in Figure for the case of C6 for two H atoms. The values of C6, C8, and C10 for the same case are reported in Table , where they are compared with reference data from ref (22). In the Supporting Information,[19] we also report even second-order coefficients up to C30, which also agree within numerical accuracy with those of ref (22). From Table we see that our wave function yields essentially the exact lowering in energy due to dispersion without any distortion of the density of each individual atom.

Figure 1

Relative error on C6 between two H atoms with increasing number of terms in eq relative to the reference result.[23]

Table 1

Dispersion Coefficients C6, C8, and C10 between Two H Atoms Computed with 30 Terms in the Variational Expression of eqs –15, Compared to the Reference Data of Thakkar[22] and Masili and Gentil[23]a

	Thakkar[22]	Masili et al.[23]	this work (30 terms)
C6	6.4990267054058405 × 100	6.4990267054058393 × 100	6.4990267054058393 × 100
C8	1.2439908358362235 × 102	n.a.	1.2439908358362234 × 102
C10	3.2858284149674217 × 103	n.a.	3.2858284149674217 × 103

See also the Supporting Information.[19]

We now turn to the many-electron case, by considering two systems with N and N electrons, for which our wave function readswhere indicates all the electronic coordinates (including spins) of the electrons in system A (or B) and, again, Ψ0 is the ground-state wave function of system A (or B) alone. We choose J such that the diagonal of the full many-body density matrix of each system is not allowed to deviate from its ground-state expression:which is enforced by choosing for J the same form as eq , i.e., by imposing thatwith ρ0 the ground-state one-electron densities of the two systems. Again, it is only the interfragment pair density and the off-diagonal elements of the first-order reduced density matrix that are allowed to change to lower the energy, implying that also in the many-electron case the interaction energy is given by eq , because the expectation value of the electron–electron interaction inside each monomer is kept unchanged by virtue of eqs and 18. The variational minimization of the interaction energy is similar to the N = 2 case. For example, C6 takes the same form as eq , where now w and L readwhere (with the same expressions for B)and where τ and S are defined in eqs and 15. We see that the dispersion coefficients C for the many-electron case become explicit functionals of the spin-summed ground-state pair densities P0(r1, r2) and P0(r1, r2) of the two monomers, C = C[P0, P0]. The variational solution for the optimal c can also be turned into a Sylvester equation, which can be solved in an efficient way.[19,24]

Being derived from a wave function, the expression for the interaction energy is variational (i.e., it provides a lower bound to the exact C6) as long as we use exact (or very accurate) pair densities for the monomers. As an example, with the same simple choice for the f(r) used for two H atoms, we have computed C6 for He–He using the accurate variational wave function for the He atom of Freund, Huxtable, and Morgan[25] (see also ref (26)). The convergence is again fast, very similar to the one of Figure , yielding C6 = 1.458440 with an error of 0.17% with respect to the accurate value 1.460978 of ref (27). For the He–H case, we have a similar error of 0.15% with respect to ref (27). Whether these residual errors are due to the reduced variational freedom of our wave function or to inaccuracies in the multiple moment integrals of the correlated He pair density of ref (25) is for now an open question, but we see that the results can be rather accurate, at a computational cost which is determined by the individual monomer calculations. This opens many possibilites, as our expressions can be coupled with any wave function method for the monomers, provided it gives access to the second-order reduced density matrices and their integrals with the f(r), whose optimal choice needs to be assessed in a systematic way.

We can of course use pair densities from lower-level theories, which opens another whole realm of possibilities, bearing in mind that the expressions will be no longer variational for the interaction energy. For example, in Figure we show the result for C6 for Ne–Ne using the Hartree–Fock (HF) pair density, including the basis set dependence. We see that the error is around 5%, as HF is a good approximation for the Ne atom, and that we are below the exact interaction energy. It is clear that a detailed study of the basis set dependence is also needed, which should be coupled to an optimal choice for the f(r). The use of Kohn–Sham orbitals from different approximate functionals to determine the exchange pair densities is also worth investigating.

Figure 2

Convergence of the C6 dispersion coefficient with the number of functions f(r) for Ne–Ne using Hartree–Fock pair densities of the individual atoms in different basis sets. HF calculations were carried out using Gaussian 16,[28] and the arbitrary order multipole moments were calculated using HORTON 3.[29] The accurate value (thick line) is from ref (30). Values for vdW-DF-04/10 and VV-09 are from ref (31), while the value for XDM was taken from ref (13). vdW-DF-09 was not included, because it is parametrized based on this dispersion coefficient.

Within DFT, one could make even simpler approximations for the exchange–correlation holes h of the monomers, transforming the dispersion coefficients into density functionalsFor example, the Becke–Roussel[32] exchange-hole model could be used to estimate higher moments of the exchange pair density,[33] in a spirit similar to the XDM model.[13] The main difference is that now the polarizability of the monomers is no longer needed, as the expression for the dispersion coefficients requires only the exchange–correlation holes.

An important ingredient in atomistic force fields models and DFT+D are the dispersion coefficients C between pairs of atoms a b, with a in system A and b in system B, which are usually computed from atomic polarizabilities. Our variational wave function can provide a framework to define and compute them. One should obtain atomic volumes Ω and Ω using one of the various definitions (e.g., Hirshfeld partitioning[34]) and then obtain the pair dispersion coefficients using, order by order, the energy density coming from our wave function with the c optimized for the AB interaction, restricting the integration on the atomic volumes. For example, for C6 (with J below being the optimized factor at leading order)and similarly for the higher-order coefficients.

In conclusion, the class of variational wave functions we have introduced here provides a neat theoretical framework to build new approximations to tackle dispersion interactions, which are reduced to a competition between kinetic energy and monomer–monomer interaction only. Although this class of variational wave functions violates the virial theorem (similarly to second-order perturbation theory), and although we know that for the exact wave function the density must be distorted, the results can be very accurate (or even exact), with a computational efficiency that is promising for the many-electron case. Moreover, we know what is being approximated, so that one could also devise a scheme in which the monomer (pair) densities are relaxed in a second step, or perturbatively, if this is needed. The realm of open possibilities is vast, ranging from the use of various wave function methods to compute the monomer pair densities or approximations thereof, to the definition and calculation of atomic-pair dispersion coefficients. The influence of the level of theory on the calculation of the monomer pair densities, the basis set, and different choices for the f(r) need all to be investigated. We have also shown that exact dispersion coefficients up to C10 between two hydrogen atoms can be obtained without introducing any density distortion, which is interesting by itself. The framework can be also generalized to many monomers and to describe the dispersive part of the interaction between charged systems.