Basics of

Molecular Dynamics Simulation

Submitted By:

Purnima Sharma,

Department of Bioinformatics

A Nobel Prize for Molecular Dynamics and QM/MM

What is Classical Molecular Dynamics?

Simulation of explicit particles (atoms, ions, .. . )

Particles interact via relatively simple analytical potential functions

Newton’s equations of motion are integrated for all particles simultaneously

100-1,000,000’s of particles depending on model Time 10 ps to 1 µs depending on model (typically ns)

Why Classical Molecular Dynamics?

Includes temperature and pressure e ects

Able to treat comparatively large systems for a relatively long time

Help interpret experiments, and to provide alternative interpretations

Test theories of intra- and intermolecular interaction Obtain detailed molecular level information

Structure and dynamics (compared to QC, MM, MC .. . )

Improve our understanding of nature through model-building (!!)

Phase Space and Time Averages

A system containing N atoms has 6 N values defining the state of the system (3 N coordinates ( r ), 3 N momenta ( p)) Each combination of coordinates and momenta define a point in the 6N-dimensional phase space : Γ _N

Consider a property that can be described as a function of the coordinates and momenta, A(Γ) , (for instance the total

A_obs = < A > _time= < A(Γ( t)) > _time= lim

1

^{τ ->∞} τ

energy or the (instantaneous) pressure)

The average value of the property can then be written as:

τ

t= 0

A(Γ( t)) dt

The ergodic hypothesis: “The long time average is equal to the ensemble average in the limit of an infinite number of members in the ensemble”

Molecular Dynamics

A MD simulation generates a sequence of points in phase space connected in time

The result is a trajectory of all particles in the system as a function of time

Time averages and other properties can be calculated from this trajectory

Motion of the system through phase space is governed by Hamiltonian equations of motion :

_r˙_i

=

∂ H ₌ p_i

∂ p_i m_i

_p˙

i

∂ r _i

= − ^{∂ H} = f _i

The Hamiltonian :

H( r , p)

= K ( p) + V ( r )

N

K ( p) =

i

_p2

i

2m_i

V is a potential energy function Newtons equations of motion :

_x˙

=

∂ H ₌ p_{x = vx}

∂ p_x m

p^˙ _x = − ^{∂ H} = − ^{∂ V ( x )} = F = ma_x

∂ x ∂ x

The Born-Oppenheimer approximation :

Electrons move much faster than the nuclei

QM: Solve Schrod^¨inger equation for static nuclei MD: Nuclei experience electrons as an average field

Steps in Performing an MD Simulation

Selection of interaction model

Selection of boundary conditions

Selection of initial conditions (positions, velocities .. . )

Selection of ensemble (NVE, NVT, NPT .. . )

Selection of target temperature, density/pressure . . .

Selection of integrator, thermostat, barostat . . .

Perform simulation until equilibration is reached (property dependent)

Perform production simulation to collect thermodynamic averages, positions, velocities

Analyze the results via post-processing

Initialize: Initial particle positions, velocities ...

Calculate forces from particle positions

Solve Newton’s eqns. of motion (integrate)

Loop over time- steps

Perform analysis / Write data to disk

The Interaction Model

The potential energy can be divided into interactions between pairs, triplets, .. . of particles :

V ( r ) = V_ij ( r _i , r _j )+ V_ijk ( r _i , r _j , r _k )+ V_ijkl ( r _i , r _j , r _k , r _l ) . . .

ij ijk ijkl

V(r_i,r_j)

V(r_i,r_k)

V(r_j,r_k)

i

k

j

V(r_i,r_j,r_k)

i

k

j

The interaction energy is (usually) largely dominated by the pair term. Truncate the series and make the assumption of

pairwise additivity :

V ( r ) ≈ V_ij ( r _i , r _j ) ij

In simulations the “true” pair potential is usually replaced by an effective pair potential, which incorporate the average many-body e ects into the first term :

ij

These potentials are usually parameterized to reproduce some experimental data

ij

V ( r ) ≈ V ^eff( r _i , r _j )

The Lennard-Jones Potential

^rij

V (kJ/mol)

Lennard-Jones interaction (Argon)

σ (3.41 Å)

ε (0.996 kJ/mol)

5

r (Å)

-1

​

-1.5

-0.5

0

1

​

0.5

1.5

2

1

2

3

4

6

7

8

9

Two dimensional simulation of argon

The Initial Configuration

For a solid it may be possible to start from the experimentally determined structure

For a liquid it is also common to start from a “solid” lattice, and letting the structure melt

For a liquid it is also possible to start from a “random” placement of molecules (this may be tricky)

If previous simulation data is available, one might modify a configuration to suit the new simulation, for instance by replacing one or a few particles

Macromolecule in liquid: Use experimental coordinates of macromolecule from solid phase structure determination, surround it by liquid

The Initial Velocities

The velocities of the particles is related to the kinetic energy :

K =

N

i

m_iv ²

i

2

The kinetic energy is related to the temperature :

2 < K > = gk_B T = 3Nk_B T

Assign velocities according to a Boltzmann distribution :

P( v_i ) =

m

i

2πk_B T

1/ 2

e

−

i

m v

2

i

k_B T

Alternatively assign particle velocities from a uniform distribution. They will obtain a Boltzmann distribution after a short time.

It is often necessary to adjust the initial velocities so the total linear momentum of the system is zero :

N N

p = p_i = m_iv _i = 0

i i

This momentum corresponds to translation of the whole system

Having p = 0 in simulations of “liquid” water, has been referred to as a “flying ice cube”

Newton’s equations of motion will preserve the linear momentum resulting in a NVE P ensemble, resulting in a system with three fewer degrees of freedom

For a system with periodic boundary conditions, the

angular momentum will not be conserved

Statistical Mechanical Ensembles

The “natural” ensemble for MD is the microcanonical ensemble: NVE (constant number of particles, constant volume, constant energy)

Typical experiments: Constant pressure and/or temperature

Canonical ensemble: NVT (constant number of particles, volume, temperature); Isobaric, isothermal (Gibbs): NPT

(constant number of particles, pressure, temperature) Some ways to include such e ects in MD

Stochastic events (change kinetic energy of single particle)

^{Periodic rescaling of velocities (2} K = gk_B T )

Ad hoc partial rescaling of velocities (Berendsen’s thermal bath)

Extended system methods: Include extra variables in the equations of motion such as the volume of the system

Periodic Boundary Conditions

A “real” system might contain about N_A = 6 · 10²³ molecules Simulation may contain about 100-1,000,000 molecules

A small cube of 10x10x10 molecules has about half of the molecules on the surface

Not appropriate to study “bulk” properties Periodic boundary conditions:

Small box replicated in all directions

A particle that leaves the box on one side is replaced by an image particle that enters from the other side

There are no walls and no surface particles

Some systems inherently contains a boundary, for instance a liquid droplet, or a surface

The Minimum Image Convention

Particles interact only with the closest periodic image of the other particles:

How to compute minimum image :

Simple and straightforward

rx = x i − x j ;

i f ( rx > 0 . 5 * lengt h ) rx = rx − l eng t h ;

i f ( rx < − 0.5 * lengt h ) rx = rx + l eng t h ;

​

​

Fortran has a function , anint, which returns the nearest integer

rx = rx - length* anint (rx/length)

In C, the floor function can be used

rx = rx − length*floor(rx/length + 0.5);

Truncating the Potential

The most time consuming part: Computing the nonbonded (vdw & electrostatic) energies and forces

For a pair additive potential there are N ² such interactions

Consider the Lennard-Jones argon potential again:

-1₀

0

-0.1

-0.2

-0.3

-0.4

-0.5

-0.6

-0.7

-0.8

-0.9

1 2 3 4 5 6 7 8 9

r (Å)

V (kJ/mol)

Spherical nonbonded c

ij

ij

ij

V_ij ( r_ij ) =

ij ij

0

^rij

V ( r ) = V ( r ) r ≤ r_c

> r_c

Disadvantage: Discontinuity in the energy and forces

Switching function ( S( r_ij ) smoothly from 1 to 0)

^rij ^{< r}l

r_l ≤ r_ij ≤ r_c

V_ij ( r_ij ) = V_ij ( r_ij ) = V_ij ( r_ij ) =

V_ij ( r_ij )

S ( r_ij ) V_ij ( r_ij )

0

^rij

> r_c

Minimum Image and Spherical Cutoff

Very commonly minimum images and spherical cuto s are used together

r_c

Searching for Neighbors

The search loop is a double loop over all particles ( n). Execution time for the loop is proportional to n²

The Verlet neighbor list : Maintain list of neighbor pairs

closer than the radiu r_c

^rskin

r_c

^rc^+rskin

∆t

r_c

^rc^+rskin

Update list from time to time ( about once every 25 timesteps, but this can be automated)

Use of Linked Lists (Cell Method)

For large systems the execution time will be dominated by the search loop (because building the Verlet neighbor lists are still O( N ² )

One solution is to divide the space into cells

a

}

x

_} b

y

Cell (0,0)

Cell (7,2)

Finding which cell a particle belongs to is an O ( N ) problem. A particle with coordinate ( x , y ) belongs to cell

( floor (x/a ) , floor (y/b ) )

All particles are sorted into so-called linked lists . There is one list for each cell .

When all particles are sorted into lists , the double loop can be done over cells instead of particles. We don’t have to consider all cells . Only cells with particles closer to each other than the c distance must be considered :

For the particles inside the central cell only particles inside the central cell itself and the shadowed cells are considered if the cuto distance is equal to, for instance, 1.2 cell lengths

Domain Decomposition

A 3x3 decomposition of a two dimensional system

For very large systems the memory requirement on each node becomes large, when using a replicated data algorithm

Also the program will become slow due to that forces from all particles must be sent to all processors each timestep

​

One solution for very large systems is to use domain decomposition . Particles are assigned to processors depending on their position in space.

When domain decomposition is used larger systems may be simulated simply by using more processors

For large systems the interaction cuto radius, r_c , is substantially smaller than the size of the system

r_c

Particles closer to the “edge” of the domain

than the cuto radius, interact with particles

on the node.

UP

DOWN

DOWN/ LEFT

LEFT

UP/ LEFT

UP/ RIGHT

RIGHT

DOWN/ RIGHT

8 neighbours in 2D , and 26 in 3D

Integrating the Equations of Motion

Motion of all particles is coupled (many-body problem), This is impossible to solve analytically . . .

Solve using finite di erences

Positions, velocities ... approximated as a Taylor series:

r ( t + ∆ t)

= r ( t) + ∆ tv ( t) + ¹ ∆ t² a ( t) + ¹ ∆ t³ b( t) + ¹ ∆ t⁴ c ( t) + . . .

2 6 24

v ( t + ∆ t)

= v ( t) + ∆ ta ( t) + ¹ ∆ t² b( t) + ¹ ∆ t³ c ( t) + . . .

2 6

a ( t + ∆ t)

2

= a ( t) + ∆ tb( t) + ¹ ∆ t² c ( t) + . . .

b( t + ∆ t)

= b( t) + ∆ tc ( t) + . . .

The simplest method truncate the series after the first-order term. This gives Euler’s method :

This is a first order method (and should not be used) The Verlet algorithm (1968, Loop Verlet) :

r ( t + ∆ t)

2

= r ( t) + ∆ tv ( t) + ¹ ∆ t² a ( t)

2

r ( t − ∆ t) = r ( t) − ∆ tv ( t) + ¹ ∆ t² a ( t)

Add the equations, and rearrange :

r ( t + ∆ t) = 2r ( t) − r ( t − ∆ t) +

∆ t² a ( t)

The Verlet algorithm is a second order algorithm Drawbacks of the Verlet algorithm :

Velocities do not appear as they have been eliminated by the addition. They must be estimated, for instance by taking another integration step :

v ( t) =

r ( t + ∆ t) − r ( t − ∆ t)

2 ∆ t

Velocities are needed, for instance, when computing the kinetic energy of the system

It is not a self-starting algorithm

A small term ( ∆ t² a( t) ) is added to the di erence of two large terms (2 r ( t) − r ( t − ∆ t) )

There are several variations of the Verlet method :

The Leap-Frog Verlet method :

2

r ( t + ∆ t ) = r ( t) + ∆ t v ( t + ¹ ∆ t)

v ( t + ¹ ∆ t) = v ( t − ¹ ∆ t) + ∆ t a ( t)

2 2

Disadvantage: Velocities and positions not available at the same time

​

The Velocity Verlet method :

2

r ( t + ∆ t ) = r ( t) + ∆ t v( t) + ¹ ∆ t² a( t)

v ( t + ∆ t ) =

2

v ( t) + ^{∆ t} _[a ( t) + a ( t + ∆ t)_]

Solved in two steps with a force evaluation in between Yields the positions, velocities and accelerations at the same time

(a) 1-dt

I t+dt

,-cSt

I

t+cSt ^t-cSt

I

r

V

a

(b)

r

V

a

(c) r

V

a

​

Fi&- 3.2 Various forms of the Verlet algorithm. (a) Verlet's original method. (b) The lcap­ fr_og form. (c) The velocity form. We show successive steps in the implementation or each algorithm. lo each case, the stored variables are in grey boxes.

MOLECULAR DYNAMICS

I t+dt

t-6t

I t+dt

Apart from numerical di erences, the three Verlet methods result in the same trajectories

The Gear predictor-corrector integrator :

First new positions, velocities and accelerations are predicted according to the Taylor expansions :

r ^p( t + ∆ t )

2

= r ( t) + ∆ t v ( t) + ¹ ∆t²

a ( t) +

6

¹ ∆ t³ b( t) +

24

¹ ∆ t⁴

c ( t)

v ^p( t + ∆ t )

2

= v ( t) + ∆ t a ( t) + ¹ ∆t²

b( t) +

6

¹ ∆ t³ c ( t)

a ^p( t + ∆ t )

2

= a ( t) + ∆ t b( t) + ¹ ∆ t² c ( t)

b^p( t + ∆ t ) = b( t) + ∆ t c ( t)

​

Then the forces are evaluated to give the true accelerations

a ( t + ∆ t)

The Gear predictor-corrector integrator (cont) :

The di erence between the predicted and calculated accelerations is computed :

∆ a ( t + ∆ t) = a ( t + ∆ t) − a ^p ( t + ∆ t)

The correction step :

The values of the coefficients c₀ , c₁ ,. . . depend on the

120

19 _{, c1 =}

4

2

³ , c₂ = 1, c₃ = ¹ , c₄ =

1

12

order of the method ( ^c₀ ^{= )}

Choosing the Timestep

A balance is needed in choosing the timestep. A too short timestep lead to a very expensive solution of the equations of motion, and leads to a limited coverage of phase space. A much too small timestep might lead to

round-o errors.

Too large of a timestep leads to instabilities

​

Example: Simulations of the 2D Lennard-Jones Argon System

-356.032

0 1 2 3 4 5 6 7 8 9 10

Time (ps)

Total energy as a function of time

-356.031

-356.03

-356.029

Energy (kJ/mol)

“Exact” trajectory of one argon atom

Velocity Verlet Integrator

∆ t = 1 fs (blue)

-356.032

0 1 2 3 4 5 6 7 8 9 10

Time (ps)

-356.031

-356.03

-356.029

Energy (kJ/mol)

Velocity Verlet Integrator

∆ t = 10 fs (blue)

-355.99

-356

-356.01

-356.02

-356.03

-356.04

-356.05

-356.06

-356.07

-356.08

-356.09

0 1 2 3 4 5 6 7 8 9 10

Time (ps)

Energy (kJ/mol)

Gear Predictor-Corrector Method

∆ t = 1 fs (blue)

-356.032

0 1 2 3 4 5 6 7 8 9 10

-356.031

-356.03

-356.029

Energy (kJ/mol)

Time (ps)

Gear Predictor-Corrector Method

∆ t = 10 fs (blue)

-347

-348

-349

-350

-351

-352

-353

-354

-355

-356

-357

0 1 2 3 4 5 6 7 8 9 10

Time (ps)

Energy (kJ/mol)

Comparison of Velocity Verlet (red) and the Gear Predictor-Corrector (blue) at ∆ t = 1 fs :

-356.022

-356.023

-356.024

-356.025

-356.026

-356.027

-356.028

-356.029

-356.03

-356.031

-356.032

0 100 200 300 400 500 600 700 800 9001000

Time (ps)

Energy (kJ/mol)

The Gear method has short time fluctuations, but it has a (substantial) long-time drift

The Verlet method uses a central di erence and is thus time-reversible

The trajectories will eventually always diverge from the “true” trajectory

The energy (Hamiltonian) of the system should be conserved Use a shorter timestep for :

Lighter particles

Higher temperatures Sti er potentials

_-:::::i

-255.280

-255.285

:::::-

⁰_e -255.290

-255.295

;;,,-.,

cu

-255.300

-255.305

C

G,)

"c;

...

-255.310

, ..

-255.315O

20

40

60

80

100

Time (ps)

90

-300

KE

80

_{-310 0},--...

----0

e

e

:::::i

:::::i

u

u

70

-320

-

-

>.

_e;;,_.,-_o.,

00

...

G,)

C

G,)

G,) 60

C

-330

G,)

"c;

·^u-

·.-..

C:

G,)

C:

G,)

...

0

50

-340

Q...

40

-350

0

20

40

60

80

100

Time (ps)

​

​

Fig. 6.3 Variation in total energy versus time for the production phase of a molecular dynamics simulation of 256 argon atoms at a temperature of 1OOK and a density of 1.396 gcm-³(top). The time step was lOfs and the equations of motion were integrated using the velocity Verlet algorithm. The variations in the kinetic and potential energies are also shown. The graphs have different scales.

)

/

I

)

)

10° r----.--.------,---

10-

1

^10° r----.---r---.---

10-

2

-

C

0

:.:;

_(lj 10-

3

>

Q)

Q)

"O

«i 10-⁴

C

"O

"-

5

0

8 10-

en

cc

J

10-⁶

J,

I

10-⁷ I

I

I

I

I

I

I

I

I

r

.I I

I

I

I

I

I

I

I

I

I

I

I

I

I

I

J

I

I

I

I

I

\.,,- ,._J\oa.•l"Ji

,,,.^Ii

I

I

I

,

't = 0.1 fs

e = + 10-⁴A E = + 10-⁶A E = ± 10-SA

^10-S0,- --:-:;;----'---

5

10 15 20

Time (ps)

10-

1

10-

2

-_C

0

:.:; -3

(lj 10

">

Q)

"O

.s -4

(lj 10

C

=o

"-

0

0

-5

o 10

en

cc

r...,.,,1-·-

,r -"",

₍ ( ,.., _,,,..., ...... - _.'.I; ,,.-

I

I

I ,

I

• I

I

I

I

I

I

I

I

J

I

,.

I

I

_.I

I

I

I

I

I

I

I

I

I

I

r

I

10-⁷ r

I

I

I

I

I

J

't = 1 fs

e = + 10-⁴A E = ± 10-⁶A e = + 10-⁸A

^10-Bol ---:;;--- ,L -

1

20

5

10 15

Time (ps)

Figure 12.3. The rapid divergence of dynamic trajectories for four water molecules differ­ ing slightly in initial conditions.