Introduction to MD

Download Report

Transcript Introduction to MD

Molecular Dynamics Simulations

An Introduction

Pingwen Zhang

Molecular Dynamics

Definitions, MotivationsForce fieldsAlgorithms and computationsAnalysis of Data

Molecular dynamics - Introduction

• Molecular dynamics (MD) simulation technique: is a computer the time evolution of interacting atoms is followed by integrating their equations of motion.

• We follow the laws of classical mechanics, and most notably Newton's law: F i = m i a i a i = d 2 r i =dt 2

Scale in Simulations

mesoscale continuum Monte Carlo 10 -6 S 10 -8 S 10 -12 S molecular dynamics quantum chemistry exp( D E/kT) F = MA Hy = Ey 10 -10 M 10 -8 M 10 -6 M Length Scale domain 10 -4 M

Why Not Quantum Mechanics?

• Modeling the motion of a complex molecule by solving the wave functions of the various subatomic particles would be accurate… • But it would also be

very

hard to program and take more computing power than anyone has!

Molecular dynamics - Introduction

• Given an initial set of positions and velocities, the subsequent time evolution is determined.

in principle completely • Atoms and molecules will ‘move’ in the computer, bumping into each other, vibrating about a mean position (if constrained), or wandering around (if the system is fluid), oscillating in waves in concert with their neighbours, perhaps evaporating away from the system if there is a free surface, and so on, in a way similar to what real atoms and molecules would do.

Molecular dynamics -Motivation

• The computer experiment .

• In a computer experiment, a model is still provided by theorists, but the calculations are carried out by the machine by following a recipe (the implemented in a suitable algorithm , programming language ). • In this way, complexity can be introduced (with caution!) and more realistic systems can be investigated, opening a road towards a better understanding of real experiments.

Molecular dynamics -Motivation

• The computer calculates a trajectory of the system • 6 N -dimensional phase space (3 momenta ). N positions and 3 N • A trajectory obtained by molecular dynamics provides a set of conformations of the molecule, • They are accessible without any great expenditure of energy (e.g. breaking bonds) • MD also used as an efficient tool for optimisation of structures ( simulated annealing ).

Molecular dynamics - Motivation

• MD allows to study the dynamics of large macromolecules • Dynamical events control processes which affect functional properties of the bio molecule (e.g. protein folding).

• Drug design is used in the pharmaceutical industry to test properties of a molecule at the computer without the need to synthesize it.

Molecular dynamics – Time Limitations

• Typical MD simulations are performed on systems containing millions of atoms • Simulation times: picoseconds to nanoseconds.

• A simulation is reliable when the simulation time is much longer than the relaxation time of the quantities we are interested in.

Historical Perspective on MD

Procedure of MD

Initialization

• Position – X-ray, NMR, simulation or analytical calculation • Velocity – The initial velocities are assigned taking them from a Maxwell distribution at a certain temperature

T

• Another possibility is to take the initial positions and velocities to be the final positions and velocities of a previous MD run

Choosing the Time Step

d

t

No MD follows the true trajectories for very many time steps – errors always accumulate

.

The time over which the MD trajectory = true trajectory is called

correlation time

.

No MD truly conserves energy since there are always errors

.

The goal is to have a constant average E with fluctuations as small as possible • Time step d t should be as large as possible to still get accurate trajectories (on the time scale needed) and conserve of energy • In general, d t should be ≈0.01 x the fastest behavior of your system (E.g., atoms oscillate about once every 10 -12 s in a solid  MD time steps are ≈ 10 -14 s in simulations of solids

Potentials

• ab initio potential – Quantum calculation, DFT, OFDFT, Tight-binding, etc.

• Empirical potential – Comes from quantum – Consistent with continuum – Fit the database: elastic moduli, surface energy, etc.

• Connection with continuum – Constitutive relation

The Lennard-Jones Potential

•The Lennard-Jones potential

u LJ

= 4     

r

12  

r

 6    •The truncated Lennard-Jones potential

u

= 

u

 

LJ

0

r r

r c

r c

•The truncated and shifted Lennard-Jones potential

u

= 

u

 

LJ

u LJ

 

c

0

r r

r c

r c

FENE potential

• FENE stands for: Finitely extensible nonlinear elastic E = 4²[( ¾ ) 12 r ¡ ( ¾ ) 6 ] ¡ 0:5· ln[1 ¡ ( r r R o ) 2 ]

EAM potential

• Embedded Atom Method works for metallic solids • Two contributions • nuclear-nuclear interaction • embedding an atom to the electron cloud E i = F ® 0 @ X 1 ½ ® r i j A + 1 X 2 © ®¯ r i j

Potential for Covalent Carbon

• The Stillinger-Weber potential • The Tersoff Potential Not only accounts for the contribution of bond lengths, but also for the bond angles

Non-Bonded Atoms

There are two potential functions we need to be concerned about between non-bonded atoms: • van der Waals Potential • Electrostatic Potential

The van der Waals Potential

• Atoms with no net electrostatic charge will still tend to attract each other at short distances • Atoms tend to repel when they get too close The Constants A and C depend on the atom types, and are derived from experimental data

The Electrostatic Potential: Coulomb’s Law

• Opposite Charges Attract • Like Charges Repel • The force of the attraction is inversely proportional to the square of the distance

Bonded Atoms

There are three types of interaction between bonded atoms: • Stretching along the bond • Bending between bonds • Rotating around bonds

Bond Length Potentials

Both the spring constant and the ideal bond length are dependent on the atoms involved.

Bond Angle Potentials

The spring constant and the ideal angle are also dependent on the chemical type of the atoms.

Torsional Potentials

Described by a dihedral angle and coefficient of symmetry (

n=1,2,3

), around the middle bond.

Effects of solvents

• Implicit models – “Generalized Born” solvent Model coarse-graining the effects of solvent by approximately solving the Poisson equation • Explicit models – Explicitly adding the water molecules which are regarded as rigid bodies

Integrator: Verlet Algorithm

Start with {r(t), v(t)}, integrate it to {r(t+

D

t), v(t+

D

t)}:

{r(t+ D t), v(t+ D t)} The new position at t+ D t:

r

(

t

 D

t

) =

r

(

t

) 

v

(

t

) D

t

 1 2 D

t

2

a

(

t

) 

O

( D

t

3 ) Similarly, the old position at t D t:

r

(

t

 D

t

) =

r

(

t

) 

v

(

t

) D

t

 1 2 D

t

2

a

(

t

) 

O

( D

t

3 ) Add (1) and (2):

r

(

t

 D

t

) = 2

r

(

t

) 

r

(

t

 D

t

)  D

t

2

a

(

t

) 

O

( D

t

4 ) {r(t), v(t)} (1) (2) (3) Thus the velocity at t is:

v

(

t

) =

r

 (

t

) = 1 2 D

t

(

r

(

t

 D

t

) 

r

(

t

 D

t

)) 

O

( D

t

2 ) (4)

Verlet Scheme

• Is time reversible • Does conserve volume in phase space • (Is symplectic) • Doe not suffer from energy drift

Velocity Verlet scheme

r

(

v

(

t t

 D

t

 D

t

) ) = =

r

(

t

)

v

(

t

) 

v

(

t

) D

t

f

(

t

)  

f

1 2 (

t

D

t

2 

f

( D

t

)

m t

) 2

m

• Velocity calculated explicitly • Possible to control the temperature • Stable in long time simulation • Most commonly used algorithm

Periodic Boundary Conditions Minimum Image

r c

Central Simulation box

Periodic Boundary Conditions Ewald Method Ewald Sum: split into two teems, real space sum and k-space sum

Cell list

Saving CPU time

Verlet list

Ensembles

• NVE – micro-canonical ensemble • NVT – canonical ensemble • NPT – grand-canonical ensemble • Temperature control – Berendsen thermostat (velocity rescaling) – Andersen thermostat (velocity resampling) – Nose-Hoover chain • Pressure control – Berendsen volume rescaling – Andersen piston 

NVT

NPT

= =  exp[  

H

(

p

,

q

)] 

dpdq

exp[  

H

(

p

,

q

)] exp[   (

H

(

p

,

q

) 

PV

)]

dpdq

exp[   (

H

(

p

,

q

) 

PV

)] 

T

= 1  D 

T

(

T B T

 1 ) 

P

= 3 1   D 

P

(

P B

P

)

MD as Optimization tool

Simulated Annealing

Cooling Schedules: • • • • Most popular global optimization algorithm Start at high T, decrease T in small steps (cooling schedule) Easy to understand & implement Drawback: might be easily trapped in local minima

Molecular dynamics – Analyses

• The simplest way of analyzing the system during (or after) its dynamic motion is looking at it.

• One can assign a radius to the atoms, represent the atoms as balls having that radius, and have a computer program construct a ‘photograph’ of the system.

• We may also color the atoms according to its properties (charge, displacement, ‘temperature’…)

Molecular dynamics – Analyses

• We also can measure instantaneous and time averages of various physically important quantities • To measure time averages: If the instantaneous values of some property A at time t is A(t) = f (r 1 (t); : : : ; r N (t); v1(t); : : : ; v N (t)) then its average is where N T T 1 < A > = A(t) N T t = 1 is the number of steps in the trajectory

Softwares

• AMBER • CHARMM • VASP (DFT) • XMD • CPMD(DFT) (Car-Parrinello MD)

References

• M. P. Allen, D. J. Tildesley (1989) Computer simulation of liquids. Oxford University Press. • Frenkel Daan; Smit, Berend [2001].

Understanding Molecular Simulation : from algorithms to applications

. Academic Press. D. C.

• Rapaport (1996) The Art of Molecular Dynamics Simulation. • Tamar Schlick (2002) Molecular Modeling and Simulation. Springer.