1. science
2. physics

D G Kanhere,
Center for Simulation and modeling
Pune University,
Mastani School- IISER Pune ,
July 2014


What is Molecular Dynamics (An Overview)
Ab Initio MD
 Born Oppenheimer Dynamics
 Car Parrinello Dynamics




General Comment on “ How to set up
dynamics run”
Thermodynamics of clusters
Multiple histogram methods
Results of clusters – As Examples

System
Many particles interacting with each other
 Solid/liquid : Box of Volume V, Temperature T and
Pressure P: Periodic BC
 Finite size systems clusters :Free BC

Given the forces acting on all the ions, and initial
state at time t=0, Compute
 Trajectories of all the particles as a function of time t,
by using Newton's laws: Essentially exact
 We have entire phase space trajectories! Therefore all
the information to compute various statistical
quantities from microscopic description.
Molecular Dynamics





Consider N particles interacting VIA Lennard
Jones Potential (N ~ 10-13)
For given co-ordinates {Ri}write down total
energy : V ( R1,R2,…)
Compute Force on atom i.
Write a subroutine which takes coordinates
as input and returns total potential energy
and forces on each of the atoms.
You are now ready to do simulation
experiments after one more step.
How to integrate Newton’s equations
xI (t )
d 2 RI
MI
 FI
2
dt
t
t
(1) RI (t  t )  RI (t )  vI (t )t 
1 FI
t 2  O(t 3 )
2 MI
(2) RI (t  t )  RI (t )  vI (t )t 
1 FI
t 2  O(t 3 )
2 MI
(1)  (2) RI (t  t )  RI (t  t )  2 RI (t ) 
FI
t 2  O(t 4 )
MI
Verlet’s algorithm
FI
RI (t  t )  2 RI (t )  RI (t  t ) 
t 2  O(t 4 )
MI
A typical molecular dynamics protocol
Initialize: select starting atomic positions and velocities as close as
possible to thermal equilibrium
Integrate: compute all forces and determine new positions using
Verlet’s algorithm
Equilibrate: let the system loose memory of the initial configurations,
i.e. let it reach thermal equilibrium
Average: accumulate averages of observables of interest (A)
T


1
A   A( R(t ), P(t )) dt
T0
 
 
A( R, P) exp(  E ) dRdP

A 
 
 exp( E )dRdP
Average in Molecular Dynamics = Average in Statistical Mechanics
Molecular Dynamics
Sum of pair interactions?
How to treat chemical
complexity?
Where are the electrons?
“Empirical” potentials
Starting point: Pair potentials (Lennard-Jones, BornMayer, Coulomb, etc)
 R0 12  R0  6 
Ze 2
4      , 
  e R , etc...
R
 R  
 R 
+ three-body corrections
 (3) ( I , J , K ) V0 cos   1 / 3
  IJˆK
I
R  RI  RJ
K
J
+ density dependent terms (embedded atom models)
+ atomic distorsion terms (includes polarization)
+ charge transfer terms
Parameters are determined from “empirical” data, such as
experimental EOS, vibrations, phase diagrams, dynamical properties,
etc.
Quantum simulations: The “standard
model”
“Molecular dynamics”
for atoms
Ma = F = -dE/dR
Schroedinger equation
for electrons
Hy  Ey
R. Cohen
“Ab-initio” molecular dynamics =
e--e- interactions:
Density Functional Theory
e--nuclei interactions:
Pseudopotentials
Classical molecular dynamics in the
potential energy surface generated by the
electrons in their quantum ground state
Where are the electrons? The adiabatic approximation
d 2 RI
dE ({R})
MI
 FI  
2
dt
dRI
Electrons respond much faster than
nuclei to external forces, because of
their lighter mass, therefore:
Electrons are always in their
instantaneous quantum mechanical
ground state, for each given {R}.
ER  y e H e Ry e
Consequence:
Every MD step requires the
calculation of the QM ground
state of the electrons
ee-
e-
e-
eNuclei
e-
The electronic Hamiltonian He
depends parametrically on the
nuclear positions {R}
Density functional theory
Electron-electron (many-body)
interaction
Density-functional theory [W. Kohn et al. 1964-1965] states that the
e-e interaction can be written as a one-electron “effective” term:
å
j
e2
ri - rj
º
ò
re (r')
ri - r'
dr'+ Vxc (re (ri ))
However, the exact functional form of Vxc is not (yet) known.
Current approximations to the exact Vxc go under different names (LDA, BLYP,
BP, GGA, “hybrid”, et al). While GGA and hybrid functionals provide slightly
better results than other approximations, the choice of the Vxc is often made in
such a way to improve agreement with exps (“ab fine”??).
“Ab-initio” forces: the Hellmann-Feynman theorem
d y 0 H e ({R})y 0
d 2 RI
dE ({R})
MI
 FI  

2
dt
dRI
dRI

d y 0 H e ({R})y 0
dRI
  y0
dH e ({R})
dy 0
dy 0
y0 
H e ({R})y 0  y 0 H e ({R})
dRI
dRI
dRI
  y0
dH e ({R})
dy 0
dy 0
y 0  E ({R})
y 0  E ({R}) y 0
dRI
dRI
dRI
d y0 y0
dH e ({R})
  y0
y 0  E ({R})
dRI
dRI
  y0
dH e ({R})
y0
dRI
Forces can be calculated without
recalculating the ground state
wave function for small atomic
displacements
The Born Oppenheimer
Dynamics
How to keep electron in the ground state (1)
d 2 RI
dE ({R})
MI

F


I
dt 2
dRI
E ({R})  y 0 H e ({R})y 0
If we expand wavefunctions in a basis set
y   cii
i
y H e y   ci c*j i H e  j
i, j
finding the ground state is equivalent to
minimizing a quadratic form in the {c}’s
(variational principle). So, standard
minimization
schemes
can be used, e.g. steepest descent:


y 
y H e y   H ey
y
Carr- Parrinello dynamics
How to keep electron in the ground state (2)
d 2 RI
dE ({R})
MI

F


I
dt 2
dRI
E ({R})  y 0 H e ({R})y 0

y   H e ({R})y
down to the minimum
for each {R}

y   H e ({R})y
???

dE ({R})
M I RI  
dRI
{R}
The Car-Parrinello algorithm
Car-Parrinello Molecular Dynamics
(i) integrate the equations of motion on the (long) time scale set by
the nuclear motion but nevertheless
(ii) take intrinsically advantage of the smooth time evolution of the
dynamically evolving electronic subsystem as much as possible.
The second point allows to circumvent explicit diagonalization or
minimization to solve the electronic structure problem for the next
molecular dynamics step as it is done in Born-Oppenheimer
Molecular dynamics. Car-Parrinello molecular dynamics is an
efficient method to satisfy requirement (ii) in a numerically stable
fashion and makes an acceptable compromise concerning the
length of the time step (i).
In CPMD a two-component quantum / classical problem is mapped onto a twocomponent purely classical problem with two separate energy scales at the
expense of loosing the explicit time dependence of the quantum subsystem
dynamics.
Now, in classical mechanics the force on the nuclei is obtained from the
derivative of a Lagrangian with respect to the nuclear positions. This suggests
that a functional derivative with respect to the orbitals, which are interpreted as
classical fields, might yield the force on the orbitals, given a suitable Lagrangian.
In addition, possible constraints within the set of orbitals have to be imposed,
such as e.g.
orthonormality (or generalized orthonormality conditions that include an overlap
matrix).
Car-Parrinello Lagrangian

y
yyy
y

Car and Parrinello (1985) have postulated the following Lagrangian
1 1
L

M
R


|

E
{
}

(

|

)



C
P
R
ii
i
D
F
T
i
i
j
i
i
i
j
2i
2
R
i
j
Kinetic Energy
Potential Energy
Constraints
where i are fictitious masses and yi are classical fields. Classical
action needs to be minimized which results in equation of motions
S   LCP (t)dt
d dLCP (t) dLCP (t)

dt dR
dR
d dLCP (t) dLCP (t)

dt dyi (r, t) dyi (r, t)
Car-Parrinello Equations of Motions
Car and Parrinello (1985) have derived equations of motions
2

E
[
{
y
}
,
{}
R
]
Z
e
D
F
T
i
R
M
R




(,)
r
t

d
r
R


R
|r
R
()
t|

E
[{
y
,{}
R]
D
F
T
i}

y
rt 


ijy
rt 

i
i(,)
j(,)

y
rt
j
i(,)

H
y
rt 
ijy
rt

D
F
T
i(,)
j(,)
j
which are obviously transformed back to Born-Oppenheimer
molecular dynamics if fictitious masses for the electrons i
At the eqilibrium, there are no forces on electrons, therefore
Ground state of density functional theory is reached with the
eigenvalues being the Kohn Sham eigenstates.
Why does the Car-Parrinello Method works?
Conserved Energy in CPMD is not a physical energy but
supplemented with a small fictitious kinetic term
1
E
M
R

E
{}
y

p
h
y
s
R
D
F
T
i
2
R
1
1
E
M
R


y
y

E
[
{}
y
R
}
]

E
T


c
o
n
s
R
i
i|
i
D
F
T
i{
p
h
y
s
f
i
c
t
2
2
R
i
1
T

y
y

f
i
c
t
i
i|
i
i 2
Various Energies extracted from CPMD for a model system.
EDFT
Tfict
The fictitious kinetic energy of the electrons is found to perform bound oscillations
around a constant, i.e. the electrons do not heat up “ systematically “ in the presence
of the nuclei; note that Tfict is a measure for deviations from the exact BornOppenheimer surface. Closer inspection shows actually two time scales of oscillations:
the one visible in the Figure stems from the drag exerted by the moving nuclei on the
electrons and is the mirror image of the EDFT fluctuations.
As a result the physical energy (the sum of the nuclear kinetic energy and the electronic
total energy which serves as the potential energy for the nuclei) is essentially constant
on the relevant energy and time scales.
Given the adiabatic separation and the stability of the propagation, the central
question remains if the forces acting on the nuclei are actually the “correct" ones
in Car-Parrinello molecular dynamics. As a reference serve the forces obtained
from full self-consistent minimizations of the electronic energy at each time step, i.e. BornOppenheimer molecular dynamics with extremely well converged wavefunctions.
How to control adiabaticity?
Since the electronic degrees of freedom are described by much heavier masses
than the electronic masses, time step to perform CPMD simulations needs not to
be too small as compared to Ehrenfest molecular dynamics.
For a system with a gap in the spectrum, the lowest possible frequency of
“fictitious” electronic oscillations
2min ~ Egap
1/2
 Egap 
min ~ 




To guarantee adiabatic separation this frequency should be
much larger than the typical phonon energy and/or the gap
in the spectrum which would make sure that the electrons
follow the nuclei adiabatically. Hence fictitious mass .
At the same time, small fictitious mass would imply
smaller and smaller time step because maximum fictitious
electronic frequency is proportional to the plane-wave
cutoff energy
 2 m ax ~ E cut
 m ax
E 
~  cut 
  
1/ 2
  
 t m ax ~ 

E
 cut 
1/ 2
As a result a compromise fictitious mass needs to be found
in CPMD simulations.
For metals gap is zero and zero frequency “fictitious” electronic
modes occur in the spectrum overlapping with the
phonon spectrum. Thus, a well-controlled Born-Oppenheimer
approach can only be recommended
CP Method as dynamical solution of DFT equations
CP Method invented a new way to solve Kohn-Sham
equations alternative to diagonalization.
Consider variational principle which can be used to find an
upper bound for the lowest eigenstates of the hamiltonian
y|H|y
y cnn
using basis set expansion
n
c |H|c 

n
n
m
n
m m
2
|
c
|
n 1
n
m
CPMD offers a way to determine the coefficients without
reduction to the eigenvalue problem.
Born-Oppenheimer –versus- Car-Parrinello AIMD
Born-Oppenheimer
Car-Parrinello
Born-Oppenheimer –versus- Car-Parrinello forces
Born-Oppenheimer –versus- Car-Parrinello: summary
from M Sprik
Ab-initio Molecular Dynamics: bibliography
• R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985)
• M. Payne, M. Teter, D. Allan, T. Arias, J. Joannopoulos, Rev. Mod.
Phys. 64, 1045 (1992).
• D. Marx, J. Hutter, "Ab Initio Molecular Dynamics: Theory and
Implementation", in "Modern Methods and Algorithms of Quantum
Chemistry" (p. 301-449), Editor: J. Grotendorst, (NIC, FZ Jülich 2000)
• http://www.theochem.ruhr-uni-bochum.de/research/marx/cprev.en.html
• R. Rousseau and S. Scandolo, “Car-Parrinello Molecular Dynamics”, in
“Encyclopedia of Condensed Matter Physics”, edited by G. Bassani, G.
Liedl, and P. Wyder, Elsevier, Amsterdam (2005)
Finite Temperature Properties
of finite size systems
•D G Kanhere
•Pune University
Lindemann Criteria
Mean Square Displacements
Extracting ionic density of
States
•Constant Temperature
•Constant Volume
•Phase space trajectories
•The Multiple Histogram Method
•To
•Extract Density of States
Sodium Clusters
●
●
●
●
●
Simplest of atomic clusters
Jellium model works ( Depends…)
Nice delocalized charge density
Magic Clusters at N=8,20,40,58,92,138,…
Icosahedra for N=13,55,147,…
Chacko et al., Phys. Rev. B 71,155407 (2005)
Lee et al., J. Chem. Phys. 123,164310 (2005)
Lee et al., Phys. Rev. B (2007)
Shahab et al., Phys Rev B ( 2007)
Gazi et al J Chem Phys ( 2008)
The heat Capacities
N=40 : peaked
N=50 : flat
N=55 : very sharp
N=58 : peaked but broad
BUT
Highest melting Point
( Tm = 375 K )
Sodium Clusters
●
●
●
●
●
Simplest of atomic clusters
Jellium model works ( Depends…)
Nice delocalized charge density
Magic Clusters at N=8,20,40,58,92,138,…
Icosahedra for N=13,55,147,…
Chacko et al., Phys. Rev. B 71,155407 (2005)
Lee et al., J. Chem. Phys. 123,164310 (2005)
Lee et al., Phys. Rev. B (2007)
Shahab et al., Phys Rev B ( 2007)
Gazi et al J Chem Phys ( 2008)
Application II
Thermodynamics
Size sensitive specific heats
Magic melters
Higher than Bulk Melting
temperatures
The Case of Gallium and
Aluminum clusters
●
●
●
Gallium and aluminum clusters show
extreme size sensitivity
Gallium clusters melt at Higher than their
Bulk Tm
Magic Melters ----- Non - Melter
Gallium Clusters: Heat Capacity
Experimental Data from
Indiana group
*
N=30 Flat
N=31 Peaked
*
Tm Shifts by 200K @ N=45
And by 350K across the
series
Calculations
•
•
•
•
Density functional with LDA/GGA
Born Oppenheimer MD
Delta T 100 au
Total simulation time per temp 100 ps or
more
• Soft pseudo potentials, plane waves
• Extensive search for equilibrium
geometries. Simulated annealing, Basin
Hopping …
The Heat Capacity : Ga30- Ga31
Joshi et al. Phys. Rev. Lett. 96,135703 (2006)
The MSD : Ga30 and Ga31
Finite Size Effect : Amorphous Clusters
• In amorphous cluster each atom may have different
environment.
• Atoms may be bonded with the rest of the cluster with different
strength.
• When heated, they will begin diffusive motion at different
temperatures
• This may result in continuous phase change.
• No sharp peak in specific heat, very broad transition.
• Cluster with large island having local order will show a peak in
the heat capacity ..Most of the atoms “melt” together.
Thermodynamics: Ga27Si3
Ghazi and Kanhere J Phys Chem (2012)
Ga13 – Ga12C and Al13 - Al12C
Ga13 is Decahedra
Ga12C is a Perfect
Icosahedra
Calculated Specific Heat Ga
Specific Heat