introduction to Ab Initio Molecular Dynamics, Notas de estudo de Engenharia de Produção

Baixe introduction to Ab Initio Molecular Dynamics e outras Notas de estudo em PDF para Engenharia de Produção, somente na Docsity!

Lecture Notes

Introduction to

Ab Initio Molecular Dynamics

J¨urg Hutter

Physical Chemistry Institute

University of Zurich

Winterthurerstrasse 190

8057 Zurich, Switzerland

[email protected]

September 14, 2002

Contents

• 1 Molecular Dynamics
  • 1.1 Introduction
  • 1.2 Equations of Motion
  • 1.3 Microcanonical Ensemble
  • 1.4 Numerical Integration
  • 1.5 Extended System Approach
    • 1.5.1 Barostats
    • 1.5.2 Thermostats
• 2 Ab initio Molecular Dynamics
  • 2.1 Born–Oppenheimer Molecular Dynamics
    • 2.1.1 Forces in BOMD
  • 2.2 Car–Parrinello Molecular Dynamics
    • 2.2.1 How to Control Adiabaticity?
    • 2.2.2 Forces in CPMD
    • 2.2.3 Velocity Verlet Equations for CPMD
  • 2.3 Comparing BOMD and CPMD
• 3 Plane Waves
  • 3.1 Unit Cell and Plane Wave Basis
  • 3.2 Kinetic Energy and Local Potentials
  • 3.3 Electrostatic Energy
  • 3.4 Exchange and Correlation Energy
  • 3.5 Car–Parrinello Equations
• 4 Pseudopotentials
  • 4.1 Why Pseudopotentials?
  • 4.2 Norm–Conserving Pseudopotentials
    • 4.2.1 Hamann–Schl¨uter–Chiang Conditions
    • 4.2.2 Bachelet-Hamann-Schl¨uter (BHS) form
    • 4.2.3 Kerker Pseudopotentials
    • 4.2.4 Trouiller–Martins Pseudopotentials
    • 4.2.5 Kinetic Energy Optimized Pseudopotentials
  • 4.3 Pseudopotentials in the Plane Wave Basis
    • 4.3.1 Gauss–Hermit Integration
    • 4.3.2 Kleinman–Bylander Scheme
  • 4.4 Dual–Space Gaussian Pseudopotentials
  • 4.5 Example: Pseudopotentials for Oxygen
• 5 Implementation
  • 5.1 Total Energy and Gradients
    • 5.1.1 Plane Wave Expansion
    • 5.1.2 Total Energy
    • 5.1.3 Wavefunction Gradient
    • 5.1.4 Nuclear Gradient
  • 5.2 Fast Fourier Transforms
  • 5.3 Density and Force Calculations in Practice
  • 5.4 Saving Computer Time
  • 5.5 Exchange and Correlation Functionals
• 6 Extensions to the Original Method
  • 6.1 Non-linear Core Correction
    • Augmented–Wave Method 6.2 Ultrasoft Pseudopotentials and the Projector
    • 6.2.1 The PAW Transformation
    • 6.2.2 Expectation Values
    • 6.2.3 Ultrasoft Pseudopotentials
    • 6.2.4 PAW Energy Expression
    • 6.2.5 Integrating the CP Equations
  • 6.3 Metals; Free Energy Functional
  • 6.4 Charged Systems
• 7 Properties from CPMD
  • 7.1 Position Operator in Periodic Systems
  • 7.2 Dipole Moments and IR Spectra
  • 7.3 Localized Orbitals, Wannier Functions
• 8 Chemical Reactions with CPMD
  • 8.1 Simulated Annealing
  • 8.2 Rare Events
    • 8.2.1 Thermodynamic Integration
    • 8.2.2 Adiabatic Free Energy Sampling
    • 8.2.3 Bias Potential Methods
• 9 Linear Response
  • 9.1 Adiabatic Density–Functional Perturbation Theory
  • 9.2 Coupled Perturbed Kohn–Sham Equations
    • 9.2.1 Exchange-Correlation Functionals
  • 9.3 Nuclear Hessian
    • 9.3.1 Selected Eigenmodes of the Hessian
  • 9.4 Polarizability
  • 9.5 NMR Chemical Shifts
    • 9.5.1 Chemical shifts and susceptibilities
    • 9.5.2 The gauge origin problem
    • 9.5.3 The position operator problem
    • 9.5.4 Density functional perturbation theory
    • 9.5.5 Pseudopotential correction
• 10 Excited States
  • 10.1 Restricted Open–shell Kohn–Sham Method
  • 10.2 Time–Dependent Density Functional Theory
  • 10.3 Excited States from Linear Response
    • 10.3.1 Tamm-Dancoff Approximation

Lecture 1

Molecular Dynamics

1.1 Introduction

The aim of molecular dynamics is to model the detailed microscopic dynamical behavior of many different types of systems as found in chemistry, physics or biology. The history of molecular dynamics goes back to the mid 1950’s when first computer simulations on simple systems were performed[1]. Excellent modern textbooks [2, 3] on this topic can be found and collections of articles from summer schools are a very good source for in depth information on more specialized aspects [4, 5, 6]. Molecular Dynamics is a technique to investigate equilibrium and transport properties of many–body systems. The nuclear motion of the particles is modeled using the laws of classical mechanics. This is a very good approximation for molecular systems as long as the properties studied are not related to the motion of light atoms (i.e. hydrogen) or vibrations with a frequency ν such that hν > kB T. Extensions to classical molecular dynamics to incorporate quantum effects or full quantum dynamics (see for example Refs. [7, 8] for a starting point) is beyond the scope of this lecture series.

1.2 Equations of Motion

We consider a system of N particles moving under the influence of a potential function U [9, 10]. Particles are described by their positions R and momenta P = M V. The union of all positions (or momenta) {R 1 , R 2 ,... , RN } will be called RN^ (PN^ ). The potential is assumed to be a function of the positions only; U (RN^ ). The Hamiltonian H of this system is

H(RN^ , PN^ ) =

∑^ N

I=

P^2 I

2 MI

+ U (RN^ ). (1.1)

properties, implying that it is not necessary to know the precise detailed motion of every particle in a system in order to predict its properties. It is sufficient to simply average over a large number of identical systems, each in a different configuration; i.e. the macroscopic observables of a system are formulated in term of ensemble averages. Statistical ensembles are usually characterized by fixed values of thermodynamic variables such as energy, E; temperature, T ; pressure, P ; volume, V ; particle number, N ; or chemical potential μ. One fundamental ensemble is called the microcanonical ensemble and is characterized by constant particle number, N ; constant volume, V ; and constant total energy, E, and is denoted the N V E ensemble. Other examples include the canonical or N V T ensemble, the isothermal–isobaric or N P T ensemble, and the grand canonical or μV T ensemble. The thermodynamic variables that characterize an ensemble can be regarded as experimental control parameters that specify the conditions under which an experiment is performed. Now consider a system of N particles occupying a container of volume V and evolving under Hamilton’s equation of motion. The Hamiltonian will be constant and equal to the total energy E of the system. In addition, the number of particles and the volume are assumed to be fixed. Therefore, a dynamical trajectory (i.e. the positions and momenta of all particles over time) will generate a series of classical states having constant N , V , and E, corresponding to a microcanonical ensemble. If the dynamics generates all possible states then an average over this trajectory will yield the same result as an average in a microcanonical ensemble. The energy conservation condition, H(RN^ , R˙N^ ) = E which imposes a restriction on the classical microscopic states accessible to the system, defines a hypersurface in the phase space called a constant energy surface. A system evolving according to Hamilton’s equation of motion will remain on this surface. The assumption that a system, given an infinite amount of time, will cover the entire constant energy hypersurface is known as the ergodic hypothesis. Thus, under the ergodic hypothesis, averages over a trajectory of a system obeying Hamilton’s equation are equivalent to averages over the microcanonical ensemble. In addition to equilibrium quantities also dynamical properties are defined through ensemble averages. Time correlation functions are important because of their relation to transport coefficients and spectra via linear response theory [11, 12]. The important points are: by integration Hamilton’s equation of motion for a number of particles in a fixed volume, we can create a trajectory; time averages and time correlation functions of the trajectory are directly related to ensemble averages of the microcanonical ensemble.

1.4 Numerical Integration

In a computer experiment we will not be able to generate the true trajectory of a sys- tem with a given set of initial positions and velocities. For all potentials U used in real applications only numerical integration techniques can be applied. These techniques are

based on a discretization of time and a repeated calculation of the forces on the particles. Many such methods have been devised [13] (look for ”Integration of Ordinary Differen- tial Equations”). However, what we are looking for is a method with special properties: long time energy conservation and short time reversibility. It turns out that symplectic methods (they conserve the phase space measure) do have these properties. Long time energy conservation ensures that we stay on (in fact close) to the constant energy hyper- surface and the short time reversibility means that the discretize equations still exhibit the time reversible symmetry of the original differential equations. Using these methods the numerical trajectory will immediately diverge from the true trajectory (the divergence is exponential) but as they stay on the correct hypersurface they still sample the same microcanonical ensemble. On the other hand, a short time accurate method will manage to stay close to the true trajectory for a longer time and ultimately will also exponentially diverge but will not stay close to the correct energy hypersurface and therefore will not give the correct ensemble averages. Our method of choice is the velocity Verlet algorithm [14, 15]. It has the advantage that it uses as basic variables positions and velocities at the same time instant t. The velocity Verlet algorithm looks like a Taylor expansion for the coordinates:

R(t + ∆t) = R(t) + V(t)∆t +

F(t) 2 M

∆t^2. (1.9)

This equation is combined with the update for the velocities

V(t + ∆t) = V(t) +

F(t + ∆t) + F(t) 2 M

∆t^2. (1.10)

The velocity Verlet algorithm can easily be cast into a symmetric update procedure that looks in pseudo code

V(:) := V(:) + dt/(2M(:))F(:) R(:) := R(:) + dtV(:) Calculate new forces F(:) V(:) := V(:) + dt/(2M(:))*F(:)

To perform a computer experiment the initial values for positions and velocities have to be chosen together with an appropriate time step (discretization length) ∆t. The choice of ∆t will be discussed in more detail in a later chapter about ab initio molecular dynamics. The first part of the simulation is the equilibration phase in which strong fluctuation may occur. Once all important quantities are sufficiently equilibrated, the actual simulation is performed. Finally, observables are calculated from the trajectory. Some quantities that can easily be calculated are (for these and other quantities see the books by Frenkel and Smit [3] and Allen and Tildesley [2] and references therein)

• The average temperature 〈

M V^2 〉 =

kBT (1.11)

possibility of non–isotropic external stress; the additional fictitious degrees of freedom in the Parrinello–Rahman approach [17] are the lattice vectors of the supercell. These variable–cell approaches make it possible to study dynamically structural phase transitions in solids at finite temperatures. The basic idea to allow for changes in the cell shape consists in constructing an extended Lagrangian where the lattice vectors a 1 , a 2 and a 3 of the simulation cell are additional dynamical variables. Using the 3 × 3 matrix h = [a 1 , a 2 , a 3 ] (which fully defines the cell with volume det h = Ω) the real–space position RI of a particle in this original cell can be expressed as

RI = hSI (1.14)

where SI is a scaled coordinate with components Si,u ∈ [0, 1] that defines the position of the ith particle in a unit cube (i.e. Ωunit = 1) which is the scaled cell [17]. The resulting metric tensor G = hth converts distances measured in scaled coordinates to distances as given by the original coordinates. The variable–cell extended Lagrangian can be postulated

L =

∑^ N

I

MI

( S˙t I G S˙I

) − U (h, SN^ ) +

W Tr h˙t^ h˙ − p Ω , (1.15)

with additional nine dynamical degrees of freedom that are associated to the lattice vectors of the supercell h. Here, p defines the externally applied hydrostatic pressure, W defines the fictitious mass or inertia parameter that controls the time–scale of the motion of the cell h. In particular, this Lagrangian allows for symmetry–breaking fluctuations – which might be necessary to drive a solid–state phase transformation – to take place spontaneously. The resulting equations of motion read

MI S¨I,u = −

∑^3

v=

∂U (h, SN^ ) ∂RI,v

( ht

)− 1 vu

− MI

∑^3

v=

∑^3

s=

G uv−^1 G˙vs S˙I,s (1.16)

W h¨uv = Ω

∑^3

s=

( Πtot us − p δus

) ( ht

)− 1 sv

where the total internal stress tensor

Πtot us =

∑

I

MI

( S˙t I G S˙I

) us + Πus^ (1.18)

is the sum of the thermal contribution due to nuclear motion at finite temperature and the internal stress tensor Π derived from the interaction potential. A modern formulation of barostats that combines the equation of motion also with thermostats (see next section) was given by Martyna et al. [19].

1.5.2 Thermostats

Standard molecular dynamics generates the microcanonical or N V E ensemble where in addition the total momentum is conserved [3]. The temperature is not a control variable and cannot be preselected and fixed. But it is evident that also within molecular dy- namics the possibility to control the average temperature (as obtained from the average kinetic energy of the nuclei and the energy equipartition theorem) is welcome for phys- ical reasons. A deterministic algorithm of achieving temperature control in the spirit of extended system dynamics [18] by a sort of dynamical friction mechanism was devised by Nos´e and Hoover [16, 20], see e.g. Refs. [2, 16, 3] for reviews of this technique. Thereby, the canonical or N V T ensemble is generated in the case of ergodic dynamics. It is well–known that the standard Nos´e–Hoover thermostat method suffers from non– ergodicity problems for certain classes of Hamiltonians, such as the harmonic oscilla- tor [20]. A closely related technique, the so–called Nos´e–Hoover–chain thermostat [21], cures that problem and assures ergodic sampling of phase space even for the patholog- ical harmonic oscillator. This is achieved by thermostatting the original thermostat by another thermostat, which in turn is thermostatted and so on. In addition to restoring ergodicity even with only a few thermostats in the chain, this technique is found to be much more efficient in imposing the desired temperature. The underlying equations of motion read

MI R¨I = −∇I EKS^ − MI ξ˙ 1 R˙I (1.19)

Qn 1 ξ¨ 1 =

[ ∑

I

MI R˙^2 I − gkBT

] − Qn 1 ξ˙ 1 ξ˙ 2

Qn k ξ¨k =

[ Qn k− 1 ξ˙ k^2 − 1 − kBT

] − Qn k ξ˙k ξ˙k+1 (1 − δkK ) where k = 2,... , K.

By inspection of Eq. (1.19) it becomes intuitively clear how the thermostat works: ξ˙ 1 can be considered as a dynamical friction coefficient. The resulting ”dissipative dynamics” leads to non–Hamiltonian flow, but the friction term can aquire positive or negative sign according to its equation of motion. This leads to damping or acceleration of the nuclei and thus to cooling or heating if the instantaneous kinetic energy of the nuclei is higher or lower than kBT which is preset. As a result, this extended system dynamics can be shown to produce a canonical ensemble in the subspace of the nuclear coordinates and momenta. In spite of being non–Hamiltonian, Nos´e–Hoover (–chain) dynamics is also distinguished by conserving an energy quantity of the extended system; see Eq. (1.21). The desired average physical temperature is given by T and g denotes the number of dynamical degrees of freedom to which the nuclear thermostat chain is coupled (i.e. con- straints imposed on the nuclei have to be subtracted). It is found that this choice requires a very accurate integration of the resulting equations of motion (for instance by using a high–order Suzuki–Yoshida integrator [22]). The integration of these equations of motion is discussed in detail in Ref. [22] using the velocity Verlet algorithm. One of the advan- tages of the velocity Verlet integrator is that it can be easily used together with higher

Lecture 2

Ab initio Molecular Dynamics

In this lecture the two most popular extension of classical molecular dynamics to include first–principles derived potential functions are discussed. The focus is on the Kohn–Sham method of density functional theory [23, 24], as this is the method almost exclusively used in applications. However, most other electronic structure theories can be used as well and it is anticipated that thanks to more powerful computers and wider acceptance of dynamical methods within the quantum chemistry community in the future a larger variety of methods will be applied. The total ground–state energy of the interacting system of electrons with classical nuclei fixed at positions {RI } can be obtained

min Ψ 0 {〈Ψ 0 |He| Ψ 0 〉} = min {φi} EKS[{φi}]

as the minimum of the Kohn–Sham energy [25, 26]

EKS[{φi}] = Ts[{φi}] +

∫ dr Vext(r) n(r)

∫ dr VH(r) n(r) + Exc[n] + Eions(RN^ ) , (2.1)

which is an explicit functional of the set of auxiliary functions (Kohn–Sham orbitals) {φi(r)} that satisfy the orthonormality relation

〈φi | φj 〉 = δij. (2.2)

This is a dramatic simplification since the minimization with respect to all possible many– body wavefunctions {Ψ} is replaced by a minimization with respect to a set of orthonormal one–particle functions. The associated charge density

n(r) =

∑^ occ

i

fi | φi(r) |^2 (2.3)

is obtained from a single Slater determinant built from the occupied orbitals, where {fi} are integer occupation numbers. The first term in the Kohn–Sham functional Eq. (2.1) is the kinetic energy of a non– interacting reference system

Ts[{φi}] =

∑^ occ i

fi

〈 φi

∣∣ ∣∣−^1 2

∇^2

∣∣ ∣∣ φi

〉 (2.4)

consisting of the same number of electrons exposed to the same external potential as in the fully interacting system. The second term comes from the fixed external potential Vext(r), in most cases the potential due to the classical nuclei, in which the electrons move. The third term is the classical electrostatic energy of the electronic density and is obtained from the Hartree potential

VH(r) =

∫ dr′^

n(r′) | r − r′^ |

which in turn is related to the density through

∇^2 VH(r) = − 4 πn(r) (2.6)

Poisson’s equation. The second last contribution in the Kohn–Sham functional, is the exchange–correlation functional Exc[n]. The electronic exchange and correlation effects are lumped together and basically define this functional as the remainder between the exact energy and its Kohn–Sham decomposition in terms of the three previous contributions. Finally, the interaction energy of the bare nuclear charges is added in the last term. The minimum of the Kohn–Sham functional is obtained by varying the energy functional Eq. (2.1) for a fixed number of electrons with respect to the orbitals subject to the orthonormality constraint. This leads to the Kohn–Sham equations { −

∇^2 + Vext(r) + VH(r) +

δExc[n] δn(r)

} φi(r) =

∑

j

Λij φj (r) (2.7) { −

∇^2 + V KS(r)

} φi(r) =

∑

j

Λij φj (r) (2.8)

HKSφi(r) =

∑

j

Λij φj (r) , (2.9)

which are one–electron equations involving an effective one–particle Hamiltonian HKS with the local potential V KS. Note that HKS^ nevertheless embodies the electronic many– body effects by virtue of the exchange–correlation potential

δExc[n] δn(r)

= Vxc(r). (2.10)

V(:) := V(:) + dt/(2M(:))F(:) R(:) := R(:) + dtV(:) Optimize Kohn-Sham Orbitals (EKS) Calculate forces F(:) = dEKS/dR(:) V(:) := V(:) + dt/(2M(:))*F(:)

Extensions to other ensembles along the ideas outlined in the last lecture are straight forward. In fact, a classical molecular dynamics program can easily be turned into a BOMD program by replacing the energy and force routines by the corresponding routines from a quantum chemistry program.

2.1.1 Forces in BOMD

The forces needed in an implementation of BOMD are

d dRI

[ min {φi}

EKS[{φi}; RN^ ]

]

. (2.16)

They can be calculated from the extended energy functional

EKS^ = EKS^ +

∑

ij

Λij (〈φi | φj 〉 − δij ) (2.17)

to be

dEKS dRI

∂EKS

∂RI

∑

ij

Λij

∂RI

〈φi | φj 〉

∑

i

  ∂E

KS ∂〈φi |

∑

j

Λij | φj 〉

  ∂〈φi^ | ∂RI

The Kohn–Sham orbitals are assumed to be optimized, i.e. the term in brackets is (almost) zero and the forces simplify to

F KS(RI ) = −

∂EKS

∂RI

∑

ij

Λij

∂RI

〈φi | φj 〉. (2.19)

The accuracy of the forces used in BOMD depends linearly on the accuracy of the min- imization (see Fig. 2.1) of the Kohn–Sham energy. This is an important point we will further investigate when we compare BOMD to the Car–Parrinello method.

Figure 2.1: Accuracy of nuclear forces for a system of 8 silicon atoms in a cubic unit cell at 10 Ry cutoff using norm–conserving pseudopotentials

2.2 Car–Parrinello Molecular Dynamics

The basic idea of the Car–Parrinello approach can be viewed to exploit the time–scale separation of fast electronic and slow nuclear motion by transforming that into classical– mechanical adiabatic energy–scale separation in the framework of dynamical systems the- ory. In order to achieve this goal the two–component quantum / classical problem is mapped onto a two–component purely classical problem with two separate energy scales at the expense of loosing the explicit time–dependence of the quantum subsystem dynam- ics. This is achieved by considering the extended Kohn–Sham energy functional EKS^ to be dependent on {Φi} and RN^. In classical mechanics the force on the nuclei is obtained from the derivative of a Lagrangian with respect to the nuclear positions. This suggests that a functional derivative with respect to the orbitals, which are interpreted as classical fields, might yield the force on the orbitals, given a suitable Lagrangian. Car and Parrinello postulated the following Lagrangian [28] using EKS.

LCP[RN^ , R˙N^ , {Φi}, { Φ˙i}] =

∑

I

MI R˙^2 I +

∑

i

μ

〈 Φ˙i

∣∣ ∣ Φ˙i

〉 − EKS^

[ {Φi}, RN^

] (2.20).

The corresponding Newtonian equations of motion are obtained from the associated Euler–Lagrange equations

d dt

∂LCP

∂ R˙I

∂LCP

∂RI

d dt

δLCP δ〈 Φ˙i |

δLCP δ〈Φi |

like in classical mechanics, but here for both the nuclear positions and the orbitals. Note that the constraints contained in EKS^ are holonomic [9]. Following this route of ideas, Car–Parrinello equations of motion are found to be of the form

MI R¨I (t) = −

∂EKS

∂RI

∑

ij

Λij

∂RI

〈Φi | Φj 〉 (2.23)

μ Φ¨i(t) = −

δEKS δ〈Φi |

∑

j

Λij | Φj 〉 (2.24)

where μ is the ”fictitious mass” or inertia parameter assigned to the orbital degrees of freedom; the units of the mass parameter μ are energy times a squared time for reasons of dimensionality. Note that the constraints within EKS^ lead to ”constraint forces” in

shows that this frequency increases like the square root of the electronic energy difference Egap between the lowest unoccupied and the highest occupied orbital. On the other hand it increases similarly for a decreasing fictitious mass parameter μ. In order to guarantee the adiabatic separation, the frequency difference ωmine −ωmaxn should be large. But both the highest phonon frequency ωmaxn and the energy gap Egap are quantities that are dictated by the physics of the system. Therefore, the only parameter to control adiabatic separation is the fictitious mass. However, decreasing μ not only shifts the electronic spectrum upwards on the frequency scale, but also stretches the entire frequency spectrum according to Eq. (2.26). This leads to an increase of the maximum frequency according to

ωemax ∝

( Ecut μ

) 1 / 2 , (2.28)

where Ecut is the largest kinetic energy in an expansion of the wavefunction in terms of a plane wave basis set. At this place a limitation to decrease μ arbitrarily kicks in due to the maximum length of the molecular dynamics time step ∆tmax^ that can be used. The time step is inversely proportional to the highest frequency in the system, which is ωmaxe and thus the relation

∆tmax^ ∝

( (^) μ

Ecut

) 1 / 2 (2.29)

governs the largest time step that is possible. As a consequence, Car–Parrinello simulators have to make a compromise on the control parameter μ; typical values for large–gap systems are μ = 500–1500 a.u. together with a time step of about 5–10 a.u. (0.12– 0.24 fs). Note that a poor man’s way to keep the time step large and still increase μ in order to satisfy adiabaticity is to choose heavier nuclear masses. That depresses the largest phonon or vibrational frequency ωmaxn of the nuclei (at the cost of renormalizing all dynamical quantities in the sense of classical isotope effects). Other advanced techniques are discussed in the literature [22].

2.2.2 Forces in CPMD

The forces needed in a CPMD calculation are the partial derivative of the Kohn–Sham energy with respect to the independent variables, i.e. the nuclear positions and the Kohn–Sham orbitals. The orbital forces are calculated as the action of the Kohn–Sham Hamiltonian on the orbitals F (Φi) = −fiHKSφi. (2.30)

The forces with respect to the nuclear positions are

F (RI ) = −

∂EKS

∂RI

These are the same forces as in BOMD (the constraint force will be considered later), but there we derived the forces under the condition that the wavefunctions were optimized

and therefore they are only correct up to the accuracy achieved in the wavefunction optimization. In CPMD these are the correct forces and calculated from analytic energy expressions are correct to machine precision. Constraint forces are

Fc(Φi) =

∑ j

Λij | Φj 〉 (2.32)

Fc(RI ) =

∑

ij

Λij

∂RI

〈Φi | Φj 〉 , (2.33)

where the second force only appears for basis sets (or metrics) with a nuclear position dependent overlap of wavefunctions.

2.2.3 Velocity Verlet Equations for CPMD

The appearance of the constraint terms complicates the velocity Verlet method slightly for CPMD. The Lagrange parameters Λij have to be calculated to be consistent with the discretization method employed. How to include constraints in the velocity Verlet algorithm has been explained by Andersen [15]. In the following we will assume that the overlap is not position dependent, as it is the case for plane wave basis sets. The more general case will be explained when needed in a later lecture. These basic equations and generalizations thereof can be found in a series of papers by Tuckerman et al. [22]. For the case of overlap matrices that are not position dependent the constraint term only appears in the equations for the orbitals.

μ | Φ¨i(t)〉 =| ϕi(t)〉 +

∑

j

Λij | Φj (t)〉 (2.34)

where the definition | ϕi(t)〉 = −fiHKS^ | Φi(t)〉 (2.35)

has been used. The velocity Verlet scheme for the wavefunctions has to incorporate the constraints by using the rattle algorithm. The explicit formulas will be derived in the framework of plane waves in a later lecture. The structure of the algorithm for the wavefunctions is given below

CV(:) := CV(:) + dt/(2m)CF(:) CRP(:) := CR(:) + dtCV(:) Calculte Lagrange multiplier L CR(:) := CRP(:) + LCR(:) Calculate forces CF(:) = HKSCR(:) CV(:) := CV(:) + dt/(2m)CF(:) Calculte Lagrange multiplier L CV(:) := CV(:) + LCR(:)