Ab Initio Molecular Dynamics with the Projector Augmented Wave

Chapter 4

Ab Initio Molecular Dynamics with the Projector Augmented Wave Method Peter E. Blöchl , Peter Margl , and Karlheinz Schwarz 1

Downloaded by MONASH UNIV on March 16, 2013 | http://pubs.acs.org Publication Date: May 5, 1996 | doi: 10.1021/bk-1996-0629.ch004

1

2

3

IBM Research Division, Zurich Research Laboratory, Säumerstrasse 4, CH-8803 Rüschlikon, Switzerland Department of Chemistry, University of Calgary, 2500 University Drive Northwest, Calgary, Alberta T2N 1N4, Canada Technical University Vienna, Getreidemarkt 9/158, A-1060 Vienna, Austria 2

3

A n introduction to the ab-initio molecular dynamics approach of Car and Parrinello and to the projector augmented wave method is given. The projector augmented wave method is an all-electron electronic structure method that allows ab-initio molecular dynamics simula tions to be performed accurately and efficiently even for first-row and transition metal elements. We describe the supercell approach and how it can be extended to isolated charged or polar molecules. A p plications to organometallic compounds, including ferrocene and the fluxional molecule beryllocene, demonstrate the capabilities of this methodology.

Ten years ago, C a r and Parrinello (1) invented the ab-initio molecular dynamics approach, which allows one to simulate the motion of the atoms from first p r i n c i ples. T h i s tool was promising for applications i n chemistry, because i t combines the accuracy of density functional calculations w i t h the ability to study the f i nite temperature dynamics of molecules. A new electronic structure method, the projector augmented wave ( P A W ) method (2), overcomes the limitations of the pseudo-potential approach, which has been employed i n most ab-initio molecular dynamics simulations. T h e P A W method makes the full all-electron wave func tions accessible and allows first-row and transition-metal elements to be studied w i t h moderate computational effort. T h i s method has proven useful i n a number of recent applications to physics (3, 4)-, chemistry (5-9), and biochemistry (10). T h e expectations that ab-initio molecular dynamics allows direct simulation of chemical processes on a picosecond time scale and a nanometer length scale are now being rapidly realized. 0097-6156/96/0629-0054$15.00Α) © 1996 American Chemical Society In Chemical Applications of Density-Functional Theory; Laird, B., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1996.

4. BLÔCHL ET AL.

Ab Initio Molecular Dynamics with PAW

55

In this article we sketch the basic ideas underlying the current methodology and describe some recent applications. We make an attempt to show the general picture at the expense of details that can be found i n more specialized publications. Furthermore this paper is restricted to our own work, and thus the reader should be aware that there are other interesting developments (11-13) i n this field that are not mentioned here.


Ab-initio Molecular Dynamics Approach Electronic structure calculations are usually performed for a static atomic structure. T h e dynamical behavior can be studied w i t h molecular dynamics techniques that describe the interaction between atoms by parameterized inter-atomic potentials w i t h parameters adjusted to experiment. Numerous such potentials have been devised and tuned to different environments. T h e i r computational simplicity allows very large systems to be studied for long time scales. However, parameterized inter-atomic potentials have difficulties describing chemical reactions, because of the complexity of the interactions during bond breaking and formation as well as the lack of experimental information required to specify the parameters of the essential part of the potential energy surface. A b - i n i t i o molecular dynamics is one way to combine the virtues of these two distinct approaches. Here the dynamics of the system is simulated, but i n each time step the forces are obtained directly from an accurate electronic structure calculation. E v e n though some simulations have been performed w i t h the H a r t r e e Fock method (14~16), most ab-initio molecular dynamics calculations are based on density functional theory (17, 18). W i t h i n the ab-initio molecular dynamics approach, it is currently possible to simulate systems having about a hundred atoms for a few to a few tens of picoseconds. Let us list here the time scales of some common processes to give a feeling for their magnitude: • Electronic transitions corresponding to visible light occur on a time scale of less than 2.5 fsec, so they can be regarded as instantaneous because the simulation is discretized i n steps of about 0.25 fsec. • A period of a carbon-hydrogen stretch frequency takes about 10 fsec-0.01 psec. Other molecular stretch vibrations and hydrogen bond-bending modes take about twice as long, i.e. 20 fsec or longer. • C h e m i c a l reactions have two distinct time scales: one is the waiting time discussed above, and the other is the reaction event. If we talk about chemical reactions we can estimate the time it takes to overcome a barrier ΕA, i.e. the waiting t i m e , as l / i / e x p ^ u ^ / f c ^ T ) , where ν is a typical molecular frequency. If we use ν = 1/20 fsec and at room temperature we find a waiting t i m e of 1.4 psec for a barrier of a 10 k J / m o l corresponding to the breaking of a hydrogen bond, but 10 psec for a barrier of 50 k J / m o l such as i n low-barrier chemical reactions. 7

In Chemical Applications of Density-Functional Theory; Laird, B., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1996.


56

C H E M I C A L APPLICATIONS O F DENSITY-FUNCTIONAL T H E O R Y

T h i s implies that at room temperature we can easily study molecular vibrations and low barrier processes such as hydrogen bond rearrangement or the fluxional motion of organometallic compounds. Note, however, that ab-initio molecular dynamics does not describe nuclear tunnel processes. C h e m i c a l reactions involving larger barriers, however, can be studied directly by simulations only at elevated temperatures. Above 1000 Κ we can access reactions involving barriers up to about 50 k J / m o l on a picosecond time scale. Whereas the waiting time is often substantially larger than the affordable simulation time, the reaction event typically occurs on the time scale of molecular vibrations, i.e. a few tenths or hundredths of a picosecond, which is well w i t h i n the reach of ab-initio molecular dynamics simulations. Thus for systems w i t h large barriers we first must locate the transition states using standard methods, and then we can study the dynamics of a particular reaction event by direct molecular dynamics simulations. W h a t is the underlying idea of ab-initio molecular dynamics? A s i n classical molecular dynamics, the nuclear trajectories can be obtained from Newton's law MR i

,

= -V E

i

Ri

(1)

where Ri is a coordinate of atom ζ, M , its mass and Ε the total energy expression. T h i s equation can be discret ized using the Ver let algorithm, which replaces the acceleration by M

t

)

=

M

±

^

m

±

M

^

l

,

( 2 )

where Δ is the time step. T h e Verlet algorithm appears to be superior to other, more complex algorithms if large time steps are chosen (19). T h e stability l i m i t of the Verlet algorithm lies at Δ = T j / 3 , where T \ is the oscillation period of the fastest vibration of the system. Hence the time step for a system w i t h C - H bonds must clearly be shorter than 3 fsec. In practice, we use a time step of 0.25 fsec. If we use the total energy of density functional theory to derive the forces, we must know the self-consistent wave functions at each time step. Considering that a typical simulation requires several thousand to ten thousand t i m e steps, it is hardly possible to perform an independent self-consistent calculation for each t i m e step. T h i s problem has been solved by C a r and Parrinello w i t h an ingeniously simple trick: T h e y introduced an additional classical equation of motion for the quantum mechanical wave functions describing the electrons m

m

n

n

™ |Φ ) = - # | Ψ „ ) + £ | * ) A m Φ

η

m

m n

,

(3)

where H is the H a m i l t o n i a n and A the m a t r i x of Lagrange multipliers, which ensure that the wave functions remain orthogonal to each other; 77ΐψ is a fictitious mass for the wave functions, a free parameter that ideally should be chosen very small. A t y p i c a l value for τηψ is 1000 a.u. It is not immediately obvious that this coupled system of equations yields mean ingful results. T h e m a i n requirement is that the electrons remain close to the electronic ground state, i.e. at the B o r n - O p p e n h e i m e r surface. T h i s requires the m

n



4. BLOCHL ET AL.


57

wave function coefficients to be " c o l d " . A t the same t i m e , however, the nuclei may be " h o t " , corresponding to a high temperature. A s w i t h any thermally coupled system the temperatures w i l l equilibrate, ultimately rendering the forces acting on the nuclei meaningless because the density functional theory requires the wave functions to be i n their ground state. T h e underlying reason why ab-initio molec ular dynamics simulations are still feasible is the adiabatic principle. It states that the thermal coupling of two subsystems w i t h well separated vibrational spectra is small. T h i s is indeed the case, at least for systems w i t h a finite H O M O - L U M O gap, and hence the heat transfer from hot nuclei to cold wave functions is minute, often not even noticeable on a picosecond time scale. In practice the simulation is made stationary using thermostats for the electronic and the atomic subsystem absorbing the remaining heat transfer (20, 21). T h e adiabatic principle can easily be demonstrated for a simple classical model: It consists of a light particle w i t h a mass m and a position r and a heavy particle with a mass M and a coordinate R. T h e two particles are coupled by a spring w i t h a given force constant c. T h e light particle is analogous to the wave functions, whereas the heavy particle is analogous to the nuclei. T h i s model reproduces many of the typical features of the coupled system of equations i n ab-initio molecular dynamics such as slow heat transfer and increased effective mass of the heavy particles (or the nuclei). We obtain a coupled system of equations MR

= F(t) - c(R - r)

mr = -c(r

- R)

,

(4)

where F(t) is a time-dependent force acting on the nuclei. T h i s system is trans formed into a center-of-mass coordinate X = (MR + mr)/(M + m) and a relative coordinate χ = r — R. If we fix the position of the heavy particle and relax the light particle, the relative coordinate vanishes. Hence χ = 0 characterizes the i n stantaneous ground state or the B o r n - O p p e n h e i m e r surface. Now our system of equations reads (M + m)X

=

F(t)

mx = -c(l

+ ^)x-^F(t)

.

(5)

For the adiabatic principle to work, the frequency ω = yjc/m should be much higher than the frequencies contributing to F(t). T h i s frequency can be controlled by choosing the mass of the light particle, which is analogous to the fictitious mass of the wave functions m $ , sufficiently small. T h a t this is indeed the case i n abinitio molecular dynamics, at least for closed-shell systems, has been demonstrated previously (22). In this l i m i t we can determine the relative coordinate analytically

where A and φ are arbitrary constants. We can now summarize the m a i n features that also occur i n ab-initio molecular dynamics.


58


• If we start on the B o r n - O p p e n h e i m e r surface, for example χ = χ = F = F = 0, the coefficient A remains zero for a l l times. T h e remaining deviation x(t) = ( l / u ; ) ( F ( < ) / m + M) from the instantaneous ground state remains small because of the small factor 1 /ω . 2


2

• T h e light particle can also embark on a free oscillation described by the sec ond term i n equation 6. T h i s is an uncontrolled deviation from the B o r n Oppenheimer surface. If such oscillations are large, the forces acting on the nuclei i n an ab-initio molecular dynamics simulation become meaningless. However, i n the adiabatic l i m i t , where ω goes to infinity, the particles de couple and there is no heat transfer. Therefore the systems remains on the B o r n - O p p e n h e i m e r surface, if that is where we started. In practice, if the adiabatic principle is only fulfilled approximately, the amplitude of the free oscillation w i l l increase, but only very slowly. Thus it can easily be controlled by a small constant friction or a thermostat acting on the small particle (21). • T h e center of mass variable follows the forces like a particle w i t h a slightly heavier mass M + m . T h i s implies that the nuclear mass needs to be renormalized, so that the effective mass of the electronic wave functions, which is fictitious, is accounted for. Ways to estimate this renormalization have been described previously (2). • T h e position of the heavy particle is, however, not identical to the center of mass variable. T h e difference is proportional to the deviation of the light par ticle from the ground state and m u l t i p l i e d by a small factor m/(M + ra). In many simulations a close look reveals a tiny j i t t e r i n the atomic trajectories, which is reminiscent of the remaining small free oscillation of the wave func tions. However, these oscillations are not random and cancel exactly as long as the electrons remain close to the B o r n - O p p e n h e i m e r surface. It should be noted that such difficulties as the necessity of renormalizing the masses are so small that they remained undetected for a long time. However, as they surfaced, remedies have been developed.

Projector A u g m e n t e d W a v e M e t h o d Before one can perform an electronic structure calculation, one must choose a basis set. Here, a tradeoff between the size of the basis set needed to obtain converged results and the effort to evaluate the total energy, the H a m i l t o n i a n , and the overlapm a t r i x elements, must be found. T h e smallest basis sets are probably obtained w i t h augmented wave methods such as the linear muffin t i n orbital ( L M T O ) method (23), which uses sophisticatedly composed basis sets. A t the other extreme are the pure plane wave methods. Simple plane waves require enormous basis set sizes, so i n practice pseudopotentials (24, 25) must be employed. O n the other hand, the integrations are t r i v i a l because the relevant operators can be cast i n a diagonal or separable (β6, 27) form using Fast Fourier transforms ( F F T ) , m a k i n g the



4. BLOCHL ET AL.


59

computational effort per basis function minute. Plane wave methods are promising, i n particular for density functional theory, because here the electron density needs to be represented on a real space grid. For plane waves this is achieved by just another F F T , whereas the same operation easily becomes costly for more complex basis sets. A further advantage of plane wave methods is that an enormous knowhow has been accumulated w i t h them for ab-initio molecular dynamics. However, even when pseudopotentials are employed, the basis set size required for describing first-row and transition-metal elements can be daunting. Further more, pseudopotentials obscure the electron density near the nucleus. T h i s puts information that depends on the local environment close to the nucleus, such as electric field gradients or hyperfine parameters, out of reach except through i n d i rect reconstruction techniques (28). These are clear disadvantages for most t y p i c a l applications i n chemistry. T h e design goal for the projector augmented wave ( P A W ) method was to com bine the advantages of the plane wave pseudopotential method w i t h that of the all-electron augmented wave methods such as the linear augmented plane wave method (23). T h e P A W method borrows the idea of augmentation from the existing a l l electron augmented wave methods and uses composed basis sets. In this way augmented wave methods can accommodate the various shapes of the wave func tions i n different regions. Close to the nucleus the kinetic energy is large, ruling out the use of a plane wave expansion. However, i n this so-called augmentation region the potential is almost spherically symmetric, suggesting an expansion of the wave function into functions similar to atomic orbitals, which are separated according to their angular momenta and treated on a radial grid. Here the potential is d o m i nated by the nucleus and the tightly bound core electrons, which are affected only slightly by the chemical environment of the atom. T h u s , a small number of basis functions is sufficient to describe the wave functions i n this region. O n the other hand, far from the nucleus the potential is shallow and the kinetic energy of the electrons is small. T h i s is the region where the chemical bonds are located, so the wave function varies drastically from one molecular environment to another. Hence plane waves, which form a complete basis set, seem ideally suited to this region. T h e P A W method describes the all-electron valence wave function as follows: far from the nucleus, i.e. farther than the covalent radius, the wave function is described by a plane wave expansion. Near the nuclei, however, we subtract and add partial waves, which are similar to atomic orbitals, to the plane wave part i n order to incorporate the proper nodal structure of the true wave functions. T h u s , the plane wave part near the nucleus can be a smooth continuation of the wave function outside the atomic regions. T h e partial wave expansions are sufficiently localized w i t h i n the augmentation regions that the overlap of partial waves centered on different nuclei need not be considered. In the P A W method an all-electron valence wave function |Φ) is therefore a superposition of three parts

|Φ) = |Φ> + ΣΙ^)(P,-|Φ)-ΣΙ^ΙΦ> •

t


(7)

60


Equation 7 is the definition of the basis set used i n the P A W method. It defines the variational degree of freedom for the all-electron wave functions. T h e plane wave expansion coefficients of |Φ) correspond to the orbital coefficients. E q u a t i o n 7 above can alternatively be read as a linear transformation from a fictitious pseudowave function to the corresponding all-electron wave function. We briefly comment on each of the three individual terms.


|Φ) : A smooth function extending over a l l space, which can be expanded i n plane waves. T h i s is called a pseudo-wave function and is denoted |Φ). One can think of the tilde as an operator that turns the all-electron wave function into the corresponding pseudo-wave function. Σ ι Ι&)(p,|Φ) ' A partial wave expansion of the all-electron wave function near an atom. P a r t i a l waves are solutions of the Schrodinger equation for a given energy and the atomic all-electron potential (note: they need not be bound states!). T h e y are denoted by the symbol |

Ab Initio Molecular Dynamics with the Projector Augmented Wave

Recommend Documents