Transcript Standard and accelerated molecular dynamics

How can we integrate Newton’s equations of motion?

..



mi ri   ri V (r1 , r2 , , rN )
Many-particle problem. We can forget about analytical solutions.
Finite differences methods: time is discretized.
Key quantity: Dt (timestep). Idea:
algorithm
Configuration at time t
Configuration at t+Dt
Nothing strange: think about a Taylor expansion:
.
1 ..
1 ...
2
x(t  Dt )  x(t )  x(t )Dt  x(t )Dt  x(t )Dt 3  
2
6
The goal is to have small truncation errors, so that even “long” time steps
allow a reliable description of the system dynamics.
A good test to judge different algorithms is to compare how accurately they
conserve the total energy (rule of thumb: dE/E<10-5) of an isolated system
for a given time step. Obviously, the time step should be much smaller than
the typical times of the fundamental processes taking place. Solids>vibrations->time step << 1ps.
Usually the energy conservation gives a more strict condition on the time
step which is typically of the order of 1-10fs. WARNING: higher
temperatures might need a smaller timestep.
Let us follow the idea of the Taylor expansion. If I want errors of the
order of Dt4 ….
1
1
2
3
x(t  Dt )  x(t )  v(t )Dt  a (t )Dt  b(t )Dt  
2
6
1
v(t  Dt )  v(t )  a (t )Dt  b(t )Dt 2  ...
2
a (t  Dt )  a (t )  b(t )Dt  ...
b(t  Dt )  b(t )
If I know the initial values of x,v,a,b I can generate the evolution at any time, but ….
1) Where are Newton equation of motion (e.g. atomistic potential)?
2) How bad is the constraint that the time derivative of the acceleration is constant?
The Verlet algorithm
Simple trick: expand
x(t  Dt ) and x(t  Dt )including third-order terms
and add together the two power expansions.
1
1
2
3
4
x(t  Dt )  x(t )  v(t )Dt  a(t )Dt  b(t )Dt  err (Dt )
2
6
1
1
2
x(t  Dt )  x(t )  v(t )Dt  a(t )Dt  b(t )Dt 3  err (Dt 4 )
2
6
2
4
x(t  Dt )  x(t  Dt )  2 x(t )  a(t )Dt  err (Dt );
f x  m ax ....
In 3N dimensions ….



1 
ri (t  Dt )  2ri (t )  ri (t  Dt )  [ f i (t )]Dt 2  error(Dt 4 )
mi
This is the standard Verlet algorithm. It is simple and it
conserves the total energy very well even for “large” time
steps. Time-reversal guaranteed. Disadvantage: velocities
can be calculated only indirectly:



r (t  Dt )  r (t  Dt )
v (t ) 
 error(Dt 2 )
2Dt
The Velocity Verlet algorithm is an equivalent (same trajectories)
algorithm which handles better the velocities.
Several other algorithms can be used (predictor-corrector, RungeKutta, etc.)
Standard Verlet algorithm
x(t  Dt)  x(t  Dt)  2x(t)  a(t)Dt 2  err(Dt 4 );
Velocity Verlet algorithm
1
2
4
x(t  Dt )  x(t )  v(t )Dt  a (t )Dt  err (Dt );
2
a (t )  a (t  Dt )
v(t  Dt )  v(t ) 
Dt  err (Dt 3 )
2
Noi useremo l’algoritmo “Velocity Verlet”
Microcanonical ensemble (N,V,E)
• Isolated system  the total energy E=K+V is
conserved. This is our situation.
• What is temperature in the microcanonical
ensemble? We definite it as a measure of the
instantaneous kinetic energy:
1
2
K (t )   mi vi (t ) ;
2 i
N dof k BT (t )
K (t ) 
; N dof  3 N
2
1
2
K (t )   mi vi (t ) ;
2 i
2 K (t )
T (t ) 
N dof k B
This “kinetic” or “istantaneous” temperature fluctuates (is not conserved as it would be in
the canonical ensemble), with fluctuations scaling as 1/sqrt(N), exactly as the kinetic energy
fluctuates. If one takes the time average over a long run, then the equipartition theorem (typical
of the canonical ensemble) is recovered in an averaged sense.
Initial conditions
In order to solve Newton equations, we need to assign initial positions and
velocities. The particles start to exchange energy, until the system
equilibrates. Typical initial conditions:
 Positions: ideal situation (perfect lattice positions)
 Velocities: extracted from a Maxwellian distribution (see after)
Corresponding to which temperature ?
Equipartition: if the system has “zero” (minimum) potential energy, and
we supply an initial kinetic energy K0, if the temperature is sufficiently
low to guarantee that only the harmonic region of the potential energy is
sampled, then half of it would transform into potential energy

1
T f  Ti
2
(valid if the potential is harmonic)
Approximating at step 0 the real Maxwellian
f (v x )dv x 
1
2
m
e
k BT

1
mv x2
2 k BT
dv x
It is not necessary to actually use a Maxwellian. Instead, we
shall use a rectangular distribution, normalized and with the
same momenta (one and two) of the Maxwellian. After
equilibration the evolution is exactly the same
f(v)
1/2C
v
C
-C
C
One time in the code:
seed=16324478;
rng(seed);
then each time I find
vx(i)=C*(2*rand(0)-1)
rand() generates a new
random number between
0 and 1
1
1 m
k BT
3k BT
2
2
mv 
v dv 
C 
;

2
2 2C C
2
m
vx(i)=C*(2*rand(0)-1)
f(v)
1/2C
v
-C
C
Since I make only a finite number of random extractions,
the average v would not be 0. Since Newton equations of
motion conserve the momentum the whole cell is going
to shift during evolution. Boring!
Also, the average T will not be exactly the one I wanted.
Boring. Easy solution: after extracting the velocity in one
direction ….
3k BT
C
;
m
xm=108x1.66x10-27/16
is the number you
should put in your code!!!!
XkbT (kB*T) [eV]=T(K)/11603.d0
Warning! First set to zero the
momentum, then rescale!!! If
you do the opposite it does not
work.
vxtot
vxtot
  v x ,i ; v x , i  v x , i 
N
i
To adjust the temperature …
tot
v
vxtot   vx ,i ; vx ,i  vx ,i  x
N
i
1
3
m vx2,i  v y2,i  vz2,i  Nk BT ;
2 i
2

v x ,i  v x ,i

T
;
T
Units, units, units …..
Say we want energies in eV, distances in Å and time in seconds …..
What is the MASS unit?
2
eV

s
E  mc 2  m 
A2
Conversion from Kg?
1eV=1.6x10-19 J; 1A2=10-20m2 
1eV x s2/A2= 1.6x10-19 J x s2x1020/m2 =16 Kg

108 is the weight in grams of an avogadro number of atoms;
Ag
108
108x10-3 is the weight in Kg of an avogadro number of atoms;
108x10-3/(6.023x1023)=108x1.66x10-27 is the weight in Kg of 1 atom
xm=108x1.66x10-27/16 is the number you should put in your code!!!!
XkbT (kB*T) [eV]=T(K)/11603.d0
Matlab implementation (build a function velinit; start using global!):
Ti=300;
mass=108*1.66e-27/16 (why? Length=Angstrom, En=eV, time=s, temp=K …)
Kb=1/11603 (eV/K)
cost=sqrt(3*Ti*Kb/(mass)); Angstrom/s: 2.631056523155901e+12 (sounds ok?)
ONLY ONE TIME IN THE CODE (MAIN/VELINIT FUNCTION):
seed=16324478;
rng(seed);
Poi, ogni volta che appare il comando «rand» il codice genera un numero casuale tra 0 e 1
2*rand=tra 0 e 2
2*rand-1 tra -1 e 1
cost*(2*rand-1) tra –cost e cost …
Estraggo vx,vy,vz per ogni atomo.
Calcolo momento lineare totale e energia cinetica totale Ekin
Calcolo la temperatura «Tnow» invertendo Ekin=3N/2 kB*Tnow; (Primo check)
Sottraggo momento lineare totale alle velocità. Ricalcolo Tnow; (Secondo check)
Aggiusto la temperatura scalando la velocità basandomi su sqrt(Ti/Tnow) (Terzo check).
Ok, siamo pronti per Verlet.
Initial and final T in a microcanonical simulation: example
call leggipos-> leggo pos. iniziali
call steepest-> minimizzo struttura iniziale (se non voglio metti call vicini)
call velinit -> assegno vel. iniziali
call forze(fxold,fyold,fzold)-> old vuol dire al tempo t, senza old al tempo t+dt
do imd=1,md!ciclo nel main di dinamica molecolare
do i=1,nat
x(i)=x(i)+vx(i)*dt+0.5d0*(fxold(i)/xm)*dt*dt
y(i)=
end do
call vicini
call forze(fx,fy,fz)
do i=1,nat
vx(i)=vx(i)+0.5d0*((fx(i)+fxold(i))/xm)*dt
fxold(i)=fx(i)! il vecchio valore di fxold non serve piu’
vy(i) …
compute here sum of v2 to get kin energy
v2=v2+…
end do
ekin=0.5d0*xm*v2
temp=…
call epot()
write(*,*)imd*dt,temp,epot+ekin
end do
Molecular dynamics in the canonical ensemble.
 Temperature (not the istantaneous one! see after) is constant, energy
fluctuates (N,V,T).
The system is in contact with a thermal bath.
How do we introduce the thermal bath in the equations of motion ?
Very simple way: Andersen thermostat
We introduce a thermostat by introducing a probability of suffering a
“thermalizing collision” with a virtual bath particle.
Probability of suffering a collision with the thermal bath in the time
interval Dt: nDt, where n is the average frequency of collisions.
P(t ,n )  n exp[nt ]
Other thermostats: Langevin, Nosé-Hoover,
Gaussian,…..
Initial and final T in a canonical simulation: example
Istantaneous T still fluctuates !!What does it
really mean to simulate in the canonical
ensemble ?
Ensemble averages in MD simulations
 A  N ,V ,T



A(q, p)e


e
H (q, p )
k BT
H (q, p)
k BT
dqdp
dqdp
t

1
 lim
A(q( ), p( ))d
t  t
0
canonical
 A  MD 
1
N steps
N steps
 A(t )
t 1
microcanonical
Propagating a trajectory  sampling the phase space G(q,p)
G(q,p)
Algorithms may become instable for very long time evolutions +
we may be interested only in the system behavior in some
portions of the phase space
 often many, shorter, trajectories with different initial
conditions are run. Typically (ex: rate of escape from a given
spatial configuration)
q
Interesting
configuration
Initial positions
This time
is lost
p
Results of a MD simulation
1) Look at what happens and learn about interesting trajectories
(Does the system change configuration? How? How fast? At what
temperature?)
2) Measure statistical quantities (Mean square displacement,
pressure, correlation functions, etc.)
Si; diamond lattice
Let us fix the attention on a particular atom: at T=0K its “neighbors’ shells” are
very easily defined: n.n. are exactly at a distance d=2.35Å; second nn at d=3.84Å,
third at 4.5Å, ....
Second
First
Third
0
d
T=0K
Second
First
Third
0
d
T>0K
Second
First
Third
0
d
T>>0K
A temperatura ambiente la situazione per il Si è simile a:
Secondi
Primi vicini
Terzi
0
Buon separatore
tra primi e secondi
vicini
d
T>0K
Let us compute the “g(r)” of Si
Stillinger-Weber potentials
(216 atoms, time-step: 1fs)
r
N (r )  4  x g ( x)dx
2
0
N

V
(radial distribution re-normalized over )
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Crystalline Si at T=300 K
Liquid Si at T=1800 K
Liquid Si at T=1800 K
Liquid Si at T=1800 K
Liquid Si at T=1800 K
Liquid Si at T=1800 K
Liquid Si at T=1800 K
Liquid Si at T=1800 K
Liquid Si at T=1800 K
Liquid Si at T=1800 K
First neighbors shell
Second neighbors shell
MD codes should be simple. If not,
they are probably badly written!
A microcanonical MD code could
look like ..
program md
implicit, dimension, etc.
call init_pos(x,y,z) [read from file, create a lattice, etc.]
call steepest(x,y,z,....) [if I want to start from a local minimum]
call init_vel(vx,vy,vz) [velocity initilialization]
call nbrs(...)
call forces(fx,fy,fz,...)
do i=1,nsteps
do j=1,natoms
x(j)=x(j)+v(j)*dt+(0.5/mass)*fx(j)*dt*dt [same for y and z]
end do
call nbrs(....)
call forces(fxnew,fynew,fznew)
do j=1,natoms
vx(j)=vx(j) + (0.5/mass)*(fxnew(j)+fx(j))*dt [same for vy and vz]
fx(j)=fxnew(j) [same for fy and fz]
end do
call temperature(ekin,temp,vx,vy,vz,...) [computes temperature and ekin]
call elennard(x,y,z,epot,...) [computes the potential energy]
write(*,*)i,temp,epot,ekin,epot+ekin
end do
stop
end
Molecular statics
Instead of looking at the time evolution of a certain system, we might want
to find minimum-energy configurations.
Many-particle systems have, in general, several local minima. Finding them
all is highly unpractical and for some systems impossible. Practical and
possible: find a local minimum starting from a configuration “not to
far” from the minimum itself
Molecular statics
Many “geometrical alghoritms” available such as
-steepest descent
-conjugate gradient
Molecular statics: quenched MD
Minimization can be also achieved through a simple
manipulation of the equation of motions:
t=t0

v
 
f v  0
Molecular statics: quenched MD
Minimization can be also achieved through a simple
manipulation of the equation of motions:
t=t0+Dt

v
 
f v  0
Molecular statics: quenched MD
Minimization can be also achieved through a simple
manipulation of the equation of motions:
t=t0+2Dt

v
 
 
f v  0  v  0
Molecular statics: quenched MD
Minimization can be also achieved through a simple
manipulation of the equation of motions:
t=t0+nDt
 
v  0  T  0K
(configuration quenched)