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 2nd 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 [11], 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 [1]. 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 [17]. Although some combinations of exchange and correlation functionals do produce overly bound rare-gas dimers, [18] 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. [24] on the correction of the Hartree-Fock model, as well as the work of Moorij et al. [25], Wu & Yang [26] 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 [26]. 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 [26].

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.

$$\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}$$

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

$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 [26]. 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 [26].

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 [26]. 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 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 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).

$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 [27] (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.

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 [26] 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.

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).

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).

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).

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, $f_s$ 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 [26]. Therefore, the B3LYP density functional potential is an obvious choice for comparison with the known, experimental [27] interaction potential in order to obtain the value of the scaling function, $C_6$ 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 [26]. For He dimer, this approach fails to give binding for any of the density functional potentials [26].

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 [30] (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. [18] as well as that of Wu & Yang [26] recently. The 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 [26]. 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 [26]. The Gaussian-03 program was used to perform all density functional calculations [38].

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.

For each of the scaling functions, the $\frac{C_6}{R^6}$ attraction is attenuated for $\frac{R}{\sigma }<1$ , where $\sigma$ 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.

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

Figure 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, f_s^p-p as a function of R/σ for Ne dimer, (b) The point-by-point scaling function, f_s^p-p as a function of R/σ for He dimer.

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

Diagram 2: Dispersion interaction between two spherically distributed dipoles.

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 [11]. 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, ${\rho }_1$ 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 $\theta$ 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 ${\lambda }_1$ 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 $\varphi$ from 0 to $2\pi$ for all values of $\theta$ between 0 and $\pi$ , and is given by Eq. (13), below. Integration was performed with Mathematica [39]. 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, $\tau$ between the center of the sphere and the point dipole, as given in Eq. (14), below. Again, integration was performed with Mathematica [39]. 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.

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 $\tau$ 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.

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 the $E_1$ 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 $\lambda$ 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, E2 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 [39]. 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 [39]. 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.

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 $\delta \left(r,R\right)$ 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)

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

Figure 7b: Error_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.

By increasing the value of the constant, $\gamma$ 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 $\sigma$ is the hard-sphere diameter, $f_s^{sim}$ can be parameterized. The radius, r, is then given by Eq.

$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)

$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$ (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 [39] 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 [27] 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.

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 [27] interaction potentials for Ne and He, 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, f_s^p-p the simulated scaling function, f_s^sim and Wu and Yang’s damping function I, fdI as functions for R>σ. n=1.35605931.
(b) The point-by-point scaling function, f_s^p-p the simulated scaling function, f_s^sim and Wu and Yang’s damping function I, fdI as functions for R>σ. n=1.16882509.

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

Figure 9b: Uinteract=UB3LYP-fC6/R6 for Ne Dimer (Bonding Region).
Figure 9: Uinteract=UB3LYP-f_s^simC6/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-f_s^simC6/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 [27] 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 [27].

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, $\rho$ 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 $\rho$  is the electron density of an atom or molecular fragment [40]. Andersson et al. [41] 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 [41]. 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 [42].