Durham Research Online
Deposited in DRO:
26 April 2011
Version of attached le:
Published Version
Peer-review status of attached le:
Peer-reviewed
Citation for published item:
Allen, M. J. and Tozer, D. J. (2002) 'Helium dimer dispersion forces and correlation potentials in density
functional theory.', Journal of chemical physics., 117 (24). pp. 11113-11120.
Further information on publisher's website:
http://dx.doi.org/10.1063/1.1522715
Publisher's copyright statement:
Copyright (2002) American Institute of Physics. This article may be downloaded for personal use only. Any other use
requires prior permission of the author and the American Institute of Physics. Allen, M. J. and Tozer, D. J. (2002)
'Helium dimer dispersion forces and correlation potentials in density functional theory.', Journal of chemical physics.,
117 (24). pp. 11113-11120. and may be found at http://dx.doi.org/10.1063/1.1522715
Additional information:
Use policy
The full-text may be used and/or reproduced, and given to third parties in any format or medium, without prior permission or charge, for
personal research or study, educational, or not-for-prot purposes provided that:
•
a full bibliographic reference is made to the original source
•
a link is made to the metadata record in DRO
•
the full-text is not changed in any way
The full-text must not be sold in any format or medium without the formal permission of the copyright holders.
Please consult the full DRO policy for further details.
Durham University Library, Stockton Road, Durham DH1 3LY, United Kingdom
Tel : +44 (0)191 334 3042 | Fax : +44 (0)191 334 2971
https://dro.dur.ac.uk
Helium dimer dispersion forces and correlation potentials in density
functional theory
Mark J. Allen and David J. Tozer
a)
Department of Chemistry, University of Durham, South Road, Durham DH1 3LE, United Kingdom
共Received 18 July 2002; accepted 30 September 2002兲
The dispersion interaction in the helium dimer is considered from the viewpoint of the force on a
nucleus. At large internuclear separations, Brueckner coupled cluster BD共T兲 forces agree well with
near-exact dispersion forces. The atomic density distortion associated with the dispersion force is
quantified by comparing the BD共T兲 dimer density with a superposition of atomic densities. For
density functional theory calculations in the Hartree–Fock–Kohn–Sham 共HFKS兲 formalism, the
accuracy of the dispersion force is governed by the correlation potential. Calculations using the
conventional Lee–Yang–Parr 关Phys. Rev. B 37, 785 共1988兲兴 potential only generate a small density
distortion, giving forces significantly smaller than BD共T兲. The BD共T兲 electron densities are
therefore used to determine improved correlation potentials using a modified Zhao–Morrison–Parr
共ZMP兲 approach 关Phys. Rev. A 50, 2138 共1994兲兴. HFKS calculations using these ZMP potentials
quantitatively reproduce the distortion, giving dispersion forces in good agreement with BD共T兲. The
dimer ZMP correlation potential is partitioned into two parts, one equal to the sum of two
unperturbed spherical atomic correlation potentials and the other representing an interaction
potential. HFKS calculations using the former do not generate the distortion; forces are close to
Hartree–Fock. Calculations using the latter do generate the distortion, giving forces essentially
identical to those from the full dimer potential. The origin of the distortion is traced to the
asymmetric structure of the interaction correlation potential in the vicinity of each nucleus. © 2002
American Institute of Physics. 关DOI: 10.1063/1.1522715兴
I. INTRODUCTION
The description of van der Waals interactions is a major
challenge for density functional theory 共DFT兲 approxima-
tions. Calculations using conventional exchange-correlation
functionals have been performed for a range of systems, in-
cluding rare gas dimers,
1–10
C
6
H
6
dimer,
3,7,9,11
CH
4
and
C
2
H
2
dimers,
9,12
He¯ CO
2
,
13,14
N
2
dimer,
15
C
6
H
6
¯ X
关
X
⫽ O
2
,N
2
,CO共Ref. 16兲, Ne, Ar 共Ref. 3兲兴, and other non-
bonded dimeric complexes.
17
Common conclusions are
reached. The local density approximation 共LDA兲 tends to
overbind
1,2,4,5,16
while the performance of generalized gradi-
ent approximation 共GGA兲 and hybrid functionals is sensitive
to the choice of exchange approximation. Functionals based
on Becke 1988 exchange
18
often predict a repulsive
interaction;
1–3,5–9,11–13,15–17
those based on PW91 共Ref. 19兲
or PBE 共Ref. 20兲 exchange do tend to bind, although quan-
titative accuracy is lacking.
4–9,14–16
This sensitivity to the
exchange functional has been attributed
5,16
to the behavior of
the exchange enhancement factor at large reduced density
gradient s; the Becke 1988 enhancement factor diverges,
whereas the PW91 and PBE factors are better behaved. In
our preliminary studies, we obtained results consistent with
this assessment. The HCTH93 共Ref. 21兲 exchange-
correlation functional, whose exchange enhancement factor
increases rapidly with s, does not bind the helium dimer. The
1/4 functional,
22
whose enhancement factor increases more
gradually,
23
does bind. For a recent review of van der Waals
studies using conventional functionals, see Ref. 24.
A DFT calculation using an appropriately chosen con-
ventional functional can therefore provide a qualitative de-
scription of van der Waals systems at intermediate separa-
tions, where there is a non-negligible overlap between the
interacting fragments and the interaction energy is composed
of several terms 共dispersion, exchange-dispersion, electro-
static, exchange-repulsion, etc.兲. At larger separations, how-
ever, where overlap is negligible, the interaction energy is
dominated by the long-range dispersion energy, arising from
correlated interactions between electrons on the separate
fragments. The local nature of conventional functionals
means they are fundamentally unable to describe this feature,
failing to recover the leading ⫺ C
6
R
⫺ 6
interaction energy.
Although this term can be introduced in an empirical
manner,
25
more advanced methods must be used to introduce
it rigorously. These include long-range
26–33
and
seamless
34–38
approaches; see Ref. 39 for an assessment of
some of these methods. Kohn–Sham orbitals have also been
used within symmetry-adapted perturbation theory.
40,41
In this study we consider the long-range dispersion in-
teraction in DFT from the viewpoint of the force on a
nucleus, rather than from the usual viewpoint of the elec-
tronic energy. Given that dispersion is a correlation effect,
we treat exchange exactly using the Hartree–Fock–Kohn–
Sham 共HFKS兲 formalism.
42,43
The electronic energy is
written
E
DFT
⫽ E
HF
关
兵
i
其
兴
⫹ E
C
关
兴
, 共1兲
a兲
Fax: 0191 384 4737; Electronic mail: D.J.Tozer@Durham.ac.uk
JOURNAL OF CHEMICAL PHYSICS VOLUME 117, NUMBER 24 22 DECEMBER 2002
111130021-9606/2002/117(24)/11113/8/$19.00 © 2002 American Institute of Physics
Downloaded 26 Apr 2011 to 129.234.252.65. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
where E
HF
关
兵
i
其
兴
is the Hartree–Fock functional, E
C
关
兴
is
an approximate correlation energy functional, and
(r) is the
electron density
共
r
兲
⫽
兺
i
i
2
共
r
兲
. 共2兲
Expansion of the orbitals
兵
i
其
in a basis set
兵

其
allows the
HFKS equations
冕
F
HF
共
r,r
⬘
兲
i
共
r
⬘
兲
dr
⬘
⫹
v
C
共
r
兲
i
共
r
兲
⫺
⑀
i
i
共
r
兲
⫽ 0 共3兲
to be recast as secular equations
兺

冕
␣
共
r
兲
冋
冕
F
HF
共
r,r
⬘
兲

共
r
⬘
兲
dr
⬘
⫹
v
C
共
r
兲

共
r
兲
⫺
⑀
i

共
r
兲
册
drC

i
⫽ 0, 共4兲
where F
HF
(r,r
⬘
) is the coordinate representation of the
Hartree–Fock operator 共nonmultiplicative due to orbital ex-
change兲 and
v
C
共
r
兲
⫽
␦
E
C
关
兴
␦
共
r
兲
共5兲
is the correlation potential. The HFKS force on nucleus A is
then
F
DFT
⫽⫺
E
DFT
R
A
⫽⫺E
HF
R
A
关
兵
i
其
兴
⫺
冕
R
A
共
r
兲
v
C
共
r
兲
dr
⫹
兺
i
⑀
i
S
ii
R
A
, 共6兲
where E
HF
R
A
关
兵
i
其
兴
,
R
A
(r), and S
ii
R
A
are the basis-function-
only derivatives of the Hartree–Fock functional, density, and
orbital overlap matrix, respectively, with respect to the
nuclear coordinate vector R
A
. Other than
v
C
(r), all the
quantities in Eq. 共6兲 can be constructed from the solutions to
Eq. 共4兲. Given that
v
C
(r) is the only approximated term in
Eq. 共4兲, it follows that the quality of this potential alone
determines the quality of the force for a given basis set. In
essence, E
C
关
兴
governs the accuracy of the total energy, so
its functional derivative governs the energy derivative. This
dependence on
v
C
(r) is particularly evident when the basis
set is complete, since Eq. 共6兲 reduces to the Hellmann–
Feynman force.
44
For a given Born–Oppenheimer configu-
ration, this depends only on the density, whose accuracy is
governed by
v
C
(r) through Eqs. 共2兲 and 共3兲.
At large separation, the force on a nucleus in a van der
Waals molecule is almost exclusively due to the dispersion
interaction. To describe this dispersion force accurately
within the HFKS formalism therefore requires an accurate
representation of
v
C
(r) at large separation; integration of this
force along the dissociation path yields the dispersion inter-
action energy. We regard
v
C
(r) as a key quantity, containing
essential physics of DFT dispersion. The aim of this study is
to use ab initio electron densities to learn about the structure
of
v
C
(r) and its relationship to dispersion forces in the he-
lium dimer He
2
.
We commence in Sec. II by providing computational de-
tails and choosing internuclear separations in He
2
where dis-
persion dominates. The physical origin of the dispersion
force—a distortion of the atomic densities—is discussed and
quantified using BD共T兲 densities. Deficiencies with conven-
tional correlation potentials are highlighted by considering
HFKS forces using the Lee–Yang–Parr 共LYP兲 potential.
45
Correlation potentials are then determined directly from the
BD共T兲 electron densities, using a modified Zhao–Morrison–
Parr 共ZMP兲共Ref. 46兲 approach. Self-consistent HFKS calcu-
lations are performed using these potentials and forces are
determined. The partitioning of the correlation potential into
atomic and interaction components is investigated. Conclu-
sions are presented in Sec. III.
II. RESULTS
A. Computational details
All calculations were performed using a modified ver-
sion of the
CADPAC program
47
with an extensive 7s5p4d
basis set on the He atoms, corresponding to the nuclear cen-
tred part of the DC
⫹
BS 共Dc147兲 basis set of Ref. 48, with
the f functions removed for technical reasons. Unless other-
wise stated, the BD共T兲, Hartree–Fock, and DFT forces were
all evaluated analytically, using conventional rigorous energy
derivative expressions. Where possible, numerical stability
was confirmed by comparing the analytic forces with nu-
merical forces determined from energies at perturbed geom-
etries. Basis set superposition errors 共BSSE兲 affect the shape
of the interaction energy curve and so also affect the calcu-
lated forces. All forces were corrected for BSSE by differen-
tiating the counterpoise energy correction. This requires the
force on a single helium atom, calculated in the presence of
additional ghost atom basis functions. For DFT calculations,
the integration grid on the ghost atom was also included, in
order to account for integration grid superposition error.
Given our extensive basis set, large internuclear separations,
and near-saturated integration grids, the BSSE corrections to
the total forces are very small. To the number of decimal
places quoted they are negligible for all methods except
BD共T兲关where it contributes 0.1⫻ 10
⫺ 6
a.u. 共about 2%兲 to
the forces at 8.0 and 8.5 a.u.兴. BSSE corrections to
Hellmann–Feynman forces were slightly larger.
All BD共T兲 densities are relaxed densities. HFKS corre-
lation potentials were determined from these densities using
the methodology of Refs. 49 and 50, which is a modification
of the ZMP approach,
46
based on the constrained search
formulation.
51
The method is as follows. The only terms in
Eq. 共1兲 that are not explicit functionals of the density are the
noninteracting kinetic and orbital exchange energies in
E
HF
关
兵
i
其
兴
. Minimization of the sum of these two terms with
respect to the orbitals, subject to the constraint that the den-
sity equals the BD共T兲 density, gives the HFKS orbitals asso-
ciated with that density. The addition of appropriate explicit
density functionals to the minimization ensures that the re-
sulting one-electron equations take the same form as Eq. 共3兲,
but does not change the solution. Comparison of the equa-
tions then allows
v
C
(r) to be identified in terms of the BD共T兲
density, the iterating density, and a Lagrange multi-
11114 J. Chem. Phys., Vol. 117, No. 24, 22 December 2002 M. J. Allen and D. J. Tozer
Downloaded 26 Apr 2011 to 129.234.252.65. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
plier associated with the density constraint. The one-
electron equations are solved within a basis set framework
and the potential is tabulated numerically on a DFT numeri-
cal integration grid.
The Lagrange multiplier is formally infinite, although
for finite basis sets the BD共T兲 densities cannot be reproduced
exactly and so this is inappropriate.
52
We considered three
approaches to circumvent this problem. We considered a fi-
nite value of ⫽ 900, as was used in Ref. 49 and in studies
of exchange-correlation potentials. We also considered two
extrapolation schemes. The first, from Ref. 53, involves an
expansion from
⫺ 3
to
⫹ 1
, where the latter term represents
basis set incompleteness. The second extrapolation, similar
to that in Ref. 46, replaces this basis set term with
⫺ 4
.To
choose an optimal scheme, we used the fact that the ZMP
iterating density should equal the BD共T兲 density, and so
Hellmann–Feynman forces from the two should be identical.
Comparison of the forces led us to conclude that the first
extrapolation scheme did not work well; a similar conclusion
was reached in Ref. 49. Forces from ⫽ 900 were close to
those from the second extrapolation scheme and so for sim-
plicity we use a finite value of ⫽ 900 throughout.
In all calculations, He
2
is oriented along the z axis. The
nuclei, which are labeled A and B to distinguish them, are
positioned at z coordinates z
A
and z
B
, respectively, where
z
B
⬎ z
A
and the internuclear separation is R⫽ z
B
⫺ z
A
. All
quoted forces correspond to the force on nucleus A. This acts
along the z axis and is equal and opposite to the force on
nucleus B. 共Forces constructed using approximate ZMP po-
tentials do not generally satisfy this translational invariance
condition because the potentials are not exact functional de-
rivatives. The homonuclear nature of He
2
ensures that this
condition is satisfied in the present study.兲 A positive force
pulls nucleus A towards nucleus B and so represents an at-
traction; a negative force is a repulsion.
B. Dispersion forces and the atomic density
distortion
Our first task is to choose a set of internuclear separa-
tions R, where the interaction is dominated by dispersion.
Korona et al.
48
fitted an accurate symmetry-adapted pertur-
bation theory 共SAPT兲 interaction energy for He
2
to the form
E
SAPT
⫽ Ae
⫺
␣
R⫹

R
2
⫺
兺
n⫽ 3
8
f
2n
共
R,b
兲
C
2n
R
2n
, 共7兲
where all parameters are defined in Ref. 48. At large R, this
approaches the long-range dispersion energy
E
disp
⫽⫺
兺
n⫽ 3
8
C
2n
R
2n
⫽⫺
C
6
R
6
⫺
C
8
R
8
⫺ ¯ ⫺
C
16
R
16
. 共8兲
The force on nucleus A is
F
SAPT
⫽⫺
E
SAPT
z
A
⫽
E
SAPT
R
, 共9兲
which at large R approaches the long-range dispersion force
F
disp
⫽⫺
E
disp
z
A
⫽
E
disp
R
⫽
兺
n⫽ 3
8
2nC
2n
R
2n⫹1
⫽
6C
6
R
7
⫹
8C
8
R
9
⫹ ¯ ⫹
16C
16
R
17
. 共10兲
Figure 1 presents E
SAPT
and F
SAPT
as a function of R. Here
E
disp
and F
disp
are also presented in order to estimate the
distance where the SAPT terms reduce to these limiting
long-range forms. The curves become indistinguishable be-
yond R⫽ 7.5 a.u. and so we shall concentrate on R⫽ 8.0, 8.5,
and 9.0 a.u. Table I presents F
disp
and F
SAPT
at the three R
values. The small difference between the two forces reflects
the nonvanishing atomic overlap. We regard F
SAPT
as near-
exact reference forces.
Before presenting forces from approximate electronic
structure methods, it is informative to consider the physical
FIG. 1. SAPT interaction energy E
SAPT
and force F
SAPT
, together with
long-range dispersion contributions E
disp
and F
disp
for He
2
.
TABLE I. The force on nucleus A in He
2
, in units of ⫻ 10
⫺ 6
a.u., for internuclear separations R⫽8.0, 8.5, and 9.0 a.u. All forces act along the He–He bond
axis. A positive force represents an attraction between the nuclei.
R
F
disp
F
SAPT
F
HF
F
BD共T兲
F
DFT
关v
C,LYP
dimer
兴
F
DFT
关v
C,ZMP
dimer
兴
F
DFT
关v
C,ZMP
atoms
兴
F
DFT
关v
C,ZMP
int
兴
8.0 5.4 5.1 ⫺ 0.2 4.9 0.9 4.8 ⫺ 0.3 4.8
8.5 3.4 3.3 ⫺ 0.1 3.1 0.3 3.2 ⫺ 0.1 3.2
9.0 2.2 2.2 ⫺ 0.0 2.1 0.1 2.1 ⫺ 0.0 2.1
11115
J. Chem. Phys., Vol. 117, No. 24, 22 December 2002 He
2
dispersion forces
Downloaded 26 Apr 2011 to 129.234.252.65. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
origin of the dispersion force in the exact He
2
. The electro-
static theorem of Feynman
44
—obtained by applying the dif-
ferential Hellmann–Feynman theorem to a nuclear
perturbation—states that the force on a nucleus is just the
classical electrostatic force exerted by the other nuclei and
the electron density. For He
2
, the exact force on nucleus A is
then
F⫽⫺
4
R
2
⫹ 2
冕
共
r
1
兲
r
A1
3
z
A1
dr
1
, 共11兲
where the first term is the repulsive force due to nucleus B
and the second is the force due to the exact density. At large
R, correlation between electrons on the two atoms causes the
atomic densities to be distorted towards one another. The
force due to the density then pulls nucleus A in the direction
of nucleus B more strongly than they repel one another, re-
sulting in a net attractive force. Feynman described this dis-
tortion:
‘‘The Schro
¨
dinger perturbation theory for two interacting
atoms at a separation R, large compared to the radii of the
atoms, leads to the result that the charge distribution of each
is distorted from central symmetry, a dipole moment of order
1/R
7
being induced in each atom. The negative charge dis-
tribution of each atom has its center of gravity moved
slightly toward the other.’’
Feynman also conjectured that the leading term in the
dispersion force on a nucleus arose entirely from the attrac-
tion of that nucleus to its own distorted density. For He
2
, this
implies that the R
⫺ 7
component of the force on nucleus A
arises from atom A’s density contribution to the second term
in Eq. 共11兲; this density is polarized towards atom B and so
pulls nucleus A in that direction. Although atom B causes the
dispersion force on nucleus A, its density and nucleus only
explicitly contribute to the higher-order forces. Feynman’s
conjecture was verified by Hirschfelder and Eliason
54
for two
hydrogen atoms. A general proof has been provided by
Hunt;
55
see Refs. 56–58 for further discussion.
In an approximate electronic structure calculation the
electrostatic theorem does not hold exactly, due to finite basis
sets or nonvariational methodology. Nevertheless, we do find
that all our Hartree–Fock, BD共T兲, and DFT forces are very
close to the Hellmann–Feynman forces 共11兲 calculated using
their respective densities, demonstrating that the forces in
these methods can be rationalized in terms of simple electro-
statics; the force reflects the density distortion produced by
the method. We shall quantify the distortion using the func-
tion ⌬
(r), defined as the density of the dimer minus the
sum of two isolated atomic densities, positioned at the dimer
nuclear coordinates; ghost atoms are not included in the
atomic calculations since the isolated atoms must be spheri-
cal. A positive ⌬
(r) corresponds to a region where the den-
sity increases upon dimer formation; a negative value corre-
sponds to a density decrease.
Table I presents forces from Hartree–Fock and BD共T兲,
denoted F
HF
and F
BD共T兲
, respectively. The former are small,
repulsive, and vanish as overlap reduces; the latter are in
reasonable quantitative agreement with the near-exact val-
ues. The absence of electron correlation in Hartree–Fock
means that the only distortion is due to overlap 共exchange兲
effects, which distorts the atomic densities away from one
another.
59,60
The attractive force due to the density is then
smaller than the nuclear repulsion, resulting in a net repul-
sive force. The small overlap at these large R values means
this distortion—and hence the force—is very small. We do
not present plots of ⌬
(r) for Hartree–Fock since it can be
difficult to distinguish the distortion from the numerical
noise in the calculations—the Hartree–Fock dimer is essen-
tially two spherical atoms at these large separations. By con-
trast, electron correlation in BD共T兲 distorts the atomic den-
sities towards one another as in the exact case, overcoming
the small exchange distortion. This is demonstrated in Fig. 2
where ⌬
(r) for BD共T兲 is plotted along the He–He bond
axis for the three R values. At the nuclei, the plots approach
⫺ 182, ⫺ 126, and ⫺ 88⫻ 10
⫺ 7
a.u. for R⫽ 8.0, 8.5, and 9.0
a.u., respectively, which are not visible on the scale. On ei-
ther side of each nucleus, ⌬
(r) exhibits a positive peak,
which is much more pronounced on the side of the nucleus
nearest to the other atom. Analogous plots have been pre-
sented for H
2
共Ref. 57兲; further discussion on the density of
rare gas dimers can be found in Ref. 61. The quality of the
BD共T兲 density is quantified by calculating the associated
Hellmann–Feynman forces. For R⫽ 8.0, 8.5, and 9.0 a.u.,
the forces are ⫹4.8, ⫹3.2, and ⫹ 2.1⫻ 10
⫺ 6
a.u., which are
close to the near-exact values of Table I.
We next consider forces from HFKS calculations. The
FIG. 2. Density differences ⌬
(r) determined from BD共T兲 densities, plotted
along the He–He bond axis for the three R values. The nuclei are symmetri-
cally placed either side of z⫽ 0.
11116 J. Chem. Phys., Vol. 117, No. 24, 22 December 2002 M. J. Allen and D. J. Tozer
Downloaded 26 Apr 2011 to 129.234.252.65. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
