Perspective: Advances and challenges in treating van der Waals dispersion forces in density functional theory

1.4 MB·English

by  Jiří Klimeš

Perspective: Advances and challenges in treating van der Waals dispersion forces in
density functional theory
Jiˇr´ı Klimeˇs and Angelos Michaelides∗
Thomas Young Centre, London Centre for Nanotechnology and Department of Chemistry,
University College London, London, WC1E 6BT, UK
(Dated: January 30, 2013)
Electron dispersion forces play a crucial role in determining the structure and properties of
biomolecules, molecular crystals and many other systems. However, an accurate description of
dispersion is highly challenging, with the most widely used electronic structure technique, density
functional theory (DFT), failing to describe them with standard approximations. Therefore, appli-
cations of DFT to systems where dispersion is important have traditionally been of questionable
3
1 accuracy. However,thelastdecadehasseenasurgeofenthusiasmintheDFTcommunitytotackle
0 this problem and in so-doing to extend the applicability of DFT-based methods. Here we discuss,
2 classify, and evaluate some of the promising schemes to emerge in recent years. A brief perspective
ontheoutstandingissuesthatremaintoberesolvedandsomedirectionsforfutureresearcharealso
n
provided.
a
J

9
2 Copyright (2012) 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. The following article appeared in J. Chem. Phys. 137, 120901 (2012) and may be found
i
c at http://jcp.aip.org/resource/1/jcpsa6/v137/i12/p120901 s1.
s
-
l PACSnumbers: 71.15.Mb31.15.-p87.15.A-
r
t
m
I. INTRODUCTION persion is becoming one of the hottest topics in compu-
.
t tational chemistry, physics, and materials science. Fig. 1
a
m The theoretical description of matter as well as of underlines this point, where it can be seen that over 800
dispersion-based DFT studies were reported 2011 com-
- many chemical, physical, and biological processes re-
d quires accurate methods for the description of atomic pared to fewer than 80 in the whole of the 1990s.
n
and molecular-scale interactions. Whilst there are many
o
quantummechanicalapproaches,inthepastfewdecades 800
c
[ Kohn-Sham density functional theory (DFT)1,2 has es- ons
tablisheditselfasthetheoreticalmethodofchoiceforthis cati600
1 task,undergoingameteoricriseinlargepartsofphysics, ubli
v chemistry, and materials science. The rise of DFT and of p400
0 er
6 its uptake in academia and industry has been widely umb200
9 discussed, and was perhaps illustrated most clearly in N
6 Burke’s recent Spotlight article on DFT.3
0
. Although DFT is in principle exact, in practice ap- 1991 1996 2001 2006 2011
1
Year
0 proximations must be made for how electrons interact
3 with each other. These interactions are approximated FIG.1. ThenumberofdispersioncorrectedDFTstudieshas
1 withso-calledexchange-correlation(XC)functionalsand increased considerably in recent years. The total number of
v: much of the success of DFT stems from the fact that studies performed is difficult to establish precisely, however,
i XC functionals with very simple forms often yield accu- anestimateismadeherebyillustratingthenumberofpapers
X
rate results. However there are situations where the ap- that cite at least one of 16 seminal works in the field (Refs.
r proximate form of the XC functional leads to problems. 4–19). (DatafromWebofKnowledge,July2012fortheyears
a
1991-2011)
One prominent example is the inability of “standard”
XC functionals to describe long-range electron correla-
tions, otherwise known as electron dispersion forces; by Dispersion can be viewed as an attractive interaction
standard XC functionals we mean the local density ap- originating from the response of electrons in one region
proximation (LDA), generalized gradient approximation to instantaneous charge density fluctuations in another.
(GGA) functionals or the hybrid XC functionals. The The leading term of such an interaction is instantaneous
“lack”ofdispersionforces–oftencolloquiallyreferredto dipole-induced dipole which gives rise to the well known
asvanderWaals(vdW)forces–isoneofthemostsignif- 1/r6 decay of the interaction energy with interatomic
−
icantproblemswithmodernDFTand,assuch,thequest separationr. StandardXCfunctionalsdon’tdescribedis-
forDFT-basedmethodswhichaccuratelyaccountfordis- persion because: (a) instantaneous density fluctuations
2
are not considered; and (b) they are “short-sighted” in
that they consider only local properties to calculate the
XCenergy. Theconsequencefortwonoblegasatoms,for
example, is that these functionals give binding or repul-
siononlywhenthereisanoverlapoftheelectrondensities
ofthetwoatoms. Sincetheoverlapdecaysexponentially
withtheinteratomicseparation,sotoodoesanybinding.
We show this in Fig. 2 for one of the most widely used
GGAs – the so-called Perdew-Burke-Ernzerhof (PBE)
functional20 –forabindingcurvebetweentwoKratoms.
FIG. 3. Two binding configurations of the DNA base pairs
adenineandthymine. Ahydrogenbondedstructureisshown
on the left (hydrogen bonds indicated by red dots) and a
“stacked” geometry on the right. For the hydrogen bonded
configuration a standard XC functional such as PBE gives a
bindingenergythatisinrathergoodagreementwiththeref-
erence value (lower part of the graph, data in kJ/mol).22 In
contrast, PBE fails to describe the binding of the “stacked”
configuration where dispersion interactions contribute signif-
FIG. 2. Binding curves for the Kr dimer obtained with the
icantly to the intermolecular binding.
PBE exchange-correlation functional and an accurate model
potential.21 Dispersion originates from fluctuations in the
electrondensitywhichpolarizedifferentatomsandmolecules,
schemes has been proposed to treat dispersion within
as shown schematically in the lower right diagram. The in-
teraction then exhibits the well-known −1/r6 decay. This DFT. Here we will discuss some of the main approaches
−1/r6 decayisnotreproducedwithPBE(orothersemi-local developed and in the process attempt to provide a use-
functionals)whichinsteadgivesanexponentialdecayforthe ful classification of them. We also highlight some of the
interaction because PBE relies on the overlap of density to obstacles that must be overcome before improved DFT-
obtain the interaction. based methods for including dispersion are available. By
itsnature,thisarticleisalimitedpersonalviewthatcan-
The binding of noble gases is a textbook dispersion not cover every development in this thriving field. For
bonded system. However, it has become increasingly the most part we have tried to keep the overview simple,
apparent that dispersion can contribute significantly to and have aimed it primarily at newcomers to the field,
the binding of many other types of materials, such as although we do go into more depth at the later stages.
biomolecules,adsorbates,liquids,andsolids. Fig.3illus- For more detailed discussions of some of the methods
trates a more “real world” example, where an accurate shown here the interested reader should consult the re-
description of dispersion is critical. The figure reports viewsofGrimme,23 Tkatchenkoet al.,24 Johnsonet al.25
bindingenergiesobtainedfromPBEandanaccurateref- or others.26–29 Some relevant developments are also dis-
erencemethodforDNAbasepairsintwodifferentconfig- cussedinthereviewsofSherrill30andRileyetal.31which
urations. In the first configuration the binding between are more focused on post Hartree-Fock (HF) methods.
the base pairs is dominated by hydrogen bonding. Hy-
drogen bonding is governed mainly by electrostatics and
a standard functional such as PBE predicts reasonable II. A CLASSIFICATION OF THE COMMON
hydrogen bond strengths and as a result the stability of DFT-BASED DISPERSION METHODS
the dimer is quite close (within 15%) to the reference
value.22 However, in the other “stacked” configuration Many DFT-based dispersion techniques have been de-
the binding is dominated by dispersion forces and for velopedandratherthansimplylistingthemall,itisuse-
this isomer PBE is hopeless, hardly predicting any bind- fultotrytoclassifythem. Anaturalwaytodothisisto
ing between the base pairs at all. This huge variability considerthelevelofapproximationeachmethodmakesin
inperformanceisfarfromidealandsincethestackedar- obtaining the long range dispersion interactions, that is,
rangement of base pairs is a common structural motif in the interactions between well separated fragments where
DNA,theresultsuggeststhatDNAsimulatedwithPBE dispersion is clearly defined. In doing this it turns out
would not be stable. that groups of methods which exploit similar approxi-
As we will see throughout this article there are many mations emerge. Here, we simply identify these group-
other examples of the importance of dispersion to the ings and rank them from the most approximate to the
bonding of materials and in recent years a plethora of more sophisticated approaches in a manner that loosely
3
resemblesthewell-known“Jacob’sladder”offunctionals persionplaysamajorrolesuchasgraphiteornoblegases
introduced by Perdew.32 Therefore, by analogy, we in- on metals. However, the results with LDA for dispersion
troduce a “stairway to heaven” for long range dispersion bonded systems have limited and inconsistent accuracy
interactionsandplaceeachgroupofdispersioncorrection and the asymptotic form of the interaction is incorrect.
schemes on a different step of the stairway. In complete More promising approaches on the ground level of our
analogy to the ladder, when climbing the stairway pro- stairway are density functionals specifically fitted to re-
gressively higher overall accuracy can be expected until produceweakinteractionsaroundminimaaswellasspe-
exact results, and thus heaven, are reached. We stress cially adapted pseudopotentials.
the point “overall accuracy” because, as with the lad- The “Minnesota functionals”34 are an example of a
der, climbing the stairway does not necessarily mean new breed of functionals that are fitted to a dataset that
higher accuracy for every particular problem but rather includes binding energies of dispersion bonded dimers,
a smaller probability to fail, i.e., a statistical improve- amongst other properties. Although these functionals
ment in performance.33 A schematic illustration of the can describe binding accurately at separations around
stairway is shown in Fig. 4 and in the following sections minima, they suffer from the same incorrect asymptotics
each step on the stairway is discussed. as the LDA does. On the plus side, the reference data
used in their construction also contains properties other
than weak interactions (e.g. reaction barriers), so that
they can be rather accurate for general chemistry prob-
Functional Heaven
lems. This is a clear advantage over “proper” dispersion
correction schemes which often utilise GGA functionals
that can be less accurate for such problems.
4+: MBD, RPA,…
For electronic structure codes that make use of pseu-
dopotentials,dispersioncanalsobemodelledbyaddinga
3: vdW-DF,…
specially constructed pseudopotential projector. Within
this class are the dispersion corrected atom-centered
Accuracy? 2: DFT-D3, vdW(TS), potentials (DCACP)11 and the local atomic potentials
Cost … (LAP) methods.35 These approaches have shown a lot
1: DFT-D,… of promise,36 however, effort is required to carefully fit
the potentials for each element and they are not easily
0: DCACP,… transferable to all-electron methods.
FIG. 4. In analogy with the Jacob’s ladder classification of
B. Step one – Simple C corrections
functionals the “stairway to heaven” is used here to clas- 6
sify and group DFT-based dispersion correction schemes. At
groundlevelaremethodswhichdon’tdescribethelongrange The basic requirement for any DFT-based dispersion
asymptotics. Simple C6 correction schemes sit on the first schemeshouldbethatityieldsreasonable 1/r6 asymp-
step, on the second step are approaches that utilise environ- −
totic behavior for the interaction of particles in the gas
mentdependentC corrections. Thelongrangedensityfunc-
6 phase, where r is the distance between the particles. A
tionals sit on step 3 and on step 4 and above are approaches
simpleapproachforachievingthisistoaddanadditional
which go beyond pairwise additive determinations of disper-
energy term which accounts for the missing long range
sion. Uponclimbingeachstepofthestairwaythelevelofap-
proximation is reduced and the overall accuracy is expected attraction. The total energy then reads
to increase.
E =E +E , (1)
tot DFT disp
where E is the DFT total energy computed with a
DFT
given XC functional. The dispersion interaction is given
A. Ground – Binding with incorrect asymptotics
by
(cid:88)
First, the ground is occupied by methods that sim- E = CAB/r6 , (2)
disp − 6 AB
ply don’t describe the long range asymptotics. These
A,B
approaches give incorrect shapes of binding curves and
underestimate the binding of well separated molecules. where the dispersion coefficients CAB depend on the el-
6
Some of these methods are, nevertheless, being used for emental pairs A and B. Within this approach dispersion
weakly bonded systems. Although this might seem odd, is assumed to be pairwise additive and can therefore be
this somewhat misguided approach is used because some calculated as a sum over all pairs of atoms A and B.
standard DFT functionals bind dispersion bonded sys- Methods on step 1 of the stairway use coefficients that
tems at short separations. A prominent example is the are tabulated, isotropic (i.e. direction independent) and
LDA which has been used to study systems where dis- constant,andthesemethodsaregenerallytermed“DFT-
4
D”. ingfunctioncanactuallyaffectthebindingenergieseven
Mostly because of the simplicity and low computa- more than the asymptotic C6 coefficients.38 The fitting
tional cost this pairwise C6/r6 correction scheme is also effectively includes the effects of C8/r8 or C10/r10
widelyused. Nonetheless,ithasatleastfourclearshort- and higher contributions, although some methods treat
comingswhichlimittheaccuracyonecanachievewithit. themexplicitly.19,39 Aninterestingdampingfunctionhas
First, the C /r6 dependence represents only the leading been proposed within the so-called DFT coupled cluster
6
term of the correction and neglects both many-body dis- approach (DFT/CC), it has a general form that can ac-
persion effects37 and faster decaying terms such as the tually force the dispersion correction to be repulsive.40
C /r8 or C /r10 interactions. Second, it is not clear
8 10
where one should obtain the C coefficients. Various
6
Energy
formulae, often involving experimental input (ionization
potentials and polarizabilities), have been proposed for
this.7–9,12 However, this reliance on experimental data
limitedthesetofelementsthatcouldbetreatedtothose
E
present in typical organic molecules. The third issue is E
DFT
that the C6 coefficients are kept constant during the cal- tot Distance
culation, and so effects of different chemical states of the
atom or the influence of its environment are neglected.
The fourth drawback, which we discuss later, is that the E
C /r6 function diverges for small separations (small r) disp
6
and this divergence must be removed.
For a more widely applicable method a more consis-
tent means of deriving the dispersion coefficients is re-
quired. In 2006 Grimme published one such scheme, re-
ferred to as DFT-D2.17 In this approach the dispersion
FIG. 5. Schematic illustration of a binding curve (E ) ob-
coefficients are calculated from a formula which couples tot
tained from a dispersion corrected DFT calculation and the
ionizationpotentialsandstaticpolarizabilitiesofisolated
contributions to it from the regular DFT energy (E ) and
DFT
atoms. Data for all elements up to Xe are available and the dispersion correction (E ). The function −1/r6 used
disp
this scheme is probably the most widely used method to model the dispersion correction diverges for small r (red
for accounting for dispersion in DFT at present. Whilst dashed curve) and must be damped (solid red curve). Since
incredibly useful, DFT-D2 is also not free from prob- the details of the damping strongly affect the position of the
lems. In particular, for some elements arbitrary choices energy minima on the binding curve, the damping function
ofdispersioncoefficientsstillhadtobemade. Forexam- needs to be fit to reference data. In addition, the damping
function will be different for different XC functionals used to
ple,alkaliandalkaliearthatomsusecoefficientsthatare
obtain E .
averages of the previous noble gas and group III atom. DFT
Furthermore, the dispersion energy, E , is scaled ac-
disp
cording to the XC functional used and as a result the Before leaving the simple C6/r6 schemes we note
interactionenergyoftwowellseparatedmonomersisnot that although the approaches developed by Grimme are
constant but sensitive to the choice of XC functional. widely used, other parameterizations, which differ in,
e.g., the combination rules used, are available.40–43 Fur-
With the simple C /r6 correction schemes the disper-
6 thermore, functionals such as B97-D17 and ωB97X-D44
sioncorrectiondivergesatshortinter-atomicseparations
havebeenspecificallydesignedtobecompatiblewiththis
and so must be “damped”. The damped dispersion cor-
level of dispersion correction.
rection is typically given by a formula like
(cid:88)
E = f(r ,A,B)CAB/r6 , (3)
disp − AB 6 AB C. Step two – Environment-dependent C
A,B 6
corrections
where the damping function f(r ,A,B) is equal to one
AB
for large r and decreases E to zero or to a constant A problem with the simple “DFT-D” schemes is that
disp
for small r. We illustrate the divergence and a possible the dispersion coefficients are predetermined and con-
dampingbytheredcurvesinFig.5. Howthedampingis stant quantities. Therefore the same coefficient will be
performed is a thorny issue because the shape of the un- assigned to an element no matter what its oxidation or
derlying binding curve is sensitive to the XC functional hybridization state. However the errors introduced by
used and so the damping functions must be adjusted so this approximation can be large, e.g. the carbon C co-
6
as to be compatible with each exchange-correlation or efficients can differ by almost 35% between the sp and
exchange functional. This fitting is also sensitive to the sp3 hybridized states.9 Thus the emergence of methods
definition of atomic size (van der Waals radii are usu- where the C coefficients vary with the environment of
6
ally used) and must be done carefully since the damp- the atom has been a very welcome development. We put
5
these methods on the second step of our stairway. They atoms,thiscreatesasymmetricelectrondensityandthus
still use Equation 3 to obtain the dispersion correction non-zero dipole and higher-order electrostatic moments,
and, as with the step 1 methods, some reference data which causes polarization in other atoms to an extent
(such as atomic polarizabilities) is used to obtain the C given by their polarizability. The result of this is an at-
6
coefficients. tractive interaction with the leading term being dipole-
We now discuss three step 2 methods: DFT-D3 induced dipole. The biggest question in the BJ model is
of Grimme, the approach of Tkatchenko and Scheffler how to quantify the XC hole and in the method it is ap-
(vdW(TS)),andtheBecke-Johnsonmodel(BJ).Theuni- proximated as the exchange-only hole (calculated using
fying concept of these methods is that the dispersion co- the Kohn-Sham orbitals). An average over the positions
efficientofanatominamoleculedependsontheeffective r1 ofthereferenceelectronthengivestheaveragesquare
volume of the atom. When the atom is “squeezed”, its of the hole dipole moment, denoted d2 which together
(cid:104) x(cid:105)
electron cloud becomes less polarizable leading to a de- with the polarizabilities α enters the formula for the dis-
crease of the C coefficients. persion coefficient
6
Grimmeetal.19capturetheenvironmentaldependence
α α d2 d2
oftheC6 coefficientsbyconsideringthenumberofneigh- C6AB = d2A αB(cid:104) +x(cid:105)Ad(cid:104)2x(cid:105)αB . (4)
bors each atom has. When an atom has more neighbors (cid:104) x(cid:105)A B (cid:104) x(cid:105)B A
it is thought of as getting squeezed and as a result the
This is, in fact, very similar to the vdW(TS) formula
C coefficient decreases. This effect is accounted for by
6 but in vdW(TS) precalculated C coefficients are used
having a range of precalculated C coefficients for vari- 6
6 insteadoftheholedipolemoments. TheBJmodelisvery
ouspairsofelementsindifferentreference(hybridization)
intriguingandseveralauthorshavestudiedhowitrelates
states. In the calculation of the full system the appro-
to formulae derived from perturbation theory.47–49
priate C coefficient is assigned to eachpair of atoms ac-
6
IntheBJmodeltheC coefficientsarealteredthrough
cordingtothecurrentnumberofneighbors. Thefunction 6
two effects. First, the polarizabilities of atoms in
calculating the number of neighbors is defined in such a
moleculesarescaledfromtheirreferenceatomvaluesac-
waythatitcontinuouslyinterpolatesbetweentheprecal-
cording to their effective atomic volumes. Second, the
culated reference values. Therefore if the hybridization
dipole moments respond to the chemical environment
stateofanatomchangesduringasimulation,theC coef-
6
throughchangesoftheexchangehole,althoughthiseffect
ficientcanalsochangecontinuously. Despitethesimplic-
seems to be difficult to quantify precisely. The details of
ity of this approach, termed “DFT-D3”, the dispersion
howtoobtainatomicvolumesandtheexchangeholeare
coefficients it produces are pretty accurate. Specifically,
knowntoaffecttheresultsobtainedtosomeextent,50–52
basedonthereferencedataofMeathandcoworkersthere
butoveralltheaccuracyoftheasymptoticC coefficients
is a mean absolute percentage deviation (MAPD) in the 6
obtained from the BJ model is quite satisfactory, with a
C coefficientsof8.4%. Moreover,theadditionalcompu-
6 MAPDof12.2%forthedataofMeathandcoworkers.19
tational cost is negligible since the number of neighbors
Compared to the two previous step 2 approaches, how-
can be obtained quickly.
ever, the computational cost is relatively high, in the
In2009TkatchenkoandScheffler18 proposedamethod
same ballpark as the cost of a hybrid functional.26
which relies on reference atomic polarizabilities and ref-
erenceatomicC coefficients45tocalculatethedispersion
6
energy. These quantities are sufficient to obtain the C
6
coefficientforapairofunlikeatoms.46Toobtainenviron- D. Step three – Long-range density functionals
ment dependent dispersion coefficients effective atomic
volumes are used. During the calculation on the system
The methods discussed so far require predetermined
of interest the electron density of a molecule is divided
input parameters to calculate the dispersion interaction,
between the individual atoms (the Hirshfeld partition-
either the C coefficients directly or the atomic polariz-
6
ing scheme is used) and for each atom its corresponding
abilities. Now we discuss approaches that do not rely on
density is compared to the density of a free atom. This
external input parameters but rather obtain the disper-
factor is then used to scale the C coefficient of a ref-
6 sion interaction directly from the electron density. This,
erence atom which changes the value of the dispersion
inprinciple,isamoregeneralstrategyandthusweplace
energy. The accuracy of the final (isotropic, averaged)
these methods on step 3 of our stairway. The methods
C coefficientsisquitehigh,withtheMAPDonthedata
6 have been termed non-local correlation functionals since
of Meath and coworkers being only 5.4%. However, it
they add non-local (i.e. long range) correlations to local
is not yet clear if the scaling of C coefficients with vol-
6 or semi-local correlation functionals.
ume will be accurate for more challenging cases such as
Thenon-localcorrelationenergyEnl iscalculatedfrom
different oxidation states in e.g. ionic materials. c
The most complex step 2 method is the BJ (cid:90)(cid:90)
model,13–16,39 which exploits the fact that around an Ecnl = dr1dr2n(r1)ϕ(r1,r2)n(r2). (5)
electron at r there will be a region of electron density
1
depletion, the so-called XC hole. Even for symmetric This is adouble space integral where n(r) is theelectron
6
densityandϕ(r ,r )issomeintegrationkernel. Theker- DF and several of its offspring are now implemented in
1 2
nelisanalogoustotheclassicalCoulombinteractionker- widely distributed DFT codes such as Siesta,68 VASP,61
nel1/r r butwithamorecomplicatedformulaused QuantumESPRESSO,69andQChem.66Otherwaystore-
1 2
| − |
for ϕ(r ,r ) with (1/r r 6) asymptotic behavior. duce the cost of vdW-DF calculations have also been re-
1 2 1 2
O | − |
The formula has a pairwise form and thus ignores the ported, for example, Silvestrelli has shown70,71 that uti-
mediumbetweenpointsr andr . VariousformsforEnl lizing Wannier functions to represent the electron den-
1 2 c
were proposed in the nineties4–6 but were restricted to sity allows for an analytic evaluation of the function-
non-overlappingfragments. Thisratherseverelimitation als proposed in the nineties.5,6 Particularly interesting
was removed by Dion et al. in 200410 who proposed a is the local response dispersion (LRD) approach of Sato
functional form which can be evaluated for overlapping andNakai72,73 whichyieldsveryaccurateC coefficients.
6
molecules and for arbitrary geometries. Within this ap- Overall, by not relying on external reference data, the
proach the XC energy E is calculated as step 3 approaches are less prone to fail for systems out-
xc
side of the reference or fitting space of step 2 methods.
Exc =ExGGA+EcLDA+Ecnl, (6) However, very precise calculations of the dispersion en-
ergy are difficult and the formulae underlying vdW-DF
where the terms on the right hand side are the exchange and similar approaches can be somewhat complicated.
energy in the revPBE approximation,53 the LDA corre-
lation energy, and the non-local correlation energy term,
respectively. This method has been termed the van der
Waals density functional (vdW-DF). vdW-DF is a very E. Higher steps – Beyond pairwise additivity
important conceptual development since it adds a de-
scription of dispersion directly within a DFT functional
Themaincharacteristicofthemethodsdiscussedsofar
and combines correlations of all ranges in a single for-
is that they consider dispersion to be pairwise additive.
mula.
As a consequence the interaction energy of two atoms or
SincethedevelopmentoftheoriginalvdW-DFanum- moleculesremainsconstantnomatterwhatmaterialsep-
ber of follow up studies have aimed at understanding arates them and all atoms interact “on their own” with
and improving the performance of the method. First, it no consideration made for collective excitations.74 While
turns out that vdW-DF tends to overestimate the long such effects don’t seem to be crucially important in the
range dispersion interactions: Vydrov and van Voorhis gasphase,especiallyforsmallmolecules,theyareimpor-
have shown that when the C coefficients are evaluated tant for adsorption or condensed matter systems where
6
the errors can be as large as 30%.54 To address this thebareinteractionisscreened.75 Wenowbrieflydiscuss
∼
these authors proposed (computationally cheaper) func- someofthemethodsbeingdevelopedwhichtreatdisper-
tionals that reduce the average errors by approximately sion beyond the pairwise approximation. The range of
50%.54–57 The developers of the original vdW-DF tried approaches is quite wide, from methods based on atom
to address its tendency to overbind at large separations centered interactions, to methods involving density, to
by proposing vdW-DF2.58,59 This functional, which in- methods using electron orbitals. Because of the fresh-
volveschangestoboththeexchangeandnon-localcorre- ness of most of the methods we discuss them together.
lationcomponentstendstoimprovethedescriptionofthe
binding around energy minima, however, the C coeffi-
6
cientsitpredictsareconsiderablyunderestimated.56 Sec-
ond, aside from the particular form of Enl, the choice of
c
the exchange functional used in Equation 6 has received
considerable attention. The original revPBE exchange
functional chosen for vdW-DF often leads to too large
intermolecular binding distances and inaccurate binding
energies. To remedy these problems alternative “less re- 1/r6 ?
pulsive” exchange functionals have been proposed60–63 ∼− ?
which when incorporated within the vdW-DF scheme
(Equation 6) lead to much improved accuracy. Of these
the “optB88” and “optPBE” exchange functionals have FIG. 6. The long-range electron correlations of two isolated
been shown to offer very good performance for a wide atoms can be described by an effective −1/r6 formula. In a
range of systems.60,61,64,65 condensed system the interaction is modified (screened) by
the presence of the other electrons.
Initially the non-local correlation functionals came, to
a lesser or greater extent,66,67 with a higher computa-
tional cost than GGAs or hybrid functionals. However, Recently, the Axilrod-Teller-Muto76,77 formula has
thanks to the work of Rom´an-P´erez and Soler the com- been used to extend the atom-centred pairwise ap-
putational cost of vdW-DF is now comparable to that of proaches to include three-body interactions.19,78 The
a GGA.68 In addition, self-consistent versions of vdW- triple-dipole interaction between three atoms A, B, and
7
C can be written as range correlations are treated by, e.g., second-order per-
turbation theory, the coupled cluster method, or RPA.
3cos(α)cos(β)cos(γ)+1
EABC =CABC , (7) Range separated functionals can be very accurate, can
9 rA3BrB3CrA3C account for many-body interactions, and can deliver rel-
ativelyfastconvergenceofthecorrelationenergywithre-
where α, β, and γ are the internal angles of the triangle
spect to basis set size.104–109 Another example from this
formed by the atoms A, B, and C. The dispersion coef-
class are the so called double hybrid functionals, which
ficient CABC is again obtained from reference data. No-
9 include Fock exchange and a second order perturbation
tably, vonLilienfeldandTkatchenkoestimatedthemag-
theory correlation contribution. Generally the coeffi-
nitudeofthethreebodytermstobe 25%(destabilis-
∼− cients for these contributions are fitted to reproduce ref-
ing)ofthetwobodytermfortwographenelayers.78Since
erence data.110,111 Since the original double hybrid func-
thethree-bodyinteractionisrepulsiveforatomsforming
tionals were more concerned with reaction energies and
an acute angled triangle, the interaction was found to
barrier heights than dispersion, dispersion interactions
be destabilising for stacked configurations of molecular
were actually underestimated. Newer functionals, e.g.,
clusters such as benzene dimers.
mPW2PLYP-D, add a dispersion correction in the sense
The atom-centered approach can be used to approxi- of step 1 or 2 methods, however, they are computation-
matethemany-bodydispersioninteractionusingamodel allyveryexpensiveandatpresentprohibitivelyexpensive
of coupled dipoles.79,80 Here quantum harmonic oscilla- for most condensed phase and surface studies.33,111
tors with characteristic frequencies occupy each atomic
position and the dispersion interaction is obtained from
the shifts of the frequencies of the oscillators upon
switching on their interaction. The model has been ap-
plied outside DFT for some time81–85 and has recently
been used to show the non-additivity of dispersion in
anisotropic materials.86 The initial applications of this F. Summary
approachtermedmany-bodydispersion(MBD)arequite
promising and suggest that the higher-order dispersion
The higher steps of the stairway with the promising
terms can be important even for systems such as solid
benzene.80 However, getting accurate relations between schemesthatgobeyondthepairwiseapproximationclose
our overview and classification of methods for treating
atoms and oscillator models seems to be rather challeng-
dispersion within DFT. In Table I we summarize some
ing, especially when describing both localized and delo-
of the key aspects of the most relevant schemes, such as
calized electrons.
what properties they use to obtain the dispersion con-
In the context of DFT, orbitals can be used to calcu-
tribution(via,e.g.,theC )andtheapproximaterelative
6
latethecorrelationenergyusingtheadiabatic-connection
computationalcost. Thecomputationalcostofamethod
fluctuation-dissipation theorem (ACFDT).6,87 The par-
is, of course, a key factor since no matter how accurate
ticular approach which has received the most attention
it is, it will not be used if the computational cost is pro-
recently is the so-called random phase approximation
hibitive. In this regard the cost of the simple pairwise
(RPA).88,89 Results from RPA have been very encour-
corrections on step 1 is essentially zero compared to the
aging and it exhibits, for example, a consistent accu-
cost of the underlying DFT calculation and these meth-
racy for solids90,91 and the correct asymptotic descrip-
ods can be recommended as a good starting point for
tion for the expansion of graphite, a feature which the
accountingfordispersion. Becauseofthecorrectasymp-
pairwise methods fail to capture.37,92 RPA is analogous totic behavior (at least for gas phase molecules)37 they
to post-HF methods and indeed direct links have been
arepreferabletomethodsfromground. Comparedtothe
established.93,94 Unfortunately, it also shares with the
step1methodstheDFT-D3andvdW(TS)schemesdon’t
post-HFmethodsahighcomputationalcost(thecostin-
add a significant computational cost, whereas the vdW-
creases approximately with the fourth power of the sys-
DF approach increases the computational time by about
temsize)95,96 andslowconvergencewithrespecttobasis
50%comparedtoaGGAcalculation. Sincetheaccuracy
setsize.97–101 ApproachesgoingbeyondRPAarealsore-
of these approaches tends to be higher than that of the
ceiving attention, for example the second-order screened
step 1 methods, they should be preferred over the step 1
exchange102 or single excitations corrections.103
corrections. However,whenappliedtoaparticularprob-
Anothermethodinvolvingorbitalstoobtainthecorre- lem, the accuracy of any method used should be tested
lation energy is to combine DFT with post-HF methods or verified against experimental or other reference data
in the hope of exploiting the benefits of each approach. since any approach can, in principle, fail for a specific
While the post-HF methods can describe long range in- system. The methods on higher steps are mostly in de-
teractions accurately (albeit expensively) using orbitals, velopment and currently the ACFDT or range separated
DFT is efficient for the effective description of the short functionals require a computational cost two to four or-
range part of the interactions. So called range separated ders of magnitude higher than a GGA calculation which
methods are an example of this approach where long limits their applicability.
8
TABLEI.SummaryofsomeofthemostwidelyusedmethodsforcapturingdispersioninteractionsinDFT,orderedaccording
tothe“Step”theysitoninthestairwaytoheaven(Fig. 4). InformationonthereferenceusedfortheC coefficients,whatthe
6
C depend on, and the additional computational cost of each approach compared to a regular GGA calculation is reported.
6
Additional
Method Step Reference for C C depend on Ref.
6 6 computational cost a
Minnesota 0 None N/A None e.g. 34
DCACP 0 None N/A Small 11
DFT-D 1 Various Constant Small e.g. 110
DFT-D3 2 TDDFT Structure Small 19
Polarizabilities
vdW(TS) 2 Atomic volume Small 18
and atomic C
6
BJ 2 Polarizabilities Atomic volume, X hole Large 14
LRD 3 C calculated Density Small 72
6
vdW-DF 3 C calculated Density 50% 10
6
≈
Dbl. hybrids 4 None or as “-D” Orbitals Large 110
a TheBJmodelandthedoublehybrids(labelled“Dbl. hybrids”)
aremorecomputationallyexpensivethanthesimpler“DFT-D”
methods26 andthecalculationofthecorrelationenergyin
vdW-DFleadstoa≈50%slow-downcomparedtoaGGA
calculationwhendoneefficiently.68
9
III. GENERAL PERFORMANCE steps of the stairway can surpass those from methods on
highersteps. Thisisnotthatsurprisingsincesomeofthe
methodsreportedinTableIIwereactuallydevelopedby
We have commented in passing on how some of the
fitting to the S22 data set itself.
methods perform, mainly with regard to the dispersion
coefficients. Now we discuss accuracy in more detail by
focussingonbindingenergiesforafewrepresentativegas
TABLEII.Meanabsolutedeviationsandmeanabsoluteper-
phase clusters, a molecular crystal, and an adsorption
centage deviations of different dispersion based DFT meth-
problem. Before doing so, it is important to emphasise
ods on the S22 data set of Jureˇcka et al..22,115 Results are
two points. First, establishing the accuracy of a method
in kJ/mol for MAD and percent for MAPD. Results from
is not as straightforward as it might seem since it must
second order Møller-Plesset perturbation theory (MP2) – a
be tested against accurate reference data (e.g. bind- widely used post-HF method – are shown for comparison as
ing energies). Experimental data can be inaccurate and are results from the LDA and PBE functionals as examples
is rarely directly comparable to theory since quantum of methods that don’t treat dispersion explicitly.
nuclear and/or thermal effects, which affect the experi-
Method MAD MAPD Ref.
mental values, are often neglected in simulation studies.
Theoretical reference data, on the other hand, is more Ground
appropriate but rather hard to come by since accurate LDA 8.66 30.5 116
reference methods such as post-HF methods (MP2, CC)
PBE 11.19 57.6 116
or quantum Monte Carlo (QMC) are so computationally
M06-2X 2.52 11.3 117
expensivethattheycangenerallybeappliedtoonlyvery
small systems (tens of atoms).30 Second, the results ob- LAP 2.51 7.8 35
tained with many dispersion based DFT approaches are Step 1
oftenstronglyaffectedbythefittingprocedureusedwhen B97-D2 1.51 7.3 44
combiningthe long rangedispersioninteractionwiththe ωB97X-D 0.76 5.6 44
underlying exchange or XC functionals.
Step 2
BLYP-D3 0.96* – 19
A. Gas phase clusters PBE-vdW(TS) 1.25** 9.2** 18
rPW86-BJ 1.50 6.1 51
Because of growing computer power and algorithmic Step 3
improvements, systems of biological importance such as LC-BOP+LRD 0.86 4.6 73
DNA and peptides can now be treated with DFT. As
vdW-DF 6.10 22.0 67
discussedintheintroduction, dispersionplays animpor-
optB88-vdW 1.18 5.7 60
tant role in stabilising the folded state and thus disper-
vdW-DF2 3.94 14.7 118
sion corrected functionals are required. However, since
accuratereferencedatacanonlybeobtainedforsystems VV10 1.35 4.5 118
with tens of atoms, it’s just not possible to obtain refer- Higher steps
encedataforDNAoraproteindirectlyandfragmentsof B2PLYP-D3 1.21* – 19
the large molecules are used to build test sets of binding
ωB97X-2 1.17 8.8 119
energiesandgeometriesinstead. Onesuchtestset112–114
PBE-MB – 5.4** 80
is the popular S22 data set of Jureˇcka et al..22 It is use-
ful since it contains 22 different dimers with a range of RPA 3.29** – 101
weak bonding types and with a wide range of interac- post-HF
tion energies. Results for some of the methods are re- MP2 3.58 19.4 115
ported in Table II; specifically we report mean absolute
* using the original reference values of Jureˇcka et al.
deviations (MAD) and mean absolute percentage devi-
(Ref. 22)
ations (MAPD). Included in Table II are results from
** using the values of Takatani et al. (Ref. 120) which
LDA and semi-local PBE. These functionals fail to re-
are very similar to those of Podeszwa et al. (Ref. 115)
produce the correct interaction energies, yielding MADs
of 10 kJ/mol and MAPDs in excess of 30%. Much im-
∼
proved performance is observed with even the simplest Small water clusters, in particular water hexamers,
pairwise corrections schemes. Indeed, on all steps of the are an interesting system where dispersion plays a sig-
stairway above ground there is at least one method with nificant role (see, e.g. Refs. 121–125). The four
a MAD below 1.5 kJ/mol (or 0.4 kcal/mol, 16 meV) relevant isomers (known as “book”, “cage”, “cyclic”,
∼
and MAPD close to 5%. This is very good agreement “prism”, Fig. 7) all have total energies that differ by
with the reference data and a clear indication of the re- <1.3 kJ/mol per molecule according to accurate post-
centimprovementsmadeintheDFT-baseddescriptionof HF methods.60,121 While the “prism” and “cage” iso-
dispersion. Ascanbeseen, theresultsobtainedonlower mers are preferred by post-HF methods, standard XC
10
Ref.130)andthis,therefore,couldserveasarichtesting
ground for DFT-based dispersion schemes.
Solid benzene is one of the most widely examined test
systems,withstudiesfocussingontheexperimentalden-
sity and the cohesive energy of the crystal.131 The lat-
ter is obtained from the experimental sublimation en-
ergy from which the effects of temperature and quan-
tum nuclear effects must be subtracted giving a value
in the 50 to 54 kJ/mol per molecule range.132 Calcula-
tions with PBE give an abysmal cohesive energy of only
10 kJ/mol per molecule and a volume 30% too large
FIG. 7. The four low energy isomers of the water hexamer comparedtoexperiment.133Thishugediff≈erencebetween
are an important example where dispersion interactions play
theoryandexperimentsuggeststhatdispersionisimpor-
asignificantrole. Inthepanelontherightbindingenergydif-
tant to the binding of the crystal and indeed the error
ferencesperwatermoleculeaspredictedbydifferentXCfunc-
is greatly reduced when even rather simple schemes are
tionals are compared to accurate reference data.60,121 While
used. For example, PBE-D2 and the DCACP poten-
thereferencemethodpredictsthe“prism”and“cage”isomers
tobemorestable,standardXCfunctionalssuchasPBEgive tials give estimates of the lattice energy (55.7 kJ/mol133
the“book”and“cyclic”clustersastheoneswithlowerenergy. and 50.6 kJ/mol,36 respectively) and densities that are
It was shown in Ref. 121 that the agreement with the refer- within 10% of experiment. While the PBE-vdW(TS)
encecalculationsisimprovedwhendispersioninteractionsare schemeoverbindsslightly(66.6kJ/mol),ithasbeensug-
accounted for, shown here by results with various dispersion gested that the many-body interactions are quite im-
correction schemes. portant for this system and their inclusion reduces the
PBE-vdW(TS) value by 12 kJ/mol.80 Non-local vdW-
DF overestimates the cell volume by 10% but gives
functionals find the more open “cyclic” and “book” iso- a rather good cohesive energy of 58.3 k≈J/mol.116 vdW-
mers to be more stable. Recent work reveals that dis- DF2 reduces both the binding energy (to 55.3 kJ/mol)
persion impacts profoundly on the relative energies of and the error in the cell volume, which turns out to be
the isomers and improved relative energies are obtained 3% larger than the experimental value.116 Using func-
when dispersion is accounted for. Indeed the water ≈tionals proposed to improve upon the overly repulsive
hexamers have become an important test system with behavior of vdW-DF, such as optB88-vdW leads to an
techniquessuchasDFT-D,121 vdW(TS),121 vdW-DF,126 improvement of the reference volume ( 3% smaller than
optB88-vdW,60 and the modified vdW-DF approach of experiment), but the cohesive energy≈is overestimated,
Silvestrelli (BLYP-vdW(Silv.))127 all having been ap- being 69.4 kJ/mol.116 Again, the overestimated cohe-
plied. The performance of these schemes in predicting sive energy may result from missing many-body inter-
the relative energies of the water hexamers is shown in actions. The RPA has also been applied to this system
Fig. 7, where one can also see the improvement over andwhilethedensityisinverygoodagreementwithex-
a standard functional such as PBE. One general point periment, the cohesive energy is slightly underestimated
to note, however, is that even without a dispersion cor- at 47 kJ/mol.132 Overall, the improvement over a semi-
rection certain functionals can already give quite accu- local functional such as PBE is clear but more applica-
rate absolute binding energies for systems held together tions and tests on a wider range of systems are needed
mainly by hydrogen bonding and adding dispersion cor- to help establish the accuracy of the methods.
rections can actually lead to too large absolute binding
energies. This is indeed the case for the water hexam-
ers where, for example, the vdW(TS) correction to PBE
C. Adsorption – benzene on Cu(111)
givesbindingenergiesforthehexamersthataretoolarge
by 4 kJ/mol per molecule.
≈ Adsorption on solid surfaces is another area where
great strides forward have recently been made with re-
gard to the role of dispersion. Dispersion is of obvious
B. Molecular crystals – solid benzene importancetothebindingofnoblegasestosurfacesbutit
can also be important to chemisorption134 and, e.g., wa-
Molecularcrystalpolymorphpredictionisanotherim- teradsorption.64 Indeedforweaklychemisorbedsystems
portant area where dispersion forces can play a criti- dispersion forces attain an increased relative importance
cal role, and even for molecular crystals comprised of andonecanestimatefromvariousstudies(e.g.,Refs.135
small molecules several polymorphs often exist within a and136)thatdispersioncontributesabout4to7kJ/mol
small energy window.128,129 Identifying the correct ener- to the adsorption energy of a carbon-sized atom; not a
geticorderingofthepolymorphscanthereforebeastrin- negligible contribution.
gent test for any method. A large number of molecular Of the many interesting classes of adsorption system,
crystalshavebeencharacterisedexperimentally(see,e.g. organicmoleculesonmetalshavebecomeahotareaofre-