## Abstract

Modelling dispersion interactions with traditional density functional theory (DFT) is a challenge that has been extensively addressed in the past decade. The exchange-dipole moment (XDM), among others, is a non-empirical add-on dispersion correction model in DFT. The functional PW86+PBE+XDM for exchange, correlation and dispersion, respectively, compromises an accurate functional for thermochemistry and for van der Waals (vdW) complexes at equilibrium and non-equilibrium geometries. To use this functional in optimizing vdW complexes, rather than computing single point energies, it is necessary to evaluate accurate forces. The purpose of this paper is to validate that, along the potential energy surface, the distance at which the energy is minimum is commensurate with the distance at which the forces vanish to zero. This test was validated for 10 rare gas diatomic molecules using various integration grids and different convergence criteria. It was found that the use of either convergence criterion, 10^{−6} or 10^{−8}, in Gaussian09, does not affect the accuracy of computed optimal distances and binding energies. An ultra-fine grid needs to be used when computing accurate energies using generalized gradient approximation functionals.

This article is part of the themed issue ‘Multiscale modelling at the physics–chemistry–biology interface’.

## 1. Introduction

As first enlightened by Kristyan & Pulay [1], traditional functionals in density functional theory (DFT) lack the apposite physics to properly capture weak dispersion interactions. Dispersion is estimated to account for one-thousandth of the correlation energy [2]. This is why conventional exchange–correlation functionals fail to describe van der Waals (vdW) complexes [1,3,4] or, at most, give erratic results. Nonetheless, this issue has been resolved with functionals derived from first principles [2,5–14], and functionals with fitted empirical parameters [15–22]. These functionals have been implemented, at the quatum level and the multiscale quantum mechanics/molecular mechanics level, in numerous applications related to biological systems, inorganic complexes and material sciences [23–32].

The exchange-dipole moment (XDM) model proposed by Becke & Johnson [2,13,14] can accurately capture dispersion interactions, yet it is a simple model. This model is applicable to any DFT-D (dispersion) feasible method. For instance, it can be combined with PW86 [33] (a generalized gradient approximation (GGA) functional for exchange) and PBE [34] (a GGA functional for correlation) to obtain PW86+PBE+XDM, which comprises a non-empirical DFT method for vdW complexes. This functional, PW86+PBE+XDM, is of high accuracy for thermochemistry [35,36], binding energies of vdW complexes at equilibrium [36] and at distorted intermonomer separations [37,38]. The outstanding performance of this functional is due to the accurate reproduction of the exact Hartree–Fock exchange energy [39,40]^{1} by the PW86 GGA functional [41], and the non-empiricism of the XDM model for dispersion.

The XDM model is based on the position-dependent dipole moment of the exchange hole given by
1.1where *ρ*_{σ} is the *σ*-spin density and *ψ*_{iσ} are occupied Hartree–Fock (HF) or Kohn–Sham (KS) orbitals. Orbitals are assumed to be real. This hole is a measure of depletion of the probability (with respect to the total *σ*-spin electron density, *ρ*_{σ}) of finding an electron of spin *σ* at a point *r*_{2} next to a reference electron of the same spin *σ* at point *r*_{1}.

The damped [5,16,42,43] dispersion energy in XDM is given by
1.2where *R*_{vdW,ij} is an effective vdW separation that linearly depends on a critical radius *R*_{c,ij}
1.3where *a*_{1} and *a*_{2} are two universal fit parameters, and the critical radius *R*_{c,ij} is the separation at which the *C*_{m}/*R*^{m} (*m*=6,8,10) terms have roughly the same value
1.4

As shown in equation (1.2), the XDM model includes not only *C*_{6} coefficients but also *C*_{8} and *C*_{10} coefficients, all of which are environment dependent. This allows the functional to be more widely applicable for its accurate description of *inter*- or even *intra*-molecular interactions, which are ubiquitous in biological systems.

The functional PW86+PBE+XDM, for exchange, correlation and dispersion, respectively, has been proved to predict very accurate thermochemistry and binding energies from single point calculations [35–38]. The ultimate goal is to use this functional for optimizing vdW complexes. For reliable optimizations, accurate forces need to be evaluated.

Within the self-consistent-field (SCF) approximation, there is an assumption made when evaluating dispersion forces. The dispersion coefficients (shown in equation (1.2)) are assumed to be constant when taking the derivative of the energy with respect to *x*, *y* and *z* coordinates. Neglecting the infinitesimal changes in the electron-nuclear part of the KS potential (*V* _{eN}) at an infinitesimal change in nuclear coordinates avoids the complexity of taking the derivative of the dispersion energy with respect to the density. The purpose of this paper is to validate this approximation made within the post-GGA approach, where XDM dispersion is added perturbatively.

The computational cost incurred by the perturbative addition of XDM scales linearly with respect to the size of the system, an advantage attributed to the computationally efficient Becke–Roussel [44] exchange-hole function in XDM. Thus, the perturbative addition of XDM to a PW86+PBE calculation scales as *N*^{3}, which is basically the scaling of a DFT GGA functional with system size *N*. The computational cost of adding XDM, with the *C*_{6} coefficients, is little regardless of whether computed perturbatively or self-consistently [45]. Nevertheless, the post-SCF approach simplifies the implementation of the non-local correction, XDM, into existing DFT codes.

Kong *et al.* [45] studied the accuracy of the perturbative addition of XDM compared with the full self-consistent computation using B3LYP+XDM. It has been found that the error in the density matrix, as a result of the perturbative addition of XDM, is of the order of 10^{−5}. They have also shown that the square of this error in the density matrix is approximately equal to the error in the total energy. Error of the order of 10^{−5} in the density matrix is a borderline precision of gradient calculations.

This paper aims to verify that forces computed with PW86+PBE+XDM vanish at the equilibrium geometry so that the assumption of constant dispersion coefficients upon small nuclear changes (as mentioned above) remains valid. This assessment is crucial as a proper evaluation of the analytical forces guarantees a self-consistent convergence of the DFT energy and densities. This is because an error in the nuclear gradient evaluation is similar in magnitude to the error in the density matrix.

## 2. Computational details

In this study, calculations are performed using PW86+PBE+XDM in a post-SCF fashion. The GGA part of this functional PW86+PBE is computed self-consistently with the G09 package [46]. The XDM part is then added as a perturbative correction to the self-consistent GGA exchange–correlation energy. The binding energy and the force curves of 10 rare gas diatomic molecules (from He, Ne, Ar and Kr noble gases) are computed. Single point calculations are performed on the diatomic molecules at the experimental equilibrium internuclear separation (obtained from [47]) with ±40 points built by increasing or decreasing the distances by increments of 0.01 Å. The forces are simple analytical forces evaluated by taking the derivative of the energy with respect to the infinitesimal changes in the nuclear coordinates (d*E*/d*r*), i.e. it is the gradient calculation with respect to the atomic movement.

Augmented Dunning double (aug-cc-pVDZ) and triple zeta (aug-cc-pVTZ) basis sets are tested because the latter gives accurate binding energies, while the former is sufficient for geometry optimizations [36]. Ahlrichs density fit W06 [42] is used as it was proved to accelerate calculations without changing the accuracy of the binding energy [38]. The use of W06 is greatly advantageous, especially with larger systems, as it saves a significant amount of time. A computation with the triple zeta basis set would be one order of magnitude less expensive when W06 is introduced, i.e. the cost with a triple zeta basis set including W06 would be comparable to the cost with a double zeta basis set without W06.

Two SCF convergence criteria were tested: 10^{−6} and the more tight criterion 10^{−8}.

For the GGA part, two different grids are considered: 099590 and 200590. The grid notation xxxyyy means that there are xxx radial shells surrounding each nucleus and yyy angular points distributed evenly on each shell. For the XDM part, 200302 grids and smaller grids, denoted by ‘default’ in this paper, are also considered in the assessment. For all calculations where the grid is not explicitly mentioned, the 200590 grid is used.

## 3. Results

Using PW86+PBE+XDM, the energies and forces will be tested on diatomic molecules of closed-shell rare gas dimers as they are the prototype of vdW interactions.

The forces are total forces (GGA+XDM forces) acting on each nucleus. As these molecules are only diatomic molecules, the force on nucleus A is equal in magnitude to the force on nucleus B, but opposite in sign. This means it is facile to obtain the magnitude of forces for an unequivocal analysis. In the case of triatomic or polyatomic molecules, the analysis might be ambiguous depending on the choice of the atom on which the force is exerted, especially if the molecules have no atoms at the centre of mass.

Figure 1 is a graph combining two plots for the binding energies and the forces in He–Kr as a function of the internuclear separations. The results shown are from calculations using aug-cc-pVTZ, a 200590 grid and a convergence criterion of 10^{−8}. This figure is representative of 10 similar graphs for the 10 rare gas diatomic molecules considered in this study. The binding energy curve and the force curve are displayed in the same graph to investigate whether the vanishing forces are commensurate with the energy minimum along the potential energy surface. Binding energies are plotted rather than absolute energies for scaling purposes. The graph in figure 1 illustrates how, at distance 3.67 Å, the binding energy of He–Kr is minimal, −0.0646 kcal mol^{−1}, and the force vanishes to 2.64×10^{−7} atomic units (a.u.), which is essentially considered to be zero. The calculated equilibrium separation is in great agreement with the experimental equilibrium separation, i.e. 3.69 Å [35], with only 0.5% difference. The functional PW86+PBE+XDM slightly overbinds He–Kr by 3.4% compared with the experimental binding energy, which is −0.0625 kcal mol^{−1} [35]. As claimed in many references [48–50], the discrepancy in the binding energy predictions can be attributed to inadequacies of the exchange functional used. However, this is unlikely to be the situation in this case, as the exchange functional PW86 used in PW86+PBE+XDM has been tested to accurately reproduce the exact HF repulsion energy at equilibrium separations [41]. The discrepancy in this case might be due to the basis set or due to the values of the parameters *a*_{1} and *a*_{2}. Nevertheless, in this paper, the accurate reproduction of the experimental values, i.e. binding energies and equilibrium separations, is not as crucial as obtaining a commensurate equilibrium distance where the binding energy is minimal and the forces vanish to zero, for a given diatomic molecule optimized at a particular basis set with its corresponding parameters *a*_{1} and *a*_{2}. For all 10 rare gas dimers considered in this study, the distance of minimal binding energy is perfectly commensurate with the distance of vanishing forces. In the worst scenarios, the distance at which the force is minimal (virtually zero) would be 0.01 Å greater than the distance at which the binding energy is at a minimum. It is thus confirmed, with no doubts, that negligible errors are introduced by excluding the electronic part of XDM from nuclear derivative calculations.

In agreement with the results shown above for the PW86+PBE+XDM functional, Kong *et al.* [45] have shown that the deviation of the density matrix between SCF XDM and the perturbative addition of XDM in B3LYP+XDM is of the order of 10^{−5}. They have also quantified the error in differentiating only the classic part (excluding the electronic part) of XDM with respect to nuclear coordinates; the errors do not exceed 10^{−3} a.u.

For the ultimate goal of a structural optimization, there are other minor, yet essential, details that should be considered when computing energies and forces. These details include the choice of a suitable SCF convergence and a fine grid for the Coulomb integrals.

### (a) Assessment of the SCF convergence criterion: binding energies

While an SCF convergence criterion of 10^{−8} ensures a tighter threshold, it might complicate the convergence process of many complexes. For example, with the 10^{−8} threshold, a stacked benzene dimer does not converge even after 129 SCF cycles, whereas it converges within eight cycles if the 10^{−6} SCF convergence criterion is specified (AA Arabi 2010, unpublished data). The purpose of this section is to investigate the effect of lowering the SCF convergence criterion from 10^{−8} to 10^{−6} on the binding energies of the 10 rare gas diatomic molecules. This assessment will be repeated with both augmented double and triple zeta Dunning basis sets (figure 2) as aug-cc-pVDZ is satisfactory for optimization purposes and aug-cc-pVTZ is accurate for computing binding energies [36].

The graphs for the rest of the rare gas diatomic molecules are not shown as the binding energy curve of Ar–Kr is similar to that of Ar–Ar (figure 2*a*) and the rest of the diatomic molecules (apart from Kr–Kr, figure 2*b*) are similar to that of Ne–Kr (figure 2*c*). It is evident from figure 2 that, for the same basis set, using a different convergence criterion would make no difference in the binding energies along the whole range of internuclear separations. Figure 2 depicts how, with a triple zeta basis set, the curves are perfectly overlapping regardless of the SCF convergence criterion used. The same observation is made with the double zeta basis set. With either basis set, and for all 10 diatomic molecules, the maximum difference in binding energy between the SCF convergence criteria of 10^{−6} and 10^{−8} is of the order of μkcal mol^{−1}, i.e. negligible. This means that a convergence criterion of 10^{−6} allows high accuracy of binding energies even for the weakest cases such as the helium dimer. These results give us confidence in using an SCF convergence criterion of 10^{−6} to appreciably facilitate the SCF convergence without being concerned about the accuracy of the binding energies.

These graphs also illustrate how a double zeta basis set binds the rare gas diatomic molecules less than a triple zeta basis set would. For Ne–Kr (figure 2*c*) and molecules with similar plots, at a certain distance beyond the equilibrium separation, this trend is flipped, i.e. the molecules are more bound with a double zeta basis set than those computed using the triple zeta basis set. This might be attributed to a basis set superposition error. Nonetheless, it is worth mentioning that the difference in binding energies between aug-cc-pVDZ and aug-cc-pVTZ, with either convergence criterion, is only 0.1 kcal mol^{−1}.

### (b) Assessment of the SCF convergence criterion: forces

The graphs in figure 3 are similar to the graphs in figure 2; however, the forces are assessed rather than the binding energies. For all rare gas diatomic molecules, the nuclear derivatives computed with aug-cc-pVDZ are smaller than those computed with aug-cc-pVTZ. As depicted in figure 3*a*,*b*, the only difference between the 10 diatomic molecules is the separation between the double and the triple zeta curves. Again the differences are very minuscule: at the convergence criterion of 10^{−6}, the maximum difference in the force between aug-cc-pVDZ and aug-cc-pVTZ is 100 microHartree. This difference drops to 90 microHartree with a convergence criterion of 10^{−8}.

The concise message to be drawn from figure 3 is that the choice of the SCF convergence criterion (10^{−6} versus 10^{−8}) does not affect the values of the forces regardless of the basis set used. There is a maximum of 0.09 microHartree difference between different convergence criteria used at either basis set.

### (c) Assessment of the grids

The integrals in the exchange–correlation potential in DFT GGA functionals are too complicated to be solved analytically. They have to be evaluated numerically on a real space grid [51]. Johnson *et al.* have shown that fine grids are needed to obtain reliable results and smooth potential energy curves with meta-GGA functionals [52]. The grids required for computation with a GGA functional are investigated in this section (figure 4).

These graphs show that changing the grid in the computation of the XDM part would make no difference in either the forces or the binding energies. However, it is apparent that changing the grid for the GGA part changes the shape of the binding energy and force curves. Using the 099590 grid gives sporadic results with oscillations in the force curve of He–Kr, but the change in the binding energy curve is much less pronounced. To remedy the situation, a refined grid, 200590, should be used. This refined grid includes more radial shells around each nucleus and gives smooth binding energy and force curves.

As to observations made with the rest of the diatomic molecules: He–Ar and He–Ne are sensitive to the grids used (just as He–Kr is); He–He, He–Ne and Ne–Ne are less sensitive; while Ar–Kr, Ar–Ar, Kr–Kr, Ne–Ar and Ne–Kr are not sensitive at all to the change of grids.

## 4. Concluding remarks

The perturbative addition of XDM to PW86+PBE is regarded as accurate and reliable, which is in good agreement with what has been shown for B3LYP+XDM in [45]. For all 10 rare gas diatomic molecules, the distance at which the binding energy is minimal along the potential energy surface is commensurate with the distance at which the forces vanish to zero. This validates the assumption of neglecting the electronic part of XDM from nuclear derivative calculations, which allows a more efficient treatment of larger systems. For a smooth potential energy surface and accurate binding energies and forces, a fine grid must be used. This paper shows that the 200590 grid is more reliable than the 099590 grid. The use of an SCF convergence criterion of 10^{−6} is sufficient. For SCF convergence criteria of 10^{−6} and 10^{−8}, the difference in the evaluated binding energies and forces of the 10 rare gas diatomic molecules considered in this study is very minuscule; therefore, it can be neglected.

## Competing interests

I declare I have no competing interests.

## Funding

I thank the Natural Sciences and Engineering Research Council of Canada (NSERC), the Killam Trust of Dalhousie University and Zayed University (Abu Dhabi) for the financial support.

## Acknowledgements

I am grateful to Professor Axel D. Becke for providing me with the code that computes the forces and for his valuable feedback on this manuscript. I would also like to thank Felix O. Kannemann for the helpful discussions about basis sets, and the Atlantic Computational Excellence Network (ACEnet) for the computer resources.

## Footnotes

One contribution of 17 to a theme issue ‘Multiscale modelling at the physics–chemistry–biology interface’.

↵1 Exact exchange repulsion energy for rare gas diatomic molecules was computed with the basis set free NUMOL code.

- Accepted July 5, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.