On the theory underlying the Car-Parrinello method and the role of the fictitious mass parameter

0.96 MB·English

by  Paul Tangney

On the theory underlying the Car-Parrinello method and the role of the fictitious
6 mass parameter.∗
0
0
2 Paul Tangney†
International School for Advanced Studies, via Beirut 2-4, 34013 Trieste, Italy
n
The Molecular Foundry, Lawrence Berkeley National Laboratory, CA 94720
a
Department of Physics, University of California, Berkeley, CA 94720
J
(Dated: February 6, 2008)
6
The theory underlying the Car-Parrinello extended-lagrangian approach to ab initio molecular

] dynamics(CPMD) is reviewed andreexamined using’heavy’ice asa test system. Itis emphasized
h
that theadiabatic decouplingin CPMD isnot a decoupling of electronic orbitals from theions but
c
e onlyadecouplingofasubsetoftheorbitalvibrationalmodesfromtherestofthenecessarily-coupled
m system of orbitals and ions. Recent work ( J. Chem. Phys. 116 ,14 (2002) ) haspointed out that,
due to the orbital-ion coupling that remains once adiabatic-decoupling has been achieved, a large
-
t valueofthefictitiousmassµcanleadtosystematicerrorsinthecomputedforcesinCPMD.These
a errorsarefurtherinvestigatedinthepresentworkwithafocusonthosepartsoftheseerrorsthatare
t
s notcorrected simplybyrescalingthemassesoftheions. Itissuggested thatanycomparison ofthe
. efficiencies of Born-Oppenheimer molecular dynamics (BOMD) and CPMD should be performed
t
a at a similar level of accuracy. If accuracy is judged according to the average magnitude of the
m systematicerrors in thecomputed forces, theefficiency of BOMD compares more favorably to that
- of CPMD than previouscomparisons havesuggested.
d
n
o I. INTRODUCTION follows,thisapproachtoAIMDwillbereferredtoasCar-
c Parrinello molecular dynamics (CPMD). In CPMD the
[
SinceitsintroductionbyCarandParrinelloin1985[1], equationsofmotionoftheions(withcoordinatesRI and
1 masses M ) and the electronic orbitals |ψ i are derived
the field of ab initio molecular dynamics (AIMD) has I i
v from the classical lagrangian
grown very rapidly. Although many of its initial suc-
0
3 cesses involved its application to questions in condensed L=µ hψ˙ |ψ˙ i+ 1 M R˙2−E[{ψ },{R }] (1)
1 matterphysicsandmaterialsscience,ithasnowbeenap- i i 2 I I i I
1 plied with success to a wide rangeof problems in diverse Xi XI
0 areasofsciencesuchaschemistry,biologyandgeophysics and the motion is further constrained to the hyper-
6
and contributed a great deal to our understanding of surface defined by the orbital orthogonality condition
0
/ these fields. The central idea of the method is that the hψi|ψji = δij;∀i,j[7]. Inspired by the ideas of Car
at Kohn-Sham orbitals[2] that describe the electronic state and Parrinello, an alternative approach to AIMD was
m of a system within density functional theory[3] (DFT) proposed[8] in which the electronic orbitals are kept as
are evolved simultaneously with the ions as classical de- close as possible to the instantaneous groundstate using
-
d grees of freedom. In the original approach of Car and a combination of extrapolation from earlier times and
n Parrinello, an inertial parameter is introduced in order explicit minimization. This general approach will be re-
o to associate a Newtonian dynamics with the electronic ferredtohereasBorn-Oppenheimermoleculardynamics
c orbitals. This parameter is an unphysical quantity that (BOMD).
:
v iscommonlyreferredtoasthefictitiousmass,µ. Whileit It has been known since its introduction that CPMD
Xi isknownthattheintroductionofthisextrainertiaintoa is an approximationto a fully convergedBOMD simula-
system affects its dynamics, the extent to which dynam- tioninwhichtheelectronicsystemisexactlyinitsground
r
a ics are altered and the way in which they are alteredare stateateveryinstantoftheionicmotion,howeverCPMD
only beginning to be studied in detail[4, 5, 6]. is frequently preferred over BOMD because it is deemed
In this paper some of the theoretical ideas that un- computationally more efficient. In CPMD simulations,
derlie the original, and still widely-used, Car-Parrinello because the orbitals are treated as classical degrees of
approach to ab initio molecular dynamics are examined freedom on the same footing as the ions, there exists an
andthe possible importance ofthe fictitious masseffects energy(whichincludesthekineticenergyassociatedwith
are investigatedby studying the example of ice. In what the fictitious classicalmotionofthe orbitals)that is per-
fectly conserved as long as the time step that is used
is small enough to adequately integrate the equations of
motion. Ontheotherhand,whiletheconservedquantity
∗The following article has been accepted by The Journal of inBOMDisthephysicallymeaningfultotalenergyofthe
Chemical Physics. After it is published, it will be found at system of electrons and ions, its conservationin practice
http://jcp.aip.org/
is always imperfect. There always exists a drift in the
†Current address: The Molecular Foundry, Lawrence Berke-
calculated total energy of a BOMD simulation because
ley National Laboratory, CA 94720; Electronic address:
[email protected] the electronic orbitals are never perfectly converged to
2
the ground state. While the magnitude of this drift can to zero on time scales that are relevant to ion dynamics.
systematically be reduced by more fully converging the The fact that a dynamics is associated with electronic
minimization of the electronic system at each time step, orbitals means that, according to the equipartition the-
any such improvement in accuracy must always be bal- orem, the subsystems of orbitals and ions should equili-
ancedbythecomputationalexpenseinvolved. Theargu- bratewithoneanotherandreachacommontemperature.
ment has frequently been made[9, 10, 11], that for many Ithas been observed,however,thatfor systems inwhich
systems CPMD is computationally superior to BOMD thereisasubstantialgapintheKohn-Shamenergyspec-
because the computational overhead required to make trumbetweenoccupiedandunoccupiedelectronicstates,
BOMD conservative within acceptable limits is greater there is no systematic net transfer of kinetic energy be-
than that required to perform a perfectly conservative tweenorbitalsandionsandthereforetheorbitalsremain
CPMD simulations with a value of µ chosen in line with at a temperature that is very low relative to the ionic
conventional wisdom[5, 9, 12]. Since BOMD can be as temperature. This low temperature means that orbitals
fast as, or faster than, CPMD depending on the size of do not have the kinetic energy required to leave the en-
the energydriftthatis tolerated,animplicit assumption ergy basin of the electronic ground state.
of this argument is that BOMD can only match the ac- The non-equilibration of orbitals and ions has been
curacy of CPMD when this drift of energy is very small. explained[12]asbeingduetoagapinthevibrationalfre-
In this paper, this assumption is questioned by pointing quency spectrum of the coupled orbital-ion system that
outthat,forcommonlyusedvaluesofthe fictitious mass separatestheultrahighfrequencyorbitalmodesfromall
µ, therecanexistlargeerrorsinthe forcesonthe ionsin the lowerfrequencymodesofthe system. Theultrahigh
a perfectly conservative and well-controlled CPMD sim- frequency oscillations of the orbitals are “adiabatically
ulation. decoupled” from the rest of the system.
Theemphasisofthispaperisonthemicroscopicdetails While the arguments presented above appear to work
of CPMD simulations, i.e. on the ability of the CPMD well for silicon[4, 12], little work has been done to ver-
approach to evolving the electronic orbitals in time to ify that CPMD forces averageto the ground state forces
produce (on average) the correct ground state forces on for other systems. Furthermore, some important details
the ions. The important question of how deviations of of the arguments of Pastore et al. have received little
forcesfromthecorrectforcesrelatetopossibledeviations attention: Adiabatic decoupling is not a decoupling of
in thermodynamic quantities from those that would be orbitals from ions. Rather, it is only a decoupling of
obtained using the correct forces is not addressed in the the quasi-independent orbital modes that consist of ul-
present work. tra high frequency oscillations of the orbitals from the
Car and Parrinello introduced their method in the lower frequency modes of the system. The lower fre-
context of coupling Kohn-Sham density functional quency modes of the system generally involve both ions
theory[2, 3] to moleculardynamics,however,it has since and orbitals. The slow (ionic time scale) component of
been applied successfully to other electronic structure the orbital motion is frequently completely neglected in
methods[18]. Ithasalsobeenappliedinapurelyclassical discussions of the CPMD method[5, 11], however,and it
context to evolve induced dipoles[19, 20] and fluctuating is often stated that the lowest frequency motion in the
charges[21] in molecular dynamics simulations. None of orbital subsystem is approximately related to the Kohn-
these applications will be discussed in the present work Sham energy gap of the system E through [5, 9, 11, 13]
g
but, clearly, many of the concepts discussed here have
analogies within these methods.
ω ≈ 2E /µ (2)
g
q
II. BACKGROUND . Inaninsulator,thisfrequencyisveryhighcomparedto
typical ionic frequencies. Clearly, this view is not com-
It has been argued, and demonstrated for the case of patible with the electronic orbitals remaining at or near
silicon[9, 12], that although there exist small instanta- the electronic ground state because the ground state by
neousdeviationsoftheforcesinaCPMDsimulationfrom definition evolves on ionic time scales.
the groundstate forcesatthesameionic positions,these The slow component of orbital motion due to the evo-
deviationsaveragetoclosetozeroonafemtosecondtime lution of the ground state is alwayspresent and, as a re-
scaleandthereforedonotresultinseriouserrorsinther- sult,thereisalwaysacontinuousexchangeofenergyand
modynamic properties[9, 11, 12, 23]. The explanation momentum between orbitals and ions. Furthermore, for
proposed for this behaviour[12] is that the electronic or- manysystems,thisionic time scalecomponentofthe or-
bitals perform small-amplitude ultra high frequency os- bital motion has been found to be appreciable[4, 14, 15].
cillations around an equilibrium that slowly evolves as It will be demonstrated in section V of the present work
the ions move. If this equilibrium corresponds to the that this ionic-time scale contribution to the motion of
electronic ground state, and the amplitudes of the oscil- the orbitals actually accounts for the vast majority of
lationsaresmall,theforcesontheionsoscillateaboutthe their classical “fictitious” kinetic energy in the simula-
true ground state forces and errors in the forces average tions of ice reported here.
3
The argument of Pastore et al.[12] that forces in ture of the individual molecules from their average elec-
CPMDoscillateaboutthe true groundstate forcesrelies tronic structure are less likely to occur in ice than in liq-
on the assumption that the ultra high frequency oscilla- uid water. As mentioned above, a rigid electronic struc-
tionsoftheorbitalsareabouttheelectronicgroundstate, ture means that errors associated with µ only affect dy-
however,for any finite µit will be shownin sectionIIIA namicalpropertieswhereaschangesinthelocalelectronic
that this is not the case. As a consequence of the ionic structureduetointeractionsbetweenmoleculesmaylead
time scale motion of the electronic orbitals in CPMD, to errors in thermodynamic properties. For these rea-
the orbitals do not oscillate aroundthe groundstate but sons,the errorsinCPMDreportedherefor ice aretaken
about averagevalues that are displaced from the ground to be indicative of similar or more serious problems in
state. Furthermore, although the resulting errors in the liquid water.
forces are small for materials in which valence electrons The fictitious mass effects described here and in refer-
aredelocalizedandthereforehavealowquantumkinetic ence 4 should not be confused with those described by
energy (such as silicon), for other materials the forces Grossmanet al.[5]. TheworkofGrossmanet al. pointed
can deviate strongly from their ground state values if outsome problems thatoccurredinsimulations ofwater
commonly used values of the fictitious mass parameter duetotheadiabaticdecouplingconditionbreakingdown
are used[4], however, these errors can systematically be and energy passing continuously into the ultra high fre-
reduced by reducing the magnitude of µ. quency orbital modes from the lower frequency modes.
It has previously been proposed that an effect of the The present work is only concerned with fictitious mass
orbital-ion coupling is simply to rescale the masses of dependent errors that are present in a well-performed
the ions[4, 15, 16], an effect that does not alter ther- Car-Parrinellosimulationwheretheadiabaticdecoupling
modynamic properties, however, theoretically, the error condition is perfectly maintainedand therefore the ultra
induced by orbital-ion coupling only reduces to a mass high frequency orbital modes are energetically isolated.
correctionin the limit in whichions do not interactwith
oneanother,i.e. inthelimitofaninfinitelydilutegas[15].
Forcondensedsystemsmost,butnotall,oftheerrorcan III. THEORY
be corrected by simply rescaling the masses of the ions.
Onegoalofthepresentpaperistoreview,toextend,and Inthissectionthedynamicsoftheorbitalsandtheions
to help clarify some of the analysis of reference 4. An- thatresultfromtheCar-Parrinellolagrangianofequation
other goal is to stimulate further investigation by using 1areanalysedinorderto gainabetter understandingof
‘heavy’ ice as a test system to highlight the importance the different contributions to the errors in the forces on
oftheeffectsdescribedinreference4withparticularem- the ions.
phasis on the part of the errors in the forces that is nei- The equations of motion of the orbitals and ions that
ther oscillatory on a very short time scale, or equivalent are derived from equation 1 are
to a simple renormalization of the masses of the ions. It
is shown that, for commonly used values of µ, the mag- µψ¨(r)=−δE[{ψi},{RI}] (3)
nitude ofthispartofthe forceerrorsismuchlargerthan i δψ∗(r)
i
errors in the forces in typical BOMD simulations.
Water is arguably the most important system studied
with ab initio molecular dynamics. CPMD has been ap- M R¨α =−∂E[{ψi},{RI}] (4)
pliedextensivelytothe studyofliquidwater[5,6,10,17] I I ∂Rα
I
and to the study of biological systems where water is
where the orbital orthonormality condition is implicit
present. The water molecules in ice have much the same
in the functional derivatives δ/δψ∗(r) and all functional
propertiesastheydoinliquidwater,butice hasalarger i
derivatives throughout this paper. The notation ∂/∂Rα
electronic energy gap making it easier to simulate with I
indicates that a partial derivative is performed with re-
CPMD.Thegoalhereistoinvestigatetheaccuracyofthe
specttoRαwithallotherionCartesiancoordinatesfixed
bare Car-Parrinello method without any further sources I
and with the orbitals ψ fixed. In the remainder of the
of error. In simulations of liquid water the smaller en- i
paper it will also be necessary to use the notation ∇α
ergy gap between occupied and unoccupied states and I
to denote a partial derivative with respect to Rα with
the high frequency of the O-H stretch vibrational mode I
allother ionCartesiancoordinatesfixedbut not with the
can result in a continuous drift of kinetic energy into
orbitals fixed.
the orbital subsystem that would further complicate our
Equation 4 produces the correct Born-Oppenheimer
analysis[5]. Theseproblemsareavoidedherebystudying
‘heavy’ ice rather than water. The larger energy gap of dynamics if the set of orbitals {ψi} at which the deriva-
ice and the lower frequency of the O-D stretch in D O tives ∂E[{ψi},{RI}] areevaluatedisthesetofgroundstate
2 ∂Rα
ensures that no such coupling occurs in the simulations orbitals {ψgI.s.}, i.e. the set of orbitals that minimize
i
reported here. theKohn-ShamenergyfunctionalE[{ψ },{R }]forfixed
i I
Thefactthatthewatermoleculesareconfinedtoalat- ionicpositions{R }. Thequalityoftheforcesontheions
I
tice means that large deviations of the electronic struc- thereforedependsontheabilityofequation3toproduce
4
a dynamics for the orbitals that keeps them close to this positions {R }, the velocities {R˙ }, and on the orbitals
I I
instantaneous ground state as the ions move. In section {δψ }. The dependence of ψ¨g.s.(r) on the {δψ } comes
i i i
IIIA,therefore,thedynamicsoftheorbitalsareanalysed in via the dependence of the ionic accelerations on the
in the regime of small deviations from the ground state {δψ } through equation 4. However, equation 4 may be
i
in order to gain a better understanding of the different rewritten as
contributions to the errors in the forces on the ions. In
section IIIB the effects of the orbital dynamics on the R¨α = − 1 ∂ E + δE δψ (r)
forces on the ions are discussed. I MI ∂RIα( (cid:12)(cid:12)g.s. Xj Z "δψj(r)(cid:12)(cid:12)g.s. j
(cid:12) (cid:12)
δE (cid:12) (cid:12)
+ δψ∗(cid:12)(r) dr+O(δψ2) (cid:12)
A. The dynamics of the orbitals δψj∗(r)(cid:12) j # )
(cid:12)g.s.
(cid:12)
Asstatedpreviously,acommonviewofCar-Parrinello = − 1 ∇(cid:12)(cid:12)αE +O(δψ2) (8)
dynamics is that the Car-Parrinello orbitals {ψi} oscil- MI I (cid:12)
late rapidly around the electronic ground state so that (cid:12)g.s.
(cid:12)
deviationsfromthe groundstate averageto closeto zero from which it is c(cid:12)(cid:12)lear that the dependence of R¨Iα, and
on the longer time scales that are relevant to ionic dy- therefore of ψ¨g.s.(r), on {δψ } is at higher than linear
i i
namics. Inthissectionitisdemonstratedthatthisisnot
order. To first order in {δψ }, therefore, equations 3, 5
i
the case if µ is large, and an expression for the average
and 6 can be combined to give
value of the electronic orbitals is derived.
For a given set of ionic coordinates {RI}, the elec- δψ¨(r)=−1 δ2E δψ∗(r′)
tvreonnieicntgrtoouinndtrostdautceeisthweeqlluadnefitintyed|.δψIitii≡s t|hψeirie−for|ψeigc.os.ni-, i µXj Z "δψj∗(r′)δψi∗(r)(cid:12)(cid:12)(cid:12)g.s. j
fai.rteo.matghiitevsedngervioniuasnttiadonntstooafftetahveCalaiutrhe-.PCaTrarrhi-nePealslrtoraitdneyelnoloafmothircbesitscayalsn|tψebmiei +δψj(rδ′2)Eδψi∗(r)(cid:12)(cid:12)g.s.δψj(r′)#dr(cid:12)(cid:12)′−ψ¨ig.s.(r) (9)
(cid:12)
completely specified by specifying the values of the dy- Byassumingthatthe(cid:12){ψ }areeitherreal,orarecomplex
namical variables {R },{δψ } and their time derivatives (cid:12) i
I i but with phase factors that do not vary as the system
{R˙ },{δψ˙ } at that instant. In what follows, the dy-
I i evolves, this equation can be evaluated as
namics of the {δψ } will be analyzed. An assumption
i
that will be made throughout the analysis is that the δψ¨(r)=−1 K(ir,j r′)δψ (r′)dr′−ψ¨g.s.(r) (10)
{δψi} are small enough that certain ψ-dependent quan- i µ j i
j Z
tities are well described by truncations at linear order X
in δψ of Taylor expansions about the electronic ground where
state.
A Taylor expansion of δE/δψi∗(r) about the ground K(ir,j r′) = fj δijδ(r−r′)H
state gives (cid:12)
(cid:12)g.s.
(cid:12)
δE δE δ2E δv(r) (cid:12)
= + δψ∗(r′) + 2 ψ∗g.s.(cid:12)(r′)ψg.s(r) (11)
δψi∗(r) δψi∗(r)(cid:12)(cid:12)g.s. Xj Z "δψj∗(r′)δψi∗(r)(cid:12)(cid:12)g.s. j δn(r′)(cid:12)(cid:12)g.s. j i !
(cid:12) (cid:12) (cid:12)
δ2(cid:12)E (cid:12) H istheKohn-Shamsin(cid:12)gleparticleHamilitonian,v(r)=
+ δψj(r′)(cid:12)δψi∗(r)(cid:12) δψj(r′)#dr′+O(δ(cid:12)ψ2) (5) vH(r)+vXC(r) is the su(cid:12)m of the Hartree and exchange-
(cid:12)g.s. correlationpotentials,andf istheoccupationnumberof
(cid:12) j
where the first term on th(cid:12)e right hand side vanishes by thejth orbital. ItisclearthatK(ir,j r′)isindependent
(cid:12)
definition of the ground state. of the {δψi} and therefore evolves on ionic timescales.
The secondtime derivative of ψi(r) may be written as If the {δψi} are sufficiently small, their motion is gov-
erned by equation 10 which describes driven coupled os-
ψ¨(r)=δψ¨(r)+ψ¨g.s.(r) (6) cillationsofthe{δψ },wherethedrivingtermisgivenby
i i i i
the acceleration of the ground state {ψ¨g.s.}. By means
where of a suitable unitary transformation, K(ii r,j r′) can be
diagonalized and equation 10 recast into a simpler form
ψ¨ig.s.(r)= R¨Iα∇αIψig.s.(r)+ R˙IαR˙Jβ∇αI∇βJψig.s.(r) involving transformations of the functions {δψi(r)} and
XI,α I,Xα,J,β {ψ¨g.s.(r)}, however, for the purposes of the present dis-
i
(7) cussion it is sufficient to notice that, by reducing the
value of the free parameter µ, one can make the magni-
Clearly, ψ¨g.s.(r) is finite unless the ions are not moving. tude of K/µ arbitrarily large thereby increasing the fre-
i
At a given instant in time ψ¨g.s.(r) depends on the ionic quencies of the oscillations to “ultra high” values, while
i
5
thedrivingtermψ¨g.s.(r)remainsconstant. Itisassumed Forlargervaluesofµ, orifthe particularsofasimula-
i
here that a value of µ has been chosen such that there tionaresuchthatultrahighfrequencyoscillationsofthe
exist three distinct time scales τ ≪ τ ≪ τ in the {δψ }havealargeamplitude,contributionstotheerrors
S I L i
Car-Parrinello dynamics. τ is a short time scale that in the forces that are of higher than linear order may
S
is comparable to the period of the ultra high frequency play a role,evenif the separationof time scales τ ≪τ
S L
oscillation of the orbitals, τ is the long time scale on can still be maintained, however,such effects will not be
L
whichtheionsmove,andτ isanintermediatetimescale. discussed in the present work.
I
When averagedovera time scale of order ∼τ , equation
I
10 becomes
B. The forces on the ions
1
δψ¨(r)=− K(ir,j r′)δψ (r′)dr′−ψ¨g.s.(r) (12) Inthis section,theimpactofthe orbitaldynamicsdis-
i µ j i cussed in the previous section on the forces on the ions
j Z
X
arediscussed. Ofprimaryconcernis the systematic part
where the averaging only significantly affects δψ¨(r) and of the error in the forces due to the displacement of the
i
δψ (r) because both K and ψ¨g.s.(r) vary on the much averagevalues of the orbitals from the groundstate, δψ.
i i This contribution is now derived.
longer time scale τ . If the orbitals are to remain close
L
Equation 4 can be written as
to the ground state, then δψ¨ ≈0, ∀i, because the short
i
∂E
time scale dynamics of the orbitals should be oscillatory Fα = −
rather than accelerating systematically in one direction. I ∂Rα
I
By inverting equation 12, the average deviation of the δE ∂ψ∗(r) δE ∂ψ (r)
wavefunctions from their ground state values can there- = −∇αE+ i + i dr
I δψ∗(r) ∂Rα δψ (r) ∂Rα
fore be written as i Z " i I i I #
X
(14)
δψ (r′)=−µ K−1(j r′,ir)ψ¨g.s.(r)dr (13)
j i Substitution of equation 3 yields
i Z
X
δψ (r′) is, in general, non-zero and dependent only on Fα =−∇αE− µ ψ¨(r)∇αψ∗(r)+ψ¨∗(r)∇αψ (r) dr
j I I i I i i I i
ionic positions and velocities. Therefore, the electronic i Z " #
X
orbitals do not average to their ground state values on (15)
time scales much longer than the time scale of their ul-
since ∇αψ (r) = ∂ψ (r)/∂Rα. The first term on the
tra high frequency oscillations but shorter than the time I i i I
scale of ionic motion. On the other hand, δψ (r′) de- right-handsideofequation15isnowexpandedinaTay-
j
lor series about the ground state :
pends linearly on the fictitious mass µ and therefore can
be madearbitrarilysmallbyreducingitsvalue. Further- δE
mthoerev,aalusepoofinµtedthoeuftrebqyuPenacsiteosreofetthael.[o1s2c],illbaytidonecsreoafstinhge −∇αIE = −∇αI(cid:26)E(cid:12)(cid:12)g.s.+Xi Z "δψi∗(r)(cid:12)(cid:12)g.s.δψi∗(r)
(cid:12) (cid:12)
δψi about their slowly evolving average values can be δE (cid:12) (cid:12)
mthaedmefsroomlartgheetlhoawtervi(ritounailcl)yfnroeqeuneenrcgyymisotdreasnosffertrheedotro- + δψi(r)(cid:12)(cid:12)g.(cid:12)s.δψi(r)#dr+O(δψi2)(cid:12)(cid:27)
bitals and ions. If the amplitudes of these oscillations = Fα +(cid:12)(cid:12)0+O(δψ2) (16)
are initially small, therefore, they will remain small un- BOI (cid:12) i
whereFα istheα-thCartesiancomponentoftheBorn-
less the system is artificially perturbed (by, for example, BOI
Oppenheimer forceonatomI. As pointedoutin section
discontinuously changing the velocities of the ions).
To summarize this discussion of the orbital dynam- IIIA, on time scales relevant to ionic motion ψ¨i(r) =
ics: it has just been shown that for small values of the δψ¨(r)+ψ¨g.s.(r) = ψ¨g.s.(r), therefore, equations 7, 15,
i i i
fictitious mass,µ, the dynamicsofthe Car-Parrinelloor- and 16 can be combined to give the error in the Car-
bitals can be described as consisting of ultra high fre- Parrinello force on timescales ( > τ ) relevant for ionic
I
quency oscillations about average values that evolve on motion as[25]
ionic timescales and that are given by equation 13. For
∆Fα = Fα−Fα
sufficiently small values of µ ( and {δψi}) the deviations I I BOI
of the forces on the ions from their values at the elec-
= −2µ ℜ R¨β∇αψ∗g.s.(r)∇βψg.s.(r)
tronic ground state can therefore be separated into an J I i J i
(
ultrahighfrequencyoscillatorypartthataveragestozero Xi Z (cid:20)XJ,β
onthe time scalerelevanttoiondynamics andasystem-
aticpartthatdoesnot. Thesystematicpartoftheerrors + R˙JβR˙Kγ ∇αIψi∗g.s.(r)∇γK∇βJψig.s.(r) dr
in the forces on the ions is the main focus of the present J,β,K,γ (cid:21) )
X
work. (17)
6
This expression is valid in the limit of small {δψ } and these increased values of the masses. In reference 4 this
i
depends only on ionic positions and velocities. It has no wasdemonstratedforthecaseofMgO:Thetruetemper-
first order dependence on the {δψ }. ature of the oxygensubsystem was that calculated using
i
ThissystematicerrorintheCar-Parrinelloforcesarises amassfortheoxygenionsofM +∆M ,whereM was
O O O
from the fact that the electronic ground state evolves the bare oxygen mass of 16 a.m.u. and ∆M a correc-
O
on ionic time scales. As ions move, they must push tionthataccountedfortheinertiaoftheorbitalsmoving
the orbitals around. The orbitals carry inertia and so rigidly with the oxygen ions. Car-Parrinello simulations
the ions must impart momentum to them in order for thatneglectthe correctionto the temperature due to in-
them to follow the ground state. This means that the creased ionic masses have effectively been performed at
ions lose momentum to the orbitals, i.e. the ions en- temperatures that are higher than those reported. Fur-
counter a resistance to their motion. It is as though the thermore, this correction is important if more than one
ions move through a viscous, inhomogeneous and time- ionic thermostat is used in a simulation. For example, if
varying medium. mass correctionsarenot accountedforin the calculation
The right hand side of equation 17 depends on the of temperature, and if a different thermostat is applied
details of the electronic ground state, as well as the ve- to each ionic species, the result can be that each species
locities of the ions, and so it is difficult to make general is effectively maintained at a different temperature. The
statementsabouthow∆F affectsthedynamicsandther- definition of temperature may be particularly important
modynamics of a simulation. In the absence of further if so-called “massive thermostating” is employed[26], in
information there is no reason to expect either dynam- which each degree of freedom is coupled to its own ther-
ical properties or thermodynamic properties to be cor- mostat.
rect if µ is large enough to make the error of equation In a system in which there is only one atomic species
17 significant. However, in the limit in which individual and every atom is in an identical chemical environment
ions are spherically-symmetric and do not interact with andiftherigid-ionapproximationworkswell,dynamical
one another, this error has been shown to reduce to the propertiescanbecorrectedsimplybyrescalingtimebya
form ∆FIα =−∆MR¨Iα, where ∆MI is a constant that is factor M/(M +∆M). If the rigid-ion approximation
proportional to the total quantum kinetic energy of the works, but there is more than one atomic species, dy-
p
electronic states on ion I [4, 15]. In other words, if ions namicalpropertiescanonlybemadecorrectbyrescaling
don’t interact with one another (i.e. in the limit of the themassofeachatomapriori. Inotherwords,thesimu-
infinitely dilute gas[15]) the systematic errorin the force lationisperformedwithanominalmassforeachionI of
due to the displacement of the average values of the or- M −∆M so that the mass of the ion once ’dressed’by
I I
bitals from the ground state, {δψi}, has the same effect theheavyCar-ParrinelloorbitalsisthecorrectmassMI.
assimplyincreasingthemassesoftheionsbyanamount If this is not done, large errors in dynamical properties
∆MI. Thisresulthadbeennoticedmuchearlier[14],and can result[4].
describedinreference16. The moregeneralresultofref-
Strictly speaking, the rigid ion approximation is not
erence 4was alsoindependently foundby Blo¨chl[15]. An
applicable whenions interactwith one another. In other
equation of motion for the ions that is an approxima-
words,whenthereexistelectronicorbitalsφthatdepend
tion to the equation of motion of the ions in CPMD is
on the positions of more than one ion, i.e. for which
therefore given by[4, 25]
∇ φ 6= 0 and ∇ φ 6= 0 for some ions I 6= J. Of course,
I J
∂E[{ψg.s.},{R }] this is generally the case for all condensed matter and
(M +∆M )R¨α =− i I =Fα (18) therefore for all systems of interest. What this means
I I I ∂Rα BOI
I is that the dynamics and thermodynamics of CPMD al-
ways differ from BOMD. Fortunately, for all condensed
The degree to whichequation18describes the motionof
systems tested to date, most of the error of equation 17
the ions in CPMD depends on the degree to which the
canbecorrectedbyapplyingmasscorrectionstotheions.
electronic structure of the condensed system in question
The concern that remains, however, is that even though
can be described in terms of a linear superposition of
itismuchsmaller,thepartoftheerrorthatremainsonce
tightly bound ionic orbitals that are transported rigidly
masscorrectionshavebeenappliedmaynotbenegligible.
(i.e. withoutchangingshapeinresponsetotheirenviron-
ment) as the ions move. This approximation to CPMD A further difficulty associated with correcting the
is known as the “rigid-ion” approximation[4]. Equation masses of the ions is that the quantum kinetic energy
18 can be understood as arising from the fact that such ofthe electronic states onagivenatomicspecies is not a
tightly bound orbitals have inertia (∝ µ) and that this well defined quantity. This is because each orbital can-
adds to the total inertia of the ions: the ions have to not, in general, be associated with only one atom. A
carry the “heavy” orbitals. method must therefore be devised to approximate the
In a dynamics described by equation 18, i.e. a Born- mass corrections. In the simulations of D2O presented
Oppenheimer dynamics but with increasedionic masses, here, the mass corrections that are used are those that
thermodynamic properties are unchanged. This is only minimize the errors in the forces on the ions.
the case, however, if thermodynamic quantities such as In order to differentiate between the part of ∆F that
temperature and thermal pressure are calculated using amounts to a simple mass correction, and the part that
7
care is taken to ensure that the velocities of the ions are
notchangeddiscontinuously,thelargeamplitudeofthese
oscillations disappears and the contribution of δF to the
instantaneous error in the Car-Parrinello force can be
neglected. At higher temperatures, however, large am-
plitude oscillations may always be present, as was the
case in simulations of MgO in reference 4.
InsectionVI,thesystematicpartoftheerror∆F,will
be examined. It has been assumed in the past[10] that
∆F≈∆Fm and that ∆Fr can be neglected. Therefore,
in section VI, the validity of this assumption will be in-
vestigated. The degree to which ∆F can be eliminated
by applying mass corrections to the ions will be tested.
Inotherwords,themagnitudeof∆Fr,whichhasthepo-
tential to alter thermodynamic properties and for which
no correction yet exists, will be examined.
IV. SIMULATION DETAILS
FIG. 1: Some samples of force components on oxygen ions
along a 4 fs segment of CPMD trajectory 1 ps after temper- Inthe CPMDsimulationsofD2Oreportedherenorm-
ature has been adjusted using velocity rescaling. The Car- conserving pseudopotentials have been used to describe
Parrinello force and the Car-Parrinello force once corrected oxygen[27] and deuterium[28] and the valence electronic
with the best constant mass correction is compared to the orbitals were expanded in plane waves with a maximum
ground state (Born-Oppenheimer) force, FBOI at the same energy of 70 Rydberg. Simulations were performed on
ionic positions. a 24×24×12 a.u. simulation box containing 32 D O
2
molecules. Only the Γ−point was used to sample the
Brillouinzone. A gradient-correctedfunctional wasused
remains once this mass correction has been applied, the
to treat the effects of exchange and correlation[30]. In a
quantities ∆Fm and ∆Fr are introduced where ∆Fm =
preliminary CPMD simulation of liquid water mass cor-
−∆MR¨ and ∆Fr = ∆F−∆Fm is the remaining error.
rections of ∆M = 6.77µ a.u. and ∆M = 0.213µ a.u.
The expression for ∆F in equation 17 has been derived O D
, were found for the oxygen and deuterium ions, respec-
bytakinganaverageovertheultrahighfrequencyorbital
tively. These mass correctionswere found by calculating
oscillations,andthereforeamorecompleteexpressionfor
the ground state forces at selected points along a trajec-
the Car-Parrinelloforce is
tory and by minimizing the error in the CPMD forces
F = F ({R })+∆F ({R },{R˙ }) relative to these ground state forces. These mass cor-
I BOI I I I I rections were then used in all of the CPMD simulations
+ δF ({δψ })+O(δψ2)
I i of ice reported here by reducing the values of the ionic
≈ FBOI +∆FmI +∆FrI +δFI (19) masses used in calculating the acceleration from New-
where δF is the oscillatory part of the force which, for ton’s equation of motion i.e.
I
smFailgluerneou1gdhemµ,oanvsetrraagteesstthoazte,riondoenetdi,mtehescdailffeserτeInc≫e bτSe-. R¨Iα =−M −1∆M ∂E[{ψ∂iR},α{RI}] (20)
tween the Born-Oppenheimer and Car-Parrinello forces I I I
consistsofbothasystematicpart∆F,andanoscillatory The values of the fictitious mass used in the simulations
part δF. The systematic part appears to be mostly cor- reported here were 900 a.u. and, for comparison, 100
rectedbythe applicationofmasscorrectionsto the ions. a.u. The equations of motionfor electrons andions were
The forces in figure 1 are those from a Car-Parrinello integratedusingamoleculardynamicstimestepof6a.u.
simulation in which velocity rescaling was employed for for the µ = 900 a.u. simulations and 2 a.u. for the µ =
0.5 ps to bring the temperature from T ∼ 120K to 100 a.u. simulations. No mass preconditioning scheme
T ∼ 220 K. The velocities of the ions were adjusted was used in any of the simulations reported here.
only 4 times, in total, during this 0.5 ps . The system Much care has been taken to minimize errors in the
was then equilibrated for one picosecond after which the simulations reported here. For example, as pointed out
forces were examined along the 4 fs segment of trajec- byRemlerandMadden[31],anddiscussedinsectionIIIB
tory and compared to the forces at the same ionic po- it is important that the velocities of the ions are not
sitions, but with the orbitals in their ground state, i.e. changed discontinuously as this can give a sudden kick
the Born-Oppenheimer forces F . The amplitude of to the electronic orbitals that results in large-amplitude
BOI
the oscillations,δF, infigure 1 appears worryinglylarge. ultra high frequency oscillations of the orbitals around
It will be shown in sections IV and VI, however, that if their equilibria. Great care has therefore been taken to
8
FIG. 2: A schematic illustrating thesequence of Car-Parrinello molecular dynamics simulations that were performed.
eliminate this source of error from the simulations re-
ported here with the result(evident in figures 7 and 8 of
sectionVI)thatinstantaneouserrorsintheforcesdue to
ultra high frequency oscillations of the orbitals (i.e. δF
) are small enough relative to the systematic errors of
equation 17 that they may safely be neglected.
TheinputcoordinatesfortheCPMDwereobtainedby
performingaverylongmoleculardynamicssimulationof
ice at low temperature (∼ 100 K ) using an ab initio
parametrizedpolarizable atomistic potentialof the same
form as that constructed for silica in reference 32. This
potential does not provide a very realistic description of
water but it was deemed preferable to using randomized
initial coordinates. From these initial conditions a se-
quence of CPMD simulations were performed as shown
schematically in figure 2. CPMD simulations were be-
gunwiththe electronsintheir groundstate andthe ions
at zero velocity. A fictitious mass of µ = 900 a.u was
used and therefore the ionic masses of oxygen and deu- FIG. 3: The temperatures of the oxygen and deuterium
subsystems during the 3.5 ps equilibration simulation with
terium were rescaledto M −∆M =12.659a.m.u and
O O
µ = 900 a.u. In (a) the temperatures have been calculated
M −∆M =1.895a.m.u,respectivelyandtemperature
D D usingmassesfortheionsthathavebeencorrectedtoaccount
was calculated using the true ionic masses of M = 16
O for the extra inertia due to the weight of the Car-Parrinello
a.m.uandM =2a.m.u. Afterhalfapicosecondofsim-
D orbitals (in this case 16 a.m.u and 2 a.m.u. due to the a
ulation in the microcanonical ensemble, the ions were at priori rescaling of the masses). The oxygen and deuterium
a temperature of 120 K. At this point two separate con- subsystems appear to be out of thermal equilibrium. In (b)
tinuations of the simulation were performed. The first the temperatures have been calculated in the standard way
continuation was performed in order to demonstrate the using the bare ionic masses (12.659 a.m.u and 1.895 a.m.u).
different contributions to the error in the forces on the Theoxygenanddeuteriumsubsystemsappearperfectlyequi-
ions and the effect of velocity rescaling on these forces, librated.
and has been discussed in section IIIB. In the second
continuation,aNos´e-Hooverthermostat[33]wasattached
to the ions and the temperature was smoothly increased corrections are used the oxygen and deuterium subsys-
to approximately250K during a further half picosecond tems are at different temperatures indicating that the
ofsimulation. Thesystemwasthenequilibratedwithout system is not well equilibrated. If mass corrections are
athermostatfor3.5picoseconds. Asshowninfigure3the not used, the system appears perfectly equilibrated. At
temperaturesofthesubsystemsofoxygenanddeuterium thispoint,thesimulationwasagaincontinuedintwosep-
ionswerecalculatedandcomparedwithandwithoutthe arate simulations. In one of these simulations the value
use of the corrections to the masses of the ions in the of µ = 900 a.u. was used as before. In another simula-
definition of temperature. As can be seen, when mass tionthe fictitious masswaschangedto µ=100a.u. The
9
is an approximation and deviations from the rigid-ion
limit always occur. In the opposite limit to the rigid-ion
approximation, i.e. in the limit of very weak orbital-ion
coupling, the temperature of the ions should be calcu-
lated with the bare ionic masses. For systems that are
not perfectly ionic, therefore, the correct definition of
temperature is unclear. While the results of section VI
suggestthataconstantmasscorrectionfortheionsisnot
perfectly appropriate for D O, the results of section V
2
suggest that the rigid-ion approximation does a remark-
ably good job of estimating the fictitious kinetic energy
oftheorbitals. Inaddition,the factthatduringthefinal
5.5 ps of simulation the oxygen and deuterium subsys-
tems remained at the same mass-corrected temperature
suggeststhat itis appropriateto calculatedtemperature
using rescaled masses.
Thefactthatreducingthevalueofµappearedtofacil-
itatethermalizationmayindicatethatthethermalization
problemwasanartifactofthefictitiousmassµ,however,
FIG. 4: The temperatures of the oxygen and deuterium sub- further work is required to clarify this issue.
systemscalculatedusingthecorrectedionicmasses(M+∆M)
in the two continuations of the simulation shown in figure 3.
In(a)afictitiousmassof100a.u. hasbeenusedandthermal
V. THE TIME SCALE OF ORBITAL MOTION
equilibration between the subsystems is observed. In (b) a
fictitiousmassof900a.u. isusedandthesubsystemsremain
out of thermal equilibrium. A common view of CPMD is that, if a large enough
gap exists between the energies of the occupied and un-
occupied electronic states, the electronic orbitals move
rescaledmassesoftheionsinthesimulationwithµ=100 on time scales that are much faster than typical ionic
a.u. were 15.629a.m.u. and1.988a.m.u. If the rigid-ion time scales. For example, it has been suggestedthat the
approximationwas perfectly applicable the fact that the motion of the orbitals in CPMD can be approximately
ionicmasseshavebeenrescaledaprioriwouldmeanthat described as a superposition of oscillations whose fre-
thetwosimulationswithdifferentvaluesofµshouldhave quenciesaregivenbyω =(2(ǫj−ǫi))1/2,whereǫ andǫ
ij µ i j
identical dynamic and thermodynamic behaviour. are the Kohn-Sham energy eigenvalues of occupied and
It was found (figure 4) that after more than a further unoccupied electronic states, respectively[5, 9, 11]. In a
2.5 ps the simulation with µ = 900 a.u. still showed no system with a substantial energy gap between occupied
sign of thermalization according to the mass-corrected andunoccupiedstates,allthesefrequenciesareveryhigh
definition of temperature. However, the simulation with compared to typical ionic frequencies. In this view of
µ = 100 a.u. equilibrated quickly and the subsystems CPMD the lowest vibrational frequency present in the
wereatthe sametemperatureafter 1.5−2ps (figure 4). orbitalsubsystemisdeterminedbythesizeoftheenergy
Inordertocomparethephononspectraoficeusingthese gap and the fictitious mass. However, it should be clear
two different values of µ for reasonably well equilibrated fromtheworkofPastoreetal. andsections IIand IIIof
simulations at the same temperature it was decided to the present work that, if the orbitals are to remain close
continue both simulations from the end of this initial 2.4 to the ground state and the ions are moving, their mo-
ps run with µ = 100 a.u. For both µ = 100 a.u. and tion must contain an ionic time scale component. Here
µ= 900 a.u., a further 5.5 ps of simulation were carried itisshownthat,infact,itmakesthedominantcontribu-
outduringwhichtimetheoxygenanddeuteriumsubsys- tion to the total orbital kinetic energy. The slow orbital
temsremainedatthesame(mass-corrected)temperature motion results from the coupling between ions and the
in both simulations. For both values of µ, the velocity orbitals and it is precisely this coupling that leads to a
autocorrelationfunctionwascalculatedonatimedomain systematicdepartureoftheaverageforcesinCPMDfrom
of 1.45 ps by averaging over this final 5.5 ps of simula- the ground state forces.
tion. These velocity autocorrelation functions were then Iftherigidionapproximationholdsperfectly,thenthe
fourier transformedto obtain the phononpower spectra. partofthekineticenergyoftheorbitalsduetotheevolu-
The reason for the inability of the µ = 900 a.u. sim- tionoftheelectronicgroundstatecanbeobtainedsimply
ulations to thermalize during the first 7 ps of simulation from the mass corrections and the velocities of the ions,
remains unclear. By using the dressed ionic masses in i.e.
thecalculationoftemperatureweareimplicitlyassuming
1 1
that orbitals move rigidly with the ions and are unper- µhψ˙ |ψ˙ i= ∆M R˙αR˙α+ ∆M R˙αR˙α
turbedbytheirenvironment. Inacondensedsystemthis i i 2 O I I 2 D I I
i IǫO IǫD
X X X
10
tion of the electrons and an accurate calculation of the
physical property of interest. It is important to know
that a calculatedphysicalproperty that agreeswell with
experimentdoessobecausethetrajectoryisrealisticdue
to good forces calculated from a good treatment of the
electronic structure. When there is a microscopic ba-
sis for agreement with experiment the method can be
applied with some confidence to situations or to prop-
erties for which experimental data is unavailable. This
point is stressed because agreement with experiment in
molecular dynamics simulations can often occur due to
a cancellation of errors. This has been observed both
in first principles simulations[5] and in simulations us-
FIG. 5: The fictitious kinetic energy (FKE) of the electronic ingempiricalpotentialswhere potentialsthatagreewith
orbitalsasafunctionoftimecomparedtotheFKEpredicted
experiment on many properties[36] have been found to
byassumingthatorbitalsareinertandrigidlyfollowtheions.
provide a very poor microscopic description of the inter-
The close agreement between the curves verifies that orbital
atomic forces[32].
motion occurs primarily on ionic time scales and that the
It is obvious, therefore, that in order to proceed with-
contributionoftheultra-highfrequencyoscillationsoftheor-
out recourse to empiricism or a posteriori experimental
bitals about their instantaneous equilibria to the total FKE
is tiny by comparison. verificationoneshouldbeabletodependontheaccuracy
of the computed forces. What is much less obvious is
howoneshouldjudge the accuracyofforces. Onewayof
(21) quantifying the average magnitude of the departure of a
setofforces{F }fromsomereferenceforces{Fref}isby
In figure 5 the fictitious kinetic energy (FKE) as pre- I I
computing the percentage root-mean-squared difference
dictedbytherigidionmodelinequation 21iscompared
between the two sets of forces relative to the root-mean-
totheFKEthatisextractedfromtheCPMDsimulation.
squared force component, i.e. :
The difference between them is also plotted and this is
made up both of ultra high frequency oscillations of the
kF −Frefk2
orbitalsabouttheirinstantaneousequilibriaandthepart I,ν I I
∆ =100× (22)
of the ionic time scale evolutionof the groundstate that rms qP kFrefk2
I,ν I
cannot be described within the rigid ion approximation.
q
The very close agreement between the predicted FKE where the sum Pis over a large number of ions I
I,ν
andthe actualFKE demonstrates clearlythat electronic and over a large number of points ν along the molecu-
orbitalsmoveonionictimescalesand,furthermore,that lar dynamics traPjectory. In general, this is an extremely
this slow motion accounts for almost all of their ficti- crude measure of the departure of the forces from the
tious kinetic energy. This also demonstrates that the reference forces, however there is frequently no alterna-
lowest vibrational frequencies of the orbital subsystem tive to using it. A number of classical potentials have
are almost[34] independent of µ and are simply given beenparametrizedbyminimizing this quantitywhile us-
by the lowestvibrationalfrequencies of the ionic subsys- ing ground state density functional theory forces as the
tem. Some previous simulations have calculated the or- reference. Ithasbeenfoundthat,forsomesimpleoxides,
bitalvibrationalpowerspectrumbyfouriertransforming values as low as ∆ = 5−20% can be achieved and
rms
theorbitalvelocityautocorrelationfunctionandhavenot that, as long as ∆ is sufficiently small (i.e. < 20% ),
rms
detected such low frequencies[5, 12]. However, in these the accuracy of these force fields for many experimental
simulationstheorbitalvibrationswereanalysedwhilethe quantities is well correlated with its value[32, 37]. Al-
ions were kept stationary and therefore the contribution thoughitiscrudeandunsatisfactory,intheabsenceofan
totheorbitalmotionfromtheevolutionoftheelectronic alternativequantitativegeneralmeasureofthequalityof
ground state was not present. In reference 35 Herbert forces in the presence of errors of unknown consequence,
and Head-Gordon have analysed the orbital vibrations this quantity is used in the present work.
during ion dynamics, and figure 7 of that paper clearly In this section the forces on the ions in the CPMD
shows a contribution to the orbital frequency spectrum simulations of ice are inspected along a segment of the
thathasalow(ionic)frequencyandthatappearsalmost trajectory. After approximately 1 ps of simulation the
independent of the fictitious mass. groundstateforceswerecalculatedalonga21fssegment
ofthe CPMD trajectory. These groundstate forceswere
thenusedasthe referenceforcesforcalculatingtheroot-
VI. FORCES mean-squared (r.m.s) relative error in the CPMD forces
according to equation 22 both before and after they had
The goal of first principles molecular dynamics is to been corrected using the rigid-ion mass correction. In
make a direct connection between an accurate descrip- other words (using the notation introduced in section