A method for incorporating dispersion forces into density functional theory through the addition of a scaled atom–atom pair wise attractive potential term to the repulsive density functional potential is presented. By considering the deviation from the $\frac{C_6}{R^6}$ potential arising from the dispersion interaction of two, non–point, spherically distributed dipoles, an analytic function was derived that when multiplied by $\frac{C_6}{R^6}$ resulted in an attractive potential capable of inducing binding. The scaling function found to give bound Ne and He dimers using Becke’s three parameter hybrid density functional with the Lee–Yang–Parr correlation functional (B3LYP) density functional potential is quite different from the damping functions used in previous studies, causing the $\frac{C_6}{R^6}$ attraction to be strengthened in the bonding region. Using the scaled attractive potential presented in this work, a dissociation energy 0.4% in error (+0.000015 kcal/mol error) was found for Ne dimer, compared to 60% error (–0.0023 kcal/mol error) when the previous damped attractive potential is used. For He dimer, use of the new scaled attractive potential yields a dissociation energy in 8% error (+0.000074 kcal/mol error); however, the previous damped attractive potential cannot induce binding at all. The interaction potential curves are quite close to the experimental potential curves for both Ne₂ and He₂. This model has been incorporated into the APF–D (Austin–Petersson–Frisch) Density Functional theory for dispersion forces within the Gaussian 09 Suite under keyword, APFD.

Density functional methods have become a popular means of assessing the electronic structure and quantum mechanical properties of molecular systems, owing to their accuracies comparable to Moller–Plesset 2^nd Order Perturbation Theory while being more cost–effective.^1–10

The computationally less demanding density functional methods would seem ideal for investigations of large systems. However, dispersion forces, which arise from the interaction of instantaneously fluxuating dipoles.¹¹ are often important for large systems, such as in the stabilization of protein structures and stacked DNA bases.^12–16 Yet, the current density functional methods in use cannot account for these long–range effects since the exchange–correlation energy depends only on the local density and at most its gradient.¹ For example, the use of gradient–corrected exchange functionals has been shown to display purely repulsive potential energy curves for the interaction of rare gases.¹⁷ Although some combinations of exchange and correlation functionals do produce overly bound rare–gas dimers,¹⁸ such as those using exchange functionals rich in local density approximation, the interaction potential does not exhibit the familiar R–6 internuclear dependence given by perturbation theory.^19–23

Following the work of Alrichs et al.²⁴ on the correction of the Hartree–Fock model, as well as the work of Moorij et al.²⁵ Wu & Yang.²⁶ introduced an empirical interaction term to the first–principles density functional theory in order to simulate the dispersion interaction energy.

This empirical term, given by Eq. (1), is completely independent of the density functional and depends only on the positions of the nuclei.

$U_{attract}\left(R\right)=\frac{-f_{d\left(R\right)C_6}}{R^6}$   (1)

For each atom pair separated a distance of $R,$ $U_{attract}$ is added to the repulsive density functional potential.²⁶ The coefficient $C_6$ controls the asymptotic behavior of the attractive term for large internuclear separation.$f_d\left(R\right)$ is a damping function that weakens $U_{attract}$ for small internuclear separations.^24,25 Thus, the damping function is zero for small R, in the region where R<0.2RM, and becomes one for large R, in the region where R>1RM, where RM is the summation of van der Waals radii. As shown in Figure 1, the effects of two different damping functions were examined, damping function I approaching unity much more gradually than damping function II.²⁶

$\begin{array}{l}{f}_{d}I\left(\frac{R}{R_M}\right)={\left(1-e^{-C_{damp{\left(\frac{R}{R_M}\right)}^3}}\right)}^2\\ f_{d}II\left(\frac{R}{R_M}\right)=\frac{1}{1+e^{-\beta \left(\frac{R}{R_M}-1\right)}}\end{array}$

Figure 1a Damping Functions I and II as a Function of Internuclear Separation.

Figure 1b Lennard Jones attractive Energy term as a function of Internuclear Separation for Ne-Ne dimer.
Figure 1 (a). fd(R/RM) I and II as a function of R/ RM, where RM is the summation of van der Waals radii; fdI approaches unity much more rapidly than fdII. (b) Uattract=- fd (R/ RM)C6/R6 for fdI and II.

For six rare gas systems (He₂, Ne₂, Ar₂, HeNe, HeAr, and NeAr), the repulsive density functional potential, given in Eq. (2).²¹ was added to the attractive term described above in order to obtain the interaction potential, $U_{\mathrm{interact}}$ , given in Eq. (3).²⁶

$U_{DFT}\left(R\right)\equiv E_{{}_{DFT}}\left(R\right)-E_{DFT}\left(20A\right),$   (2)

$U_{\mathrm{interact}}\left(R\right)=U_{DFT}\left(R\right)+U_{attract}\left(R\right).$   (3)

Binding was obtained in only three systems (Ne₂, Ar₂, and NeAr) using the B3LYP hybrid density functional while no binding was produced using the BLYP or BPW91 functionals for any of the six systems.²⁶ Figures 2 & 3 illustrate the effect of adding the empirical term to the repulsive density functional potentials with damping function I for Ne₂ and He₂, respectively. ²⁶

Figure 2a UDFT as a function of Internuclear Separation for Ne Dimer.

Figure 2b Uinteract=UDFT-fdIC6/R6 for Ne dimer.

Figure 2c Uinteract=UDFT-fdIC6/R6 (Bonding Region) for Ne dimer.
Figure 2 (a) UDFT as a function of internuclear separation for Ne dimer. All density functional potentials are purely repulsive; (b) Uinteract= UDFT+Uattract as a function of internuclear separation for Ne dimer. Only the B3LYP potential yields binding. De=0.0015eV (experimental*=0.0038eV) Uattract=-fdI(R)C6/R6 (c) Enlargement of bonding region in (b).

Figure 3a UDFT as a function of internuclear separation for He dimer.

Figure 3b Uinteract=UDFT-fdIC6/R6 for He dimer.

Figure 3c Uinteract=UDFT-fdIC6/R6 for He dimer.
Figure 3 (a) UDFT as a function of internuclear separation for He dimer. All density functional potentials are purely repulsive; (b) Uinteract=UDFT+Uattract as a function of internuclear separation for He dimer. experimental* curve also plotted. None of the density functional potentials yield binding, Uattract=-fdI(R)C6/R6 (c) Enlargement of bonding region in (b).

For Ne dimer, Figure 2 shows that only the B3LYP functional potential has sufficiently little repulsiveness that binding can occur when the damped attractive term is incorporated.²⁶ However, the addition of the damped attractive term to the more repulsive BLYP and BPW91 functional potentials cannot yield binding. Indeed, as Figure 4 illustrates, even when the full, undamped $R^{-6}$ attraction, given in Eq. (4), is added to the density functional potentials for Ne₂, the more repulsive BLYP and BPW91 functionals cannot produce binding.

Figure 4a Uinteract=UDFT-C6/R6 for Ne Dimer.

Figure 4b Uinteract=UDFT-C6/R6 for Ne Dimer (Bonding Region).
Figure 4 Uinteract=UDFT+Uattract as a function of internuclear separation for Ne dimer. Even using the full, undamped attractive potential, only the B3LYP potential yields binding.; De=0.0020eV (experimental*=0.0038eV) Uattract=-C6/R6, (c) Enlargement of bonding region in (b).

$U_{{}_{attract}}\left(R\right)=-\frac{C_6}{R^6}$   (4)

For the least repulsive B3LYP potential, addition of the undamped $R^{-6}$ attraction actually gives a slight improvement in the dissociation energy compared to the highly accurate values of Oglivie & Wang.²⁷ (1.5 versus 2.0 compared to 3.5 meV). In the case of He dimer, Figure 3 shows that the B3LYP, BLYP, and BPW91 functionals exhibit potential curves that are so repulsive, the addition of the damped attractive term cannot induce binding. Again, as shown in Figure 5, even when the undamped $R^{-6}$ attraction is added to these density functional potentials, binding cannot be induced.

Figure 5a Uinteract=UDFT-C6/R6 for He dimer.

Figure 5b Uinteract=UDFT-C6/R6 for He Dimer
Figure 5 Uinteract=UDFT+Uattract as a function of internuclear separation for He dimer. Even using the full, undamped attractive potential, none of the density functional potentials yield binding. Experimental* potential also plotted. Uattract=-C6/R6, (c) Enlargement of bonding region in (b).

Since the addition of the full, undamped $R^{-6}$ attraction cannot induce binding for those density functional potentials which are more repulsive than $U_{B3LYP}$ of Ne dimer.²⁶ and the addition of the undamped $R^{-6}$ attraction improves the dissociation energy of the bound cases compared to that produced by the damped attractive term, then damping the attractive energy cannot provide a general solution to the problem of simulating dispersive interaction potentials. Two apparent solutions thus present themselves: we may either (a) strengthen the attractive term, or (b) damp the repulsive term.

It is not desirable to damp the repulsive density functional potentials, since the energy lowering would cause a reduction in the force constant and hence frequencies of bound vibrational modes. Indeed, if this energy were damped, the dispersion forces responsible for weak interactions could be introduced only at the expense of the covalent chemistry, the latter of which current density functional methods describe quite well. Therefore, the goal of this work will be to simulate the dispersive interaction potential, typified by the rare gases, through the addition of the density functional potential to a scaled attractive potential. Since the density functional methods studied in the present work are shown to be too repulsive to induce binding, as discussed above, the effect of this scaling will be to strengthen the attractive term.

If the objective is to scale the attractive potential, $U_{attract}$ in order to induce binding, one plausible approach is to simulate the behavior of a scaling function, that when multiplied by $\frac{C_6}{R^6}$ will result in an attractive potential capable of inducing binding, as given by Eq. (5), below.

$U_{\mathrm{interact}}=U_{DFT}-f_s\frac{C_6}{R^6}$   (5)

In general, the errors associated with the use of various density functional potentials could be quite different, resulting in perhaps drastically dissimilar scaling functions. It has been shown that interaction potentials obtained using the B3LYP density functional method are closer to the experimental potentials than those obtained with either the BLYP or BPW91 methods.²⁶ Therefore, the B3LYP density functional potential is an obvious choice for comparison with the known, experimental.²⁷ interaction potential in order to obtain the value of the scaling function, at any particular internuclear separation, as shown in Eq. (6)

$f_s{}^{p-p}\left(R\right)=\frac{U_{B3LYP}\left(R\right)-U_{exact}\left(R\right)}{\frac{C_6}{R^6}}$   (6)

Once the point–by–point values of the scaling function, f $f_s{}^{p-p}$ have been determined, it may be possible to derive an analytic function that will simulate its shape by considering the deviation from the $\frac{C_6}{R^6}$ potential arising from the dispersion interaction of two non–point, spherically distributed dipoles. This analytic function would depend on the position of the nuclei and a fitting parameter to the experimental potential. If the shape of the scaling function is relatively constant among atom pairs of varying size, then the simulated scaling function could potentially provide a viable means of obtaining the interaction potentials for other systems (ie: rare gas dimers) using the B3LYP hybrid density functional method. Since the scaling function is fit to the experimental potential through parameterization, its predictive utility is compromised; however, considering the interaction between two spherically distributed dipoles may prove insightful for developing a more general, predictive theory.

Two examples were selected from the literature in order to determine the appropriate behavior for the scaling function, f_s, using the B3LYP density functional potential. He and Ne dimers were obvious choices. For Ne dimer, addition of the damped attractive potential to the density functional potential binds only with the least repulsive B3LYP functional potential of Ne dimer.²⁶ For He dimer, this approach fails to give binding for any of the density functional potentials.²⁶

For both He₂ and Ne₂, single point energy calculations were performed using Becke’s three parameter hybrid density functional.^28,29 with the Lee–Yang–Parr correlation functional.³⁰ (B3LYP) density functional method between 2.10 to 4.80 angstroms at 0.02 angstrom intervals. The molecular energy at 20 angstroms was subtracted from each of these total molecular energies in order to obtain the B3LYP density functional potential, given in equation (2). This method is similar to that of Perez–Jorda et al.¹⁸ as well as that of Wu & Yang.²⁶ recently. The $C_6$ coefficients of Ne and He used in this work have been taken from the work of Kumar and Meath.^31–37 Although Wu and Yang use the cc–pV5Z basis set (augmented with diffuse functions for Ne₂) in their single point calculations, use of the cc–pVQZ basis set yielded density functional potentials reported to be in good agreement with those found with the pentuple zeta basis set.²⁶ As a timesaving measure, therefore, the quadruple zeta basis set (again, augmented with diffuse functions for Ne₂. Diffuse functions for the correlation consistent basis sets are not available at this time for He) was used in all single point energy calculations. Tight convergence criteria were used in the SCF calculations of Ne, since diffuse functions were used. A pruned (99590) grid was used for integration.²⁶ The Gaussian–03 program was used to perform all density functional calculations.³⁸

The point–by–point scaling functions, $f_s{}^{p-p}$ necessary to induce binding in Ne and He, given by equations (5) and (6) above, are presented in Figures 6a & 6b for Ne and He dimers, respectively.

Figure 6a Point-by-point Scaling Function (fs, p-p) as a function of internuclear distance for Ne dimer

Figure 6b 6b Point-by-Point Scaling Fcn (fs, p-p) as a function of Internuclear Separation for He dimer
Figure 6 (a) The point-by-point scaling function, fsp-p as a function of R/σ for Ne dimer, (b) The point-by-point scaling function, fsp-p as a function of R/σ for He dimer.

For each of the scaling functions, the $\frac{C_6}{R^6}$ attraction is attenuated for $\frac{R}{\sigma }<1$ , where is the atomic “hard–sphere” diameter. As the internuclear separation approaches the hard–sphere diameter, the scaling functions gradually rise to unity such that the $\frac{C_6}{R^6}$ term is no longer being attenuated when $R=\sigma$ . Thus, in this region, the scaling functions behave quite similarly to the damping functions used in previous studies.^24–26 However, in the region where $R>\sigma$ , the scaling functions each become larger than one, such that the $\frac{C_6}{R^6}$ attractive term is being strengthened; yet, purely damping functions cannot account for these effects. Finally, in the limit of infinite internuclear separation, the scaling functions for both Ne and He dimers smoothly decline to unity. Although Ne and He differ by eight electrons, the shape of the scaling functions necessary to induce binding remains roughly constant in these systems.

The oscillatory behavior of the scaling functions at large internuclear separations may be mere artifacts of the size of the grid used for integration. For example, when a smaller pruned (75302) grid was used for numerical integration, more oscillatory behaviour was observed. It seems reasonable that such unphysical behaviour could be eliminated with the use of larger grid sizes. Therefore, these oscillations were ignored when deriving an analytical function that simulates the curvature of the scaling function.

Unlike Gauss’s law, where Coulomb’s inverse square force is at work, the dispersion interaction between two atoms cannot be modeled as point dipoles centered at the nuclei. Indeed, the point–by–point scaling functions, $f_s{}^{p-p}$ , necessary to induce binding in Ne and He show that the $\frac{C_6}{R^6}$ attraction between two point charges must be strengthened in order to induce binding, indicating that the interaction distance is shorter than the internuclear distance. This strengthening of the $\frac{C_6}{R^6}$ attraction may be more important than the expansion in higher–order atomic polarizabilities used in traditional treatments.¹¹ Therefore, the potential energy function describing the dispersion interaction between two spherically distributed dipoles must be derived.

First consider the dispersion interaction between a spherically distributed dipole and a point dipole, as depicted in Diagram 1, below.

Here, the point dipole is located a distance τ away from the center of the sphere, and ${\rho }_1$ is the interaction distance between a point on the surface of the sphere, which is described in spherical–polar coordinates $\left(r_1,{\theta }_1,{\varphi }_1\right)$ and the point dipole. The potential energy between any point on the surface of the sphere and the point dipole is just the familiar $\frac{1}{R^6}$ attraction between two point dipoles. The component of this interaction that contributes to the dispersion attraction between the sphere and the point dipole must be in the direction of the axis adjoining them. However, since the potential energy itself is a scalar quantity, we cannot merely sum up the potential energies from all individual interactions between points on the sphere and the point dipole to obtain the total potential. Rather, we must sum up the electric fields, vector quantities, from all the individual interactions in order that the perpendicular components will cancel. To this end, consider Diagram 1, above. The electric field from the interaction of any point on the sphere’s surface and the point dipole has a magnitude that is the derivative of the potential between two point dipoles, as given by Eq. 7,

$\overrightarrow{\Vert E_1\Vert }=\frac{\partial V}{\partial {\rho }_1}=\frac{\partial \left(\frac{{\gamma }_6}{{\rho }^6}\right)}{\partial {\rho }_1}=\frac{-6{\gamma }_6}{{\rho }_1}$   (7)

Where ${\gamma }_6$ is a proportionality constant between the surface area of a sphere and the $C_6$ coefficient. Note here that Eq. (7) only holds for a point dipole that is outside the spherically distributed dipole. Eq. 8 gives the interaction distance, between the point on the surface of the sphere and the point dipole

${\rho }_{1^2}=r_1^{2}Sin{\theta }_1^2+{\left(\tau -r_{1}Cos{\theta }_1\right)}^2=r_1^2+{\tau }^2-2\tau r_{1}Cos{\theta }_1$   (8)

Note that the only point on the sphere’s surface that will produce an electric field totally in the direction of the axis adjoining the center of the sphere and the point dipole is the one directly in between them, at $(r_1,\theta =0,\varphi =1)$ . All other points on the sphere’s surface will produce an electric field that has a perpendicular component which does not contribute to the dispersion attraction between the sphere and the point dipole. However, rotating $\varphi$ from 0 to $2\pi$ for any fixed value of sweeps out a circle (outlined in yellow in Diagram 1 above) that is perpendicular to the axis adjoining the centre of the sphere and the point dipole. Due to the cylindrical symmetry, for every point on the circle that produces an electric field given by Eq. (9),

$${\overrightarrow{E}}_a=\left(\frac{-6{\gamma }_{6}Cos{\lambda }_1}{{\rho }_{{}_1^7}}\right){\overrightarrow{e}}_1+\left(\frac{-6{\gamma }_{6}Sin{\lambda }_1}{{\rho }_{{}_1^7}}\right){\overrightarrow{e}}_2$$   (9)

there exists a second point on the circle that produces an electric field given by Eq. (10), which cancels the component of the field perpendicular to the axis adjoining the sphere’s center and the point dipole.

${\overrightarrow{E}}_b=\left(\frac{-6{\gamma }_{6}Cos{\lambda }_1}{{\rho }_{{}_1^7}}\right){\overrightarrow{e}}_1+\left(\frac{-6{\gamma }_{6}Sin{\lambda }_1}{{\rho }_{{}_1^7}}\right){\overrightarrow{e}}_2$   (10)

Therefore, the average electric field felt by the point dipole from these two points is given by Eq. (11)

${\overrightarrow{E}}_1=\left(\frac{-6{\gamma }_{6}Cos{\lambda }_1}{{\rho }_1^7}\right)\overrightarrow{e_1}$   (11)

Substituting into Eq. (11) for the cosine of angle as presented in Diagram 1, and for the interaction distance ${\rho }_1$ , gives Eq. (12), below.

$\overrightarrow{E_1}=\frac{-6{\gamma }_{{}_6\left(\tau -r_{1}Cos\theta \right)}}{{\rho }_1^8}\overrightarrow{e_1}=\frac{-6{\gamma }_{{}_6\left(\tau -r_{1}Cos\theta \right)}}{{\left(r_1^2+{\tau }^2-2\tau r_{1}Cos{\theta }_1\right)}^4}\overrightarrow{e_1}$   (12)

Assuming a hollow sphere with an infinitely thin shell (such that there is no need to integrate radially), the total electric field produced by the entire sphere can be obtained by rotating from 0 to for all values of between 0 and , and is given by Eq. (13), below. Integration was performed with Mathematica.³⁹ The output is appended to this work (Appendix I).

$E_{total}=-6{\gamma }_6{\displaystyle \underset{{\varphi }_1=0}{\overset{2\pi }{\int }}{\displaystyle \underset{{\theta }_1=0}{\overset{\pi }{\int }}\frac{\left(\tau -r_{1}Cos{\theta }_1\right)}{{\left(r_1^2+{\tau }^2-2\tau r_{1}Cos{\theta }_1\right)}^4}}}\left(r_1^{2}Sin{\theta }_1\right)d{\varphi }_{1d{\theta }_1}=\frac{-{\gamma }_{6}8\pi r_1^2\tau \left(5r_1^2+3{\tau }^2\right)}{{\left(r_1-\tau \right)}^5{\left(r_1+\tau \right)}^5}$   (13)

The potential energy from the interaction between a spherically distributed dipole and a point dipole, then, is the integral of the total electric field with respect to the distance, between the center of the sphere and the point dipole, as given in Eq. (14), below. Again, integration was performed with Mathematica.³⁹ The output is appended to this work (see Appendix I).

$\begin{array}{l}{U}_{attract}(r_1,\tau )=-8\pi {\gamma }_6{\displaystyle \int \frac{r_1^2\tau (5r_1^2+3{\tau }^2)}{{(r_1-\tau )}^5{\left(r_1+\tau \right)}^5}}d\tau =\frac{-{\gamma }_{6}4\pi r_1^2}{{\left({\tau }^2-r_1^2\right)}^3}\left(\frac{{\tau }^2+r_1^2}{{\tau }^2-r_1^2}\right)\\ =\frac{-C_6}{{\left({\tau }^2-r_1^2\right)}^3}\left(\frac{{\tau }^2+r_1^2}{{\tau }^2-r_1^2}\right)\end{array}$   (14)

It is important to not that this potential energy function is only valid for a point dipole that is outside of a spherically distributed dipole, and does not apply to the region where $R\le r_1$ .

In order to derive the potential energy function for the dispersion interaction of two spherically distributed dipoles, let the point dipole of Diagram 1 be a point on the surface of a second sphere, as presented in Diagram 2, below.

Diagram 1 Dispersion interaction between a spherically distributed dipole and a point dipole.

Here, the internuclear distance is given by R, and the second interaction distance, ${\rho }_2$ is the distance between the center of one sphere and a point on the surface of the second sphere (note that ${\rho }_2$ is the same distance as in Diagram 1). Again, the points on the second sphere are described in spherical, polar coordinates $\left(r_2,{\theta }_2,{\varphi }_2\right)$ . From the previous discussion concerning the sphere and the point dipole, it was found that the electric field, labeled $E_1$ in Diagram 2, felt by any point on the surface of the second sphere due to interaction with the other sphere is given by Eq (13) above.

Diagram 2 Dispersion interaction between two spherically distributed dipoles.

Substituting the second interaction distance for , we have Eq. (15)

$\overrightarrow{E_1}=\frac{-{\gamma }_{6}8\pi r_1^2{\rho }_2\left(5r_1^2+3{\rho }_2^2\right)}{{\left(r_1-{\rho }_2\right)}^5{\left(r_1+{\rho }_2\right)}^5}{\overrightarrow{e}}_1$   (15)

This electric field vector can be written as the sum of two components, one parallel to the internuclear axis and the other perpendicular to it, as given by Eq. (16),

$\overrightarrow{E_1}=\overrightarrow{E_2}+\overrightarrow{E_{\perp }}$   (16)

where $E_2$ and $E_{\perp }$ are given by Eq. (17) and Eq. (18), respectively.

$\overrightarrow{E_2}=\left(Cos{\lambda }_2\right)\overrightarrow{E_1}$   (17)

$\overrightarrow{E_{\perp }}=(Sin{\lambda }_2)\overrightarrow{E_1}$   (18)

The only component of $E_1$ the that contributes to the dispersion attraction is the component in the direction of the internuclear axis $E_2$ . However, rotating ${\varphi }_1$ from 0 to $2\pi$ or any fixed value of ${\theta }_1$ sweeps out a circle (outlined in yellow in Diagram 2 above) that is perpendicular to the internuclear axis. As illustrated in Diagram 2, for any point on this circle that produces an electric field, $E_1$ given by Eq. (16), on the other sphere there exists another point on the circle that produces a field given by Eq. (19), where $E_2$ and $E_{\perp }$ are given by Eq. (17) and Eq. (18), respectively.

$\overrightarrow{E_1}=\overrightarrow{E_2}-\overrightarrow{E_{\perp }}$   (19)

Thus, the average electric field,$E_2$ from these two points is given by Eq. (17), above. As shown in Diagram 2, the cosine of angle is given by Eq. (20), below.

$Cos{\lambda }_2=\frac{R-r_{2}Cos{\theta }_2}{{\rho }_2}$   (20)

Substituting Eq. (20) in to Eq. (17) for the component of the electric field at any point on the surface of the second sphere contributing to the dispersion interaction, E₂ gives Eq. (21), below.

$\overrightarrow{E_2}=\frac{-{\gamma }_{6}8\pi r_1^2{\rho }_2\left(5r_1^2+3{\rho }_2^2\right)}{{\left(r_1-{\rho }_2\right)}^5{\left(r_1+{\rho }_2\right)}^5}\left(\frac{R-r_{2}Cos{\theta }_2}{{\rho }_2}\right)\overrightarrow{e_k}$   (21)

As before, assuming a hollow sphere with an infinitely thin shell (such that there is no need to integrate radially), the total electric field felt due to the interaction of the two spherically distributed dipoles can be obtained by rotating ${\varphi }_2$ from 0 to $2\pi$ for all values of ${\theta }_2$ between 0 and $\pi$ and is given by Eq. (22), below. Integration was performed with Mathematica.³⁹ The output is appended to this work (see Appendix I).

$\begin{array}{l}{E}_{total}=-{\gamma }_{6}8\pi {\displaystyle \underset{{\varphi }_2=0}{\overset{2\pi }{\int }}{\displaystyle \underset{{\theta }_2=0}{\overset{\pi }{\int }}\frac{r_1^2{\rho }_2\left(5r_1^2+3{\rho }_2^2\right)}{{\left(r_{1-{\rho }_2}\right)}^5{\left(r_{1+{\rho }_2}\right)}^5}}}\left(\frac{R-r_{2}Cos{\theta }_2}{{\rho }_2}\right)\left(r_2^{2}Sin{\theta }_2\right)d\varphi d{\varphi }_2\\ =16{\pi }^2{r}_1^2{r}_2^2\left(\begin{array}{l}\frac{15R^4+10R^2{r}_1^2-r_1^4+48R^3{r}_2+8R{r}_1^2{r}_2+54R^2{r}_2^2-2r_1^2{r}_2^2+24R{r}_2^3+3r_2^4}{(24R^2{r}_2{(R^2-r_1^2-2R{r}_2+r_2^2)}^4}-\\ \frac{15R^4+10R^2{r}_1^2-r_1^4-48R^3{r}_2-8R{r}_1^2{r}_2+54R^2{r}_2^2-2r_1^2{r}_2^2-24R{r}_2^3+3r_2^4}{(24R^2{r}_2{(R^2-r_1^2+2R{r}_2+r_2^2)}^4}\end{array}\right)\\ =C_6\left(\begin{array}{l}\frac{15R^4+10R^2{r}_1^2-r_1^4+48R^3{r}_2+8R{r}_1^2{r}_2+54R^2{r}_2^2-2r_1^2{r}_2^2+24R{r}_2^3+3r_2^4}{24R^2{r}_2{(R^2-r_1^2-2R{r}_2+r_2^2)}^4}-\\ \frac{15R^4+10R^2{r}_1^2-r_1^4-48R^3{r}_2-8R{r}_1^2{r}_2+54R^2{r}_2^2-2r_1^2{r}_2^2-24R{r}_2^3+3r_2^4}{(24R^2{r}_2{(R^2-r_1^2+2R{r}_2+r_2^2)}^4}\end{array}\right)\end{array}$   (22)

Finally, the potential energy from the interaction between two spherically distributed dipoles is the integral of the total electric field with respect to the internuclear distance, R, as given in Eq. (23), below. Again, this attractive potential is not valid when the spheres touch or as they unite. Integration was performed with Mathematica.³⁹ The output is appended to this work (Appendix I).

$U_{attract}=(r_1,r_2,R)={\displaystyle \int E_{total}}dR=\frac{C_6\left({\left[R^2-{(r_1-r_2)}^2\right]}^3+\frac{2}{3}\delta (r_1,r_2,R)\right)}{{\left[R^2-{(r_1-r_2)}^2\right]}^3{\left[R^2-{(r_1+r_2)}^2\right]}^3}$   (23)

This potential energy function, Eq. (23) is written as the sum of two terms, the first involving a simplified deviation from the $\frac{C_6}{R^6}$ attraction of two point–dipoles, and the second containing the error terms arising from this simplification, where δ(r₁,r₂,R) is given by Eq. (24), below.

$\begin{array}{l}\delta (r_1,r_2,R)=3R^4(r_1^2-3r_1{r}_2+r_2^2)+3{(}^{r_1}(r_1^4-r_1^3{r}_2+4r_1^2{r}_2^2-r_2^3{r}_1+r_2^4)-\\ -2R^2(3r_1^4-9r_1^3{r}_2+10r_1^2{r}_2^2-9r_2^3{r}_1+3r_2^4)\end{array}$   (24)

In the case of homonuclear diatomics, like Ne₂ and He₂, Eq. (23) reduces to Eq. (25), where δ(r,R) is given by Eq. (26)

$U_{attract}=\frac{C_6}{{\left[R^2-{(2r)}^2\right]}^3}+\frac{\frac{2}{3}{C}_6\delta (r,R)}{R^6{\left[R^2-{(2r)}^2\right]}^3}$   (25)

$\delta (r,R)=4R^2{r}^4-3R^4{r}^2$   (26)

The corresponding scaling function, $f_s^{sim}$ for strengthening the $\frac{C_6}{R^6}$ attraction of two point dipoles, is given by Eq. (25), where $\delta (r,R)$ is given by Eq. (26)

$U_{attract}=f_s^{sim}\frac{C_6}{R^6}\Rightarrow f_s^{sim}=\frac{R^6}{{\left[R^2-{\left(2r\right)}^2\right]}^3}+\frac{2/3\delta \left(r,R\right)}{{\left[R^2-{\left(2r\right)}^2\right]}^3}$   (27)

The behavior of the second, error term, above, is presented in Figure 7a for Ne dimer. Here, the radius of Ne is taken to be half its hard–sphere diameter.

Figure 7a Error_Term of Scaling Function of Eq. (27).

Figure 7bError_Term of Scaling Function of Eq. (27).
Figure 7 (a) Second term (Error Term) of Eq. (27), unphysical divergence as spheres touch; (b) Error term of Eq. (28), preventing unphysical divergence.

Note that the error term approaches negative infinity as the internuclear distance, R, approaches the sum of the radii for the two interacting Ne atoms. However, the potential energy should not diverge when the spheres touch or as they unite. This divergence can be prevented by adding a constant to the denominator of the error term in Eq. (27), to give Eq. (28), again where is given by Eq. (26)

$Error\_term=\frac{2/3\delta \left(r,R\right)}{{\left[R^2-{\left(2r\right)}^2\right]}^3+\gamma }$   (28)

By increasing the value of the constant, this error term becomes negligible for values of the internuclear distance surrounding the bonding region $\left(R=1.12\sigma \right)$ as shown in Figure 7b, below.

Therefore, it is reasonable to ignore this second, error term in the scaling function for the interaction of two spherically distributed dipoles. The simplified scaling function is given by Eq. (29), below.

$f_s^{sim}=\frac{R^6}{{\left[R^2-{\left(2r\right)}^2\right]}^3}$   (29)

Note that this simplified scaling function will diverge when the two spheres touch. To prevent this divergence, a constant can be added to the denominator, as shown in Eq. (30)

$f_s^{sim}=\frac{R^6}{{\left[R^2-{\left(2r\right)}^2\right]}^3+\gamma }$  (30)

By allowing the radius, r to be a multiple (1/n) of $\frac{\sigma }{2}$ , where is the hard–sphere diameter,$f_s^{sim}$ can be parameterized. The radius, r, is then given by Eq. (31)

$r=\frac{\sigma }{2n}$  (31)

The constant, Υ is parameterized by setting f_s^sim at R=σ equal to f_s^p–p at R=σ and solving Υ as shown in Eq. (32)

$\begin{array}{l}{f}_s^{sim}\left(R=\sigma \right)=f_s^{p-p}\left(R=\sigma \right)=\frac{{\sigma }^6}{{\left[{\sigma }^2-{\left(\frac{\sigma }{n}\right)}^2\right]}^3+\gamma }\Rightarrow {\sigma }^6\\ =f_s^{p-p}\left(R=\sigma \right){\left[{\sigma }^2-{\left(\frac{\sigma }{n}\right)}^2\right]}^3+\gamma \Rightarrow \gamma =\left[\frac{1}{f_s^{p-p}\left(R=\sigma \right)}-{\left(\frac{n^2-1}{n^2}\right)}^3\right]{\sigma }^6\end{array}$   (32)

Substituting Eq. (31) for the radius, r, and Eq. (32) for the constant, Υ in Eq. (30), gives Eq. (33), below.

$f_s^{sim}=\frac{R^6}{{\left[R^2-{\left(\frac{\sigma }{n}\right)}^2\right]}^3+\left[\frac{1}{f_s^{p-p}\left(R=\sigma \right)}-{\left(\frac{n^2-1}{n^2}\right)}^3\right]{\sigma }^6}$   (33)

The Levenberg–Marquardt method.³⁹ of non–linear least squares (please refer to the Fortran code presented in Appendix I) was used to obtain the optimized value for the parameter, n, by fitting Eq. (34), below, for R>σ to the experimental potentials.²⁷ for both He and Ne. These values for n were found to be just above unity, indicating that the optimized radii for He and Ne, as given by Eq. (31), are smaller than their corresponding van der Waals radii, as expected since most of the electron density is contained therein.

$U_{\mathrm{interact}}=U_{B3LYP}-f_s^{sim}\frac{C_6}{R_6}=U_{B3LYP}-\frac{C_6}{{\left[R^2-{\left(\frac{\sigma }{n}\right)}^2\right]}^3+\left[\frac{1}{f_s^{p-p}\left(R=\sigma \right)}-{\left(\frac{n^2-1}{n^2}\right)}^3\right]{\sigma }^6}$   (34)

The point–by point scaling functions, simulated scaling functions, as given by Eq. (33), above, and the damping functions employed in previous studies are presented in Figures 8a & 8b for Ne and He dimers, respectively.

Figure 8a Point-by-Point Scaling fcn (fs,p-p), Simulated Scaling Function (fs,sim), and Damping fcn (fdI) for Ne Dimer

Figure 8b Point-by-Point Scaling Fcn (fs,p-p), Simulated Scaling Fcn (fs,sim), and Damping Function (fdI) for He Dimer
Figure 8 (a) The point-by-point scaling function, fsp-p the simulated scaling function, fssim and Wu and Yang’s damping function I, fdI as functions for R>σ. n=1.35605931. (b) The point-by-point scaling function, fsp-p the simulated scaling function, fssim and Wu and Yang’s damping function I, fdI as functions for R>σ. n=1.16882509.

Figures 8a & 8b, above, reveal that Ne requires roughly 1.5 times more attraction in the bonding region $\left(R=1.12\sigma \right)$ in order to induce binding than purely damping functions can provide, while He requires nearly 2 times this amount of attraction. Figures 9 & 10 below, show the interaction potentials obtained using the simulated scaling functions, given by Eq. (33), compared to the experimental.²⁷ interaction potentials for Ne and He, respectively.

Figure 9a Uinteract=UB3LYP-fC6/R6 for Ne Dimer.

Figure 9b Uinteract=UB3LYP-fC6/R6 for Ne Dimer (Bonding Region).
Figure 9 Uinteract=UB3LYP-fssimC6/R6, Uinteract=UB3LYP-fdIC6/R6 and experimental* potentials for Ne dimer, (b) Enlargement of bonding region in (a) for R>σ. n=1.35605931.

Figure 10a Uinteract=UB3LYP-fC6/R6 for He Dimer.

Figure 10b Uinteract=UB3LYP-fC6/R6 for He Dimer (Bonding Region).
Figure 10 Uinteract=UB3LYP-fssimC6/R6, Uinteract=UB3LYP-fdIC6/R6, and experimental* potentials for He dimer, (b) Enlargement of bonding region in (a) for R>σ. n=1.16882509.

The equilibrium internuclear distances and dissociation energies obtained using the simulated scaling function, Wu and Yang’s damping function I, and the experimental.²⁷ values are summarized in Table 1 below.

 Although the dissociation energies obtained using the simulated scaling functions for Ne and He dimers are 0.4% (+0.000015 kcal/mol error) and 8% (+0.000074 kcal/mol error) in error, respectively, the dissociation energy of Ne dimer obtained with Wu and Yang’s damping function I is over 60% in error and this damping function cannot give binding for He dimer at all. The equilibrium geometries for Ne and He dimers obtained using the simulated scaling function are within 0.3% (10 mA) of the experimental bond lengths. Therefore, a scaling function of the type presented in Eq. (33) shows promise for improving interaction potentials compared to those obtained using a purely damping function for the B3LYP density functional method.

Although use of the simulated scaling function given in Eq. (33) has been shown to give significantly improved dissociation energies and equilibrium bond lengths (Table 1) compared to the use of a purely repulsive damping function, it is important to note that the model used to derive it not completely general. For example, while the optimized radii for He and Ne dimers are less than their respective van der Waals radii, as is expected since most of the electron density is contained therein, the optimized value of the radius for He, 1.095 A (n=1.16882509), is greater than that of Ne, 1.0136 A (n=1.35605931). Despite this discrepancy, however, the simplicity of the model presented herein affords field and potential integrals that are analytically tractable, and do provide improved dispersion potentials compared to experiment.²⁷

If the shape of the scaling function required to give binding is invariant to differing atom–atom pairs, as is suggested by the systems studied here, then it seems reasonable that an attractive scaling function like that presented in Eq. (33) could be useful for extending these ideas to other systems. In general, however, the exact interaction potentials are only known for the most trivial of chemical systems and since the potential energy function derived in this work must be fit to a known potential, its predictive value is compromised. The potential energy function derived herein is the solution to the six–dimensional volume integral over electron densities, as presented in Eq. (35),

$U_{attract}={\displaystyle \iint \frac{{\chi }_1\left(r_1\right){\chi }_2\left(r_2\right){\rho }_1\left(r_1\right){\rho }_2\left(r_2\right)}{R_6}}d{V}_{1}d{V}_2$   (35)

 for the simplified case in which the electron density, is distributed over the surface of a sphere. Here, $\chi \left(r_i\right)$ is the electric polarization response of the electron density at position i to an applied electric field, such that the atomic polarizability is given as Eq. (36),

${\alpha }_i\left(r_i\right)={\displaystyle \int {\chi }_i}\left(r_i\right)d{V}_i$   (36)

and the $C_6$ term is defined by Eq. (37)

$C_6=\frac{-3}{\pi }{\displaystyle \underset{0}{\overset{\infty }{\int }}{\alpha }_1}\left(r_{{}_1}\right){\alpha }_2\left(r_2\right)dr$   (37)

The next logical step towards extending these concepts to a complete, predictive theory is to solve the six–dimensional volume integral over electron densities for the general case where  is the electron density of an atom or molecular fragment.⁴⁰ Andersson et al.⁴¹ have recently shown that the long range $\frac{C_6}{R^6}$ attractive potential becomes Eq. (38),

$U_{attract}=\frac{6e}{4{\left(4\pi \right)}^{3/2}{m}^{1/2}}\frac{1}{{\left|r_1-r_2\right|}^6}{\displaystyle \underset{V_1}{\int }{d}^3}{r}_1{\displaystyle \underset{V_2}{\int }\frac{\sqrt{{\rho }_1\left(r_1\right){\rho }_2}\left(r_2\right)}{\sqrt{{\rho }_1\left(r_1\right)}+\sqrt{{\rho }_2\left(r_2\right)}}}{d}^3{r}_1$   (38)

where the electron densities, ${\rho }_i\left(r_i\right)$ have been extracted from Hartree–Fock wave functions with Slater basis sets, and leads to the calculation of significantly improved C₆ coefficients.⁴¹ The results of the simplified case presented in this work, in which the electron density is distributed over the surface of a sphere, indicate that the strengthening of the $\frac{C_6}{R^6}$ dispersion potential may be achieved by including the $\frac{1}{r_1-r_2}$ separation term in the integrand, in accord with Eq. (35). Thus, further studies toward this end are worthwhile.

A method for incorporating dispersion forces into density functional theory through the addition of a scaled atom–atom pairwise attractive potential term to the repulsive density functional potential is presented. By considering the deviation from the $\frac{C_6}{R^6}$ potential arising from the interaction of two spherical, non–point, dipoles, an analytic function was derived that when multiplied by $\frac{C_6}{R^6}$ resulted in an attractive potential capable of inducing binding. The scaling function found to give bound Ne and He dimers using the repulsive B3LYP density functional potential is quite different from the damping functions used in previous studies, causing the $\frac{C_6}{R^6}$ attraction to be strengthened in the bonding region. The attraction needed to give bound Ne and He dimers was found to be 1.5 and 2 times that which the damped attractive potential could provide, respectively. When the interaction potential for Ne dimer is constructed through the addition of the B3LYP density functional potential to the scaled attractive term, a dissociation energy 0.4% (+0.000015 kcal/mol error) in error is obtained, compared with 60% (–0.0023 kcal/mol error) error using the damped attraction. Similarly, for He dimer, use of the scaled attractive term results in a dissociation energy 8% (0.000074 kcal/mol error) in error; however, binding could not be induced using the damped attraction. These results, obtained by considering the dispersion interaction of two spherically distributed electron densities, are have proven insightful for developing a completely general density functional theory for predicting van der Waals interactions involving the dispersive interaction of electron distributions in atoms and molecular fragments. Through the introduction of dispersion forces, density functional theory, which has become an important tool for studies of covalent chemistry, may be extended to studies of weakly interacting systems. This model has been incorporated into the APF–D (Austin–Petersson–Frisch) Density Functional theory for dispersion forces within the Gaussian 09 Suite under keyword, APFD.⁴²

None.

None.

1. Koch W, Holthausen MC A Chemist’s Guide to Density Functional Theory. (2nd edn), Wiley–VCH Federal Republic of Germany, Germany, pp. 2001:313.
2. Hohenberg P, Kohn W. Inhomogeneous Electron Gas. Phys Rev. 1964;136:B864.
3. Kohn W, Sham L. Self–Consistent Equations Including Exchange and Correlation Effects. Phys Rev. 1965;140:A1133.
4. Kohn W, Becke AD, Parr RG. Density Functional Theory of Electronic Structure. J Phys Chem. 1996;100(31):12974–12980.
5. Ziegler T. Approximate density functional theory as a practical tool in molecular energetics and dynamics. Chem Rev. 1991;91(5):651–667.
6. Johnson BG, Gill PMW, Pople JA. The performance of a family of density functional methods. J Chem Phys. 1993;98:5612.
7. Perdew JP, Kurth S, Zupan A, et al. Accurate Density Functional with Correct Formal Properties: A Step Beyond the Generalized Gradient Approximation. Phys Rev Lett. 1999;82:2544.
8. Liangyou F, Ziegler T. Application of density functional theory to infrared absorption intensity calculations on main group molecules. J Chem Phys. 1992;96:9005.
9. Ernzenhof M, Scuseria GE. Kinetic energy density dependent approximations to the exchange energy. J Chem Phys. 1999;111:911–915.
10. Sim F, St Amant A, Papai I, et al. Gaussian density functional calculations on hydrogen–bonded systems. J Am Chem Soc. 1992;114(11):4391–4400.
11. Kauzmann Quantum Chemistry: An Introduction. Academic Press, Inc. New York, UK 1957.
12. Liu H, Elstner M, Kaxiras E, et al. Quantum mechanics simulation of protein dynamics on long timescale. Proteins. 2001;44(4):484–489.
13. Elstner M, Hobza P, Frauenheim T, et al. Hydrogen bonding and stacking interactions of nucleic acid base pairs: A density–functional–theory based treatment. J Chem Phys. 2001;114:5149.
14. Sponer J, Hobza P Nonplanar geometries of DNA bases. Ab initio second–order Moeller–Plesset study. J Phys Chem. 1994;98(12):3161–3164.
15. Hobza P, Sponer J, Polasek. H–Bonded and Stacked DNA Base Pairs: Cytosine Dimer. An Ab Initio Second–Order Moeller–Plesset Study. M J Am Chem Soc. 1995;117(2):792–798.
16. Sponer J, Leszczynski J, Hobza P. Nature of Nucleic Acid−Base Stacking: Nonempirical ab Initio and Empirical Potential Characterization of 10 Stacked Base Dimers. Comparison of Stacked and H–Bonded Base Pairs. J Phys Chem. 1996;100(13):5590–5596.
17. Kristyan S, Pulay P. Can (semi)local density functional theory account for the London dispersion forces? Chem Phys Lett. 1994;229(3):175–180.
18. Perez–Jorda JM, Becke AD. A density–functional study of van der Waals forces: rare gas diatomics. Chem Phys Lett. 1995;233(1–2):134–137.
19. Zhan Y, Pan W, Yang W. Describing van der Waals Interaction in diatomic molecules with generalized gradient approximations: The role of the exchange functional. J Chem Phys. 1997;107:7921.
20. Patton DC, Pederson MR. Application of the generalized–gradient approximation to rare–gas dimmers. Phys Rev. A56; R2495 1997.
21. Kohn W, Meir Y, Makarov DE. van der Waals Energies in Density Functional Theory. Phys Rev Lett. 80; 4153 1998:.
22. Wesolowski TA, Parisel O, Ellinger Y, et al. Comparative Study of Benzene···X (X = O2, N2, CO) Complexes Using Density Functional Theory:  The Importance of an Accurate Exchange−Correlation Energy Density at High Reduced Density Gradients. J Phys Chem. 1997;101:7818.
23. Dobson JF. Prospects for a van der Waals density functional. Int J Quantum Chem. 1998;69:615.
24. Ahlrichs R, Penco R, Scoles G. Intermolecular forces in simple systems. Chem Phys. 1977;19(2):119–130.
25. Mooij WTM, van Duijneveldt FB, van Duijneveld–van de Rijdt JGCM, van Eijck BP Effect of Methyl Ketones in the Belousov−Zhabotinskii Reaction. J Phys Chem A. 1999;103(48):9780–9872.
26. Wu Q, Yang W. Empirical correction to density functional theory for van der Waals interactions. J Chem Phys. 2002;116: 515.
27. Ogilvie JF, Wang FYH. Potential–energy functions of diatomic molecules of the noble gases: II. Unlike nuclear species. J Mol Struct. 1993;291(2–3):313–322.
28. Becke AD. Density–functional thermochemistry. III. The role of exact exchange. J Chem Phys. 1993;98(7):5648.
29. Becke AD. Density–functional exchange–energy approximation with correct asymptotic behaviour. Phys Rev A. 1988;38(6):3098.
30. Lee C, Yang W, Parr RG. Development of the Colle–Salvetti correlation–energy formula into a functional of the electron density. Phys Rev B. 1988;37(2):785.
31. Margoliash DJ, Meath WJ. Pseudospectral dipole oscillator strength distributions and some related two body interaction coefficients for H, He, Li, N, O, H2, N2, O2, NO, N2O, H2O, NH3, and CH4. J Chem Phys. 1978;68(4):1426.
32. Jhanwar BL, Meath W. Pseudo–spectral dipole oscillator strength distributions for the normal alkanes through octane and the evaluation of some related dipole–dipole and triple–dipole dispersion interaction energy coefficients. J Mol Phys. 1980;41(5):1061–1070.
33. Jhanwar BL, Meath W. Dipole oscillator strength distributions, sums, and dispersion energy coefficients for CO and CO 2. Chem Phys. 1982;67:185–199.
34. Kumar A, Meath W. Pseudo–spectral dipole oscillator strengths and dipole–dipole and triple–dipole dispersion energy coefficients for HF, HCl, HBr, He, Ne, Ar, Kr and Xe. J Mol Phys. 1985;54(4):823.
35. Kumar A, Meath W. Dipole oscillator strength properties and dispersion energies for acetylene and benzene. J Mol Phys. 1992;75(2):311–324.
36. Burton GR, Can WF, Cooper G, et al. Alence shell absolute photoabsorption oscillator–strengths, constrained dipole oscillator strength distributions, and dipole properties for CH3NH2, (CH3)(2)NH, AND (CH3)(3)N. J Chem. 1994;72(3):529–546.
37. Kumar A, Meath W. Isotropic dipole properties for acetone, acetaldehyde and formaldehyde. J Mol Phys. 1997;90(3):389–398.
38. Ayala K, Morokuma GA, Voth P, et al. Gaussian, Inc., Pittsburgh 2003.
39. Mathematica 4.1 for Windows
40. Brian P Flannery, Saul A Teukolsky, William T Vetterling. Numerical Receipts in Fortran (2nd Ed.), Cambridge University Press, England 1992.
41. Andersson Y, Langreth DC, Lundvist BI. van der Waals Interactions in Density–Functional Theory. Phys Rev Lett. 1992;76(1):102.
42. Austin A, Petersson G, Frisch. et al. A Density Functional with Spherical Atom Dispersion Terms. M J Chem Theory Comput. 2012;8(12):4989–5007.