PHYSICAL REVIEW B

VOLUME 48, NUMBER 17

NOVEMBER 1993-I

Ab initio molecular

dynamics

for open-shell transition

metals

G. Kresse and
Institut fur Theoretische Physik Technische

J. Hafner
8-10, A-10/0 Wien, Austria

Universitat Wien, Wiedner IIauptstrasse (Received 2 July 1993)

We show that quantum-mechanical simulations in a finite-temperature localmolecular-dynamics density approximation based on the calculation of the electronic ground state and of the HellmannFeynman 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 problem with classical (LDA) to the many-electron molecular dynamics (MD) has marked an important progress in the attempts to put the atomistic simulation of the structural and electronic properties of materials on a firm quantum-mechanical basis. The approach origis based on the inally proposed by Car and Parrinello introduction of a pseudo-Newtonian dynamics describing the evolution of the electronic degrees of freedom close to the adiabatic ground state. An essential condition for the practicability of the Car-Parrinello method is that the transfer of energy between the atomic and electronic subsystems be small in order to prevent the electrons from drifting away &om the ground state. Metals are notoriously difIicult to handle in the Car-Parrinello method because in the absence of a gap the resonance between the atomic and electronic dynamics drives the system rapidly into nonadiabaticity. For transition metals an additional difIiculty arises from the fact that even modern soft-core pseudopotentials require a large energy cutofF for complete total energy convergence and hence a plane-wave basis that is too large to allow for a'6 initio MD simulations. The aim of the present paper is to demonstrate that these problems can be solved and that ab initio MD calculations for open-shell transition metals are possible. Our approach is based on a finite-temperature density-functional for the electrons and approximation an exact determination of the LDA ground state after each molecular-dynamics step using conjugate-gradient this completely avoids problems connected techniques with nonadiabaticity and with the crossing of energy levels close to the Fermi edge. The feasibility of dynamical finite-temperature LDA simulations for s, p bonded elements using these techniques has been demonstrated in our recent paper. The problem of the large energy cutofI' for transition metals is solved by adopting the concept of "ultrasoft" pseudopotentials introduced by Vanderbilt. The success of the technique is demonstrated in the example of simulations for liquid copper and vanadium. We chose Cu, because convergence with respect to the energy cutofF will be most diFicult to achieve for the narrow d band and the localized d state of the noble metals, and V as the other example because the problems of nonadiabaticity and convergence to the LDA ground state are most serious for metals with a high density of states at

the Fermi level. In the limit of zero temperature the total energy of the ground state, E[n(r)), is stationary with respect to variations of the electron density. At finite temperature the free energy
o E[n(r), f;, p] = E[n(r)] S[f;] + p

W, f; i

(1)

is the proper variational functional which has to be minimized with respect to variations of the electron density n(r), the orbital occupancy f, , and the chemical potential p. At a given temperature the occupation f, of a state with energy e,. is given by the Fermi-Dirac distribution and the entropy S[f,] corresponds to noninteracting electrons (in that case o = Ic~T) Gaussian .broadening of the one-electron levels can be computationally more convenient than Fermi-Dirac broadening. In this case o is the width of the Gaussian; the corresponding form for the entropy has been worked out in Ref. 7. The variational derivatives of the free energy with respect to the atomic positions are the Hellmann-Feynman forces. ' The advantage of using a finite-temperature LDA is the smooth variation of the occupation of the orbitals; the crossing of levels close to the Fermi energy does not lead to instabilities in the equation of motion. The LDA ground state is calculated after each molecular dynamics step using a method similar to the conjugate gradient techniques pioneered by Teter, Payne, and Allan. ' ' The method is a doubly iterative one: In the inner loop the wave functions at each k point in the Brillouin zone (BZ) and each band are improved by a preconditioned conjugate-gradient method as described in Ref. 12 until the change in the energy eigenvalue is smaller than 10 eV (or smaller then 40% of the change in the first step). After running over all bands (including a sufficient number of empty bands), a subspace diagonalization is performed; the new fractional occupancies f; are calculated using a Gaussian broadening of 0 = 0.1 eV, and the charge density is updated. To prevent charge sloshing, the mixing scheme proposed by Kerker is used. The minimization of the electronic free energy [Eq. (1)] is terminated when the change in the eV per atom. energy becomes smaller than 10
13 115

0163-1829/93/48(17)/13115(4)/$06. 00

1993

The American Physical Society

13 116

BRIEF REPORTS

The atomic motion is described by Nose dynamics generating a canonical ensemble at fixed temperature; the equations of motion are integrated using a fourthorder predictor-corrector algorithm which allows one to use time steps as large as 1.5 x 10 5 s with good energy conservation. After moving the atoms, the wave functions for the new configuration are estimated by a subspace alignment procedure proposed by Arias, Payne, and Joannopoulos. To obtain reasonable predictions for the wave functions of the highest occupied orbitals it is absolutely essential to calculate the wave functions for a sufBcient number of empty states, because the new wave functions are a linear combination of old wave functions over a certain energy range. The electron-ion interaction is-described by ultrasoft pseudopotentials recently proposed by Vanderbilt. In this scheme norm conservation is given up, and the pseudization is done at several reference energies. To correct the charge deficit (i.e. , the difFerence of the pseudocharge density and the exact charge density in the core region) augmentation charges centered around each ion must be added; i.e. , the total valence electron density is represented in the form
(2)
22

a.u. for V,

for Cu), whereas the using the new scheme (R"& 2. 0 " " a.u. , B"q 2.6 a.u. for V, B,q 2.0 a. u. , B,q 2.7 a.u. for Cu; two reference energies were used). To describe the wave functions, a relatively small energy cutoff E,t 6 G, ,/(2m, ) is necessary. The action of the local potential on the wave function and the smooth part of the charge density can be calculated using a relatively coarse fast-Fourier-transform (FFT) grid which must contain all wave vectors up to G = 2G, t. Nonlocality is handled in the real space projection scheme. A finer grid is necessary to represent the augmentation charges and for the calculation of the Hartree and exchange-correlation potentials. All operations involving Q, ~. (r) are performed on this second grid in real space. The total time spent on the fine grid scales linearly with the number of ions and is negligible compared with the computational costs necessary for one conjugate-gradient step on the wave functions. One of the consequences of the introduction of ultrasoft pseudopotentials is that the gradients of the free energy with respect to the orbitals P) are now given by
d part was described

R, = R,= 2.7 a.u.

where

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

pseudo wave functions with sufficient small pseudizing radii to mimic the exact all". electron wave function with good accuracy, P, ' are ultrasoft pseudo wave functions constructed with large pseudizing radii and without the constraint of norm conservation, and i = lme is shorthand for angular momentum and magnetic quantum number and the reference energy e where the logarithmic derivatives of the pseudo and all electron wave functions are adjusted. The weights charges are given by p;~ of the augmentation

"' P, are normconserving

).

0 [H(V-, R~) eS(R~)]

BR~
(7)

(4)
where the P s are localized functions spanning the core region, defined by (P;~gP. ') = 8;~. The norm conserving pseudo wave functions P, are written as a sum of three Bessel functions, whose wave vectors q,' are chosen so that their logarithmic derivatives match that of the allelectron wave function, at the cutofF radius

"'

P;-,

R,

The coefFicients o, are chosen so that the second deriva"' tive of P, is continuous at R, . For the ultrasoft wave function a set of two Bessel functions is suKcient. We found that the optimization scheme proposed by Rappe et al. does not improve convergence properties without compromising transferability (details will be given elsewhere, Bef. 20). For the s- and p- parts a conventional norm conserving potential was used (R, = R, 2.6

Details of the modified conjugate-gradient algorithm will be published elsewhere. We tested the ultrasoft pseudopotentials for Cu and V by performing total energy calculations for the crystalline phases. With an energy cutoff of 150 eV for Cu and V the equilibrium lattice constant, bulk modulus, and the structural energy differences for different phases are fully converged and in good agreement with the fullpotential linear mufBn-tin-orbital calculations of Paxif the same exchangeton, Methfessel, and Polatoglou, correlation functional is used (Wigner interpolation ). For the simulation of the liquid metals, we used the exchange-correlations functional of Ceperley and Alder which reduces the parametrized by Perdew and Zunger, 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 repeated cubic cell of side L = 16.46 a. u. for copper (corresponding to the experimental density p = 7.98 g cm = 17.46 a. u. for vanadium (correat T = 1500 K) and at T = 2200 K). sponding to a density of p = 5.36 g cm The simulation is started for a configuration of the liquid generated using classical molecular dynamics and tightbinding-bond pair forces. For the ab initio moleculardynamics calculations the wave functions at the I' point were calculated for 300 bands and 175 bands for Cu and

BRIEF REPORTS
V respectively, i.e. , 25 and 50 bands more than necessary to accommodate the 48 and 3d valence electrons. After switching to the quantum many-body forces, the With a system approached equilibrium very rapidly. 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 energy conservation: For vanadium the total free energy of the coupled ion-electron-thermostat system remains constant 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.1% of the cohesive energy. In Fig. 2 we report the pair correlation functions of l-Cu and l-V as obtained by our simulation, and we For I-Cu the compare them with experimental data. agreement between theory and experiment is excellent, even for the second and third peaks in g(R). This is remarkable in view of the small size of the simulation cell. The position of the first peak agrees well with nearestneighbor 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 integrating over the first peak in the radial distribution function up to the first minimum is N = 13.6+0.5. The agreement between theory and experiment is less spectacular for I-V, for two reasons. First, it is known that the LDA and especially the Ceperley-Alder LDA functional lead to an underestimate of the equilibrium density compared with experiment. Second, the diraction data for the liquid "early" 3d metals Ti, V, Cr, Mn are of limited accuracy, due to the difhculties arising from high meltNonetheless, ing temperatures and high vapor pressure. the results are interesting: The first peak in g(R) is much broader and more symmetric; it is centered around the theoretical first and second neighbor distances in bodycentered cubic vanadium (dI 52 A, d2 93 A), and 2. 2. the coordination number is % = 12.6 + 0.4. Evidently this arises from much softer interatomic potentials. The important thing is that the local order in the liquid metal reflects the local order in the stable crystalline phase and that this is correctly predicted by ab initio simulations. As a brief check of the dynamical properties we have calculated the self-diAusion coefBcients. For L-Cu we find D = 5.6 + 0.4 x 10 m /s at T = 1500 K in excellent
3

13 117

(a)

4
R

I)()
(b3

I I

0
R

()(I

FIG. 2. Pair correlation function g(B) for liquid Cu (a) and V (b). Solid line, ab initio MD; squares, experimental data taken from Ref. 21.
I

I
I I I I I I I I I

(a)

I I I I I l
I I

Is' mII

-10
E

-5
E (eV1

3, 0

2. 5--

(b)

-8. 1 8. 2 -8. 3 -8. 4 8. 5-8. 6 a

1. 0
0

5.4
E

-2
(eV)

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

0. 5
time

1. 0
)

(ps

FIG. 1. Variation of the total free energy (upper curve) and of the potential energy (lower curve) of liquid V at T = 2200 K along an ab initio MD run of 1 ps (see text).

FIG. 3. Electronic density of states for fcc and l-Cu (a) and for bcc and I-V (b). Dotted curve, crystalline DOS; dashed curve, liquid DOS, calculated from the eigenvalues at the I'-point and Gaussian broadening; solid curve, liquid DOS, calculated on a 6 x 6 x 6 grid of k points. For Cu the dashed curve is aligned with the solid curve; the corresponding Fermi 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 m/s at T = 1413 K and D = 6. m/s at For vanadium we find D = T = 1533 K, respectively. 9.3 x 10 m /s at T = 2200 K; no experimental data are available. Finally we show in Fig. 3 the electronic densities of states (DOS). A calculation of the DOS for the I' point only leads to an s-electron DOS with a series of pseudogaps which are a consequence of the rather small sample. This could lead to a misestimate of the position of the d-band maximum relative to the Fermi edge. Therefore, the calculation of the DOS was repeated for a small number of configurations on a 6 x 6 x 6 grid in the irreducible wedge of the simple cubic BZ. The resulting DOS for the liquid metals is smooth and astonishingly close to that of the stable crystalline phase. For l-Cu the similarity of the electronic spectra of /- and fcc-Cu is conFor l-V firmed by the photoemission spectra of Norris. the calculated DOS shows a broad minimum just above the Fermi level, reHecting the deep DOS-minimum characteristic for the bcc 3d metals. Hence the differences in the local atomic arrangement in the two liquid metals is reBected in their electronic structure. For L-Cu, our results may be compared with recent ab initio MD simulations by Pasquarello et at. These calculations are based on ultrasoft pseudopotentials similar to ours, a pseudo-Newtonian dynamics for the electronic degrees of freedom, and a two-thermostat technique for the control of adiabaticity. The function of the electronic Nose thermostat is to drain the energy transferred from the ions to the electrons and to put a limit to the deviation from the Born-Oppenheimer surface. Both results

0x10

are in excellent agreement. We think that our method is computationally more eKcient, because of the possibility to use a much larger time step (typically a factor of 5 10; a time step of 0.26 fs was used in Ref. 25) outweighs the increased computational efFort arising for the calculation of the LDA ground state after each time step (due to the success of the prediction of the new wave function by subspace alignment only three to four conjugategradient iterations are needed; one conjugate gradient step takes about twice the time as a conventional CarParrinello time step). However, the advantages of the method should be most important for metals with an open d shell. In this case the high DOS at the Fermi level makes the nonadiabaticity problem particularly severe, and atomic displacement will induce a strong mixing of the electron states around the Fermi level. With our technique, these problems are perfectly under con-

trol.
To conclude, we have shown that ab initio finitetemperature local-density-functional MD simulations are possible for elements with closed and open d shells. This will allow us to study the structural and electronic properties in liquids, glasses, and clusters containing transition metals and at transition metal surfaces and interfaces using parameter-free ab initio techniques.
Our work has been supported by Siemens Nixdorf Aus-

tria within the contract of cooperation with the Technische Universitat Wien. Most calculations have been done on the Fujitsu-SNI S100 computer of the Computer Center of the Technische Universitat Wien.

Vita and M. J. Gillan (unpublished). M. Weinert and J.W. Davenport, Phys. Rev. B 45, 13709

R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). G. Pastore, E. Smargiassi, and I'. Buda, Phys. Rev. A 44, 6334 (1991). N. D. Mermin, Phys. Rev. 137, A1441 (1965). M. P. Teter, M. C. Payne, and D. C. Allan, Phys. Rev. B 40, 12255 (1989). G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993). D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). A. De Vita, Ph. D. thesis, Keele University, 1992; A. De

(1992). R.M. Wentzcovitch, J.L. Martins, and P.B. Allen, Phys. Rev. B 45, 11372 (1992). R.D. King-Smith, M. C. Payne, and J.S. Lin, Phys. Rev. B 44, 13063 (1991). T.A. Arias, M. C. Payne, and J.D. Joannopoulos, Phys. Rev. B 45, 1538 (1992). D.M. Bylander, L. Kleinman, and S. Lee, Phys. Rev. B 42, 1394 (1990). G. P. Kerker, Phys. Rev. B 23, 3082 (1981).

S. Nose, J. Chem. Phys. 81, 511 (1984). A. Arnold, N. Mauser, and J. Hafner, J. Phys. Condens. Matter 1, 965 (1989). A. M. Rappe, K.M. Rabe, E. Kaxiras, and J.D. Joannopoulos, Phys. Rev. B 41, 1227 (1990). A. T. Paxton, M. Methfessel, and H. M. Polatoglou, Phys. Rev. B 41, 8127 (1990). E. Wigner, Phys. Rev. 133, 1002 (1934). J. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). G. Kresse (unpublished). Ch. Hausleitner, G. Kahl, and J. Hafner, J. Phys. Condens. Matter 3, 1589 (1991). Y'. Waseda, The Structure of Non Crystalline Mat-erialsLiquids
and Amorphous

Solids (McGraw-Hill,

New York,

1980).

T. Iida

and R.I.L. Guthrie, The Physical Properties of Liq uid Metals (Clarendon, Oxford, 1988). C. Norris, in Liquid Metals 2976, edited by R. Evans and D. A. Greenwood (Institute of Physics, Bristol, 1977). A. Pasquarello, K. Laasonen, R. Car, C. Lee, and D. Van-

derbilt, Phys. Rev. Lett. 69, 1982 (1992).