introduction to Ab Initio Molecular Dynamics

introduction to Ab Initio Molecular Dynamics

(Parte 1 de 7)

Lecture Notes

Introduction to Ab Initio Molecular Dynamics

Jurg Hutter

Physical Chemistry Institute

University of Zurich

Winterthurerstrasse 190 8057 Zurich, Switzerland


1.1 Introduction4
1.2 Equations of Motion4
1.3 Microcanonical Ensemble5
1.4 Numerical Integration6
1.5 Extended System Approach8
1.5.1 Barostats8
1.5.2 Thermostats10

1 Molecular Dynamics 4

2.1 Born–Oppenheimer Molecular Dynamics14
2.1.1 Forces in BOMD15
2.2 Car–Parrinello Molecular Dynamics16
2.2.1 How to Control Adiabaticity ?18
2.2.2 Forces in CPMD19
2.2.3 Velocity Verlet Equations for CPMD19
2.3 Comparing BOMD and CPMD20

2 Ab initio Molecular Dynamics 12

3.1 Unit Cell and Plane Wave Basis24
3.2 Kinetic Energy and Local Potentials26
3.3 Electrostatic Energy27
3.4 Exchange and Correlation Energy29
3.5 Car–Parrinello Equations30

3 Plane Waves 24

4.1 Why Pseudopotentials ?3
4.2 Norm–Conserving Pseudopotentials35
4.2.1 Hamann–Schluter–Chiang Conditions35
4.2.2 Bachelet-Hamann-Schluter (BHS) form37
4.2.3 Kerker Pseudopotentials37

4 Pseudopotentials 3 1

4.2.5 Kinetic Energy Optimized Pseudopotentials38
4.3 Pseudopotentials in the Plane Wave Basis39
4.3.1 Gauss–Hermit Integration40
4.3.2 Kleinman–Bylander Scheme41
4.4 Dual–Space Gaussian Pseudopotentials43
4.5 Example: Pseudopotentials for Oxygen4
5.1 Total Energy and Gradients47
5.1.1 Plane Wave Expansion47
5.1.2 Total Energy47
5.1.3 Wavefunction Gradient48
5.1.4 Nuclear Gradient48
5.2 Fast Fourier Transforms49
5.3 Density and Force Calculations in Practice51
5.4 Saving Computer Time51
5.5 Exchange and Correlation Functionals53

5 Implementation 47

6.1 Non-linear Core Correction5
Augmented–Wave Method56
6.2.1 The PAW Transformation57
6.2.2 Expectation Values59
6.2.3 Ultrasoft Pseudopotentials61
6.2.4 PAW Energy Expression63
6.2.5 Integrating the CP Equations64
6.3 Metals; Free Energy Functional67
6.4 Charged Systems71

6 Extensions to the Original Method 5 6.2 Ultrasoft Pseudopotentials and the Projector

7.1 Position Operator in Periodic Systems75
7.2 Dipole Moments and IR Spectra76
7.3 Localized Orbitals, Wannier Functions78

7 Properties from CPMD 75

8.1 Simulated Annealing83
8.2 Rare Events84
8.2.1 Thermodynamic Integration85
8.2.2 Adiabatic Free Energy Sampling8
8.2.3 Bias Potential Methods89

8 Chemical Reactions with CPMD 83 2

9.1 Adiabatic Density–Functional Perturbation Theory93
9.2 Coupled Perturbed Kohn–Sham Equations93
9.2.1 Exchange-Correlation Functionals96
9.3 Nuclear Hessian97
9.3.1 Selected Eigenmodes of the Hessian98
9.4 Polarizability9
9.5 NMR Chemical Shifts101
9.5.1 Chemical shifts and susceptibilities101
9.5.2 The gauge origin problem103
9.5.3 The position operator problem105
9.5.4 Density functional perturbation theory107
9.5.5 Pseudopotential correction108
10.1 Restricted Open–shell Kohn–Sham Method109
10.2 Time–Dependent Density Functional Theory112
10.3 Excited States from Linear Response113
10.3.1 Tamm-Dancoff Approximation118

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ν > kBT. 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 = MV. The union of all positions (or momenta) {R1,R2,...,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

The forces on the particle are derived from the potential

The equations of motion are according to Hamilton’s equation

from which we get Newton’s second law (using PI = MI RI)

The equations of motion can also be derived using the Lagrange formalism. The Lagrange function is

and the associated Euler–Lagrange equation ddt ∂L

leads to the same final result. The two formulations are equivalent, but the ab initio molecular dynamics literature almost exclusively uses the Lagrangian techniques.

1.3 Microcanonical Ensemble

The following section is slightly shortened from the very clear and concise feature article by Tuckerman and Martyna [10]. The equations of motion are time reversible (invariant to the transformation t → −t) and the total energy is a constant of motion

These properties are important to establish a link between molecular dynamics and statistical mechanics. The latter connects the microscopic details of a system the physical observables such as equilibrium thermodynamic properties, transport coefficients, and spectra. Statistical mechanics is based on the Gibbs’ ensemble concept. That is, many individual microscopic configurations of a very large system lead to the same macroscopic 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 NV E ensemble. Other examples include the canonical or NV T ensemble, the isothermal–isobaric or NPT 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, RN) = 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 [1, 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 system 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 Differential 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 hypersurface 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:

This equation is combined with the update for the velocities

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

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

• The diffusion constant

for large times τ. • The pair correlation function

• The temporal Fourier transform of the velocity autocorrelation function 〈V(t)·V(0)〉 is proportional to the density of normal modes (in a purely harmonic system).

1.5 Extended System Approach

In the framework of statistical mechanics all ensembles can be formally obtained from the microcanonical NV E ensemble – where particle number, volume and energy are the external thermodynamic control variables – by suitable Laplace transforms of its partition function. Thermodynamically this corresponds to Legendre transforms of the associated thermodynamic potentials where intensive and extensive conjugate variables are interchanged. In thermodynamics, this task is achieved by a ”sufficiently weak” coupling of the original system to an appropriate infinitely large bath or reservoir via a link that establishes thermodynamic equilibrium. The same basic idea is instrumental in generating distribution functions of such ensembles by computer simulation. Additional degrees of freedom that control the quantity under consideration are added to the system. The simulation is then performed in the extended microcanonical ensemble with a modified total energy as a constant of motion. This system has the property that after the correct integration over the additional degrees of freedom has been performed the distribution function of the targeted ensemble is recovered. Two important special cases are: thermostats and barostats, which are used to impose temperature instead of energy and / or pressure instead of volume as external control parameters [2, 3, 7, 16, 17, 18, 19].

Keeping the pressure constant is a desirable feature for many applications of molecular dynamics. The concept of barostats and thus constant–pressure molecular dynamics was introduced in the framework of extended system dynamics by Andersen [18]. His method was devised to allow for isotropic fluctuations in the volume of the supercell. An extension of Andersen’s method consists in allowing for changes of the shape of the computational cell to occur as a result of applying external pressure [17], including the 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 a1, a2 and a3 of the simulation cell are additional dynamical variables. Using the 3 × 3 matrix h = [a1,a2,a3] (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 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

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 vu −MI

( Πtotus − p δus where the total internal stress tensor

Πtotus = 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].

Standard molecular dynamics generates the microcanonical or NV 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 dynamics 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 physical 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 Nose and Hoover [16, 20], see e.g. Refs. [2, 16, 3] for reviews of this technique. Thereby, the canonical or NV T ensemble is generated in the case of ergodic dynamics. It is well–known that the standard Nose–Hoover thermostat method suffers from non– ergodicity problems for certain classes of Hamiltonians, such as the harmonic oscillator [20]. A closely related technique, the so–called Nose–Hoover–chain thermostat [21], cures that problem and assures ergodic sampling of phase space even for the pathological 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

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, Nose–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. constraints 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 [2]). The integration of these equations of motion is discussed in detail in Ref. [2] using the velocity Verlet algorithm. One of the advantages of the velocity Verlet integrator is that it can be easily used together with higher order schemes for the thermostats. Multiple time step techniques can be used and time reversibility of the overall algorithm is preserved.

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

(Parte 1 de 7)