Universal Pairwise Interatomic van der Waals Potentials Based on Quantum Drude Oscillators

Repulsive short-range and attractive long-range van der Waals (vdW) forces play an appreciable role in the behavior of extended molecular systems. When using empirical force fields, the most popular computational methods applied to such systems, vdW forces are typically described by Lennard-Jones-like potentials, which unfortunately have a limited predictive power. Here, we present a universal parameterization of a quantum-mechanical vdW potential, which requires only two free-atom properties—the static dipole polarizability α1 and the dipole–dipole C6 dispersion coefficient. This is achieved by deriving the functional form of the potential from the quantum Drude oscillator (QDO) model, employing scaling laws for the equilibrium distance and the binding energy, and applying the microscopic law of corresponding states. The vdW–QDO potential is shown to be accurate for vdW binding energy curves, as demonstrated by comparing to the ab initio binding curves of 21 noble-gas dimers. The functional form of the vdW–QDO potential has the correct asymptotic behavior at both zero and infinite distances. In addition, it is shown that the damped vdW–QDO potential can accurately describe vdW interactions in dimers consisting of group II elements. Finally, we demonstrate the applicability of the atom-in-molecule vdW–QDO model for predicting accurate dispersion energies for molecular systems. The present work makes an important step toward constructing universal vdW potentials, which could benefit (bio)molecular computational studies.


Introduction
2][3][4][5] The accurate description of vdW interactions requires sophisticated quantum-mechanical treatment, using the adiabatic-connection fluctuation-dissipation theorem (ACFDT) in density-functional theory or high-level quantum chemistry methods, such as coupled cluster or quantum Monte-Carlo. 4,5However, the prohibitive computational cost of these methods precludes their applicability to extended (bio)molecular systems.Therefore, practical simulations of large and complex systems are often done using classical force fields such as AMBER, 6 CHARMM, 7 or GROMACS. 8r the description of vdW forces, these popular force fields resort to the seminal Lennard-Jones (LJ) 9 (or an improved Buckingham 10 ) potential as a practical shortcut.Two parameters, the well depth D e and the equilibrium position R e , fully specify the LJ potential.
However, these parameters can be determined unambiguously only for relatively simple vdW-bonded systems, such as noble gas dimers or crystals.Moreover, the LJ potential is notorious for its lack of flexibility and very limited quantitative accuracy. 11,12On the other hand, the celebrated Tang-Toennies potentials [13][14][15][16] are derived from first principles and yield high accuracy for dimers including noble gases and group II elements.To achieve such an accuracy, the Tang-Toennies potentials employ from 5 to 9 parameters, depending on the exact flavor. 16Setting these parameters requires knowledge of R e and D e for each vdW bonded dimer, 13,15 which prevents a generalization of the Tang-Toennies models to the whole periodic table.Moreover, like the LJ potential, the most recent conformal Tang-Toennies (TTS) potential 15 is prone to large errors for dispersion coefficients (see Fig. 1a) despite its high accuracy close to equilibrium distances.Hence, a vdW potential combining wide transferability across the periodic table, high accuracy and minimal parametrization is not yet available.
Here, we develop a universal conformal pairwise vdW potential, which can be parametrized for all chemical elements based solely on two non-bonded atomic properties -static dipole polarizability α 1 and dipole-dipole dispersion coefficient C 6 .Our potential is consistently derived within the framework of the quantum Drude oscillator (QDO) model 17 using the Heitler-London perturbation theory, 18,19 and it is devoid of adjustable parameters.This is achieved by building connections between atomic scaling laws, [20][21][22] the microscopic law of corresponding states, [23][24][25][26] and symmetry-adapted perturbation theory (SAPT) 27,28 for intermolecular interactions.The derived exchange repulsion term in our potential obeys correct physical limits both at R → 0 and R → ∞, and the predicted C 6 dispersion coefficients are significantly more accurate compared to the other conformal Lennard-Jones and Tang-Toennies 15 potentials.The designed vdW-QDO potentials are twice as accurate as the LJ potentials, when averaged over 15 noble-gas dimers.In addition, the vdW-QDO potential augmented by a damping function can accurately describe binding curves of dimers consisting of (closed-shell) group II atoms.Moreover, the vdW-QDO potential can be applied to molecular systems, when coupled with atom-in-molecule (AIM) approach. 29We demonstrate this by accurately reproducing the dispersion energy for dispersion-dominated molecular dimers from the S66×8 dataset. 307][38][39][40][41][42][43][44][45] Within the QDO model, the response of valence electrons is described via a quasi-particle (drudon or Drude particle) with a negative charge −q and mass µ, harmonically bound to a positively-charged pseudo-nucleus of charge q with a characteristic frequency ω.Coupled QDOs are also extensively used in the development of vdW density functionals, 29,40,46 quantum mechanical 17,36 and polarizable force fields 39,[47][48][49][50] as well as recent machine learning force fields. 51,52e QDO model has been already used to build interatomic vdW potentials for water or noble-gas dimers and crystals. 17,36,38,39,53However, within the corresponding studies, the repulsive term was added in ad hoc manner, either by fitting Born-Mayer 10,54 exponents to ab initio repulsive walls 17,36,38,39 or by directly adding the Hartree-Fock exchange energy. 53erefore, such potentials cannot be generalized beyond the systems for which direct firstprinciples simulations are possible.In contrast, here we suggest a consistent treatment of both Pauli (exchange) repulsion and vdW dispersion within the QDO framework.To our knowledge, this is the first vdW potential of such a type, which does not directly utilize the reference binding energy of dimers or the Hartree-Fock exchange energy curve, but nevertheless provides relatively good accuracy.

Model construction
The long-range vdW dispersion energy for two identical QDOs is given by the usual multipolar series 17,35 where the dispersion coefficients are related to the oscillator parameters via the closed-form expressions 17 where α 1 = q 2 /µω 2 is the QDO dipole polarizability and k e = 1/4πε 0 .Tang and Toennies showed 13,14 that including the three leading dispersion terms is sufficient to obtain the accurate vdW potential.Therefore, we also truncate the series of Eq. ( 1) at the C 10 term.
The exchange repulsion is introduced into the model according to Refs.20,21, where multipole contributions to the exchange energy of a homonuclear dimer were derived by considering two identical drudons as bosons, assuming that they represent closed valence shells of atoms with zero total spin.Consequently, the total wave function of a dimer is  1.(b) The vdW-QDO potential for neon dimer benchmarked to the TTS potential and the reference CCSD(T) potential. 55Vertical dotted lines indicate the equilibrium distance as predicted by CCSD(T) (black) and vdW-QDO potential (red).For the comparison, the LJ potentials in two different parametrizations (see our discussion in the text) are also displayed in blue and green.
represented by a permanent where ψ A (r 1 ) = (µω/πℏ) 3 /4 e −(µω/2ℏ)r 2 1 and ψ B (r 2 ) = (µω/πℏ) 3 /4 e −(µω/2ℏ)(r 2 −R) 2 are, respectively, the ground-state wave functions of drudons centered at nuclei A and B separated by R. Within the Heitler-London perturbation theory, 18,19 the exchange energy of two identical vdW-bonded QDOs for distances near equilibrium and larger is well approximated by the exchange integral 20 The evaluation of Eq. ( 4) with the multipole expansion of Coulomb coupling VC between the two QDOs results in multipole contributions to the exchange energy. 20,21In dipole approxi-mation, this yields where S is the overlap integral.Higher-order multipole contributions (l > 1) to the exchange repulsion energy J (l) ex have the same leading-term dependence on internuclear distance R, with the only difference in a proportionality coefficient, i.e.J (l) ex ∝ k e q 2 S/R. 21Therefore, we introduce an effective exchange repulsion energy as with the proportionality coefficient A to be determined self-consistently, in what follows.In this way, we effectively include multipole contributions to all orders.Importantly, our E eff ex has 1/R dependence, which properly describes the infinite repulsive wall at short distances.
Thus, in contrast to the Born-Mayer or Duman-Smirnov 56-58 functional forms for exchange repulsion possessing a finite value of E ex at R → 0, Eq. ( 6) is in agreement with the orbital overlap model for Pauli repulsion. 59,60[62] In order to determine the coefficient A in Eq. ( 6), we employ the force balance condition at the equilibrium distance, To evaluate the equilibrium distance R e in our model, we use the quantum-mechanical relation between the atomic (static) dipole polarizability and vdW radius 20 where the proportionality coefficient Φ is given by 22 with α fsc = e 2 /4πε 0 ℏc ≈ 1/137.036as the fine-structure constant.The relation given by Eqs. ( 8)-( 9) turned out to be valid for real atoms.Especially, it is very accurate for noble gases, where the mean absolute relative error (MARE) 20,22 Since by definition R vdW is a half of the equilibrium distance R e in a homonuclear vdW bonded dimer, 20,63 accurate equilibrium distances can be obtained via With α 1 and C 6 being fixed, there are two unknown quantities in Eq. ( 7), A and µω, since C 8 and C 10 are solely expressed in terms of C 6 and µω via Eq.( 2).
As shown in Ref. 45, the product µω can be obtained from the force balance in the dipole with R e substituted from Eq. (10).Solution of this transcendental equation allows to determine the three oscillator parameters {q, µ, ω} given only {α 1 , C 6 }.Let us call this parametrization scheme vdW-OQDO, similar to the recently suggested OQDO scheme. 45e details of the procedure and the corresponding values of {q, µ, ω} can be found in the Supporting Information.
Solving Eqs. ( 11) and ( 7) together, one can obtain and the total vdW potential The vdW-QDO potential for neon is displayed by the red curve in Fig. 1b, which shows excellent agreement with the TTS potential 15 as well as with the CCSD(T) calculations 55 across the whole range of distances from 0.7R e (∼4 bohr) to infinity.Inclusion of C 8 and C 10 dispersion coefficients together with the suggested approach to treat exchange repulsion energy allows us to predict the correct depth and shape of the potential without losing the accuracy in predicting the equilibrium distance, which is inherited from the dipole approximation.In addition, we compare our potential to the LJ potential, for which we use two different parametrizations: LJ1 derived from thermodynamical properties and LJ2 designed to reproduce reference R e and D e (see the Supporting Information for more details).We note that the present vdW-QDO potential (13) performs accurately in the whole range of distances, whereas the LJ1 potential (blue curve in Fig. 1b) underestimates the energy in potential minimum region and the LJ2 potential overestimates the long-range energy (green curve), although both being reasonably accurate in the repulsive region.This imbalance and lack of flexibility of the LJ potential, which is observed for all noble gases, is one of the main issues limiting its quantitative predictive power. 11,12Moreover, the LJ potential severely overestimates C 6 coefficient (Fig. 1a), which is responsible for the correct long-range energy.The proposed vdW-QDO potential overcomes these difficulties without increasing the number of parameters.Moreover, our potential recovers correct bonding behavior using only a free atom property α 1 and asymptotic interaction parameter C 6 , which do not contain information about the interaction between atoms at short distances.

Application to noble-gas dimers
4][25] Namely, for the vdW potential of other noble-gas dimers we write where with the numerical values of the starred (unitless) parameters and their definitions presented in Table 2. Thus, only R e and D e for every dimer are required to obtain their vdW potential.For R e , the accurate scaling law ( 10) is already established, whereas an analogous scaling law for D e of noble-gas dimers is not yet known.Substituting R = R e to Eq. ( 13) and using Eq. ( 7) to eliminate Ak e q 2 e − (γR) 2 2 /R yields with β = µω ℏ R 2 e = (γ * ) 2 .Analyzing reference CCSD(T) data for D e from Ref. 15, we found that Eq. ( 16) truncated at first two terms can accurately predict D e for all noble-gas dimers In Fig. 2, D e by Eq. ( 17) are compared to the reference CCSD(T) data.The bar chart shows that Eq. ( 17) is accurate for homo-and heteronuclear dimers of He-Xe with all errors below 1 meV.For dimers with Rn the error is larger, with Rn 2 being underbound by 4.5 meV, or 13%.The larger errors for Rn dimers likely stem from the fact that the reference coupled-cluster calculation 74 is less reliable than the corresponding calculations for the lighter dimers He 2 -Xe 2 .larger overestimation of D e should be expected for Rn 2 , 74 where relativistic effects are more pronounced.Accounting for that, the estimated error of Eq. ( 17) for Rn 2 would not exceed 5.5%.We conclude that the suggested scaling law (17) allows one to accurately evaluate D e for all noble-gas dimers given only {α 1 , C 6 } without involving any adjustable parameters.
To extend the developed potential to heteronuclear dimers, combination rules for potential parameters can be used.The simplest ones are given by 6][77] Therefore, instead of mixing R e and D e , we use mixing rules for α 1 and C 6 , since our potential is fully parametrized by these two quantities.With the effective mixed values {α AB 1 , C AB 6 }, we can set three oscillator parameters {q, µ, ω} through the same vdW-OQDO parametrization procedure as for homonuclear dimers, which is described in the Supporting Information.By doing so, even the heteronuclear dimer AB is effectively represented by two identical oscillators, which still allows us to apply the formalism for exchange repulsion, Eqs. ( 3)-( 4), developed for homonuclear dimers.For C 6 dispersion coefficient, the very accurate combination rule arising from the London formula is already well known 13,78 To combine polarizabilities, we employ the robust mixing rule for vdW radii which was established in Ref. 20, where it was shown that accurate equilibrium distances in noble-gas dimers (MARE = 1%) are delivered by similar to the homonuclear case (8).Thus, the effective polarizability α AB 1 can be simply represented by Thereby, combining Eqs. ( 14), ( 15), ( 17), ( 19) and ( 21) with the vdW-OQDO parametrization scheme (see the Supporting Information), we obtain vdW-QDO potentials for all 21 noble-gas dimers.They are shown in Fig. 3  with 5.35 bohr predicted by Eq. ( 10) against the reference value of 5.608 bohr. 68In addition, the actual error in potential for He dimers is small in magnitude (despite being seemingly large visually due to the scale of y-axis in Fig. 3b).
To evaluate the accuracy of our potential quantitatively, we introduce the normalized area difference metric ∆ S between tested and reference potentials as The essential physical meaning of ∆ S is illustrated in Fig. 4b which shows that this single unitless number represents a measure of difference between tested and reference potentials.
The integration limits are set to 0.8R e and 2.0R e in order to evaluate accuracy close to the minimum region, whereas the long-range accuracy can be evaluated separately in terms of dispersion coefficients.Previously, ∆-gauge was used in benchmarks of various densityfunctional codes in calculations of equations of state for solids. 79,80The TTS potential 15 was chosen as the reference potential, and we benchmark vdW-QDO and LJ1 potentials with respect to it (a similar benchmark for LJ2 can be found in the Supporting Information).
The computed ∆ S -matrices for 15 He -Xe dimers are displayed in Fig. 4c-d (Rn dimers are omitted since there are no LJ parameters for Rn available in literature).We note that the vdW-QDO potential has twice better accuracy than the LJ one with ⟨∆ QDO S ⟩ = 9.0% compared to ⟨∆ LJ1 S ⟩ = 18.4% and ⟨∆ LJ2 S ⟩ = 17.3% when averaged over 15 dimers.Helium is the obvious outlier for vdW-QDO with the highest ∆ S = 31.1%,whereas for all other dimers ∆ S is below 13.6%.In contrast, LJ potential shows much broader variations in ∆ S spanning from 3.7% for Ar-Kr to 43.2% for Xe 2 .Although the LJ potential shows better ∆ S values than vdW-QDO for some of the dimers (He 2 , Ne 2 , He-Ne, Ne-Ar, Ar-Kr), overall the performance of vdW-QDO potential is more accurate and robust.Generally, Fig. 4 supports the above conclusions about the accuracy of vdW-QDO potential based on Fig. 3.We note that the predictions of LJ potential become worse for heteronuclear dimers composed of small and large atoms (e.g.He-Ar, He-Kr, Ne-Kr) than for atoms with a relatively close size (Ne-Ar, Ar-Kr) (Fig. 4c).In contrast, the evenly accurate predictions of the vdW-QDO model (Fig. 4d) suggest that the combination rules (19) and (21) employed in this work are more accurate and robust than the Lorentz-Berthelot rules (18).
To evaluate the quality of the potentials in the long-range limit, in Fig. 4a we compare the dispersion coefficients predicted by vdW-QDO and TTS potentials to the reference ab initio values (Table 1).Such a comparison is fair, since both potentials are built as conformal ones, unlike the earlier TT-2003 potentials, 13 which directly utilize reference C 2n coefficients for every element being therefore not strictly conformal.To recover the TTS dispersion coefficients, we used the reported scaling law C 2n = C * 2n × D e R 2n e with C * 6 = 1.3499,C * 8 = 0.4147, C * 10 = 0.1716. 15For vdW-QDO dispersion coefficients, from Eqs. ( 14), (15) and ( 17) one can obtain where C * 6 = 1.1779,C * 8 = 0.3848, C * 10 = 0.1540 (see Table 2).We found that both potentials severely underestimate C 8 and C 10 and demonstrate similar magnitude of these errors increasing with the atomic number.However, for C 6 vdW-QDO potential performs much better, showing homogeneous overestimate of 12-13%, whereas the TTS potential again possesses increasing magnitude of error, reaching its maximum of 41% for Rn.While C 8 and C 10 are important to deliver accurate potential near the equilibrium, in the asymptotic limit the quality of the potential is fully determined by the leading C 6 coefficient.Therefore, we can conclude that our conformal vdW-QDO potential shows physically more reasonable long-range behavior than the conformal TTS potential.
In contrast to the Tang-Toennies potentials, 13,15 the above vdW-QDO model does not employ any damping of the dispersion energy.However, for noble-gas dimers the damping of the dispersion energy is not essential, and interatomic vdW potential can be effectively described even without a damping function, as was shown above.This provides additional reasoning why the scaling law for vdW radius (8), which was originally derived without considering dispersion damping, 20 works so well.To check the effect of damping function on the results, we derived the damped vdW-QDO potential, where the QDO damping function reads It was found that for noble-gas dimers the obtained damped vdW-QDO potential curves are practically the same as the undamped ones within the considered range of distances (see Fig. S2 in the Supporting Information).Interestingly, the damping function (24) differs from the Tang-Toennies damping function 81 just in the upper summation limit (for the TT it is 2n) and in the physical meaning of the unitless variable (z = bR for the TT damping function, with b stemming from the Born-Mayer repulsion term Ae −bR ).Note that, due to the distinct form of the Pauli repulsion in the vdW-QDO and TT models, the QDO damping function contains only even powers of R up to 2n (see the derivation of the QDO damping function and more detailed discussion in the Supporting Infomation).
We also observe that for both the TT and vdW-QDO models exchange repulsion and dispersion terms separately do not agree with their SAPT counterparts, whereas the total potentials show very close agreement with the sum of the SAPT terms (see the Supporting Information).This finding is in line with the statement of Ref. 58 that the generalized Heitler-London theory delivers a more compact expansion of interaction energy than the SAPT.

Application to group II dimers
Another class of vdW systems where Tang-Toennies potentials work well consists of Me 2 dimers, with Me = Mg, Ca, Sr, Ba, Zn, Cd, Hg.5][86] The only exception is the Be-Be dimer, which has been being a long-standing puzzle for quantum chemistry.The potential curve of Be 2 has a remarkably different shape in the long-range region, compared to other alkaline-earth elements. 85Since this dimer does not obey the principle of corresponding states, it is excluded from our consideration here.In what follows, we show that the vdW-QDO potential is also capable to describe the potential curves of the group II dimers upon several modifications.
First, in contrast to noble gases, for Me 2 dimers the damping function ( 24) cannot be omitted due to much larger polarizabilities and hence dispersion coefficients than those of noble gases (see Table 3).As a result, without damping function a pronounced divergence of the vdW-QDO potential would already occur at near-equilibrium distances.Second, the scaling laws of Eqs. ( 10) and ( 17) are not so accurate for the group II dimers, since the bonding in Me 2 is not purely of vdW type. 84,85Therefore, the reference values of R e and D e (Table 3) were used for our vdW-QDO potential.
Following Tang et al., 85 we choose Sr 2 dimer as the reference system to get the shape of the potential curve and then rescale it onto other dimers.Similar to the case of noble gases, the vdW-QDO potential shape for Sr 2 dimer was derived as (see the Supporting Information) with the numerical values of its parameters presented in Table 4.The QDO parameters {q, µ, ω} for the Sr 2 dimer were set using α 1 , C 6 and R e following the damped vdW-OQDO procedure, as explained in the Supporting Information.Altogether, the three physical quantities are employed to parametrize the vdW-QDO potential for Sr 2 , compared to the five {R e , D e , C 6 , C 8 , C 10 } in case of the Tang-Toennies potential. 85For other Me 2 dimers, we need only R e and D e from Table 3 to perform the rescaling [95] both ab initio and experimentally derived potentials (when they are available).This is a remarkable result, since the binding energies of group II dimers are up to 5 times larger than those of the heaviest noble gases.Furthermore, the shape of their potentials is distinct to the one of noble gases (Fig. 6).Thus, the vdW-QDO functional form is robust and well-suited to describe vdW potentials across various types of systems.In contrast, the LJ potential cannot be employed to describe group II dimers with any combination of parameters, since its energy well is too narrow for such binding curves (Fig. 6).

Application to molecular dimers
Finally, we can show that the developed vdW-QDO potential is also applicable to molecular systems, with an example of eight dispersion-dominated molecular dimers from the S66×8 benchmark dataset. 30They are hydrocarbons, including benzene as well as aliphatic and cyclic molecules (see the list in the Supporting Information).Such systems were chosen to diminish an influence of electrostatic term not included in the current vdW-QDO potential.
We compute the energy of vdW interaction between molecules A and B at given intermolecular separation as where summation goes over the atoms i and j of the molecules A and B, respectively, and R ij is the interatomic distance.Interaction energy between a pair (i, j) is given by the damped vdW-QDO potential To set the QDO parameters {q, µ, ω}, we apply the vdW-OQDO procedure (see the Supporting Information) coupled with the atom-in-molecule (AIM) approach to each pair of atoms (i, j).Following the Tkatchenko-Scheffler (TS) method, 29 the AIM polarizabilities and dispersion coefficients are defined by where V free i and V AIM i are the corresponding Hirshfeld volumes.To compute them, singlepoint DFT-PBE0 96 calculations for every dimer were performed at their equilibrium geometry.Then, the effective polarizability α ij 1 and dispersion coefficient C ij 6 for a pair (i, j) were defined using the combination rules of Eqs. ( 19) and (21).Finally, the vdW-OQDO parametrization procedure (see the Supporting Information) was applied to map {α ij 1 , C ij 6 } onto {q, µ, ω}.After performing pairwise summation in Eq. ( 27) and repeating the whole procedure for all 8 intermolecular separations, the vdW-QDO interaction curves for dimers are obtained and compared to the reference CCSD(T) interaction curves. 30For comparison, interaction energies of dimers at PBE0+TS, 29 PBE0+MBD, 40 and DFTB3 97 +MBD levels of theory were also calculated.All DFT calculations in this work were performed using FHI-aims code 98 with the 'tight' atomic basis sets.
We showcase the obtained results with an example of neopentane dimer (Fig. 7b).One can see that PBE0+TS overbinds neopentane dimer significantly and underestimates the equilibrium separation by 5%.Inclusion of many-body effects at the PBE0+MBD level improves the results of PBE0+TS for energy, but the predicted equilibrium separation is still 95% of the reference value.On the other hand, DFTB3+MBD method clearly lacks repulsion and attraction at short-range and long-range distances, respectively, although at the equilibrium distance the two errors largely cancel each other.As for the vdW-QDO potential, we note that the minimum of intermolecular interaction curve is very close to the reference CCSD(T) energy, although our method overestimates the equilibrium separation by 10%.When considering the overall interaction curve, however, the vdW-QDO potential deviates significantly from the reference, being too repulsive at shorter distances.These conclusions remain valid also for other seven dimers, as supported by the corresponding results in the Supporting Information.
To get insights into such a behavior, we compared the exchange and dispersion terms in vdW-QDO potential (27) to their counterparts from SAPT-DFT calculations 99 (see Fig. 7a).

The dispersion part of vdW-QDO potential provides excellent agreement with E
(2) disp from SAPT-DFT, which illustrates the well-known fact that vdW dispersion interaction between small molecules can be effectively described by a pairwise potential, despite being essentially many-body in its nature. 4,40,100In contrast, for the vdW-QDO exchange repulsion term we observe a noticeable deviation from the E (1) exc contribution of SAPT-DFT.This might indicate that exchange repulsion between molecules is not accurately described by the commonly used pairwise approach and hence exchange repulsion requires many-body treatment as well. 101PT-DFT calculations 99 also indicate that electrostatic contributions are small but not negligible even for dispersion-dominated dimers (Fig. 7a and Figs.S6-S8 in the Supporting Information).To check a possible effect of inclusion of accurate many-body exchange and electrostatic interactions on our results, we consider the corrected vdW-QDO potential, where the first-order SAPT-DFT energy was added to the dispersion energy from the vdW-QDO model elst + E (1) exc + V disp QDO .
This approach is similar to the HFDc (1) scheme of Podeszwa et al. 102 with the important difference that here the dispersion energy is calculated from the QDO model, whereas in the HFDc (1) method this energy is computed from SAPT.The corrected vdW-QDO curve in Fig. 7b delivers much better description of the interaction energy than the original vdW-QDO method.Notably, at larger distances V corr QDO also shows a good agreement with PBE0+MBD energy.
The overall statistics for the eight molecular dimers is shown in Fig. 7c in terms of the error in the equilibrium energy obtained by five methods with respect to the reference CCSD(T) values. 30Purely analytical vdW-QDO predictions are scattered within 2 kcal/mol range, which is remarkable considering the approximations made in the model.Including the first-order SAPT energy reduces errors roughly to 1 kcal/mol (chemical accuracy).This test demonstrates that the vdW-QDO approach is not limited to atomic dimers and can be generalized to molecular complexes, although in that case the accuracy is currently lower.Nevertheless, considering the approximations made, the fact that the vdW-QDO method (even without SAPT corrections) properly predicts binding of the considered weaklybound molecular complexes from a set of AIM quantities (α 1 and C 6 ) is already reassuring.

Summary and outlook
We developed a universal pairwise vdW potential devoid of empiricism and parameterized by only two atomic non-bonded parameters.The developed vdW-QDO potential combines the strengths of the Lennard-Jones and Tang-Toennies models.Similarly to the former, the vdW-QDO potential is fully determined by only two parameters.At the same time, our potential is comparable in accuracy to the Tang-Toennies potential for noble-gas dimers, as QM7-X, 104 which include information about atom-in-molecule volumes.This enables the direct applicability of the vdW-QDO potential to arbitrary molecular systems.Finally, the vdW-QDO potential can be generalized to include polarization and electrostatic contributions, which is the subject of our ongoing studies.Eventually, this effort should deliver a general coarse-grained force field for all non-covalent interactions entirely based on the model system of coupled QDOs.

Figure 1 :
Figure 1: (a) Errors in dispersion coefficients arising from the LJ and TTS potentials.The TTS dispersion coefficients are obtained as C 2n = C * 2n × D e R 2n e , 15 and the reference dispersion coefficients C ref 2n are given in Table1.(b) The vdW-QDO potential for neon dimer benchmarked to the TTS potential and the reference CCSD(T) potential.55Vertical dotted lines indicate the equilibrium distance as predicted by CCSD(T) (black) and vdW-QDO potential (red).For the comparison, the LJ potentials in two different parametrizations (see our discussion in the text) are also displayed in blue and green.

Figure 2 :
Figure 2: Potential well depth D e by Eq. (17) compared to the reference CCSD(T) values 15 for 21 noble-gas dimers.Mean error (ME) and mean absolute error (MAE) are displayed.

Figure 4 :
Figure 4: (a) Errors in dispersion coefficients predicted by vdW-QDO (dark colors) and TTS (light colors) potentials.(b) Schematic illustration of the ∆ S metric calculation.(c) Heatmaps showing ∆ S (in %) for LJ1 (left) and vdW-QDO (right) potentials with respect to the reference TTS potential.The left and right colorbars have the same scale.

Figure 6 :
Figure 6: Dimensionless shapes of the vdW-QDO potential curves for Ne 2 (red) and Sr 2 (green) compared to the shapes of the LJ (dashed blue) and TTS (dashed black) potentials.

Table 1 :
The reference values of static dipole polarizability α 1 (in a.u.), dispersion coefficients C 6 , C 8 , C 10 (in a.u.), and dimer potential well parameters R e (in bohr) and D e (in meV) for noble-gas dimers.For R e and D e , the values in Å and Kelvin, respectively, are extra given in parentheses.Also, we compare their reference values (columns 5 and 7) to the predictions of Eqs.(10)and (17) (columns 6 and 8).The values for C 6 , C 8 , and C 10 labeled with the star ( * ) are taken from Ref. 13 instead of Refs.64,65.

Table 3 :
The reference ab initio values of the dipole polarizability α 1 (in a.u.), dipolar dispersion coefficient C 6 (in a.u.), and dimer potential well parameters R e (in bohr) and D e (in meV) of the group II dimers.For R e and D e , the values in Å and Kelvin, respectively, are extra given in parentheses.The used data sources are the following: a Ref. 88, b Ref. 89, c Ref. 90, d Ref. 86, e Ref. 91, f Ref. 92.