History, Development and Basics of Molecular Dynamics Simulation - - PowerPoint PPT Presentation
History, Development and Basics of Molecular Dynamics Simulation Technique A way to do experiments in computers (Odourless) Creating virtual matter and then studying them in computers S. Yashonath Solid State and Structural Chemistry Unit
History, Development and Basics of Molecular Dynamics Simulation Technique
A way to do experiments in computers (Odourless) Creating virtual matter and then studying them in computers
Solid State and Structural Chemistry Unit Indian Institute of Science Bangalore 560 012
National Workshop on Atomistic Simulation Techniques, IIT, Guahati
Dedication
Hyderabad, India
Acknowledgements
All my students :
Sanjoy Bandyopadhyay (IIT, Kgp),
A.V. Anil Kumar(NISER, Bhubhaneswar), S.Y. Bhide(Dow Chemicals, Mumbai), C.R. Kamala(GE, USA),
Manju Singh (Bangalore), B.J. Borah (private univ., Baroda),
Plan of the talk : Introduction Classical mechanics : basics Statistical mechanics : basics Intermolecular potentials Numerical integration Algorithms Computational Tricks Trajectory Equilibrium and time dependent properties Microcanonical ensemble Other ensembles
What is molecular dynamics ?
A way of solving equations of motion numerically (hence you need a computer). Equations of motion are coupled differential equations and hence can not be solved analytically. You can get all properties about the system being simulated (from statistical mechanical relationships)
Subjects involved in MD
History and evolution of MD
Rosenbluth, Rosenbluth, Teller & Teller –simulation of a dense liquid of 2D spheres
–simulation of anharmonic 1D crystal
–molecular dynamics (MD) simulation of hard spheres (1958: First X‐ray structure of a protein)
– Simulation of damaged Cu crystal
–MD simulation of liquid Ar
–Monte Carlo simulation of water
–MD simulation of water * 1972: I.R. McDonald
Aneesur Rahman in his younger days
Post 1980 developments : the extended Hamiltonian methods
structure transformation with corrections from S. Yashonath
Ab initio MD (includes electronic degrees of freedom)
Alder & Wainright Hard sphere 1956 Metropolis et al Monte Carlo 1953 Fermi, Pasta, Ulam 1D crystal 1955 Simulation of damaged Cu crystal Vineyard 1960 Rahman MD of argon 1964 Barker and Watts MC of water 1969 Rahman and Stillinger MD simulation of water, 1971 NPT of argon McDonald 1972
A chart of development of simulation methods
Feynman vs. Einstein
in terms of jiggling and wiggling of atoms.” R. Feynman.
Possible, But Not Simpler : A. Einstein
MD jargon : Terms often used
single atom
interactions
Simple microcanonical ensemble MD
How does one get force Fi ? From the intermolecular potential : = grad of potential We need : or potential φi.
A knowledge of intermolecular interaction potential is therefore central to molecular dynamics. Without it, you can not perform a MD simulation. Before we go into the various potential functions, we need to understand the
Move to MD
Molecular Dynamics : the nitty gritty details !
through the slides of Prof. P. K. Padmanabhan
Two Excellent Books
=
i j ij i
F F
i i i
m F a / =
3
F12 F23 F13
F1 = F12 + F13 F3 = F31 + F32 F2 = F21 + F23
Newton’s IInd Law:
Three Atoms 1
ij ij
U F −∇ =
2 3 Progress in time…
x(t) → x(t+Δt) y(t) → y(t+Δt) z(t) → z(t+Δt)
1 2
Taylor Expansion: too crude to use it as such!! Δt ~ 1-5 fs (10-15sec) Advance positions & velocities of each atom:
A good Integrator …for example!
Verlet Scheme: Newton’s equations are time reversible, Summing the two equations, Now we have advanced our atoms to time t+Δt ! ! Velocity of the atoms:
1 2
1 2
F12 F23 F13
…Atoms move forward in time! t+Δt t+2Δt
Calculate 2. Update 3. Update The main O/P of MD is the trajectory.
Δt ~ 1-5 fs (10-15sec)
The missing ingredient… Forces?
Gravitational Potential:
r m Gm U
2 1− =
too weak, Neglect it!! Force is the gradient of potential:
!
r q q U
2 14 1 πε =
However, this pure monopole interaction need not be present!
Interatomic forces for simple systems
Gives an accurate description of inert gases (Ar, Xe, Kr etc.) Faithful in describing pure ionic solids (NaCl, KCl, NaBr etc.)
(non-bonded interactions) Instantaneous dipoles
The Lennard-Jones Potential
for Ar :
(kB)
Pauli’s repulsion “dispersion” interaction
r(nm)
) )( 6 12 ( 4
8 6 14 12 j i x ij
x x r r F − − = σ σ ε
ij ij
U F −∇ =
i i i x ij
x r r U x U F ∂ ∂ ∂ ∂ − = ∂ ∂ − =
Å
Should be known apriori.
Length and Times of MD simulation
Typical experiment sample contains ~ 1023 atoms! Typical MD simulations (on a single CPU) a) Can include 1000 – 10,000 atoms (~20-40 Å in size)! b) run length ~ 1 –10 ns (10-9 seconds)! Consequence of system size: Larger fraction of atoms are on the surface,
r dr m r m dr r N Ns 3 / 3 4 / 4
3 2
= = ρ π ρ π
7 8
10 ~ 10 ) 3 ( 3 ~ .) (
−
Α Α Expt N Ns 45 . ~ 20 ) 3 ( 3 ~ ) ( Α Α MD N Ns
Surface atoms have different environment than bulk atoms!
L The Simulation Cell
Insert the atoms in a perfectly porous box – simulation super-cell. The length of the box is determined as, L3 = M/Dexp = N*m/Dexp Dexp = Expt. density; m = At. mass; N= No. of atoms; If crystal structure/unit cell parameters are unknown (eg., liquids) Assign positions and velocities(=0) for each atom.
~20 A Periodic Boundary Condition
In 3-D the simulation simulation- super-cell is surrounded by 26++ image cells!
X L~20 A Y
n1, n2, n3 ∈
∈ ∈ ∈ -1,0, 1,
x’ = x + n1L y’ = y + n2 L z’ = z + n3 L Image coordinates: Now, there are no surface atoms! Construct Periodic Images:
~20 A Minimum Image Convention ~20 A Rc
A large enough system (ie, bigger sim.-cell) is chosen such that Rc ≤ ≤ ≤ ≤ L/2. Thus particle i interact either with particle j or one of its images, but not both! Interactions between atoms separated by a chosen cut-off distance (Rc) or larger (ie, rij > Rc) are neglected. Rc is chosen such that U(Rc) ~ 0
t = 10 pico-sec As Time Progress…
How to find the image of j that is nearest to i ?
i j
) )( 6 12 ( 4
8 6 14 12 j i ij ij x ij
x x r r F − − = σ σ ε Force between i & j :
i i j
L L
folding
If you have coordinates which are anywhere between (-infinity,infinity), and if you want to bring the coordinates between (0,L) then this procedure is called folding. Assuming L is 10, we see that : (234,-546) (4,-6) (4,4). If we have L = 10, we get the same answer. However, if L = 12, 12 x 19 = 228 and 12 x 45 = 540. Therefore, you get (234-228, -546+540) (6,-6) (6,6).
Unfolding
time steps then can you unfold it ? That is, can you map it to the range (-infinity, infinity) ?
4:(3.09,3.9), 5:(2.2,4.2), 6: (1.4,4.3), 7: (0.3, 5.0), 8: (9.4, 6.6)
have large difference (of the order of L), then we can guess that they have been folded. Then we can see that the unfolded coordinates at step 8 are (-0.6,6.6)
The three lines of code…
dx = x(j) – x(i) dy = y(j) – y(i) dz = z(j) – z(i) dx = dx – L*ANINT(dx/L) dy = dy – L*ANINT(dy/L) dz = dz – L*ANINT(dz/L)
Define,
ANINT(1.51) = 2 ANINT(3.49) = 3 ANINT(-3.51) = -4 ANINT(-11.39) = -11.0
rij_2 = (dx**2 + dy**2 + dz**2)
) )( 6 12 ( 4
8 6 14 12 j i ij ij x ij
x x r r F − − = σ σ ε
rij_8 = rij_2**4 rij_14 = rij_2**7
j i
L=10 A L=10 A (2,8) (0,0) (-1,6) (19,-4)
dx = dx – L*AINT(dx/L) dy = dy – L*AINT(dy/L) xi-xj=(2-19)-10*ANINT{(2-19)/10}= -17 - (-10*2)=3
yi-yj=(8+4)-10*ANINT{12/10}= 12 - (10*1)=2
(9,6)
The Lennard-Jones Potential
for Ar :
(kB)
Pauli’s repulsion “dispersion” interaction
r(nm)
) )( 6 12 ( 4
8 6 14 12 j i x ij
x x r r F − − = σ σ ε
ij ij
U F −∇ =
i i i x ij
x r r U x U F ∂ ∂ ∂ ∂ − = ∂ ∂ − =
Å
Should be known apriori.
The structure of a simple MD code
Remarks on Statistical Ensemble
There is no energy coming in Or going out of our system of atoms: micro-canonical (NVT) ensemble. Thus the total energy (E) and total linear momentum of the system should be conserved – through out our simulation!
= N i i v i m K 1 2 2 1
=
N i i
U U
1
2 1
U K E + =
remain const. (fluctuates) (fluctuates) Transient phase Equilibrated state
Calculating Temperature
Equipartition theorem: Average temperature :
1
1 ( )
M m
m
T T t M
=
>=
M – no. of MD steps performed
2 3
Even if we start with vi =0, the system picks up non-zero v (hence some T) as time progress! This some T(!) need not be what we want! So how do we control T?
N i i iv
m t T Nk
2
2 1 ) ( 2 3
Instantaneous temperature:
N i i iv
m Nk t T
2
3 1 ) (
Controlling the Temperature ?
Velocity rescaling:
N i i iv
m Nk T
2
3 1
Actual temp. at some instant.
1 2 Tr v v r T
T T T T T
r r
∆ + > < ∆ −
If T is out side the fluctuation window around Tr: Then scale all velocities: This phase of the simulation should not be used for averaging! However to sustain the temp. around Tr we will need to do this procedure several times at intervals.
Calculating thermodynamic quantities Average Pressure, Heat Capacity, Cv: Average Energy,
>= <
M i i
U M U
1
1
Remarks on Energy Conservation
(Routines available in standard packages. Or, do an MD with constant velocity scaling.)
6 ~ 10 / E ps E ∆ −
system.
Remarks on Interatomic Forces
intact irrespective of the environment around the atom! This can be a poor assumption for highly polarizable atoms/ions! Solution? Develop a shell model of atoms/ions! Or DFT-based ab-initio (Car-Parrinello) MD calculations ! FF’s are developed by empirical methods or ab-initio calculations.
Comments on Classical MD
Extensively employed to understand Physical processes at atomic resolution Phase Transitions, Diffusion and transport properties, Local structural and short-time relaxation of crystalline and amorphous solids liquids solid-fluid interfaces nano-clusters Very powerful in studying a variety of physical phenomena and under several external conditions (T & P). And, serves a very useful bridge between experiment and theory! Not useful in the study of electronic properties! Not powerful enough to describe chemical reactions!
Structural Characterization
Radial Distribution Function (rdf)
Rdf of fcc-solid
Running Coordination Numbers:
Dynamical Properties: Diffusion Coefficient Fick’s Law: Continuity Eq.: Diffusion Eq.: Einstein’s relation:
T fk D Nq
B
/
2
= σ
Nernst-Einstein’s relation:
CaF2
MSD time
…Diffusion Coefficient: CaF2 F- Ca2+
T fk D Nq
B
/
2
= σ
PRB, 43, 3180,1991.
…Dynamical Properties: Vibrational Spectrum Velocity Autocorrelation Spectrum, Water @300 K Power Spectrum, Helps interpreting IR spectrum.
…Structural Properties: Site Occupancies PRL 97, 166401 (2006)
Ion Channel’s In NASICON’s
Supriya Roy and P.P. Kumar, PCCP (2014).
Na3Zr2Si2PO12
Simulation of Live-Virus!!