Electronic charge density distortions due to dispersion: physically meaningful DMA multipoles for H₂, HeH, and He⋯He

Abstract

Feynman attributed long-range dispersion forces to the attraction of each nucleus to the local dipolar distortion of the electronic charge distribution. Here we take a step toward the first demonstration of Feynman's statement with full configuration-interaction wave functions. We have used Stone's distributed multipole analysis (DMA) to obtain the local multipoles in H₂ in the b³Σ⁺_u and X¹Σ⁺_g states and the local dipoles in HeH and He⋯He in their ground states. Except for the H₂ singlet, these states have repulsive potentials with shallow wells due to van der Waals dispersion. For H₂, the DMA dispersion dipole on each nucleus, computed ab initio with the d-aug-cc-pV6Z basis, shows excellent agreement with the sum of the R⁻⁷ and R⁻⁹ terms predicted by perturbation theory. The DMA dipoles of HeH and He⋯He also agree quite well with the prediction of perturbation theory. The signs and the R-dependence of the DMA dispersion dipoles are fully consistent with Feynman's statement. For H₂, we also find strong agreement between the results of perturbation theory and the dispersion terms in the DMA quadrupoles, DMA octopoles, DMA hexadecapoles, the total quadrupoles, and the total hexadecapoles. The dynamic correlation effects on the multipoles have physical meaning when computed with sufficiently large basis sets.

---

Introduction

Richard Feynman offered a simple physical explanation for dispersion forces, in the case of two S-state atoms interacting at long range, where exchange effects are negligible.¹ He reasoned that the identical distance dependence of the van der Waals dispersion force and the dispersion dipole suggests that the two are connected. He observed that a dipole moment proportional to R⁻⁷ develops at each of the nuclear centers due to dispersion. The electronic charge builds up between the nuclei. Feynman stated that the van der Waals dispersion force is not produced by interactions between the dispersion dipoles; instead, it is “the attraction of each nucleus for the distorted charge distribution of its own electrons that gives the attractive 1/R⁷ force.”¹

In this work, we have applied Stone's distributed multipole analysis (DMA)^2–5 to find the dispersion dipoles as the difference between the results from full configuration interaction (FCI) wave functions and the results from wave functions that do not include dynamic correlation. We make a connection between Feynman's explanation of dispersion and calculations at the FCI level for the first time. We also compare the DMA results with predictions from perturbation theory.

Hirschfelder and Eliason offered the first numerical test of Feynman's explanation of the van der Waals dispersion forces, which they termed Feynman's “conjecture.” They considered two H atoms initially in 1s states, interacting at long range.⁶ Using perturbation theory, they showed that Feynman's statement is correct to at least four figures, if exchange is neglected.⁶ At that time, existing ab initio methods were not sufficiently accurate to test the statement for heavier pairs.

Hunt observed two paradoxical aspects of Feynman's explanation.⁷ The van der Waals dispersion energy can be obtained entirely within linear response, but the dispersion dipoles depend on nonlinear response.^8–21 If at least one of the interacting molecules lacks a center of symmetry, then the dispersion dipole varies as R⁻⁶ in the separation between the molecules,¹⁶ but the van der Waals force varies as R⁻⁷, so Feynman's rationale does not hold in these cases.

Working within the polarization approximation, Hunt gave an analytical proof of Feynman's statement that resolves both issues. She expressed the dispersion energy in terms of the nonlocal polarizability density α(r, r′; iω) on each center and then differentiated the result to obtain the dispersion force.⁷ The necessary connection between linear and nonlinear response arises when α(r, r′; iω) is differentiated with respect to a nuclear coordinate, because the derivative depends on the hyperpolarizability density.^22,23

Hunt derived a generalization of the Thomas–Reiche–Kuhn summation rule,^24–26 which she used to show that the forces on the individual nuclei in a non-centrosymmetric molecule vary as R⁻⁶ (like the dispersion dipole in that case), but when the forces are summed over the nuclei, the R⁻⁶ term vanishes, leaving an R⁻⁷ dependence of the force on the center of mass. Thus, Hunt also proved the generalization of Feynman's explanation to molecules of arbitrary symmetry.⁷ This work shows very close agreement between the DMA dispersion dipoles obtained from ab initio calculations and the local dispersion dipoles from perturbation theory. The level of agreement indicates that the polarization approximation is successful for the local dipoles.

There have been no previous tests of Feynman's statement at full CI level, although relevant advances have been made by Allen and Tozer²⁷ and by Thonhauser et al.²⁸ Allen and Tozer have determined the dispersion forces due to interactions between two helium atoms, by differentiating the interaction energy at Bruckner coupled-cluster level.²⁷ They compared the results with the Coulomb forces calculated using Δρ(r), defined as the difference between the charge density of the interacting atoms and the isolated atomic charge densities. The forces computed in these two ways agreed well with each other, for He–He separations of 8.0 a.u., 8.5 a.u. and 9.0 a.u.²⁷ They differed from the force obtained with the long-range dispersion coefficients C_n (with even values of n from 6 to 16)^29–31 by ∼4.5–11%, with the largest differences at 8.0 a.u. in their work. Allen and Tozer derived a correlation potential for density functional theory, which yielded the forces with similar accuracy.²⁷

Thonhauser and coworkers derived a nonlocal correlation potential for density functional theory and used it in their calculations on multiple systems.²⁸ For Ar–Ar and Kr–Kr, they obtained a high level of agreement between the dispersion force calculated from the derivative of the energy and the Coulomb forces calculated using Δρ(r).²⁸ The calculated binding distances were ∼5% larger than the experimental values for Ar–Ar and Kr–Kr, and the calculated potential wells were too deep, by ∼10 meV for Ar–Ar and slightly less for Kr–Kr.²⁸

Neither of these references^27,28 provided a full test of Feynman's statement, because the forces were calculated with the changes in the total charge density of the atom pairs relative to the densities of the isolated atoms, rather than from the electron densities assigned to each nucleus. For He–He²⁷ and for Ar–Ar,²⁸ substantial depletions of the electron density were found along the bond, immediately adjacent to each nucleus, but with an increase in the electron density further along the bond. The distributed dipoles were not computed.

In a recent collaborative effort, Cheng et al. have developed a method to provide accurate real-space electron densities.³² They have applied the method to examine the forces on H nuclei in the ground ¹Σ⁺_g state of H₂, the Li nuclei in Li₂ in its ground state, and the N nuclei in N₂ in its ground state, but only in the vicinity of equilibrium in each case. For H₂ in the ¹Σ⁺_g state, Cheng et al. evaluated the forces between R = 1.2 a.u. and R = 1.6 a.u.³² In the current work on the states of H₂, we focus on interactions at distances between 10 a.u. and 22 a.u.

Odbadrakh and Jordan have examined the dispersion forces between coupled Drude oscillators and the dipoles induced in each oscillator.³³ For two 3D Drude oscillators, with the wave function correct to second order, they confirmed that the dispersion force varying at R⁻⁷ can be obtained from the dispersion dipole, due to special relationships among the polarizability α, the dipole–dipole–quadrupole hyperpolarizability B, and the dispersion energy coefficient C₆ that hold for Drude oscillators, but not necessarily for atoms or molecules.³³ They did not calculate the force directly from the change in charge density.

Results obtained by Kooi and Gori-Giorgi might appear to be inconsistent with Feynman's analysis. In 2019 and later work, they developed a fixed diagonal matrices formalism that gives highly accurate results for the dispersion energy coefficients C₆, C₈, and C₁₀ for H₂, yet predicts no change in the electron density at any order.^34–36 This feature is shared by the theory of dispersion forces presented by Eisenschitz and London in 1930.^37–39 Neither approach can show the connection between the dispersion forces and the dispersion dipoles. Commentary by Hirschfelder and Eliason⁶ suggests an explanation: A wave function correct to first order in a perturbation suffices to give the energy correct to second order, so a first-order theory can correctly capture the dispersion energy coefficients C₆, C₈, C₁₀, and even higher coefficients. But a wave function correct to first order yields no change in the electron density, since that arises only at second and higher order. The formalism developed in ref. 34 corresponds to a wave function that is correct to first order. We note that the approach used by Kooi and Gori-Giorgi is more efficient than calculations based on charge densities, in terms of providing numerical results for the dispersion energy.^34–36

Grimme et al. have commented that dispersion effects “are rooted in instantaneous electron correlations.”⁴⁰ We note that the electron correlations permit a static build-up of electron density in the region between the nuclei, and that the forces on the nuclei are obtained from the electronic charge density by Coulomb's law. Thus, the van der Waals forces on the nuclei derive entirely from classical electrostatic effects, due to the dispersion charge density.^1,7 The forces on the electrons vanish.

In this work, we examine results from Stone's distributed multipole analysis (DMA),^2–5 applied to the hydrogen molecule in the b³Σ⁺_u and X¹Σ⁺_g states, with the two nuclei taken as the multipole centers. For the b³Σ⁺_u state, we have found the correlation contributions to the distributed multipoles up through the hexadecapoles by subtracting the Hartree–Fock DMA values from the full configuration-interaction (FCI) DMA values. At long range, with a d-aug-cc-pV6Z basis,⁴¹ we find an excellent fit to the local dipoles D₇R⁻⁷ + D₉R⁻⁹ obtained from perturbation theory.^19–21 The higher multipoles fit the forms predicted with perturbation theory, as well.^11,15,19,21

For the X¹Σ⁺_g state, we could not obtain the dispersion effects by subtracting the Hartree–Fock results from the FCI results, because the Hartree–Fock wave function for this state does not dissociate properly into two neutral H atoms; instead, it includes ionic terms with both electrons localized on a single center. The Hartree–Fock energy is erroneously high, even for R as short as 5 a.u. We replaced the Hartree–Fock wave function with a CAS(2,2) wave function, which does not allow for dynamic correlation in H₂.⁴² Then we obtained the dispersion contributions to the energy and the DMA multipoles by subtracting the CAS(2,2) values from the FCI values. We have also determined the DMA dispersion dipoles in the ground states of HeH and He⋯He from ab initio calculations.

Stone has remarked that the DMA method “is not stable with respect to changes of the basis set” and that the distributed multipoles obtained with large basis sets “may not correspond to physical expectations.”⁴ In our calculations, we observed shifts in the FCI DMA multipoles and the Hartree–Fock DMA multipoles when we changed the basis set—we have even seen sign changes in several of the multipoles. However, we have found that the differences between the FCI and Hartree–Fock values are quite stable with respect to changes in the basis. We have found that physically meaningful correlation terms in the distributed dipoles are obtained with sufficiently large basis sets for H₂, HeH and He⋯He.

The high level of agreement between the ab initio DMA dipoles and the dipoles obtained from perturbation theory without exchange is significant: The results are consistent with Feynman's 1939 statement about the origin of long-range van der Waals forces.¹ Direct calculations of the dispersion forces from the dispersion charge densities are beyond the scope of the current work, but these calculations are in progress.

A complementary method in widespread use to treat van der Waals dispersion interactions is based on the exchange-hole model of Becke and Johnson.^43–46 Ángyán has provided a first-principles derivation of this model, based on the charge-density autocorrelation function, which gives the dispersion energy in terms of the exchange–correlation hole.⁴⁷ He has also identified the added assumptions that must be made in order to obtain the dispersion energy directly from the exchange hole.⁴⁷ Heßelmann has shown that the dispersion energy can be written in terms of the densities and exchange holes of the monomers.⁴⁸

Langbein has derived the dispersion energy by considering the fluctuating point multipoles on each center.⁴⁹ Hunt recast Langbein's approach in terms of nonlocal polarizability densities and used it when she proved Feynman's statement about the origin of dispersion forces.⁷ Reviews of the theory of van der Waals dispersion forces have been given by Dykstra and Lisy;⁵⁰ Clark, Politzer, and Murray;⁵¹ Sherrill and Merz;⁵² Stone;⁵³ and Ángyán, Dobson, Jansen, and Gould.⁵⁴

To our knowledge, the current work is the first examination of the local multipoles specifically due to dispersion effects for H–H pairs or other interacting atoms, with a full representation of correlation and exchange.

Atomic charges and distributed multipoles

The DMA method has been applied to a very wide variety of van der Waals complexes.^55–60 Faerman and Price found that atomic multipole moments obtained from DMA analyses for peptides and amides were “reasonably transferable to other molecules, provided that at least the directly bonded functional groups are the same.”⁶¹ Price, Harrison, and Guest provided an example of the application of the DMA method to a relatively large molecule with their ab initio work on the potential of an undecapeptide cyclosporin derivative.⁶² In a review of the DMA method, Buckingham, Fowler, and Hutson commented that this method accurately describes electrostatic interactions that are energetically accessible, as close as the “surface of a molecule as defined by the van der Waals radii. DMA also gives a detailed picture of the qualitative features of the charge distribution.”⁶³ The DMA method typically offers an advantage over the point-multipole expansion, because the DMA electrostatic potential converges up to the boundary formed by the van der Waals radii of the atoms,^2–5 including regions where the point-multipole series diverges. However, the DMA method offers little advantage over the point-multipole series, in terms of the range of convergence when there are just two centers,⁵ as in the current work.

Beginning with the work of Mulliken,⁶⁴ various methods have been proposed to assign partial charges to atomic centers. Mulliken's approach splits the contributions to the density matrix from orbitals on two different centers equally between the centers.⁶⁴ In contrast, Hirshfeld's approach allows for a weighted distribution of the charge between centers.^65,66 The method of iterated stockholder atoms (ISA) suggested by Lillestolen and Wheatley takes this a step further by iterating to self-consistency in the weighted assignment of partial charges.^67,68 The results have been analyzed for conformational stability⁶⁹ and for use in producing force fields;⁷⁰ they have also been adjusted by use of multipole constraints.⁷¹

Other methods include the natural atomic population analysis developed by Reed et al.⁷² as an atom-localized variant of the natural orbitals that diagonalize the full density matrix.⁷³ Bader has derived atomic charges by analyzing the electron density topologically, with a separation into regions bounded by zero flux in the gradient vector field of the charge density.⁷⁴ Other means of assigning atomic charges or atomic charge densities include the method of deformed atoms in molecules;^75–78 the introduction of fuzzy atoms,⁷⁹ an efficient partitioning that makes the atomic charges as small as possible while preserving the overall molecular multipoles,⁸⁰ and the use of point charges derived from the molecular electrostatic potential.^62,81–84

For H₂ if the total electronic charge is split between the two H centers, the atomic charges are zero. The distributed dipole and higher multipoles are non-zero, however. Misquitta, Stone, and Fazeli have adapted the ISA to obtain distributed multipoles,⁸⁵ with very useful results. Numerous other methods of determining the values of distributed multipoles have been suggested.^86–96 Glick et al. have used neural networks to obtain atomic multipoles quickly, with transferable values between molecules,⁹⁵ and Heindel et al. have used distributed multipoles in the computation of many-body interactions.⁹⁶ It is not yet known whether these methods yield dispersion multipoles that agree with the predictions of long-range perturbation theory. Additional methods allow for distributed polarizabilities, as well as distributed multipoles.^97–107 and for the effects of charge overlap.^108–113

Perturbation theory of the dispersion multipoles

In this section, we summarize the perturbation theory for the long-range dispersion multipoles of interacting atoms. We test those results in this work, by comparison with the ab initio values for the DMA dispersion multipoles.

Byers Brown and Whisnant published the first accurate expression for the leading term in the dipole moment of two unlike atoms in S states;^8,9 this term varies as R⁻⁷ in the internuclear separation.^8,9 They obtained the dispersion dipole as the sum of two integrals over frequencies, a single integral D^I₇ and a double integral D^II₇. Four atomic response properties appear in the integrals, α(iω), η(ω), Γ(u, v) and Λ(u, v). Of these properties, only α(iω) is well known in other contexts.

Hunt later developed an approximation for the dispersion dipole of interacting atoms in S states, as a function of the static polarizability and the static dipole–dipole–quadrupole hyper-polarizability B of each atom, together with the van der Waals interaction energy coefficient C₆.¹¹ Shortly afterward, Galatry and Gharbi^12–14 used the fluctuation-dissipation theorem^114,115 to derive the dispersion dipole as a single integral containing the polarizability α(iω) and the hyperpolarizability B(0, iω), both evaluated at imaginary frequencies.^12–14 Galatry and Hardisson extended the analysis to include molecules with lower symmetry.¹⁶ Craig and Thirunamachandran derived the dispersion dipole induced in coupled centrosymmetric systems, directly from perturbation theory.¹⁵ For atoms in S states, their result agrees with the expression found by Galatry and Gharbi.^12–14

Hunt and Bohr developed the first theory of the dispersion dipole that explicitly includes the changes in the correlations of quantum charge-density fluctuations due to an applied field F.¹⁷ Their approach is based on Langbein's theory of van der Waals interactions.⁴⁹ The analysis is summarized here, to show that the same approach also gives the distributed dispersion dipoles. By the fluctuation-dissipation theorem,^114,115 the correlations of the fluctuating electronic dipoles of an atom or molecule in the presence of an external field F are determined by the imaginary component of the dipole polarizability , according to¹⁷

Here, μ_α(ω) and μ_β(ω′) refer to the fluctuating dipoles on the same center. For the electronic dipoles, the intensity of the fluctuations at room temperature T is virtually identical to their intensity at zero Kelvin. For centrosymmetric species, there is no change in the correlations of the fluctuating dipoles to first-order in F. On the other hand, for molecules that lack a center of symmetry, the fluctuating dipole correlations are altered at first order in the applied field F. This gives rise to a term in the dispersion dipole that varies as R⁻⁶ in the intermolecular separation, as in the work of Galatry and Hardisson.¹⁶

The correlations of the fluctuating electronic dipole and quadrupole in the presence of the field F are determined by the imaginary component of the dipole–quadrupole polarizability ,¹⁷

For centrosymmetric species, the fluctuating dipole and quadrupole are uncorrelated in the absence of an applied field, but an applied field induces correlations that are linear in the field. This produces a contribution to the dispersion dipole that varies as R⁻⁷ in the intermolecular separation.

The analysis is presented in detail in ref. 17. Changes in the real and imaginary components of the susceptibilities are both included. So, the theory accounts for hyperpolarization of a molecule by the applied field acting together with the nonuniform fluctuating field due to the interaction partner, as well as the field-induced fluctuation correlations characterized by eqn (1) and (2). By use of complex contour integration, the integrals over the real frequencies of fluctuations from −∞ to ∞ are transformed into integrals along the positive imaginary-frequency axis. The dispersion dipole for a pair of atoms A and B, with vector R pointing from A to B along the z axis is^12,15,17

Note that ref. 17 uses the opposite convention for the direction of R. In this equation, B_αβ,γδ characterizes the change in the dipole–quadrupole polarizability A_β,γδ to first order in the applied field F_α. The frequency 0 is associated with the dipole operator μ_α and frequencies iω and −iω are associated with the dipole operator μ_β and the quadrupole operator Θ_γδ. For an isotropic system, the B tensor has the form¹¹

B_αβ,γδ(0, iω) = (1/4)B(0, iω)[3(δ_αγδ_βδ + δ_αδδ_βγ) – 2δ_αβδ_γδ]. (4)

So, the dispersion dipole to order R⁻⁷ is

This is identical to the result derived by Galatry and Gharbi^12–14 and by Craig and Thirunamachandran.¹⁵ Based on the work by Hunt and Bohr,¹⁷ if μ^AB_z,disp is positive, the sign of the pair dispersion dipole is A⁻B⁺.

The hyperpolarization of atom A can be viewed as localized to atom A. Likewise, the correlations of the fluctuating dipole and quadrupole induced by the field acting on atom A are localized to A. Therefore, we can obtain the dispersion dipole induced in atom A in the presence of atom B as

Fowler has carried the perturbation analysis to higher order to obtain the coefficient D₉ of R⁻⁹ in the series for the dispersion dipole;²⁰ for an AB pair, D^AB₉ is

where L(0, iω, −iω) is the isotropic dipole–quadrupole–octopole hyperpolarizability and C(iω) is the isotropic quadrupole polarizability, as defined by Fowler.²⁰ The quantity D^AB₉ can also be separated into terms that are localized on center A or on center B, D^A(B)₉ and D^B(A)₉, respectively.

For the quadrupole due to dispersion, which varies as R⁻⁶ in the separation of a pair of atoms, Hunt provided an approximation in terms of the static susceptibilities α and B, and C₆.¹¹ Craig and Thirunamachandran derived an expression for the dispersion quadrupole using second-order perturbation theory within the polarization approximation.¹⁵

The result given by Craig and Thirunamachandran appears to give the sum of the dispersion quadrupoles on each center, but not the total quadrupole. The total quadrupole includes contributions from the distributed dipoles, as well as from the local quadrupoles, because dipoles of opposite signs located at (0, 0, −R/2) and (0, 0 R/2) contribute to the quadrupole at the origin of an AB pair. Hence, both the local dipoles and the local quadrupoles contribute terms to the total quadrupole that vary as R⁻⁶ to leading order.

To order R⁻⁶, the local dispersion quadrupole Θ^A(B)_zz at center A in the presence of B is^15,19,21

where center A is located at (0, 0, −R/2). In the integrand in eqn (8), the frequency 0 in B^A(iω, 0) is associated with the quadrupole operator, whereas the frequency 0 in B^A(0, iω) in eqn (6) for the dispersion dipole is associated with a dipole operator.^15,17,19,21 The local quadrupole Θ^A(B)_zz is negative, because the electronic charge becomes elongated along the z axis due to dispersion.

To leading order, the total dispersion quadrupole of H₂ is given by^15,19,21

Distributed octopoles Ω^A(B)_zzz and Ω^B(A)_zzz are produced by hyperpolarization due to the field and field-gradient of the fluctuating dipole at the other center, and by fluctuation correlations induced by a nonuniform field. The distributed octopole Ω^A(B)_zzz at center A at (0, 0, −R/2) varies as R⁻⁷. It is given by an imaginary frequency integral that contains the susceptibilities L^A(iω, −iω, 0) and α^B(iω),^20,21

In the integrand, the frequency 0 in L^A(iω, −iω, 0) is associated with the octopole operator, while the frequencies iω and −iω are associated with the dipole and quadrupole operators.^20,21 The total octopole Ω_zzz vanishes, since dispersion interactions produce equal and opposite distributed octopoles in the pair.

The distributed hexadecapole scales as R⁻⁸ to leading order. It is produced by hyperpolarization at each center, and by the field-induced fluctuation correlations. We have not found a value for the coefficient of the leading term in the distributed hexadecapole. This complicates fitting the hexadecapole, compared with the moments of lower orders.

Computational methods

For the b³Σ⁺_u state of H₂, we have obtained the correlation effects on the energy and molecular properties by subtracting the Hartree–Fock values from the full configuration-interaction values,^116,117 obtained with Molpro.^118–120 The basis sets come from the basis set exchange for quantum chemistry¹²¹ (see also ref. 41). The basis sets are cc-pV6Z, aug-cc-pV5Z, d-aug-cc-pV5Z, aug-cc-pV6Z, and d-aug-cc-pV6Z.† ¹²¹ We used Stone's DMA analysis² as implemented in Molpro.^118–120 To identify the “long range” R values, we have examined the energies of H₂ and the charges associated with the basis functions on individual centers in the singlet and lowest triplet state at long range.

The aug-cc-pV6Z results for the interaction energy ΔE at FCI level are shown in Fig. 1. The cyan curve shows the dispersion energy computed with the C₆, C₈, and C₁₀ coefficients from the work of Kooi and Gori-Giorgi,³⁴

E_disp = –C₆R⁻⁶ − C₈R⁻⁸ – C₁₀R⁻¹⁰, (11)

Fig. 1 Interaction energy ΔE in the triplet state (blue), in the singlet state (red), and energy from perturbation theory (cyan) with the C₆, C₈, and C₁₀ coefficients of Kooi and Gori-Giorgi.³⁴

The highly accurate values³⁴ are C₆ = 6.4990267054058393 a.u., C₈ = 124.39908358362234 a.u., and C₁₀ = 3285.8284149674217 a.u. These results agree closely with Thakkar's values for C₆, C₈, and C₁₀,¹²² the value of C₆ given by Masili and Gentil,¹²³ results given to four figures by Maeder and Kutzelnigg,¹²⁴ and results given to eight figures by Bishop and Pipin.²⁹

Fig. S1 in the supplement shows a close-up of the FCI energies of the triplet and singlet states for internuclear distances R ranging from 10 a.u. to 22 a.u. Fig. S2 shows the net charge q that is localized within the orbitals on each H center for the triplet state, along with the overlap charge. For internuclear separations R larger than 10 a.u., the overlap charge and the charge on a single center are both nearly zero. Fig. S3 shows our results for the potential of the b³Σ⁺_u state, compared with the Kołos–Wolniewicz potential¹²⁵ and the potential obtained by Kurokawa et al.¹²⁶ Fig. S4 shows our potential for the singlet state with other potentials.^125–128 For values of R > 11–12 a.u., the FCI energies of the singlet and triplet are very similar, and the HF values for the triplet DMA dipoles, which include exchange energies, are small compared to the FCI values. On that basis, we consider R > 11–12 a.u. to be long range for the H₂ molecule in the X¹Σ⁺_g and b³Σ⁺_u states.

We investigated the dependence of the DMA multipoles on the basis set. Comparisons are provided in Fig. 2, which shows the DMA dispersion dipoles on the H center at (0, 0, −R/2), computed with various basis sets: cc-pV6Z, aug-cc-pV5Z, d-aug-cc-pV5Z, aug-cc-pV6Z, and d-aug-cc-pV6Z. The DMA dipole values generally approach each other as the basis set size increases, although the results at cc-pV6Z level clearly differ from the others, due to the specific exponents in the basis functions. Results from the augmented and doubly-augmented basis sets are similar. Results obtained with the d-aug-cc-pV5Z basis virtually superimpose on the results in the d-aug-cc-pV6Z basis, except at the shortest internuclear distances.

Fig. 2 DMA dispersion dipole of the triplet state at the H center located at (0, 0, −R/2). The DMA dispersion dipole on the H center at (0, 0, R/2) is equal in magnitude, but opposite in sign.

Results and discussion

In this section, we present results for the dispersion multipoles of H₂ in the b³Σ⁺_u and X¹Σ⁺_g states, and for the dispersion dipoles of HeH and He⋯He. We compare the ab initio results with the predictions from perturbation theory, where available.

Distributed multipoles of H₂ in the b³Σ⁺_u triplet state

We have evaluated the integral in eqn (6) for the dispersion dipole at hydrogen nucleus A in the presence of H atom B by 64-point Gaussian quadrature, using values of α(iω) and B(0, iω) tabulated by Bishop and Pipin.¹²⁹ The result is μ^H(H)_z,disp = −394.51R⁻⁷. We believe this to be the most accurate value currently available for the R⁻⁷ term in the DMA dipole for H₂ in the b³Σ⁺_u state. Fowler has found a very similar value, μ^H(H)_z,disp = −393.5R⁻⁷, with α(iω) and B(0, iω) calculated for the H atom in a 10s 8p 5d basis; he evaluated the integral by 16-point Gaussian quadrature.¹⁹ For comparison, a later calculation by Fowler with a 10s 8p 5d 4f basis for the H atom yielded D^H(H)₇ = −394.9 a.u., also with 16-point Gaussian quadrature.²⁰

Fowler has evaluated the coefficient of R⁻⁹ for a single H center in the 10s 8p 5d 4f basis, with the result D^H(H)₉ = −12803 a.u.²⁰ Using Mathematica,¹³⁰ we have found that the best fit to the d-aug-cc-pV6Z DMA dipoles gives D^H(H)₉ = −13 099 a.u., if D^H(H)₇ is fixed at D^H(H)₇ = −394.51 a.u. The difference between the two values of D^H(H)₉ is ∼2.26%.

Fig. 3 shows the results for the DMA dispersion dipole from the aug-cc-pV6Z basis and from the d-aug-cc-pV6Z basis, compared with the D^H(H)₇R⁻⁷ + D^H(H)₉R⁻⁹ dipole, calculated using our value for D^H(H)₇ and Fowler's value for D^H(H)₉.

Fig. 3 Long-range DMA dispersion dipole in red with the R⁻⁷ and R⁻⁹ terms from perturbation theory; results from the aug-cc-pV6Z basis (cyan) and the d-aug-cc-pV6Z basis (blue).

The results from the d-aug-cc-pV6Z basis show a remarkable level of agreement with the analytic form from perturbation theory, over the entire range of R values greater than or equal to 12 a.u. Interestingly, the aug-cc-pV6Z basis gives a better fit than the d-aug-cc-pV6Z basis, for the R⁻⁷ term in the dipole taken alone. At distances below 12 a.u., the fit is not as good for either basis, since exchange effects contribute to the dispersion dipoles when the internuclear separation is smaller.

Next, we consider the long-range dispersion quadrupole from eqn (8). Fowler's result is Θ^H(H)_zz = −52.2R⁻⁶ (ref. 19). The subscript disp is omitted here and below for simplicity. Fowler and Steiner evaluated the dispersion multipoles in a more extensive calculation, by expanding the wave function in complete set of atomic orbitals and then truncating the basis to 20 radial functions R_nl(r) for l from 0 to 3, with the associated spherical harmonics.²¹ Their result for the leading term in the local dispersion quadrupole is Θ^H(H)_zz = −52.300432R⁻⁶. We have not found a good fit of the ab initio values of Θ^H(H)_zz at long range with an R⁻⁶ term alone, but after adding an optimized term that varies as R⁻⁸, we have fit Θ^H(H)_zz well. The coefficient Q₈ for the R⁻⁸ term (obtained as the best fit with the R⁻⁶ term kept fixed at the value obtained by Fowler and Steiner²¹) is Q₈ = −6908.07 a.u. Fig. 4 shows the DMA dispersion quadrupoles obtained in the aug-cc-pV6Z basis and the d-aug-cc-pV6Z basis.

Fig. 4 Local dispersion quadrupole Θ^H(H)_zz of the H₂ triplet state in the aug-cc-pV6Z basis (cyan) and in the d-aug-cc-pV6Z basis (blue). The curve in red shows the function derived from long-range perturbation theory, including R⁻⁶ and R⁻⁸ terms.

The results for the DMA quadrupole obtained with the d-aug-cc-pV6Z basis are in excellent agreement with the long-range function from perturbation theory. The dispersion quadrupoles Θ^H(H)_zz and Θ^H(H)_zz are identical.

We have also examined the total dispersion quadrupole in the two basis sets. We have used the result of Fowler and Steiner²¹ to fix the coefficient of the R⁻⁶ term and then found the best fit coefficient for the R⁻⁸ contribution. That gives the dispersion quadrupole as Θ_zz = 684.42026R⁻⁶ + 13048.3R⁻⁸. In Fig. 5, the results from perturbation theory are compared with the ab initio results for the total dispersion quadrupole of the b³Σ⁺_u state of hydrogen in the aug-cc-pV6Z and d-aug-cc-pV6Z basis sets.

Fig. 5 Total dispersion quadrupole Θ_zz of the H₂ triplet state with respect to an origin at the center of symmetry, in the aug-cc-pV6Z basis (purple) and the d-aug-cc-pV6Z basis (blue). Long-range function (red), distributed dipole contribution (green).

As noted, both the DMA dispersion dipoles and the DMA dispersion quadrupoles contribute to the total dispersion quadrupole. The total quadrupole is heavily influenced by the dispersion dipoles, as shown in Fig. 5, where the points in green represent the contribution to Θ_zz from the local dispersion dipoles alone. Due to this contribution, the total quadrupole is opposite in sign to the distributed quadrupoles. For Θ_zz, the ratio of the distributed dipole contribution to the distributed quadrupole contribution increases monotonically from 4.2 to 8.6 over the range of R values from 15 a.u. to 22 a.u. The total dispersion quadrupole of H₂ in the triplet state in the d-aug-cc-pV6Z basis is in excellent agreement with the long-range form from perturbation theory.

Next, we consider the results for the local dispersion octopoles in the triplet state. In Fowler's 10s 8p 5d 4f basis, eqn (10) gives Ω^H(H)_zzz = −1090R⁻⁷ for the H center at (0, 0, −R/2).²⁰ The highly accurate value computed by Fowler and Steiner is Ω^H(H)_zzz = −1090.1104R⁻⁷ (ref. 21), with both values in a.u. The DMA dispersion octopoles at the two H centers are equal in magnitude but opposite in sign, and the total octopole is zero. We did not obtain a good fit to the ab initio values with the R⁻⁷ term alone. We added an R⁻⁹ term, and fit the DMA octopoles from R = 12 a.u. to R = 22 a.u. to the function

Ω^H(H)_zzz = −1090.1104R⁻⁷ + bR⁻⁹. (12)

The best fit value that we have obtained with Mathematica¹³⁰ is b = −169 880 a.u. Fig. 6 shows the ab initio points and the best fit function for the DMA octopole of H₂ in the triplet b³Σ⁺_u state.

Fig. 6 DMA dispersion octopole Ω^H(H)_zzz for the H center at (0, 0, −R/2) in the b³Σ⁺_u state. Ab initio results from the d-aug-cc-pV6Z basis (blue), results from the aug-cc-pV6Z basis (cyan), and the long-range form (red).

Based on the plot in Fig. 6, the DMA octopole appears to be physically meaningful, although it does show some sensitivity to the basis.

The DMA dispersion hexadecapole Φ^H(H)_zzzz scales as R⁻⁸ to leading order. We have not located a value for the coefficient of the leading term in Φ^H(H)_zzzz. However, we have determined the best fit of Φ^H(H)_zzzz to a form that includes both R⁻⁸ and R⁻¹⁰ terms using Mathematica.¹³⁰ The result is

Φ^H(H)_zzzz = −36563.7R⁻⁸ − 5.7294 × 10⁶R⁻¹⁰. (13)

This curve and the ab initio values are shown in Fig. 7. The DMA dispersion hexadecapoles on the two H centers are equal.

Fig. 7 DMA dispersion hexadecapole Φ^H(H)_zzzz of H₂ in the triplet state in the aug-cc-pV6Z basis (cyan), the d-aug-cc-pV6Z basis (blue), and the best long-range fit to the ab initio d-aug-cc-pV6Z results using R⁻⁸ and R⁻¹⁰ terms (red).

We have also examined the total dispersion hexadecapole. No long-range coefficients are known for the total dispersion hexadecapole. We have obtained an initial fitting function from a log–log plot of Φ_zzzzvs. R, from R = 10 a.u. to R = 20 a.u., using the d-aug-cc-pV6Z results. The results for larger R values were erratic. The slope of the log–log plot is −4.12771, and the intercept is 5.802265. We used this result as a starting point to find a good fit to a function with integer powers of R. Starting from Φ_zzzz = 331.049R^−4.12771 based on the log–log plot, by numerical experimentation we have fit Φ_zzzz to a function that has the form h₄R⁻⁴ + h₆R⁻⁶. We varied h₄; then with the h₄ value fixed, we found the value of h₆ that best fits the data.

The total hexadecapole Φ_zzzz = 220R⁻⁴ + 2898.09R⁻⁶ fits the data quite well, as shown in Fig. 8. The results for the total dispersion hexadecapole obtained with the aug-cc-pV6Z basis and with the d-aug-ccc-pV6Z basis are quite similar. As we found for the total dispersion quadrupole versus the distributed quadrupole, the sign of the total hexadecapole is opposite to that of the distributed hexadecapoles. Both the distributed quadrupoles and the total quadrupole fit functions containing R⁻⁶ and R⁻⁸ with appropriate coefficients. In contrast, the distributed hexadecapoles Φ^H(H)_zzzz fit eqn (13), which contains R⁻⁸ and R⁻¹⁰ terms, while the total hexadecapole falls off with an R dependence very close to R⁻⁴.

Fig. 8 Total dispersion hexadecapole Φ_zzzz of the triplet state, with respect to the center of symmetry as the origin. Results in the aug-cc-pV6Z basis (purple), the d-aug-cc-pV6Z basis (blue), and from the best long-range fit (red). Points in green show the contribution to the total dispersion hexadecapole that comes solely from the distributed dispersion dipoles.

The lower-order distributed multipoles contribute to the dispersion hexadecapole with respect to an origin at the center of symmetry. The contribution of the distributed dipoles is particularly significant, because these dipoles vary as R⁻⁷ to leading order and they are weighted by R³ in order to obtain the dispersion hexadecapole at the shifted origin. The distributed quadrupoles, which vary as R⁻⁶ to leading order, are weighted by R² to obtain the hexadecapole at the shifted origin. These terms most likely account for the R⁻⁴ dependence of Φ_zzzz to leading order. The contributions of the DMA dispersion dipoles to the total dispersion hexadecapole are shown in green in Fig. 8. They do not explain the magnitude of Φ_zzzz, but their contribution is clearly important.

Distributed dispersion multipoles of H₂ in the X¹Σ⁺_g state

As mentioned above, we have obtained the dispersion effects on H₂ in the ground singlet state by subtracting the CAS(2,2) values from the full CI results. The CAS(2,2) values for H₂ in the ground singlet state do not account for dynamical correlation,⁴² and therefore they do not include the dispersion contributions. In Fig. S5 of the SI, we compare the dispersion energy of the singlet with the dispersion energy computed from the dispersion coefficients C₆, C₈, and C₁₀ calculated by Bishop and Pipin.²⁹ The fit between the two results is excellent.

Fig. 9 shows the DMA dispersion dipole for the singlet state of H₂, for comparison with the perturbation theory result that was obtained for the triplet state.

Fig. 9 DMA dispersion dipole on the H center at (0, 0, −R/2) for the ground singlet state of H₂ in the d-aug-cc-pV6Z basis (blue). The ab initio results are compared with the form derived from perturbation theory for the triplet state, including both R⁻⁷ and R⁻⁹ terms (red curve). The fit is virtually identical in quality for the triplet and singlet states.

We have found a similar pattern for Θ^H(H)_zz, the DMA dispersion quadrupole in the ground singlet state, as shown in Fig. 10. The ab initio DMA dispersion quadrupole for the singlet state fits the Θ^H(H)_zz function obtained from perturbation theory for the triplet state very well. The DMA quadrupoles on the two centers are identical.

Fig. 10 DMA dispersion quadrupole on the H centers in the singlet state, in the d-aug-cc-pV6Z basis (blue), compared with the result from perturbation theory that fits the DMA dispersion quadrupoles in the triplet state, including R⁻⁶ and R⁻⁸ terms.

This pattern of a high-quality fit to the results obtained for triplet H₂ also holds for the total dispersion quadrupole of the singlet, as plotted in Fig. S6 in the SI. The results for Ω^H(H)_zzz and Φ^H(H)_zzzz are good, but not quite as close as we found for the lower-order multipoles. In both cases, the ab initio values lie slightly above the curve for the triplet state from perturbation theory, as shown in Fig. S7 and S8.

To illustrate the extent of agreement that we have found even for higher-order dispersion multipoles, we have plotted the total dispersion hexadecapole Φ_zzzz in Fig. 11. The high level of agreement reflects contributions to Φ^H(H)_zzzz from the DMA dipole μ^H(H)_z and DMA quadrupole Θ^H(H)_zz.

Fig. 11 Total dispersion hexadecapole for the singlet state of H₂, obtained as the difference between the FCI and CAS(2,2) values in the d-aug-cc-pV6Z basis (blue) and the long-range function for the triplet state (red).

Distributed dispersion dipoles in HeH

The extension of the results for H₂ to systems with more electrons poses a challenge. The number of configuration state functions that can be produced from n basis functions rises sharply when the number of electrons increases. When the FCI wave function is represented in terms of n orbitals for a molecule with N electrons and spin S, the number of configuration state functions D is given by¹³¹

D = (2S + 1)(n + 1)⁻¹C(n + 1, N/2 − S)C(n + 1, N/2 + S + 1). (14)

For the aug- and d-aug-cc-pV5Z and for the aug- and d-aug-cc-pV6Z basis sets used in this work, D is larger for HeH than for the H₂ states by factors of ∼100–200, increasing for the larger basis sets. The D values are tabulated in the SI.

As a test of our FCI wave function in the d-aug-cc-pV6Z basis we have compared the dispersion force derived from the ab initio calculation with the long-range dispersion force calculated with the dispersion energy coefficients provided by Bishop and Pipin for HeH,²⁹C₆ = 2.8213439 a.u., C₈ = 41.836374 a.u., and C₁₀ = 871.54066 a.u. We fit the ab initio energies to a fifth-order interpolation function and then differentiated to find the force. Fig. 12 shows that the results for the force drawing H toward He as obtained from the two calculations virtually superimpose.

Fig. 12 Dispersion force on the H nucleus in HeH. Ab initio results (blue), results from dispersion energy coefficients (red).

An additional numerical challenge arises for He-H because the DMA dipoles are smaller for systems that contain He than for systems that contain only H atoms. At ħω = 0.007 Hartree, B_zz,zz(0, iω) = −106.453978 a.u. for the H atom, while B_zz,zz(0, iω) = −7.326008 a.u. for the He atom.¹²⁹ At that same frequency, the polarizability of the H atom is 4.498717 a.u., while the polarizability of the He atom is 1.383118 a.u.¹³²

We have obtained a long-range fit to the DMA dipole at the H center in the d-aug-cc-pV6Z basis, starting with the values of D₇ and D₉ provided by Fowler,²⁰ and then adding an optimized D₁₁ term. This approach gives

μ^H(He)_z,disp = −153.4R⁻⁷ − 4220R⁻⁹ − 169 692R⁻¹¹. (15)

Fig. 13 shows the DMA dispersion dipole on H in HeH computed in this way (red curve) and the values obtained ab initio with various basis sets. The agreement is again excellent.

Fig. 13 DMA dispersion dipole for the H center in HeH from several different basis sets, and the prediction of long-range perturbation theory including R⁻⁷, R⁻⁹, and R⁻¹¹ terms (red).

Fowler²⁰ has given the DMA dispersion dipole for the He center in HeH as

μ^He(H)_z,disp = 35.65R⁻⁷ + 832.8R⁻⁹. (16)

This DMA dispersion dipole is positive, unlike the DMA dispersion dipole on H, but in both cases, charge is drawn toward the center of HeH. The ab initio results for μ^He(H)_z,disp do not fit eqn (16) as well the results for μ^H(He)_z,disp fit eqn (15), but we have found general agreement, as shown in Fig. S9.

Fig. 14 shows the results for the total dispersion dipole of HeH. The ab initio results agree very well with the long-range form obtained with the D₇ and D₉ coefficients given by Fowler,²⁰ with a D₁₁ term that we optimized by fitting the d-aug-cc-pV6Z results. The long-range form of the total dispersion dipole for HeH is

μ_z,disp = −117.75R⁻⁷ − 3387.2R⁻⁹ − 249 379R⁻¹¹. (17)

Fig. 14 Total dispersion dipole of in HeH from various basis sets, and the long-range form including R⁻⁷, R⁻⁹, and R⁻¹¹ terms.

Distributed dispersion dipoles in He⋯He

In the calculations on He⋯He, we were not able to work with the d-aug-cc-pV6Z basis. The number of configuration state functions (CSF) in this basis for He⋯He is 941 205 825 compared with only 52 975 for H₂ in the triplet state, from eqn (14).¹³¹ After taking symmetry in the D_2h group into account, the number of CSFs is still ∼117 650 280. Computer time limitations prevented our use of the d-aug-cc-pV6Z basis for He⋯He; three days were required for one iteration at one R value.

We have investigated the DMA dispersion dipole of He⋯He in the smaller basis sets, aug-cc-pV5Z, aug-cc-pV6Z, d-aug-cc-pV5Z, and t-aug-cc-pV5Z (which has 260 basis functions, so the work is tractable). The ab initio results for the DMA dispersion dipole on the He center at (0, 0, −R/2) are shown in Fig. 15 for comparison with the long-range form,

μ^He(He)_z,disp = −16.36R⁻⁷ – 278.2R⁻⁹ – 22341.5R⁻¹¹. (18)

Fig. 15 DMA dipole on the He center at (0, 0, −R/2) in He⋯He, compared with the long-range form including R⁻⁷, R⁻⁹, and R⁻¹¹ terms.

As above, this form was obtained with Fowler's results for D₇ and D₉ for the DMA dispersion dipole of He in He⋯He,²⁰ then adding an R⁻¹¹ term optimized to fit the results from the d-aug-cc-pV5Z basis. The results are plotted in Fig. 15. A close-up of the results from R = 12–20 a.u. is provided in Fig. S10, which shows a very good fit to the aug-cc-pV6Z results over this range.

Conclusions

The ab initio DMA dispersion dipoles for H₂ in the X¹Σ⁺_g and b³Σ⁺_u states agree remarkably well at long range with the known result from perturbation theory, as shown in Fig. 3 and 9. A single function containing the known R⁻⁷ and R⁻⁹ terms matches the ab initio DMA dipoles obtained in the d-aug-cc-pV6Z basis, for both the singlet and triplet states. The DMA dipole on the H center in HeH matches the long-range predictions very closely, as shown In Fig. 13. The ab initio values of the DMA dipole on the He center in HeH do not agree as well with the long-range predictions, as shown by Fig. S9 in the supplementary material, but we have still found a reasonable level of agreement. For He in He⋯He, we have found good agreement between the DMA dipoles calculated in all basis sets except for the t-aug-cc-pV5Z basis, as shown in Fig. 15. The close-up of the results in the range from R = 12 a.u to R = 20 a.u. in Fig. S10 shows that the best fit is obtained with the aug-cc-pV6Z basis. Generally, dependence on the choice of basis is detectable, but it is not large, except for the cc-pV6Z basis set for H₂ in the b³Σ⁺_u state and the t-aug-cc-pV5Z basis for He⋯He. We conclude that the DMA dispersion dipoles have physical meaning in themselves for H₂, He, and He⋯He.

We have also analyzed the correlation contributions to the DMA dispersion quadrupoles, octopoles, and hexadecapoles of H₂ in the X¹Σ⁺_g and b³Σ⁺_u states. Good matches of the ab initio results to the sums of the leading long-range terms from perturbation theory and an estimated term of the next higher order are shown in Fig. 4 and 10 for the quadrupoles and in Fig. 6 and S7 for the octopoles. The correlation contributions to the DMA dispersion hexadecapole are not currently known from perturbation theory, but the ab initio results are fit well by the sum of the terms of the two leading orders suggested by the perturbation analysis, as shown in Fig. 7. Overall, the fits of the DMA dispersion dipoles, quadrupoles, octopoles, and hexa-decapoles to the expected long-range functional forms for H₂ are striking.

The DMA dispersion quadrupoles and hexadecapoles are opposite in sign to the total properties computed with respect to the center of the H₂ molecule as an origin, as shown in Fig. 4, 5, 7, 8, 10 and 11 and Fig. S6, and S8. The difference arises because the DMA dipoles contribute significantly to the total quadrupoles and hexadecapoles. For the total quadrupole, the ratio of the contributions from the DMA dipoles to the contributions from the DMA quadrupoles ranges from 3.638 to 8.833 (in absolute value), for internuclear separations between 10 a.u. and 22 a.u. For the total hexadecapole, the relative contributions from the DMA dipoles are even greater, ranging from 18.29 to 5217 times the contribution from the DMA hexadecapoles. The total quadrupoles and hexadecapoles match the expected long-range forms quite well, as shown in Fig. 5, 8 and 11 and Fig. S6.

The energies of the ground singlet state and the lowest triplet state of H₂ obtained from full CI calculations agree very well at long range, as shown in Fig. 1 and Fig. S1. The energies also agree very well with the long-range dispersion energy calculated using E_disp = −C₆R⁻⁶ − C₈R⁻⁸ – C₁₀R⁻¹⁰, with the coefficients reported by Kooi and Gori-Giorgi.³⁴ On the scales of Fig. 1 and Fig. S1, the dispersion energies computed with the coefficients obtained earlier by Maeder and Kutzelnigg,¹²⁴ by Bishop and Pipin,²⁹ and by Thakkar¹²² are indistinguishable from the curves shown.

This work has shown that van der Waals dispersion effects produce dipolar distortions of the electronic charge distribution toward the center of the H₂ molecule at long range in both the singlet and the triplet states, and toward the centers of HeH and He⋯He in their ground states. This result is fully consistent with Feynman's statement about the connection of dispersion forces to the local dispersion dipoles.¹ Our ab initio results have been obtained from large-basis calculations that include the antisymmetrization of the wave function and the full exchange effects. To our knowledge, this is the first time that agreement between the ab initio DMA dispersion dipoles and the results from perturbation theory has been demonstrated.

Author contributions

Both Nathan D. Jansen and Katharine L. C Hunt: conceptualization, methodology, software, validation, investigation, data curation, and writing – review and editing. Katharine L. C. Hunt: writing – original draft, visualization, supervision, project administration, and funding acquisition. Hua-Kuang Lee: investigation of HeH.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included in the article and in the supplementary information (SI). Supplementary information: tables of the energies, DMA dipoles, quadrupoles, octopoles, and hexadecapoles on the nucleus at (0, 0, −R/2), and the total quadrupoles and total hexadecapoles for both the triplet and singlet states of H₂. Results for the DMA dipoles of HeH and He⋯He are tabulated, together with the total dipole of HeH. The DMA dipoles in the triplet state obtained with various basis sets are listed, as well as the charges in the orbitals centered on the nuclei and the overlap charges. Results from eqn (14) for the number of configuration state functions in various basis sets used in the work are also provided. The SI includes plots of the singlet and triplet energies at long range, the charges as noted above, and the energies of the H₂ triplet and singlet states obtained in this work, compared with results from ref. 125–128. Plots showing the dispersion energy, the total dispersion quadrupole, the DMA dispersion octopole, and DMA dispersion hexadecapole for H₂ in the singlet state are also provided. For HeH, the DMA dispersion dipole on the He center is plotted, and for He⋯He, a close-up of the DMA dispersion dipole at long range is shown.

Acknowledgements

This work has been supported in part by US National Science Foundation Grant CHE-2154028.

References

1. R. P. Feynman, Phys. Rev., 1939, 56, 340.
2. A. J. Stone, Chem. Phys. Lett., 1981, 83, 233.
3. A. J. Stone and M. Alderton, Mol. Phys., 1985, 56, 1047.
4. A. J. Stone, J. Chem. Theory Comp., 2005, 1, 1128.
5. A. J. Stone, Physical basis of intermolecular interactions, in Non-Covalent Interactions in Quantum Chemistry and Physics: Theory and Applications, ed. A. O. de la Roza and Gino A. DiLabio, Elsevier, Amsterdam, Netherlands, 2017, pp. 3–26.
6. J. O. Hirschfelder and M. A. Eliason, J. Chem. Phys., 1967, 47, 1164.
7. K. L. C. Hunt, J. Chem. Phys., 1990, 92, 1180.
8. W. Byers Brown and D. M. Whisnant, Chem. Phys. Lett., 1970, 7, 329.
9. W. Byers Brown and D. M. Whisnant, Mol. Phys., 1973, 25, 1385.
10. D. M. Whisnant and W. Byers Brown, Mol. Phys., 1973, 26, 1105.
11. K. L. C. Hunt, Chem. Phys. Lett., 1980, 70, 336.
12. L. Galatry and T. Gharbi, Chem. Phys. Lett., 1980, 75, 427.
13. L. Galatry and T. Gharbi, C. R. Seances Acad. Sci., Ser. B, 1980, 290, 401.
14. L. Galatry, C. R. Seances Acad. Sci., Ser. B, 1980, 291, 113.
15. D. P. Craig and T. Thirunamachandran, Chem. Phys. Lett., 1981, 80, 14.
16. L. Galatry and A. Hardisson, J. Chem. Phys., 1983, 79, 158.
17. K. L. C. Hunt and J. E. Bohr, J. Chem. Phys., 1985, 83, 5198.
18. B. Linder and R. A. Kromhout, J. Chem. Phys., 1986, 84, 2753.
19. P. W. Fowler, Chem. Phys., 1990, 143, 447.
20. P. W. Fowler, Chem. Phys. Lett., 1990, 171, 277.
21. P. W. Fowler and E. Steiner, Mol. Phys., 1990, 70, 377.
22. K. L. C. Hunt, J. Chem. Phys., 1989, 90, 4909.
23. K. L. C. Hunt, Y. Q. Liang, R. Nimalakirthi and R. A. Harris, J. Chem. Phys., 1989, 91, 5251.
24. W. Thomas, Naturwissenschaften, 1925, 13, 627.
25. F. Reiche and W. Thomas, Z. Phys., 1925, 34, 510.
26. W. Kuhn, Z. Phys., 1925, 33, 408.
27. M. J. Allen and D. J. Tozer, J. Chem. Phys., 2002, 117, 11113.
28. T. Thonhauser, V. R. Cooper, S. Li, A. Puzder, P. Hyldgaard and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 125112.
29. D. M. Bishop and J. Pipin, Int. J. Quantum Chem., 1993, 45, 349.
30. A. J. Thakkar, J. Chem. Phys., 1988, 89, 2092.
31. T. Korona, H. L. Williams, R. Bukowski, B. Jeziorski and K. Szalewicz, J. Chem. Phys., 1997, 106, 5109.
32. L. Cheng, P. B. Szabó, Z. Schätzle, D. P. Kooi, J. Köhler, K. J. H. Giesbertz, F. Noé, J. Hermann, P. Gori-Giorgi and A. Foster, J. Chem. Phys., 2025, 162, 034120.
33. T. T. Odbadrakh and K. D. Jordan, J. Chem. Phys., 2016, 144, 034111.
34. D. P. Kooi and P. Gori-Giorgi, J. Phys. Chem. Lett., 2019, 10, 1537.
35. D. P. Kooi and P. Gori-Giorgi, Faraday Discuss., 2020, 224, 145.
36. D. P. Kooi, T. Weckman and P. Gori-Giorgi, J. Chem. Theory Comp., 2021, 17, 2283.
37. R. Eisenschitz and F. London, Z. Phys., 1930, 60, 491.
38. F. London, Z. Phys., 1930, 63, 245.
39. F. London, Trans. Faraday Soc., 1937, 33, 8.
40. S. Grimme, A. Hansen, J. G. Brandenburg and C. Bannwarth, Chem. Rev., 2016, 116, 5105.
41. D. E. Woon and T. H. Dunning, Jr., J. Chem. Phys., 1994, 100, 2975.
42. E. Ramos-Cordoba, P. Salvador and E. Matito, Phys. Chem. Chem. Phys., 2016, 18, 24015.
43. D. Becke and E. R. Johnson, J. Chem. Phys., 2005, 122, 154104.
44. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2005, 123, 154101.
45. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2006, 124, 014104.
46. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2007, 127, 154108.
47. J. G. Ángyán, J. Chem. Phys., 2007, 127, 024108.
48. A. Heßelmann, J. Chem. Phys., 2009, 130, 084104.
49. D. Langbein, Theory of van der Waals Attraction, Springer, New York, 1974.
50. C. E. Dykstra and J. M. Lisy, J. Mol. Structure: THEOCHEM, 2000, 500, 375.
51. T. Clark, P. Politzer and J. S. Murray, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2015, 5, 169.
52. C. D. Sherrill and K. M. Merz, Quantum mechanical methods for quantifying and analyzing non-covalent interactions and for force-field development, in Many-Body Effects and Electrostatics in Biomolecules, ed. Q. Cui, M. Meuwly and P. Ren, Pan Stanford Publishing PTE Ltd., Singapore, 2016, pp. 65–120.
53. A. J. Stone, Physical basis of intermolecular interactions, in Non-Covalent Interactions in Quantum Chemistry and Physics: Theory and Applications, ed. A. O. DeLaRoza and G. A. DiLabio, Elsevier Science, Amsterdam, the Netherlands, 2017, pp. 3–26.
54. J. Angyán, J. Dobson, G. Jansen and T. Gould, Intermolecular perturbation theory, in London Dispersion Forces in Molecules, Solids, and Nano-Structures: An Introduction to Physical Models and Computational Methods, ed. A. Angyán, J. Dobson, G. Jansen, and T. Gould, Royal Society of Chemistry, Cambridge, England, 2020.
55. S. L. Price and A. J. Stone, J. Chem. Phys., 1987, 86, 2859.
56. P. W. Fowler and A. J. Stone, J. Phys. Chem., 1987, 91, 509.
57. P. L. A. Popelier, A. J. Stone and D. J. Wales, Faraday Discuss., 1994, 97, 243.
58. D. J. Wales, P. L. A. Popelier and A. J. Stone, J. Chem. Phys., 1995, 102, 5551.
59. D. J. Wales, A. J. Stone and P. L. A. Popelier, Chem. Phys. Lett., 1995, 240, 89.
60. G. J. B. Hurst, P. W. Fowler, A. J. Stone and A. D. Buckingham, Int. J. Quant. Chem., 1986, 29, 1223.
61. C. H. Faerman and S. L. Price, J. Am. Chem. Soc., 1990, 112, 4915.
62. S. L. Price, R. J. Harrison and M. F. Guest, J. Comput. Chem., 1989, 10, 552.
63. A. D. Buckingham, P. W. Fowler and J. M. Hutson, Chem. Rev., 1988, 88, 815.
64. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833.
65. F. L. Hirshfeld, Theor. Chim. Acta, 1977, 44, 129.
66. J. F. Harrison, Mol. Phys., 2005, 103, 1099.
67. T. C. Lillestolen and R. J. Wheatley, Chem. Commun., 2008, 5909.
68. T. C. Lillestolen and R. J. Wheatley, J. Chem. Phys., 2009, 131, 1441001.
69. T. Verstraelen, P. W. Ayers, V. Van Speybroeck and M. Waroquier, Chem. Phys. Lett., 2012, 545, 138.
70. T. Verstraelen, S. Vandenbrande, F. Heidar-Zadeh, L. Vanduyfhuys, V. Van Speybroeck, M. Waroquier and P. W. Ayers, J. Chem. Theory Comp., 2016, 12, 3894.
71. J. E. S. Mikkelsen and F. Jensen, J. Chem. Theory Comp., 2025, 21, 1179.
72. Alan E. Reed, Robert B. Weinstock and Frank Weinhold, J. Chem. Phys., 1985, 83, 735.
73. P.-O. Löwdin, Phys. Rev., 1955, 97, 1474.
74. R. F. W. Bader, Chem. Rev., 1991, 91, 893.
75. T. Arai, J. Chem. Phys., 1957, 26, 435.
76. J. Fernández Rico, R. López, I. Ema and G. Ramírez, J. Comput. Chem., 2004, 25, 1347.
77. J. Fernández Rico, R. López, I. Ema, G. Ramírez and E. Ludeña, J. Comput. Chem., 2004, 25, 1355.
78. R. Lopez, F. Martinez, I. Ema, J. M. G. de la Varga and G. Ramirez, Computation, 2019, 7, 64.
79. I. Mayer and P. Salvador, Chem. Phys. Lett., 2004, 383, 368.
80. M. Swart, P. Th Van Duijnen and J. G. Snijders, J. Comput. Chem., 2001, 22, 79.
81. N. K. Ray, M. Shibata, G. Bolis and R. Rein, Int. J. Quantum Chem., 1984, 27, 427.
82. U. C. Singh and P. A. Kollman, J. Comput. Chem., 1984, 5, 129.
83. S. J. Weiner, P. A. Kollman, D. A. Case, U. C. Singh, C. Ghio, G. Alagona, S. Profeta, Jr. and P. Weiner, J. Am. Chem. Soc., 1984, 106, 765.
84. C. Chipot, J. G. Ángyán and C. Millot, Mol. Phys., 1998, 94, 881.
85. A. J. Misquitta, A. J. Stone and F. Fazeli, J. Chem. Theory Comp., 2014, 10, 5405.
86. W. A. Sokalski, M. Shibata, R. Rein and R. L. Ornstein, J. Comput. Chem., 1992, 13, 883.
87. C. E. Whitehead, C. M. Brenerman, N. Sukumar and M. D. Ryan, J. Comput. Chem., 2003, 24, 512.
88. P. L. A. Popelier and M. Rafat, Chem. Phys. Lett., 2003, 376, 148.
89. E. V. Tsiper and K. Burke, J. Chem. Phys., 2004, 120, 1153.
90. J. F. Harrison, J. Phys. Chem. A, 2005, 109, 5492.
91. J. Pilmé and J.-P. Piquemal, J. Comput. Chem., 2008, 29, 1440.
92. A. Gramada and P. E. Bourne, Phys. Rev. E, 2008, 78, 0666601.
93. T. Bereau and M. Meuwly, Multipolar force fields for atomistic simulations, in Many-Body Effects and Electrostatics in Biomolecules, ed. Q. Cui, M. Meuwly, and P. Ren, Pan Stanford Publishing PTE Ltd., Singapore, 2016, pp. 233–268.
94. O. Loboda and C. Millot, J. Chem. Phys., 2017, 147, 161718.
95. Z. L. Glick, A. Koutsoukas, D. L. Cheney and C. D. Sherrill, J. Chem. Phys., 2021, 154, 224103.
96. J. P. Heindel, S. Sami and T. Head-Gordon, J. Chem. Theory Comp., 2024, 20, 8594.
97. U. Dinur, J. Comput. Chem., 1991, 12, 91.
98. A. J. Misquitta and A. J. Stone, J. Chem. Phys., 2006, 124, 024111.
99. D. Elking, T. Darden and R. J. Woods, J. Comput. Chem., 2007, 28, 1261.
100. A. Defusco, D. P. Schofield and K. D. Jordan, Mol. Phys., 2007, 105, 2681.
101. A. J. Stone and A. J. Misquitta, Int. Rev. Phys. Chem., 2007, 26, 193.
102. S. L. Price, M. Leslie, G. W. A. Welch, M. Habgood, L. S. Price, P. G. Karamertzanis and G. M. Day, Phys. Chem. Chem. Phys., 2010, 12, 8478.
103. R. Kumar, F. F. Wang, G. R. Jenness and K. D. Jordan, J. Chem. Phys., 2010, 132, 014309.
104. T. Jankowski, K. Wolinski and P. Pulay, Chem. Phys. Lett., 2012, 530, 1.
105. H. W. Kim and Y. M. Rhee, J. Comput. Chem., 1662, 2012, 33.
106. F. Rob and K. Szalewicz, Molec. Phys., 2013, 111, 1430.
107. A. J. Misquitta and A. J. Stone, Theor. Chem. Acc., 2018, 137, 153.
108. A. S. Werneck, T. M. Rocha Filho and L. E. Dardenne, J. Phys. Chem. A, 2008, 112, 268.
109. Q. T. Wang, J. A. Rackers, C. He, R. Qi, C. Narth, L. Lagardere, N. Gresh, J. W. Ponder, J.-P. Piquemal and P. Y. Ren, J. Chem. Theory Comp., 2015, 11, 2609.
110. A. Öhrn, J. M. Hermida-Ramon and G. Karlström, J. Chem. Theory Comp., 2016, 12, 2298.
111. S. A. Bojarowski, P. Kumar and P. M. Dominiak, Chem. Phys. Chem., 2016, 17, 2455.
112. S. A. Bojarowski, P. Kumar and P. M. Dominiak, Acta Cryst. B, 2017, 73, 589.
113. F. Jiménez-Grávalos and D. Suárez, J. Chem. Theory Comp., 2021, 17, 4981.
114. H. B. Callen and T. A. Welton, Phys. Rev., 1951, 82, 296.
115. H. B. Callen and T. A. Welton, Phys. Rev., 1951, 83, 34.
116. P. J. Knowles and N. C. Handy, Chem. Phys. Lett., 1984, 111, 315.
117. P. J. Knowles and N. C. Handy, Comput. Phys. Commun., 1989, 54, 75.
118. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, WIREs Comput. Mol. Sci., 2012, 2, 242.
119. H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, T. F. Miller, III, A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut and M. Sibaev, J. Chem. Phys., 2020, 152, 144107.
120. H.-J. Werner, P. J. Knowles, P. Celani, W. Györffy, A. Hesselmann, D. Kats, G. Knizia, A. Köhn, T. Korona, D. Kreplin, R. Lindh, Q. Ma, F. R. Manby, A. Mitrushenkov, G. Rauhut, M. Schütz, K. R. Shamasundar, T. B. Adler, R. D. Amos, J. Baker, S. J. Bennie, A. Bernhaardsson, A. Berning, J. A. Black, P. J. Bygrave, R. Cimiraglia, D. L. Cooper, D. Coughtrie, M. J. O. Deegan, A. J. Dobbyn, K. Doll, M. Dornbach, F. Eckert, S. Erfort, E. Goll, C. Hampel, G. Hetzer, J. G. Hill, M. Hodges, T. May, B. Mussard, S. J. McNicholas, W. Meyer, T. F. Miller, III, M. E. Mura, A. Nicklass, D. P. O’Neill, P. Palmeiri, D. Peng, K. A. Peterson, K. Pflüger, R. Pitzer, I. Polyak, P. Pulay, M. Reiher, J. O. Richardson, J. B. Robinson, B. Schröder, M. Schwilk, T. Shiozaki, M. Sibaev, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, J. Toulouse, M. Wang, M. Welborn and B. Ziegler, MOLPRO, version 2022.3, a package of ab initio programs, https://www.molpro.net.
121. https://www.basissetexchange.org/ .
122. A. J. Thakkar, J. Chem. Phys., 1988, 89, 2092.
123. M. Masili and R. J. Gentil, Phys. Rev. A, 2008, 78, 034701.
124. F. Maeder and W. Kutzelnigg, Chem. Phys., 1979, 42, 95.
125. W. Kołos and L. Wolniewicz, J. Chem. Phys., 1965, 43, 2429.
126. Y. I. Kurokawa, H. Nakashima and H. Nakatsuji, Phys. Chem. Chem. Phys., 2019, 21, 6327.
127. W. Kołos and L. Wolniewicz, J. Mol. Spectrosc., 1975, 54, 303.
128. K. Pachucki, Phys. Rev. A, 2010, 82, 032509.
129. D. M. Bishop and J. Pipin, J. Chem. Phys., 1993, 98, 4003.
130. Wolfram Research, Inc. Mathematica, Version 14.2, Champaign, IL, 2024.
131. J. Paldus, Phys. Rev. A, 1976, 14, 1620.
132. D. M. Bishop and J. Pipin, J. Chem. Phys., 1992, 97, 3375.

---

Footnote

† Basis sets used in this work:¹²¹, cc-pV6Z, [10s 5p 4d 3f 2g 1h/6s 5p 4d 3f 2g 1h], with 182 contracted functions, aug-cc-pV5Z, [9s 5p 4d 3f 2g/6s 5p 4d 3f 2g] with 160 contracted functions, d-aug-cc-pV5Z, [10s 6p 5d 4f 3g/7s 6p 5d 4f 3g] with 210 contracted functions, aug-cc-pV6Z, [11s 6p 5d 4f 3g 2h/7s 6p 5d 4f 3g 2h] with 254 contracted functions, and d-aug-cc-pV6Z, [12s 7p 6d 5f 4g 3h/8s 7p 6d 5f 4g 3h] with 326 contracted functions.¹²¹

---