Car-Parrinello molecular dynamics with Vanderbilt ultrasoft
pseudopotentials PRB 1993 47 10142-10153
Some derivation and comments are described in the footnote by Masahiro Yamamoto

Kari Laasonen (IRRMA), Alfredo Pasquarello (AT& T Bell Laboratories)
Roberto Car (IRRMA), Changyol Lee (Harvard University),David Vanderbilt (Rutgers University)
Modiﬁed on January 18, 2006
We show how the ultrasoft pseudopotentials which have recently been proposed by Vanderbilt can be implemented
eﬃciently in the context of Car-Parrinello molecular-dynamics simulations. We address the diﬀerences with respect tothe conventional norm-conserving schemes, identify certain problems which arise, and indicate how these problems can beovercome. This new scheme extends the possibility of performing ﬁrst-principles molecular dynamics to systems includingﬁrst-row elements and transition metals. URL: http://link.aps.org/abstract/PRB/v47/p10142

1

INTRODUCTION

The treatment of the electronic structure causes a considerable increase in the computational eﬀort, such that thesize of the system which can be aﬀorded is generally muchsmaller than for classical simulations. The number of planewaves depends on the size of the system and on the energy
cutoﬀ required for a sulﬁciently accurate description of theelectronic structure. This energy cutoﬀ is a property of thePP, and can sometimes aﬀect the feasibility of a CP simulation.
In the original CP version norm-conserving PP’s (Refs.
6 and 7) have been used in their fully separable form[8].
In such a PP scheme, the pseudo-wave-function matches
the all-electron wave function beyond a cutoﬀ radius which
deﬁnes the core region. Within the core region, the pseudowave-function has no nodes and is related to the all-electron
wave function by the norm- conserving condition which ensures that both wave func- tions carry the same charge. In
the spirit of reduc- ing the energy cutoﬀ, several improvements have been proposed.[9-12] However, despite these improvements, the energy cutoﬀ needed to describe the localized valence orbitals of frst-row elements or transitions metals is still frequently too high to allow MD simulations ofextended systems.
Recently, Vanderbilt has proposed a new PP scheme inwhich the norm-conserving condition has been relaxed[13].
In this scheme, the pseudo-wave-functions are allowed to
be as soft as possible within the core region, yield- ing
a dramatic reduction of the cutoﬀ energy required to describe them. Technically, this is accomplished by introducing a generalized orthonormality condition, which modiﬁes
the conventional approach signiﬁcantly In order to recover
the full electronic charge, the electron density given by thesquared moduli of the wave functions is augmented in thecore regions. Thus, the electron density can be subdivided
in a smooth part extending throughout the unit cell, anda hard part localized in the core regions. Note that theaugmented part appears only in the electron density; this
distinguishes the current scheme from others, such as thelinearized augmented plane wave (LAPW), in which similarideas have been applied to the wave functions.
In this paper we address the consequences of the gen-Molecular-dynamics (MD) simulations have been widelyused in many diﬀerent ﬁelds, which range from physical chemistry to solid-state physics, to study the microscopic behavior of temperature- and time-dependent phenomena[1]. A great deal of work has been carried out usingclassical potentials to account for the interactions between
atoms. The main drawback of these interaction potentials
is that they are usually derived semiempirically from particular bonding situations, and can fail to give an appropriate description when the chemical environment is changed.
Because of this transferability problem, it is sometimes desirable to have a MD scheme in which the variations of theelectronic structure are accounted for during the simulation.
In an important paper, Car and Parrinello[2] proposed amethod to perform molecular dynamics in which the electronic structure is described in the density-functional localdensity approximation (LDA).[3,4] In this scheme, the ionicforces are determined directly from the electronic structureof the system independently of any empirical parameter,
and are therefore highly accurate over a wide range of bonding situations. The Car-Parrinello (CP) method has beensuccessfully applied to a large variety of systems providing detailed information on electronic as well as structuralproperties.
In spite of the generality of the underlying idea, the CarParrinello method has primarily been applied in its original version, i.e., using plane-wave basis sets with periodic
boundary conditions in conjunction with pseudopotentials(PP’s). There are many advantages of such a basis set. The
mathematical formulation is particularly simple. The basis
set is independent of the ionic positions, giving an unbiased
uniform description of the simulation cell and preventing
undesirable Pulay terms[5] from appearing in the calculation of the ionic forces. Plane waves easily allow the use of
fast Fourier transforms (FFT’S) to transfer quantities fromreal space to Fourier space and vice versa. Another advantage is the possibility of testing the accuracy of the results
by increasing the energy cutoﬀ, which deﬁnes the highest
kinetic energy of the plane waves in the basis set.
1

byeralized orthonormality condition in the context of
Car-Parrinello simulations. The calculations are aﬀected
in several ways. First, a new term appears in the
Kohn- Sham equations which is dependent on the
wave functions and must thus be updated at every time step. Second, the orthonormality condition depends upon the ionic positions. As a consequence, the manner in which this condition is imposed during the ionic motion and the expressions
for the ionic forces are substantiahy modiﬁed withrespect to the norm-conserving case. Finally, thehard contribution to the electron density must be
accounted for without losing the advantage of the
low cutoﬀ energy required for the wave functions.
In some cases this problem can be ameliorated bya careful generation of a pseudodensity to represent
the hard contribution. Otherwise it can be solved
using a real-space approach which takes advantageof the fact that the augmented parts of the electron
density are localized. We will show that, in spite of
the more complex formulation, the present scheme
is well suited to handle extended systems containingﬁrst-row elements or transition metals.
The paper is organized as follows. In Sec. II we recall the general properties of the Vanderbilt ultrasoft PP
scheme. In Sec. III we implement the PP in a Car- Parrinello molecular-dynamics scheme, addressing various consequences of the generalized orthonormality condition. Section IV explains how the implementation is divided between
real and reciprocal space. We give a brief summary in Sec.
V of various applications that the method already has allowed, and conclude in Sec. VI with a discussion of perspectives and future directions. The Appendix contains analternative derivation of the constraint contribution to theionic forces. Although the method has never been describedin full, various accounts of speciﬁc applications have beenpresented elsewhere. [14-16]

2

2.1

VNL =

X

(0)
I
Dnm
|βnI 〉〈βm|(0)where the functions βnI as well as the coeﬂicients Dnm characterize the PP and diﬀer for diﬀerent atomic species. In
the following we will for simplicity consider only one atomic
species . The βnI functions are centered on site I, and thus
depend on the ionic positions throughβnI (r) = βn (r − RI )

i

nm,I

where the augmentation functions QInm (r) = Qnm (r −RI ) are also provided by the PP and are strictly
localized in the core regions. Thus, while the electron
density in Eq. (4) is still quadratic in the wave functions, it
is now separated into a soft delocalized contribution givenby the squared moduli of the wave functions, and a new
hard contribution localized at the cores. The ultrasoft PP(0)ion
is fully determined by the quantities Vloc(r), Dnm , Qnm (r),and βn (r). The algorithm used to generate these quantities
is described in Ref. 13 and is brieﬂy reviewed in Sec. II B.
The relaxation of the norm-conserving condition isachieved by introducing a generalized orthonormality condition〈φi |S({RI })|φj 〉 = δij(5)R
where qn m =drQnm (r). The orthonormality condition(5)
isconsistent
with the conservation of the chargeRdrn(r) = Nv . Note that the overlap operator S is dependent on the ionic positions through the |βnI 〉, eq(3).
The ground-state orbitals φi are those which minimize
the total energy (1) under condition (5), 1

Z Zn(r)n(r′ )drdr′|r − r′ |Z
ion+Exc [n] + drVloc(r)n(r)+U ({RI })

(6)

nm,I

In Vanderbilt’s ultrasoft PP scheme,[13] the total energy of
Nv valence electrons, described by the wave functions φi , isgiven byX
Etot [{φi }, {Ri }] =〈φi | − ∇2 + VNL |φi 〉12

I|
qnm |βnI 〉〈βmS =1+

Kohn Sham equations+(3)

Here βn is an angular momentum eigenfunction in the angular variables, times a radial function which vanishes outside
the core region; the indices n and m in Eq. (2) run over
the total number Nβ of such functions. In the ultrasoft PP
case, often two reference energies, and therfore two radialfunctions, are required for each included angular momentum channel. This leads to a number Nβ which is generally
twice as large as for a corresponding Kleinman-Bylander[8]norm-conserving PP.
The electron density in Eq. (1) is given by


X
X
I
|φi (r)|2 +QInm (r)〈φi |βnI 〉〈βm|φi 〉 (4)n(r) =VANDERBILT’SULTRASOFTPSEUDOPOTENTIAL
where S is a Hermitian overlap operator given by
SCHEME
X

i

(2)

nm,I

δEtot
= ϵi Sφi (r)δφ∗i (r)(7)

where ϵi have been introduced as Lagrange multipliers. Because of the fact that the augmentation part of the charge(1)where n(r)is the electron density, Exc is the exchange
and correlation energy, and U ({RI }) is the ion-ion interion
action energy. The PP contains a local part Vloc
(r) =
P
ion
V
(|r
−
R
|)
(Ref.
17)
and
a
fully
nonlocal
part
given
I
I loc

1

δ(Etot − ϵi 〈φi |S({RI })|φj 〉)=0δφ∗i (r)2part of the PP, Eq.(2), by deﬁning new coeﬃcients
Z
I
(0)
Dnm
= Dnm
+ drVeﬀ (r)QInm (r)density depends on the wave functions,Xδn(r′ )
I
= φi (r′ )δ(r′ − r) +QInm (r′ )βnI (r)〈βm|φi 〉 (8)∗δφi (r)

(13)

nm,I

(0)

Note, however, that the Dnm are just parameters whichIcharacterize the PP, whereas the new Dnm
depend on the
wave functions through Veﬀ , Eq.(12), and have to be updated in the self-consistency procedure. (The treatment
given here of the screening of the D matrices is consis- tentwith that of Ref. 14, but diﬀers from, and should be considered to supersede,[18] that of Ref. 13.)
At this stage, the diﬀerence with respect to the normconserving case resides in the presence of the S operator,
the wave-function dependence of the Dnm and the fact thatthe number N of β functions is twice as large. The calculaItion of the Dnm, Eq. (13), can be carried out in real spaceand produces only a modest overhead (see Sec. IV B). The
presence of the S operator requires the calculation of the〈βnI |φi 〉, which, however, are also needed for the nonlocal
PP, and thus do not require any additional computation.
The number of operations needed to calculate one scalar
product of this type scales like the number of plane waves
Npw , and the number of these scalar products is given by
Nat Nβ Nband where Nrmat is the number of atoms, Nβ is thenumber of βn functions per atom, and Nband is the number
of ﬁlled states. Since Npw and Nband both scale like Nat ,3
this part scales like Nat
and for large systems it is typicahy the most time consuming part. (When also the ionicforces [see Eq. (43)] are calculated, similar scalar products〈dβnI /dRI |φi 〉 are required, just as in conventional normconserving schemes.) Thus , for large systems in which
the computational cost related to these scalar products isdominant, one can deduce that the ultrasoft PP scheme
becomes advantageous com- pared to the norm-conserving
NC(NC) scheme when Npw> 2Npw , where we have takenNC
Nβ = 2Nβ . In terms of energy cutoﬀ Ecwf , the criterion
is Ecwf,N C > 1.6Ecwf , which is easily satisﬁed for ﬁrst-row
elements and transition metals.
Very recently King-Smith, Payne, and Lin[23] have shownthat it is possible to evaluate the scalar products between
the wave functions and the βnl functions in real space by taking advantage of the fact that the βnI are localized. In thisway, this part of the calculation scales like Nat Nband , i.e.,
Iike Nt2 . Although in principle the full calculation would
3
still scale like Nat
because of the orthonormalization proce-

additional terms appear in the Kohn-Sham equations fromthe density-dependent terms in the total energy (1). As an
example, we consider the exchange and correlation energy.
Using Eq.(8) we obtain
ZδExc [n]δExc [n] δn(r′ )=dr′∗δφi (r)δn(r′ ) δφ∗i (r)= µxc (r)φi (r)
Z
X
I
I
+
βn (r)〈βm |φi 〉 dr′ µxc (r′ )QInm (r′ )
nm,I

(9)
where µxc (r) = δExc /δn(r). The other terms can be calculated similarly. We obtain
H|φi 〉 = ϵi S|φi 〉where
H = −∇2 + Veﬀ +

X

(10)

I
I
|
|βnI 〉〈βm
Dnm

(11)

nm,I

Here Veﬀ is the screened eﬀective local potential,ZδEtotn(r′ )Veﬀ == Vloc (r) + dr′+ µxc (r)δn(r)|r − r′ |(12)2

All the terms arising from the augnented part of the electron density have been grouped together with the nonlocal
2

Etot

=

X

〈φi | − ∇2 + VNL |φi 〉

i

+

1
2

Z Z
dr′ dr′′Z

ion
drVloc(r)n(r) + U ({RI })+Exc [n] +δEtot
δφ∗i (r)n(r′ )n(r′′ )|r′ − r′′ |=(−∇2 + VNL )|φi 〉

+

1
2

+

1
2

+

δExc [n]+δφ∗i (r)=(−∇2 + VNL )|φi 〉

+

1
2

+

1
2

+

1
2

+

Z Z
dr′ dr′′δn(r′ ) n(r′′ )δφ∗i (r) |r′ − r′′ |dr′ dr′′n(r′ ) δn(r′′ )|r′ − r′′ | δφ∗i (r)

Z Z

Z

dr′′

Z Z
Z Z

Z

ion ′
dr′ Vloc(r )δn(r′ )
δφ∗i (r)n(r′′ )1
φi (r) +|r − r′′ |
2

dr′

ion ′
dr′ (VH (r′ ) + µxc (r′ ) + Vloc
(r ))

X

n(r′ )φi (r)|r′ − r|nm,I

VH (r)=n(r′ ) X I
I
I
|φi 〉
(r)〈βmdr′ dr′′ ′Qnm (r′′ )βn|r − r′′ |δEtot
δφ∗i (r)

=

dr

nm,I

X

nm,I

Z

ion
Vloc
(r)φi (r) +ZII(r)〈βm
|φi 〉
βn

ion ′
dr′ Vloc
(r )

X

I
I
|φi 〉
(r)〈βmQInm (r′ )βnZn(r′′ ) X I
I
I
|φi 〉
(r)〈βmdr′ dr′′ ′Qnm (r′ )βn|r − r′′ |µxc (r)φi (r) +ion(−∇2 + VNL + VH (r) + µxc (r) + Vloc(r))|φi 〉

Z

+

Z

nm,I

+

=

+
dr′ µxc (r′ )QInm (r′ )n(r)|r − r′ |(−∇2 + Veﬀ )|φi 〉

Xµ

(0)
Dnm +

nm,I
I
I
|φi 〉
〉〈βm
|βn

I
I
QInm (r′ )βn(r)〈βm
|φi 〉

nm,I

3

¶

Z

dr Veﬀ (r )QInm (r′ )′′dure, the cost of a previously dominant part of the compu- which are dual to the |φn 〉 and are also local; they form thetation would be considerably reduced, allowing the study of ’projectors’ of the nonlocal operators. 4 It follows from Eq.
(6) that
still larger systems .
qnm = 〈ψn |ψm 〉R − 〈φn |φm 〉R
(20)

2.2

Pseudopotential generation algorithm

5
It is straightforward to verify that the |φn 〉 obey the secularequationIn this section we give a concise description of the speciﬁc
PP generation algorithITL introduced in Ref. 13.
!
ÃXAs in other PP methods, an all-electron (AE) calculaDnm |βn 〉〈βm | |φn 〉 =T + Vloc +tion is ﬁrst carried out on a free atom in some referencenmÃ!
conﬁguration, Ieading to a screened potential VAE (r). ThenX
for each angular momentum l , a set of reference en- ergies
ϵn 1 +
qnm |βn 〉〈βm | |φn 〉(21)
ϵlτ is chosen (τ = 1, Nτ ,, typically 1 ≤ Nτ ≤ 3 ) to covernm
the energy range over which good scattering properties aredesired (usually the range of occupied bulk valence states). where
Dnm = Bnm + ϵm qnm(22)
At each ϵlτ the solution of the Schrödinger equation whichis regular at the origin is obtained (for ﬁxed VAE ),
(0)
6
ion
Finally the quantities Vloc
and Dnm are obtained from a[T + VAE (r)]ψn (r) = ϵn ψn (r)(14) ”descreening” procedure,

Here n is a composite index, n = Imτ , and T is the kineticenergy operator (−1/2)∇2 . (The amplitude of ψn is assumed to have been ﬁxed in some deﬁnite way, e.g., by its
value at an arbitrary radius R.) Despite the fact that ψn isin general non-normalizable, a bra-ket notation(T + VAE − ϵn )|ψn 〉 = 0

(0)
Dnm

4 From the deﬁnition of |β 〉iXBik |βi 〉ji

Bik |χj 〉

i,j

|χj 〉

X

δjk |χj 〉 = |χk 〉X=〈φi |χj 〉 |βi 〉| {z }i〈φk |χj 〉X

Bij〈φi |χj 〉〈φk |βi 〉

=

i

then
5q

nm =

R

〈φk |βi 〉=δki

drQnm and Qnm are localized in the core region.

6

(16)

VNL |φi 〉=X(Bkl + ϵl qkl )|βk 〉 〈βl |φi 〉| {z }

kl

X

=

P |

(17)

+

δli

Bki |βk 〉

{z

k

=

X

kn

ϵl qkl |βk 〉〈βl |φi 〉|χi 〉 + ϵi

X

|{z}
δli

=

(19)

qkl |βk 〉〈βl |φi 〉kl(ϵi − T − Vloc )|φi 〉 + ϵiX

qkl |βk 〉〈βl |φi 〉

kl

"
=

}

Bki B −1 |χn 〉=|χi 〉

kl

3

=

¢

j

m

VNL |φi 〉B −1=

is well deﬁned.
(0)ion
We now form the quantities Vloc
, Dnm , Qnm (r), andβn (r) needed to specify the PP. Qnm (r), and βn (r) are givenbyQnm (r) = ψn∗ (r)ψm (r) − φ∗n (r)φm (r)
(18)
X¡
¢
|βn 〉 =
B −1 mn |χm 〉

X¡

=

i

kn

and

(24)

This PP has been generated in such a way that the following features are reproduced in the reference con-ﬁguration.

(15)

are formed. These are local (they vanish at and beyond R,
where Vloc = VAE and φn = ψn ). 3
Thus the matrix of inner products
Bnm = 〈φn |χm 〉

(23)

7

is adopted as a stand-in for the previous equation. Quantities such as 〈ψn |ψn 〉 are ill-deﬁned, but the special notation〈ψn |ψm 〉R will be used to denote the integral ofψn∗ (r)ψm (r)inside the sphere of radius R.
Next, cutoﬀ radii rcl are chosen, and for each ψn obtained
above, a pseudo wavefunction φn is constructed, subject
only to the constraint that it join smoothly to ψn at rcl .
No norm-conservation constraint is explicitly imposed. A
smooth local potential Vloc (r) is also generated in such a
way that it matches smoothly to VAE (r) at a cutoﬀ radiusrcloc , and a diagnostic radius R is chosen such that R isslightly larger than the maximum of the rcl and rcloc . (Rmust be outside the rc ’s in order that scattering properties
calculated there will be meaningful.) Then the orbitals|χn 〉 = (ϵn − T − Vloc )|φn 〉= Vloc (r) − VH (r) − µxc (r)Z= Dnm − dr′ Vloc (r′ )Qnm (r′ )

ion
Vloc

ϵi (1 +

X

qkl |βk 〉〈βl |) − T − Vloc |φi 〉kl|χi 〉 = (ϵi − T − Vloc )|φi 〉=then φi is an eigenvector of T + Vloc + VNL . Then it is shown thatthe fully nonlocal KB-type PP can be obtained directly from VNL =|χi 〉〈χi |/〈χi |φi 〉.

(ϵi S − T − Vloc )|φi 〉7 n(r′ ) in the paper may be replaced by Q′
nm (r )

4#

(Here we assume that the eigenvalues of the occupied valence orbitals in the reference conﬁguration are included
among the chosen ϵτ l ; this is our standard practice.) (i)The pseudoeigenvalues are equal to the AE ones, and thecorresponding orbitals agree exactly out-side rcl . (ii) The
scattering properties are correct at each ϵτ l in the sense thatthe logarithnric derivative and its en-ergy derivative match
the AE one at that energy[13]. Thus, the transferability, as
measured by scattering properties in the reference conﬁguration, can be systematically im-proved by increasing thenumber of such energies. (iii) The valence charge density is
precisely equal to the AE valence density in the reference
conﬁguration (except in-sofar as it is aﬀected by replacement of the Qnm by pseudo densities as discussed in Sec. II
C). This helps improve the transferability with respect tochanges in the screen-ing environment.
We illustrate the result of such a procedure for the caseof the 3d orbital of Cu. In Fig. 1 the all-electron wave function and the pseudo-wave-function obtained with a cutoﬀ
radius of rcd = 2.0 a.u. are compared. The main advantageof the present scheme over previous ones is that no normconservation constraint is imposed during the construction
ψn → φn of the pseudowavefunction. Thus, the construction can be done in such a way as to make the ψn (r) as
smooth as possible. This is the reason we refer to thisscheme as an ”ultrasoft” PP scheme.
FIG, 1. All-electron (solid) and pseudo (dashed) radial
wave functions of the 3d orbital of Cu. A cutoﬀ radius of 2a.u. has been used.

by construction[13]. The number of possible L channels in
Eq. (26) is ﬂnite because of the fact that a nonlocal PP isrequired only for the lowest angular momentum channels.
The Qradnm in Eq. (26) are then replaced by L-dependent
counterparts QL
nm ,X
L
Qnm (r) =cnm(27)
LM YLM (r̂)Qnm (r)LMwhich satisfy the condition that for each L-component theLth moment of the electron charge density be conserved,
L
Z rin

2

r drr
0

L

QL
nm (r) =LZ rin

r2 drrL Qrad
nm (r)(28)

0

L
where L-dependent inner cutoﬀ radii rin
have been introduced. Since the high Fourier components of Qnm aremainly related to the high-L components, it is convenient,Lfor low-L comfor a given cutoﬀ Ecdens , to use smaller rin
ponents. In this way, a relatively better description of thelowest moments of the electron charge density is maintained.
In Fig. 2 we show the pseudo QL
nm for L = 0, L = 2, andL = 4 obtained from Qrad
nm for the case in which n and mcorrespond to the same reference energy in the d channel for
Cu (corresponding to the wave functions in Fig. 1). The
core cutoﬀ radius in this case is 2 a.u., whereas the innerL
cutoﬀ radii rin
range from 0.6 a.u. (L = 0) to 1.2 a.u.
L
(L = 4). With this choice of rin
it is possible to give a goodLdescription of all the Qnm using an energy cutoﬀ of about
200 Ry, as can be inferred from Fig. 3 where we show the
Fourier transforms of the functions obtained in Fig. 2.
In order to construct optimal pseudo QL (r) (where we
2.3 Construction of the pseudoelectron
have dropped the indices n and m), we expand them incharge densityLpolynomials of r inside rinas
In norm-conserving PP schemes the electron density is deL
QL (r) = rL ρL (r), r < rin(29)
ﬁned as in Eq. (4) where only the ﬁrst term is kept on thedens
right-hand side. Thus, the energy cutoﬀ Ec
re-quired todescribe fully the electron density is four times the energy whereρL (r) = C1 + C2 r2 + C3 r4 + ...
(30)
cutoﬀ Ecwf of the wave functions
suﬃcient smoothness
Ecdens = 4Ecwf(25) with the number of terms ensuringof the polynomial. We want QL (r), i.e., ρL (r) as smoothThis direct relationship between the cutoﬀ for electron den- as possible. Therefore, following lines similar to Rappe etsity and wave functions does not hold in the ultrasoft PP al.,[11] we insist that the Fourier coeﬃcients above a certain
scheme because of the presence of the augmentation func- cutoﬀ wave vector Gc should be as small as possible. Thus,
tions Qnm in the electron density (4). In this scheme, it we minimize
Z ∞is therefore appropriate to introduce two independent enG2 Q2L (G)dG(31)
I=ergy cutoﬀs: one for the soft part of the electron den-sity,
Gc
Ecsof t = 4Ecwf , and a second and generally much higher
FIG. 2. L-dependent Q rad (r) for the case of the 3d
one Ecdens to describe the augmentation functions Qnm .
orbitalof Cu. The original Q rad is given by the solid curve.
It is often possible to reduce the charge cutoﬀ Ecdens by
The pseudo QL for L = 0, 2, 4 have been obtained with inner
replacing the functions Qnm by pseudocounterparts [13,16].
cutoﬀ radii rinL of 0.6, 0.8, and 1.2 a.u., respectively.
In this construction, the Qnm are modiﬁed within an inFIG. 3. Fourier transforms of the QL given in Fig. 2.
ner core region (deﬁned by rin ). The charge density dewhereZ ∞scribed by the pseudo Qnm preserves all the charge moQ(G)=
drr2 QL (r)jL (Gr)(32)
L
ments, so that the electrostatic potential beyond rin remains
0
unchanged. This is achieved by decomposing the Qnm acand jL is the spherical Bessel function of order L. The
cording to angular momentum L as
minimization should be done subject to the constraint (27)X
radQnm (r) =cnm(26) and subject to the following continuity requirements:
LM YLM (r̂)Qnm (r)
LM

L
L
ρ(rin
) = ρAE (rin
)
′ L
′
L
ρ (rin ) = ρAE (rin )LLρ′′ (rin) = ρ′′AE (rin
)

cnm
LM

where
are Clebsch-Gordan coeﬁicients, YLM (r̂) are
spherical harmonics, and Qrad
nm (r) gives the all-electron radial dependence of Qnm and is independent of L and M
5(33)

the Euler equations of motion, 8

where the primes indicate radial derivatives. This treatment gives us a very smooth charge density with smoothand continuous ﬁrst and second derivatives, features which
are especially important in the context of gradient-corrected
LDA schemes[19,20] In Fig. 4 we demonstrate the successof this approach for the case of oxygen, speciﬁ-cally for QL
nmwith L = 0 and n, m = 2p. Figure 4 shows how the Fouriercoeﬃcients of the pseudo QL
nm behave as a function of G
when optimally generated according to the above scheme, inLcomparison to the original QL
nm and pseudo Qnm generated
in a fourth-order polynomial just to give the correct charge
and dipole moment and the continuous ﬁrst derivative of
the charge. It is clear that the reciprocal-space convergenceis much improved using the present approach.
We have found that in many cases, including that of oxygen, replacement of the Qnm by pseudodensities using theabove two reﬁnements (use of L-dependent rin and/or optimally soft construction of the QL ) allows us to obtain an excellent description of the charge density with Ecdens = Ecsof t .
In these cases, condition (25) is restored, and the solutionof the Kohn-Sham equations (10) can be obtained in theusual way by fast Fourier transforming between a single
pair of real-space and reciprocal-space meshes, just as isdone with conventional norm-conserving PP’s. However, insome cases, and in particular for transition metals such as
Cu, even the pseudo Qnm require Ecdens > 4Ec wf . In thesecases, the problem of having a large number of plane waves
can still be circumvented using the fact that the Qnm are
localized in real space. The details of our method for accomplishing this will be discussed later in Sec. IV.
FIG. 4. Fourier components of optimal pseudo (solid) ,for nonoptimal pseudo (dashed), and original (dotted) QLnm
oxygen with L = 0 and n, m = 2p. Vertical scale is arbitrary.

δEtot X
Λij Sφj
+
δφ∗i
j

µφ̈i

= −

FI

= MI R̈I = −

∂Etot
+
∂RI
X
∂S
Λij 〈φi |
|φj 〉
∂R
I
ij

(36)

(37)

During a MD run, the Lagrange multipliers are to be determined dynamically using the method of Ryckaert et al.,[21]as discussed in Sec. III B. For the special case of equilibrium, the vanishing of Eq. (36) reduces (with Λij = ϵi δij )tothe electronic Kohn-Sham equations (7) or (10), and thevanishing of Eq. (37) corresponds to the vanishing of theionic forces FI . An alternative deriva-tion of the force expression in Eq. (37), particularly applicable to direct minimization approaches, is given in the Appendix.
We now derive explicit expressions for the two terms on
the right-hand side of Eq. (37) , which correspond to thecontributions from the change in the Hamiltonian and fromthe change in the orthonormality constraints, respectively
The latter contribution appears because of the Rl dependence of the overlap operator S; note that it is absent in
the case of norm-conserving schemes, in which the orthonormality condition does not depend in any way on the ionicpositions. The ﬁrst term on the right-hand side of Eq. (37)must be obtained keeping in mind that the electron densitylalso depends on RI through the QL
nm and βn . Introducing
the quantities
XI〈φi |βnI 〉〈βm
|φi 〉
(38)
ρInm =
i

I
ωnm

=

X

I
Λij 〈φj |βnI 〉〈βm
|φi 〉

(39)

ij

3
3.1

MOLECULAR DYNAMICS
Lagrangianforcesformulation

andwith derivatives

ionic

Following the CP approach[2], the electronic wave functions
and the ionic coordinates evolve according to a classical Lagrangian
L

= µ

XZ

dr|φ̇i |2 +

i

X·
∂β I
I
〈φi | n 〉〈βm
=
|φi 〉
∂R
I
i
¸
I
I ∂βm
|φi 〉
+〈φi |βn 〉〈
∂RI
·
X
∂β I
I
|φi 〉
=
Λij 〈φj | n 〉〈βm
∂R
I
i
¸
I
I ∂βm
+〈φj |βn 〉〈
|φi 〉
∂RI

∂ρInm
∂RI

I
∂ωnm
∂RI

1X
MI Ṙ2I
2

(40)

(41)

I

−Etot ({φi , RI })and noting that

(34)

∂n(r)
∂RI

subject to a set of constraints
Nij ({φi , RI }) = 〈φi |S|φj 〉 − δij = 0

=

¸
X·
∂ρI
dQInm (r) I
QInm (r) nm +
ρnm (42)
∂RI
dRI
nm

we arrive at the expression
ZdU
dV ionFI = −− dr loc n(r)
dRI
dRI

(35)

Here µ is a ﬁctitious mass parameter for the electronic degrees of freedom, MI is the mass of the atoms, and Etot
and S are as given in Eqs. (1) and (6), respectively. The
orthonormality constraints (35) are holonomic, and do notcause energy dissipation during the MD run. They may beincorporated by introducing Lagrange multipliers Λij into

8

d
dt

6

³

∂L
∂ q̇i

´

X ∂gk
∂L
λk
=
∂qi
∂qi
K

−

k

Z
−
−

drVeﬀ (r)
X

I
Dnm

nm

X dQI (r)
nm
ρInm
dR
I
nm

with special attention to the treatment of the RI dependence of the orthonormality constraints. For the electronic
wave functions we obtain∂ρInm X∂ω I
qnm nm
+
∂RI
∂RI
nm

(43)
φi (t + ∆t)

=

(∆t)2
2φi (t) − φi (t − ∆t) −
×
µ


X
δE
 tot −
Λij (t + ∆t)S(t)φj (t)
δφ∗i
j

I
where Dnm
and Veﬀ have been deﬁned in Eqs. (13) and(12), respectively. 9 The last term of Eq. (43) corresponds
to the constraint contribution [last term of Eq. (37)]. 10
Note that because the basis set consists of plane waves, the(44)wave functions do not depend on the ionic positions and no
additional Pulay-type corrections are needed.[5]where S(t) means that the operator S is evaluated for ionicpositions RI (t). Similarly for the ionic coordinates,

3.2

Evolutionstraintsoforthonormality

con-

(δt)
2RI (t) − RI (t − ∆t) −
×
MI


X
∂S(t)
 ∂Etot −Λij (t + ∆t)〈φi (t)||φj (t)〉
∂RI
∂RI
ij

RI (t + ∆t) =Here we discuss in some detail the discretization of the equations of motion (36) and (37) using the Verlet algorithm, 11
9

Etot

=

X

〈φi | − ∇2 + VNL |φi 〉 +

Z

i

∂Etot
∂RI

=

−

X

〈φi |

1
−
2
−

=

drdr′

Z

dr

∂RI

Z

dr

ion (r)
∂Vloc

∂RI

{z

nm

I
(0) ∂ρ
Dnm nm
∂RI

nm

Z

∂ρI
dQInm (r) I
QInm (r) nm +
ρnm
∂RI
dRI

∂U ({RI })
−
∂RI

Z

−

drVeﬀ (r)

X ∂ρI

X dQI

·

nm

|

∂RI

φ̄i = 2φi (t) − φi (t − ∆t) −¸Z(∆t)2 δEtot (t)
µ
δφ∗i

(49)

In norm-conserving schemes the identity operator isfound in place of S, which leads to a simpler form of Eq.
The velocities do not appear in the above formula and they have been
eliminated by addition of the equations obtained by Taylor expansion
about ⃗ri (t):¸⃗ri (t + δt)=⃗ri (t − δt)=

drVeﬀ (r)QInm (r)

{z

(48)

12

n(r)

nm (r) I
ρnm
dRI

(0)
Dnm +

nm

∂RI

nm

dr

ion (r)∂Vloc= 〈φ̄i |S(t + ∆t)|φ¯j 〉,= 〈S(t)φi (t)|S(t + ∆t)|φ¯j 〉,= 〈S(t)φi (t)|S(t + ∆t)|S(t)φj 〉

with

}

X·

−

n(r) −

X

Aij
Bij
Cij

ion
dr [VH (r) + µxc + Vloc(r)] ×|(47)

where a dagger indicates the Hermitian conjugate (be-causeof the Hermiticity of S, λ = λ† ) and wheren(r)∂U ({RI })
∂n(r)
ion
drVloc
−
(r)
∂RI
∂RI

=Veff (r)−A + λB + B † λ† + λCλ† = 1,

ion (r)
∂Vloc

(46)

Fulﬁlling this constraint leads to the following matrix equation for the Lagrange multipliers λ = (∆t)2 Λ∗ (t + ∆t)/µ:∂n(r′ ) n(r)∂RI |r − r′ |Z=

〈φi (t + ∆t)|S(t + ∆t)|φj (t + ∆t)〉 = δijn(r′ )
drdr∂RI |r − r′ |∂U ({RI })
−
−
∂RI
−

The orthonormality condition has to be imposed at eachtime step[21]∂ X (0) I ′I′
Dnm |βn 〉〈βm
|φi 〉
∂RI

dExc ∂n(r)−dn ∂RI−(45)n(r)n(r′ )|r − r′ |′ ∂n(r)

Z Z

Z

drdr′

nmI ′

Z Z

i

1
−
2

Z Z

ion
drVloc(r)n(r) + U ({RI })+Exc [n] +
−

1
2

}

(δt)2
(δt)3 ⃗⃗ai (t) +bi (t)...
2
6
2
3
(δt)
(δt) ⃗⃗ri (t) − δt⃗⃗ai (t) −vi (t) +bi (t)...
2
6⃗ri (t) + δt⃗vi (t) +

I
=Dnm

12
10

last term

X

=

ij

X

=

nm

2

X
∂
I′
I′
Λij 〈φi |(1 +qnm |βn〉〈βm|)|φj 〉∂RI

φi (t + ∆t)

nmI ′

δij

∂ω I
qnm nm∂RIX

Λij (t + ∆t)S(t)φj (t)=〈φi (t + ∆t)|S(t + ∆t)|φj (t + ∆t)〉=〈φ¯i +(∆t)2 XΛij (t + ∆t)S(t)φj (t)|S(t + ∆t)|µj(∆t)2 Xφ¯j +Λjk (t + ∆t)S(t)φk (t)〉µtions at t + δt may be given by
=

(∆t)
φ¯i +
µ

j

11 If the classical trajectory is continuous, the estimate of the posi-⃗ri (t + δt)=k2⃗ri (t) − ⃗ri (t − δt) + (δt)2⃗ai (t) + O[(δt)4 ]7

(47) presented in Ref. 22. In the ultrasoft PP case, the solu- 4
tion of Eq. (47) is somewhat problematic because Eq. (45)is not a closed expression for RI (t + ∆t). The problem isthat Λ(t+∆t) appearing in Eq. (45) depends upon S(t+∆t)through Eqs. (47)-(49), and thus upon RI (t + ∆t). Con- 4.1
sequently, in principle it is necessary to solve iteratively for
4.2
RI (t+∆t). We do this by ﬁrst estimating the new Λ(t+∆t)using two previous values,(0)Λij (t + ∆t) = 2Λij (t) − Λij (t − ∆t)

(50)

(0)

In order to solve Eq. (47) , we generalize the iterative
procedure used in Ref. 22 because the unmodiﬁed proceduredoes not always converge. In the norm-conserving case, thematrix B converges to the identity matrix for vanishing δt.
This is not the case in the ultrasoft PP case. However,when the matrix B is decomposed into Hermitian (Bh ) andanti-Hermitian (Ba ) parts,

Fourier interpolation scheme

5

APPLICATIONS

6

CONCLUSIONS

7

APPENDIX:
ALTERNATIVEDERIVATIONOFIONICFORCES

Here we brieﬂy give an alternative derivation of the expression for the ionic forces FI given in Eqs. (37) or (43) of Sec.
111A, and in particular the last term which corresponds tothe contribution from the orthonormality constraints. The
present derivation makes no reference to the Lagrangian formulation, and is therefore more natural in the context of
direct minimization approaches such as conjugate-gradientschemes.
In the main body of this manuscript, we took the pointof view that Eq. (1) speciﬁes the total-energy functional
even when the φi depart from the constraints (5). Here we
introduce a diﬀerent generalized functional(51)

it is straightforward to see that Ba vanishes in the limit ofsmall ∆t. The ﬁrst iteration λ(0) can now be obtained from

λ(0) Bh + Bh λ(0) = 1 − A

Double-grid technique

Because of the reduced number of plane waves required,
the present scheme based on ultrasoft PP’s makes it possible to extend ﬁrst-principles molecular-dynamics simulations to systems containing ﬁrst-row elements and transition
metals, which could not be aﬀorded using norm-conserving
pseudopotentials. It should allow the study of relaxation
processes in transition-metal molecules or defects, oxidation processes, surface recon-structions, and the propertiesof liquid and amorphous materials. The present formulation
also allows the com-bined use of ultrasoft pseudopotentials
for some of the atoms and standard norm-conserving pseudopotentials for the other ones. Such an approach would
be particu-larly useful for the study of transition-metal or
ﬁrst-row defects in semiconductors.

and using this to ﬁnd the new RI (t + ∆t), which is correct to O(∆t4 ). Then Eq. (47) is solved (see below) ina similar way as in the norm-conserving case[22] giving a(1)
new set of Λij (t + ∆t), with which the whole procedure isrepeated, and so on until convergence is achieved. Fortunately, it turns out in practice that the ionic positions arevery well determined by Eq. (50), so that the procedure
almost always converges on the very ﬁrst iteration. Thus,
the operations described above are generally executed only
once per time step.

B = Bh + Ba

RECIPROCAL- AND REALSPACE IMPLEMENTATION(52)

where the C-dependent term has been neglected because ofhigher order in ∆t. Equation (52) can be solved exactly
introducing the unitary matrix U , which diagonalizes Bh ,Ẽtot ({φi }, {RI }) = Etot ({φ̄i ({φi , RI })}, {Rl })(54)
U † Bh U = D, where Dij = di δij . The solution to Eq. (47)can be obtained by iteratingwhere Etot is as given in Eq. (1), and the φ̄i are a new setof S-orthonormal functions
X(A−1/2 )ji |φj 〉(55)|φ̄i 〉 =λ(n+1) Bh +Bh λ(n+1) = 1−A−λ(n) Ba −Ba† λ(n) −λ(n) Cλ(n)j(53)
where at every step the new λ(n+1) are obtained in the same
constructed from the φi via
way as λ(0) had been obtained from Eq. (52).
Aij = 〈φi |S|φj 〉When imposing the orthonormality condition (46), Eqs.
(47)-(49) require the calculation of an additional set of
scalar products of the type 〈βnI |φi 〉 as compared to the normconserving case. However, the additional computational
overhead for these additional products is modest, since theoperations of the form 〈∂βnI /∂RI |φi 〉 are more numerous.

(56)

We imagine φ̄i = φi before the virtual displacement, andcalculate the partial derivative ∂ Ẽtot /∂RI holding φi (butnot φ̄) constant. While the φi do not generally continue
to satisfy the constraints (5) after the virtual displacement,8

8

the φ̄i do by construction, and therefore we can associate
the above derivative with the physical force,
FI = −∂ Ẽtot∂RI

1. Present address: Laboratory of Atomic and Solid State Physics,
Cornell University, Ithaca, NY 14853 2501.
2. Simulations of Liquids and Solids, edited by G. Ciccotti, D.
Frenkel, and I. R. McDonald (North Holland, Amsterdam, 1986).
3. R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
4. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
5. W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
6. P. Pulay, Mol. Phys. 17, 197 (1969) ; C. Satoko, Phys. Rev. B
30, 1754 (1984).
7. D. R. Hamann, M. Schlu”ter and C. Chiang, Phys. Rev. Lett.
43, 1494 (1979).
8. G. B. Bachelet, D. R. Hamann and M. Schlu”ter, Phys. Rev. B
26, 4199 (1982).
9. L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425
(1982).
10. G. Kerker, J. Phys. C 13, L189 (1980) [CAS][dot INSPEC].
11. D. Vanderbilt, Phys. Rev. B 32, 8412 (1985).
12. A. M. Rappe, K. M. Rabe, E. Kaxiras and J. D. Joannopoulos,
Phys. Rev. B 41, 1227 (1990).
13. N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
14. D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
15. K. Laasonen, R. Car, C. Lee and D. Vanderbilt, Phys. Rev. B
43, 6796 (1991).
16. C. Lee, D. Vanderbilt, K. Laasonen, R. Car and M. Parrinello,Phys. Rev. Lett. 69, 462 (1992).
17. A. Pasquarello, K. Laasonen, R. Car, C. Lee and D. Vanderbilt,
Phys. Rev. Lett. 69, 1982 (1992).
18. The interpretation of Vilocon as an atomic or crystalline quantity is to be determined from context.
19. We thank X. Gonze for ﬁrst pointing out that the original
treatment of the D matrices in Ref. 14 was inconsistent except in thefree atom case.
20. A. D. Becke, Phys. Rev. A 38, 3098 (1988).
21. J. P. Perdew, Phys. Rev. B 33, 8822 (1986).
22. J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comput.
Phys. 23, 327 (1977) ; G. Ciccotti and J. P. Ryckaert, Comput. Phys.
Rep. 4, 345 (1986).
23. R. Car and M. Parrinello, Simple Molecular Systems at VeryHigh Density (Plenum, New York, 1989), p. 455.
24. R. D. King Smith, M. C. Payne and J. S. Lin, Phys. Rev. B
44, 13063 (1991).
25. A. Satta, A. Pasquarello, and K. Laasonen (unpublished).
26. K. Laasonen, M. Parrinello, R. Car, C. Lee, and D. Vanderbilt,
Chem. Phys. Lett. (to be published).
27. P. Ballone and G. Galli, Phys. Rev. B 42, 1112 (1990).
28. P. Blo”chl and M. Parrinello, Phys. Rev. B 45, 9413 (1992).
29. S. Nose’, Mol. Phys. 52, 255 (1984) ; W. G. Hoover, Phys.
Rev. A 31, 1695 (1985).
30. R. D. King Smith and D. Vanderbilt, Ferroelectrics 136, 85(1992) .
31. K. W. Kwak, R. D. King Smith, and D. Vanderbilt, in Proceedings of the 7th Trieste Semiconductor Symposium on Wide Band Gap
Semiconductors [Physica B (to be published) ].

(57)

We rewrite Eqs. and (1) and (4) in the form
Etot =

X〈φ̄i |H ion |φ̄i 〉 + EHxc [n]

(58)

i

and
n(r) =

X〈φ̄i |K(r)|φ̄i 〉

(59)

i

Here
H ion

ion
= −∇2 + Vloc
+
¸
Z
X·
(0)
I
ion
I
Dnm + drQnm (r)Vloc (r) |βnI 〉〈βm
|
nm,I

K(r)

X

= |r〉〈r| +(60)I
QInm (r)|βnI 〉〈βm
|

(61)

nm,I

The notation is essentially that of Ref. 14. The minimization of (A5) Ieads to the secular equation
H|φi 〉 =X

Λij S|φi 〉

(62)

j

where Λij = 〈φi |H|φj 〉 and the screened H is given in Eq.
(11).
Using the fact that Aij = δij before the virtual displacement, we ﬁnd∂ φ̄i (r)
∂S
1X
〈φj |
=−
|φi 〉 φj (r)
∂RI
2 j
∂RI

(63)

Now the contributions to ∂Etot /∂RI which do not involve
the RI dependence of the φ̄i areÃ!
ZX∂ ẼtotdH ionδ Ẽtot dn(r)〈φi ||φi 〉 + dr
=
∂RI
dR
δn(r)
dRI
I
i
H
Z
X
dK(r)
dH ion=+ drµHxc
|φi 〉
〈φi |
dR
dRI
I
i
(64)
while Eq.(63) leads to additional terms
Ã!
X∂ Ẽtot
∂S
= −
〈φi |
|φj 〉 〈φj |H ion |φi 〉
∂RI
∂R
I
ij
S
Z
X
∂S
− drµHxc (r)〈φi ||φj 〉〈φj |K(r)|φi 〉
∂RI
ij
= −

X
ij

〈φi |

∂S
|φj 〉 Λji
∂RI

References(65)

After a few lines of algebra, it becomes evident that Eq.(64)
generates the 2nd, 3rd, and 4th terms of Eq.(43), while
Eq.(65) corresponds to the last term of Eqs.(37) or (43).
9