PHYSICAL REVIEW B VOLUME 48, NUMBER 17 1 NOVEMBER 1993-I

Ab initio molecular dynamics for open-shell transition metals

G. Kresse and J. Hafner
Institut fur Theoretische Physik Technische Universitat Wien, Wiedner IIauptstrasse 8-10, A-10/0 Wien, Austria
(Received 2 July 1993)
We show that quantum-mechanical molecular-dynamics simulations in a finite-temperature local-
density approximation based on the calculation of the electronic ground state and of the Hellmann-
Feynman forces after each time step are feasible for liquid noble and transition metals. This is
possible with the use of Vanderbilt-type "ultrasoft" pseudopotentials and eKcient conjugate-gradient
techniques for the determination of the electronic ground state. Results for liquid copper and
vanadium are presented.

The unification of the local-density approximation the Fermi level.

(LDA) to the many-electron problem with classical In the limit of zero temperature the total energy of
molecular dynamics (MD) has marked an important the ground state, E[n(r)), is stationary with respect to
progress in the attempts to put the atomistic simulation variations of the electron density. At finite temperature
of the structural and electronic properties of materials the free energy
on a firm quantum-mechanical basis. The approach orig-
inally proposed by Car and Parrinello is based on the
introduction of a pseudo-Newtonian dynamics describing
the evolution of the electronic degrees of freedom close
E[n(r), f;, p] = E[n(r)] —o S[f;] + p ) f; —W, i (1)
to the adiabatic ground state. An essential condition for
the practicability of the Car-Parrinello method is that the is the proper variational functional which has to be min-
transfer of energy between the atomic and electronic sub- imized with respect to variations of the electron density
systems be small in order to prevent the electrons from n(r), the orbital occupancy f, , and the chemical potential
drifting away &om the ground state. Metals are noto- p. At a given temperature the occupation f, of a state
riously difIicult to handle in the Car-Parrinello method with energy e,. is given by the Fermi-Dirac distribution
because in the absence of a gap the resonance between the and the entropy S[f,] corresponds to noninteracting elec-
atomic and electronic dynamics drives the system rapidly trons (in that case o = Ic~T) Gaussian .broadening of
into nonadiabaticity. For transition metals an additional the one-electron levels can be computationally more con-
difIiculty arises from the fact that even modern soft-core venient than Fermi-Dirac broadening. In this case o is
pseudopotentials require a large energy cutofF for com- the width of the Gaussian; the corresponding form for the
plete total energy convergence and hence a plane-wave entropy has been worked out in Ref. 7. The variational
basis that is too large to allow for a'6 initio MD sim- derivatives of the free energy with respect to the atomic
ulations. The aim of the present paper is to demon- positions are the Hellmann-Feynman forces. ' The ad-
strate that these problems can be solved and that ab ini- vantage of using a finite-temperature LDA is the smooth
tio MD calculations for open-shell transition metals are variation of the occupation of the orbitals; the crossing
possible. Our approach is based on a finite-temperature of levels close to the Fermi energy does not lead to insta-
density-functional approximation for the electrons and bilities in the equation of motion.
an exact determination of the LDA ground state after The LDA ground state is calculated after each molec-
each molecular-dynamics step using conjugate-gradient ular dynamics step using a method similar to the conju-
techniques — this completely avoids problems connected gate gradient techniques pioneered by Teter, Payne, and
with nonadiabaticity and with the crossing of energy lev- Allan. ' ' The method is a doubly iterative one: In
els close to the Fermi edge. The feasibility of dynamical the inner loop the wave functions at each k point in the
finite-temperature LDA simulations for s, p bonded ele- Brillouin zone (BZ) and each band are improved by a
ments using these techniques has been demonstrated in preconditioned conjugate-gradient method as described
our recent paper. The problem of the large energy cutofI' in Ref. 12 until the change in the energy eigenvalue
for transition metals is solved by adopting the concept of is smaller than 10 eV (or smaller then 40% of the
"ultrasoft" pseudopotentials introduced by Vanderbilt. change in the first step). After running over all bands
The success of the technique is demonstrated in the ex- (including a sufficient number of empty bands), a sub-
ample of simulations for liquid copper and vanadium. We space diagonalization is performed; the new fractional
chose Cu, because convergence with respect to the energy occupancies f; are calculated using a Gaussian broaden-
cutofF will be most diFicult to achieve for the narrow d ing of 0 = 0.1 eV, and the charge density is updated. To
band and the localized d state of the noble metals, and prevent charge sloshing, the mixing scheme proposed by
V as the other example because the problems of nonadi- Kerker is used. The minimization of the electronic free
abaticity and convergence to the LDA ground state are energy [Eq. (1)] is terminated when the change in the
most serious for metals with a high density of states at energy becomes smaller than 10 eV per atom.

0163-1829/93/48(17)/13115(4)/$06. 00 13 115 1993 The American Physical Society

13 116 BRIEF REPORTS

The atomic motion is described by Nose dynamics a.u. for V, R, = R,„= 2.7 a.u. for Cu), whereas the
generating a canonical ensemble at fixed temperature; d part was described using the new scheme (R"& —2. 0
the equations of motion are integrated using a fourth- a.u. , B"q —2.6 a.u. for V, B,"q —2.0 a. u. , B,"q —2.7
order predictor-corrector algorithm which allows one a.u. for Cu; two reference energies were used). To de-
to use time steps as large as 1.5 x 10 5 s with good scribe the wave functions, a relatively small energy cutoff
energy conservation. After moving the atoms, the wave E,„t —6 G, „,/(2m, ) is necessary. The action of the local
,
functions for the new configuration are estimated by a potential on the wave function and the smooth part of the
subspace alignment procedure proposed by Arias, Payne, charge density can be calculated using a relatively coarse
and Joannopoulos. To obtain reasonable predictions for fast-Fourier-transform (FFT) grid which must contain all
the wave functions of the highest occupied orbitals it is wave vectors up to G = 2G, „t. Nonlocality is handled in
absolutely essential to calculate the wave functions for a the real space projection scheme. A finer grid is neces-
sufBcient number of empty states, because the new wave sary to represent the augmentation charges and for the
functions are a linear combination of old wave functions calculation of the Hartree and exchange-correlation po-
over a certain energy range. tentials. All operations involving Q, ~. (r) are performed
The electron-ion interaction is-described by ultrasoft on this second grid in real space. The total time spent
pseudopotentials recently proposed by Vanderbilt. In on the fine grid scales linearly with the number of ions
this scheme norm conservation is given up, and the and is negligible compared with the computational costs
pseudization is done at several reference energies. To necessary for one conjugate-gradient step on the wave
correct the charge deficit (i.e. , the difFerence of the pseu- functions. One of the consequences of the introduction
docharge density and the exact charge density in the core of ultrasoft pseudopotentials is that the gradients of the
region) augmentation charges centered around each ion free energy with respect to the orbitals P„) are now given
~

must be added; i.e. , the total valence electron density is by

represented in the form

(2)
22
where S = 1++, Q, ~~P;)(P~~ is an overlap matrix [note
that (6) holds only if the Hamiltonian is diagonal in the
where subspace spanned by ~P )]. Similarly the forces on the
ions at R~ are given by (same assumption)

"'
P, are normconserving pseudo wave functions with suf-
ficient small pseudizing radii to mimic the exact all-
d F
d RN. ). 0 [H(V-, R~) —e„S(R~)]
BR~
electron wave function with good accuracy, P,". ' are ul-
trasoft pseudo wave functions constructed with large (7)
pseudizing radii and without the constraint of norm con- Details of the modified conjugate-gradient algorithm will
servation, and i = lme is shorthand for angular momen- be published elsewhere.
tum and magnetic quantum number and the reference
We tested the ultrasoft pseudopotentials for Cu and
energy e where the logarithmic derivatives of the pseudo V by performing total energy calculations for the crys-
and all electron wave functions are adjusted. The weights
talline phases. With an energy cutoff of 150 eV for Cu
p;~ of the augmentation charges are given by and V the equilibrium lattice constant, bulk modulus,
and the structural energy differences for different phases
(4) are fully converged and in good agreement with the full-
potential linear mufBn-tin-orbital calculations of Pax-
where the P s are localized functions spanning the core ton, Methfessel, and Polatoglou, if the same exchange-
region, defined by (P;~gP. ') = 8;~. The norm conserving correlation functional is used (Wigner interpolation ).
"'
pseudo wave functions P, are written as a sum of three For the simulation of the liquid metals, we used the
Bessel functions, whose wave vectors q,' are chosen so exchange-correlations functional of Ceperley and Alder
that their logarithmic derivatives match that of the all- parametrized by Perdew and Zunger, which reduces the
electron wave function, P;-, at the cutofF radius R, equilibrium volume by 3% (copper) to 5% (vanadium).
Details will be published elsewhere.
For the simulation of liquid copper and vanadium we
considered ensembles of 50 atoms in a periodically re-
peated cubic cell of side L = 16.46 a. u. for copper (cor-
The coefFicients o, are chosen so that the second deriva- responding to the experimental density p = 7.98 g cm
tive of P,"' is continuous at R, . For the ultrasoft wave at T = 1500 K) and I
= 17.46 a. u. for vanadium (corre-
function a set of two Bessel functions is suKcient. We sponding to a density of p = 5.36 g cm at T = 2200 K).
found that the optimization scheme proposed by Rappe The simulation is started for a configuration of the liquid
et al. does not improve convergence properties without generated using classical molecular dynamics and tight-
compromising transferability (details will be given else- binding-bond pair forces. For the ab initio molecular-
where, Bef. 20). For the s- and p- parts a conventional dynamics calculations the wave functions at the I' point
norm conserving potential was used (R, = R, „=2.6 were calculated for 300 bands and 175 bands for Cu and
 BRIEF REPORTS 13 117

V respectively, i.e. , 25 and 50 bands more than neces- 3

sary to accommodate the 48 and 3d valence electrons. (a)
After switching to the quantum many-body forces, the
system approached equilibrium very rapidly. With a
time step of Lt = 1.5 x 10 ' s and a Nose mass of
Q = (0.5 x 10s)m, a. u. the system shows excellent en-
ergy conservation: For vanadium the total free energy of
the coupled ion-electron-thermostat system remains con-
stant within 0.008 eV per atom over a run of 1.0 ps (650
time steps) (see Fig. 1); i.e. , the drift is smaller than 0 I I I I I « I I I I I I I I I I I I I I I I I I I I I

0.1% of the cohesive energy. 4 6

In Fig. 2 we report the pair correlation functions of R I)()
l-Cu and l-V as obtained by our simulation, and we
compare them with experimental data. For I-Cu the 3
(b3
agreement between theory and experiment is excellent,
even for the second and third peaks in g(R). This is re-
markable in view of the small size of the simulation cell.
The position of the first peak agrees well with nearest-
neighbor distance in crystalline face-centered-cubic Cu
(d = 2.55 A); the somewhat asymmetric shape of the
peak corresponds to a close random packing of rather
hard spheres. The coordination number evaluated by in- H

tegrating over the first peak in the radial distribution 0 I I I I I I I I I I I

I
I I I I I I I I I I I I I I

function up to the first minimum is N = 13.6+0.5. The 0 6

agreement between theory and experiment is less spectac- R ()(I
ular for I-V, for two reasons. First, it is known that the FIG. 2. Pair correlation function g(B) for liquid Cu (a) and
LDA and especially the Ceperley-Alder LDA functional V (b). Solid line, ab initio MD; squares, experimental data
lead to an underestimate of the equilibrium density com- taken from Ref. 21.
pared with experiment. Second, the diÃraction data for
the liquid "early" 3d metals Ti, V, Cr, Mn are of limited I

accuracy, due to the difhculties arising from high melt- I

(a)
ing temperatures and high vapor pressure. Nonetheless, I
I

the results are interesting: The first peak in g(R) is much I

I
I

broader and more symmetric; it is centered around the I

I

theoretical first and second neighbor distances in body- I

I

centered cubic vanadium (dI —2. 52 A, d2 —2.93 A), and

I
I
I

the coordination number is % = 12.6 + 0.4. Evidently

I

this arises from much softer interatomic potentials. The l

important thing is that the local order in the liquid metal

I
I

Is' mII
reflects the local order in the stable crystalline phase and
that this is correctly predicted by ab initio simulations. -10 -5
As a brief check of the dynamical properties we have E (eV1
calculated the self-diAusion coefBcients. For L-Cu we find 3, 0
D = 5.6 + 0.4 x 10 m /s at T = 1500 K in excellent E (b)
2. 5--

-8. 1
—8. 2
1. 0
-8. 3 0 5.-
-8. 4
—8. 5-
4 -2
E (eV)
-8. 6 a I I I I I I I I I I I I I I I I I I I I « I I I I I I I I I I I I I I ( I I I I I I I I I I I I I I
FIG. 3. Electronic density of states for fcc and l-Cu (a) and
0. 5 1. 0 for bcc and I-V (b). Dotted curve, crystalline DOS; dashed
curve, liquid DOS, calculated from the eigenvalues at the
time (ps )
I'-point and Gaussian broadening; solid curve, liquid DOS,
FIG. 1. Variation of the total free energy (upper curve) and calculated on a 6 x 6 x 6 grid of k points. For Cu the dashed
of the potential energy (lower curve) of liquid V at T = 2200 curve is aligned with the solid curve; the corresponding Fermi
K along an ab initio MD run of 1 ps (see text). level is shown by a dashed vertical line (see text).
13 118 BRIEF REPORTS 48

agreement with the experimental values of D = 4.6 x l0 are in excellent agreement. We think that our method is
m/s at T = 1413 K and D = 6. 0x10 m/s at computationally more eKcient, because of the possibility
T = 1533 K, respectively. For vanadium we find D = to use a much larger time step (typically a factor of 5—
9.3 x 10 m /s at T = 2200 K; no experimental data 10; a time step of 0.26 fs was used in Ref. 25) outweighs
are available. Finally we show in Fig. 3 the electronic the increased computational efFort arising for the calcu-
densities of states (DOS). A calculation of the DOS for lation of the LDA ground state after each time step (due
the I' point only leads to an s-electron DOS with a se- to the success of the prediction of the new wave func-
ries of pseudogaps which are a consequence of the rather tion by subspace alignment only three to four conjugate-
small sample. This could lead to a misestimate of the po- gradient iterations are needed; one conjugate gradient
sition of the d-band maximum relative to the Fermi edge. step takes about twice the time as a conventional Car-
Therefore, the calculation of the DOS was repeated for a Parrinello time step). However, the advantages of the
small number of configurations on a 6 x 6 x 6 grid in the method should be most important for metals with an
irreducible wedge of the simple cubic BZ. The resulting open d shell. In this case the high DOS at the Fermi
DOS for the liquid metals is smooth and astonishingly level makes the nonadiabaticity problem particularly se-
close to that of the stable crystalline phase. For l-Cu the vere, and atomic displacement will induce a strong mix-
similarity of the electronic spectra of /- and fcc-Cu is con- ing of the electron states around the Fermi level. With
firmed by the photoemission spectra of Norris. For l-V our technique, these problems are perfectly under con-
the calculated DOS shows a broad minimum just above trol.
the Fermi level, reHecting the deep DOS-minimum char- To conclude, we have shown that ab initio finite-
acteristic for the bcc 3d metals. Hence the differences in temperature local-density-functional MD simulations are
the local atomic arrangement in the two liquid metals is possible for elements with closed and open d shells. This
reBected in their electronic structure. will allow us to study the structural and electronic prop-
For L-Cu, our results may be compared with recent ab erties in liquids, glasses, and clusters containing transi-
initio MD simulations by Pasquarello et at. These cal- tion metals and at transition metal surfaces and inter-
culations are based on ultrasoft pseudopotentials similar faces using parameter-free ab initio techniques.
to ours, a pseudo-Newtonian dynamics for the electronic
degrees of freedom, and a two-thermostat technique for Our work has been supported by Siemens Nixdorf Aus-
the control of adiabaticity. The function of the electronic tria within the contract of cooperation with the Tech-
Nose thermostat is to drain the energy transferred from nische Universitat Wien. Most calculations have been
the ions to the electrons and to put a limit to the devi- done on the Fujitsu-SNI S100 computer of the Computer
ation from the Born-Oppenheimer surface. Both results Center of the Technische Universitat Wien.

R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). S. Nose, J. Chem. Phys. 81, 511 (1984).
G. Pastore, E. Smargiassi, and I'. Buda, Phys. Rev. A 44, A. Arnold, N. Mauser, and J. Hafner, J. Phys. Condens.
6334 (1991). Matter 1, 965 (1989).
N. D. Mermin, Phys. Rev. 137, A1441 (1965). A. M. Rappe, K.M. Rabe, E. Kaxiras, and J.D. Joannopou-
M. P. Teter, M. C. Payne, and D. C. Allan, Phys. Rev. B 40, los, Phys. Rev. B 41, 1227 (1990).
12255 (1989). A. T. Paxton, M. Methfessel, and H. M. Polatoglou, Phys.
G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993). Rev. B 41, 8127 (1990).
D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). E. Wigner, Phys. Rev. 133, 1002 (1934).
A. De Vita, Ph. D. thesis, Keele University, 1992; A. De J. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
Vita and M. J. Gillan (unpublished). G. Kresse (unpublished).
M. Weinert and J.W. Davenport, Phys. Rev. B 45, 13709 Ch. Hausleitner, G. Kahl, and J. Hafner, J. Phys. Condens.
(1992). Matter 3, 1589 (1991).
R.M. Wentzcovitch, J.L. Martins, and P.B. Allen, Phys. Y'. Waseda, The Structure of Non Crystalline Mat-erials-
Rev. B 45, 11372 (1992). Liquids and Amorphous Solids (McGraw-Hill, New York,
R.D. King-Smith, M. C. Payne, and J.S. Lin, Phys. Rev. B 1980).
44, 13063 (1991). T. Iida and R.I.L. Guthrie, The Physical Properties of Liq
T.A. Arias, M. C. Payne, and J.D. Joannopoulos, Phys. uid Metals (Clarendon, Oxford, 1988).
Rev. B 45, 1538 (1992). C. Norris, in Liquid Metals 2976, edited by R. Evans and
D.M. Bylander, L. Kleinman, and S. Lee, Phys. Rev. B 42, D. A. Greenwood (Institute of Physics, Bristol, 1977).
1394 (1990). A. Pasquarello, K. Laasonen, R. Car, C. Lee, and D. Van-
G. P. Kerker, Phys. Rev. B 23, 3082 (1981). derbilt, Phys. Rev. Lett. 69, 1982 (1992).