Following the work of Alrichs et al. [24] on the correction of the HartreeFock model, as well as the work of Moorij et al. [25], Wu & Yang [26] introduced an empirical interaction term to the firstprinciples 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 NeNe 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_{2}, Ne_{2}, Ar_{2}, 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_{2}, Ar_{2}, 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_{2} and He_{2}, 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_{2}, 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=UDFTfdIC6/R6 for Ne dimer.
Figure 2c: Uinteract=UDFTfdIC6/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=UDFTfdIC6/R6 for He dimer.
Figure 3c: Uinteract=UDFTfdIC6/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=UDFTC6/R6 for Ne Dimer.
Figure 4b: Uinteract=UDFTC6/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=UDFTC6/R6 for He dimer.
Figure 5b: Uinteract=UDFTC6/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.
The pointbypoint scaling functions,
${f}_{s}{}^{pp}$
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 “hardsphere” diameter. As the internuclear separation approaches the hardsphere 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 [2426]. 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: Pointbypoint Scaling Function (fs, pp) as a function of internuclear distance for Ne dimer
Figure 6b: PointbyPoint Scaling Fcn (fs, pp) as a function of Internuclear Separation for He dimer
Figure 6: (a) The pointbypoint scaling function, f_{s}^{pp} as a function of R/σ for Ne dimer, (b) The pointbypoint scaling function, f_{s}^{pp} 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 pointbypoint scaling functions,
${f}_{s}{}^{pp}$
, 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 higherorder 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 sphericalpolar 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(5{r}_{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 (5{r}_{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(5{r}_{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(5{r}_{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(5{r}_{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{15{R}^{4}+10{R}^{2}{r}_{1}^{2}{r}_{1}^{4}+48{R}^{3}{r}_{2}+8R{r}_{1}^{2}{r}_{2}+54{R}^{2}{r}_{2}^{2}2{r}_{1}^{2}{r}_{2}^{2}+24R{r}_{2}^{3}+3{r}_{2}^{4}}{(24{R}^{2}{r}_{2}{({R}^{2}{r}_{1}^{2}2R{r}_{2}+{r}_{2}^{2})}^{4}}\\ \frac{15{R}^{4}+10{R}^{2}{r}_{1}^{2}{r}_{1}^{4}48{R}^{3}{r}_{2}8R{r}_{1}^{2}{r}_{2}+54{R}^{2}{r}_{2}^{2}2{r}_{1}^{2}{r}_{2}^{2}24R{r}_{2}^{3}+3{r}_{2}^{4}}{(24{R}^{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{15{R}^{4}+10{R}^{2}{r}_{1}^{2}{r}_{1}^{4}+48{R}^{3}{r}_{2}+8R{r}_{1}^{2}{r}_{2}+54{R}^{2}{r}_{2}^{2}2{r}_{1}^{2}{r}_{2}^{2}+24R{r}_{2}^{3}+3{r}_{2}^{4}}{24{R}^{2}{r}_{2}{({R}^{2}{r}_{1}^{2}2R{r}_{2}+{r}_{2}^{2})}^{4}}\\ \frac{15{R}^{4}+10{R}^{2}{r}_{1}^{2}{r}_{1}^{4}48{R}^{3}{r}_{2}8R{r}_{1}^{2}{r}_{2}+54{R}^{2}{r}_{2}^{2}2{r}_{1}^{2}{r}_{2}^{2}24R{r}_{2}^{3}+3{r}_{2}^{4}}{(24{R}^{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 pointdipoles, and the second containing the error terms arising from this simplification, where δ(r_{1},r_{2},R) is given by Eq. (24), below.
$\begin{array}{l}\delta ({r}_{1},{r}_{2},R)=3{R}^{4}({r}_{1}^{2}3{r}_{1}{r}_{2}+{r}_{2}^{2})+3{(}^{{r}_{1}}({r}_{1}^{4}{r}_{1}^{3}{r}_{2}+4{r}_{1}^{2}{r}_{2}^{2}{r}_{2}^{3}{r}_{1}+{r}_{2}^{4})\\ 2{R}^{2}(3{r}_{1}^{4}9{r}_{1}^{3}{r}_{2}+10{r}_{1}^{2}{r}_{2}^{2}9{r}_{2}^{3}{r}_{1}+3{r}_{2}^{4})\end{array}$
(24)
In the case of homonuclear diatomics, like Ne_{2} and He_{2}, 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)=4{R}^{2}{r}^{4}3{R}^{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 hardsphere 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 hardsphere 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}^{pp} at R=σ and solving Υ as shown in Eq. (32)
${f}_{s}^{sim}\left(R=\sigma \right)={f}_{s}^{pp}\left(R=\sigma \right)=\frac{{\sigma}^{6}}{{\left[{\sigma}^{2}{\left(\frac{\sigma}{n}\right)}^{2}\right]}^{3}+\gamma}\Rightarrow {\sigma}^{6}={f}_{s}^{pp}\left(R=\sigma \right){\left[{\sigma}^{2}{\left(\frac{\sigma}{n}\right)}^{2}\right]}^{3}+\gamma \Rightarrow \gamma =\left[\frac{1}{{f}_{s}^{pp}\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}^{pp}\left(R=\sigma \right)}{\left(\frac{{n}^{2}1}{{n}^{2}}\right)}^{3}\right]{\sigma}^{6}}$
(33)
The LevenbergMarquardt method [39] of nonlinear 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}^{pp}\left(R=\sigma \right)}{\left(\frac{{n}^{2}1}{{n}^{2}}\right)}^{3}\right]{\sigma}^{6}}$
(34)
The pointby 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: PointbyPoint Scaling fcn (fs,pp), Simulated Scaling Function (fs,sim), and Damping fcn (fdI) for Ne Dimer
Figure 8b: PointbyPoint Scaling Fcn (fs,pp), Simulated Scaling Fcn (fs,sim), and Damping Function (fdI) for He Dimer
Figure 8: (a) The pointbypoint scaling function, f_{s}^{pp} 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 pointbypoint scaling function, f_{s}^{pp} 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=UB3LYPfC6/R6 for Ne Dimer.
Figure 9b: Uinteract=UB3LYPfC6/R6 for Ne Dimer (Bonding Region).
Figure 9: Uinteract=UB3LYPf_{s}^{sim}C6/R6, Uinteract=UB3LYPfdIC6/R6 and experimental* potentials for Ne dimer, (b) Enlargement of bonding region in (a) for R>σ. n=1.35605931.
Figure 10a: Uinteract=UB3LYPfC6/R6 for He Dimer.
Figure 10b: Uinteract=UB3LYPfC6/R6 for He Dimer (Bonding Region).
Figure 10: Uinteract=UB3LYPf_{s}^{sim}C6/R6, Uinteract=UB3LYPfdIC6/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].


Re (A) 
De (kcal/mol) 
% Error in Re 
% Error in De 
Ne 
Scaling attraction 
3.08 
0.003815 
0.3 
0.4 
Damping attraction 
3.00 
0.0015 
2.9 
60.5 
Experimental 
3.09 
0.0038 


He 
Scaling attraction 
2.96 
0.001019 
0.0 
7.8 
Damping attraction 
No Binding 
No Binding 
No Binding 
No Binding 
Experimental 
2.96 
0.000945 


Table 1: Equilibrium geometry, Re and dissociation energy, De values obtained from addition of the scaled attractive potential to the repulsive B3LYP density functional potential for Ne and He dimers, according to the equation:
${U}_{\mathrm{interact}}={U}_{B3LYP}{f}_{s}^{sim}\frac{{C}_{6}}{{R}^{6}}$
(Figures 9 &10 for graphs).
If the shape of the scaling function required to give binding is invariant to differing atomatom 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 sixdimensional 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 sixdimensional 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 HartreeFock wave functions with Slater basis sets, and leads to the calculation of significantly improved C_{6} 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.