Transcript Document 7168657

Nano Mechanics and Materials:
Theory, Multiscale Methods and Applications
by
Wing Kam Liu, Eduard G. Karpov, Harold S. Park
2. Classical Molecular Dynamics

A computer simulation technique

Time evolution of interacting atoms pursued by integrating the
corresponding equations of motion

Based on the Newtonian classical dynamics

Method received widespread attention in the 1970`s

Digital computer become powerful and affordable
Molecular Dynamics Today

Liquids



Defects in crystals


Provides insight on ways and speeds of fracture process
Surfaces


Improved realism due to better potentials constitutes a driving force
Fracture


Allows the study of new systems, elemental and multicomponent
Investigation of transport phenomena i.e. viscosity and heat flow
Helps understand surface reconstructing, surface melting,
faceting, surface diffusion, roughening, etc.
Friction


Atomic force microscope
Investigates adhesion and friction between two solids
Molecular Dynamics Today

Clusters




Biomolecules


Corporation of many atoms (few-several thousand)
Comprise bridge among molecular systems and solids
Many atoms have comparable energies producing a difficulty in
establishing stable structures
Dynamics of large macromolecules
(proteins, nucleic acids (DNA, RNA), membranes)
Electronic properties and dynamics



Development of Car-Parrinello method
Forces on atoms gained through solving electronic structure problem
Allows study of electronic properties of materials and their dynamics
MD System


Subdomain of a macroscale object
Manipulated and controlled form the environment via
interactions

Various kinds of boundaries and interactions are possible

We consider:


Adiabatically isolated systems that can exchange neither
matter nor energy with their surroundings
Non-isolated systems that can exchange heat with the
surrounding media (the heat bath, and the multiscale boundary
conditions developed at Northwestern)
2.1 Mechanics of a System of Particles
Particle, or material point, or mass point is a mathematical model of a body whose
dimensions can be neglected in describing its motion.
Particle is a dimensionless object having a non-zero mass.
Particle is indestructible; it has no internal structure and no internal degrees of freedoms.
Classical mechanics studies “slow” (v << c) and “heavy” (m >> me) particles.
Examples:
1) Planets of the solar system in their motion about the sun
2) Atoms of a gas in a macroscopic vessel
Note: Spherical objects are typically treated as material points, e.g. atoms comprising a
molecule. The material point points are associated with the centers of the spheres.
Characteristic physical dimensions of the spheres are modeled through particle-particle
interaction.
Generalized Coordinates
Generalized coordinates are given by a minimum set of
independent parameters (distances and angles) that determine
any given state of the system.
(q1 , q2 ,..., qs )
Standard coordinate systems are Cartesian, polar, elliptic, cylindrical and spherical.
Other systems of coordinates can also be chosen.
We will be are looking for a basic form of the equation, which invariant for all coordinate systems.
Examples:
Material point in xy-plane
Pendulum
Sliding suspension pendulum
y
x

x
q1  x, q2  y
q1  

q1  x, q2  
We are looking for a basic form of equation motion, which is valid for all coordinate systems.
Generalized Coordinates
One distinctive feature of the classical Lagrangian mechanics, as compared with special
theories (e.g., Newtonian dynamics and continuum mechanics), is that the choice of
coordinates is left arbitrary.
(q1 , q2 ,..., qs )
Number s is called the number of degrees of freedom in the system.
The mathematical formulation will be valid for any set of s independent parameters that
determine the state of the system at any given time.
Solution of a problem begins with the search for such a set of parameters.
Least Action (Hamilton’s) Principle
Lagrange function (Lagrangian)
L(q1 , q2 ,..., qs , q1 , q2 ,..., qs , t )
t2
Action integral
S   L(q, q, t )dt
t1
Trajectory variation
q(t )   q(t ),  q(t1 )   q(t2 )  0
Least action, or Hamilton’s, principle
q (t )
q(t2 )
Sb
Action is minimum along the true trajectory:
Sa  Sb , Sa  Sc
t2
 S    L(q, q, t )dt  0
t1
q(t1 )
t1
The main task of classical dynamics is to find the true trajectories
(laws of motion) for all degrees of freedom in the system.
Sa
Sc
t2
t
Remark: in ADMD this
principle is used to obtain the
trajectories.
Lagrangian Equation of Motion
Lagrange function (Lagrangian)
L(q1 , q2 ,..., qs ; q1 , q2 ,..., qs , t )
t2
Least action principle
 S    L(q, q, t )dt  0
t1
Substitution of the Lagrange function into the action integral with further application of
the least action principle yields the Lagrangian equation motion:
d L L

 0,   1, 2,..., s
dt q q
Lagrangian equation is based on the least action principle only, and it is valid for all
coordinate systems.
Lagrangian Equation of Motion: Derivation
Using the most general form of the Lagrange function,
the least action principle gives
t2
q (t )
q(t1 )
Sa
t2
t2
t1
Variation of the coordinates cannot change the observable values q(t1) and q(t2):
 q : q(t )   q(t )   q(t1 )   q(t2 )
Therefore,
Sc
t1
 L
 L d L 
L 
L
    q   q  dt 
q   
  q dt  0

q

q

q

q
dt

q


t1 
t1 
t1
t2
Sb
t2
 S    L(q, q, t )dt    L(q, q, t )dt
t1
q(t2 )
t2
L
q  0
q
t1
and the second term finally gives
d L L

 0,   1, 2,..., s
dt q q
t2
t
Lagrange Function in Inertial Coordinate Systems
General form of the Lagrangian function is obtained based on these arguments:
In inertial coordinate systems, equations of motion are:
1) Invariant as to the choice of a coordinate system (frame invariance), and
2) Compliant with the basic time-space symmetries.
Frame invariance:
Galilean coordinate transformation
r  r ' Vt , t  t '
Galilean relativity principle
L(r, r, t )  L '(r ', r ', t )
Time-space symmetries
Homogeneity and isotropy of space
L(r, r, t )  L '(r  a, r, t )
Homogeneity of time
L(r, r, t )  L '(r, r, t   )
Lagrange Function of a Material Point
Free material point
(generalized and
Cartesian coordinates)
Interacting material point
Conservative systems
m 2
m 2 m 2
L  q , L  r   x  y2  z2 
2
2
2
m 2
L  q  U (q )
2
L  T ( q)  U ( q)
E  2T  L  Const
For conservative systems, kinetic energy depends only on velocities, and potential
energy depends only on coordinates.
Lagrange Function: Examples
Pendulum
m 2 2
L  l   mgl cos( )
2
l

q1  
1) Kinetic energy; 2) Potential energy due to
external gravitational field, where g is the
acceleration of gravity).
m 2
L  y  mgy   e   y
2
Bouncing ball
q1  y
Particle in a
circular cavity
y
1) Kinetic energy; 2) Potential energy due to
external gravitational field; 3) Potential energy
of repulsion between the ball and floor.
y
x
q1  x, q2  y
   R
m 2
2
L  (x  y ) e
2
x2  y 2

1) Kinetic energy; 2) Potential energy of
repulsion between the particle and cavity wall.
Pendulum: Lagrange Function Derivation
General form
L  T ( )  U ( )
l
Kinetic energy
tangential velocity
Potential energy
m
m
T  v2  l 2 2
2
2
l 
v  lim
 l
t  0 t
l cosφ

vτ
U  mgh  mgl cos 
Potential energy depends on the height h with respect to a zero-energy level.
Here, such a level is chosen at the suspension point, i.e. below the suspension level, the
height is negative.
The total Lagrangian
m 2 2
L  l   mgl cos( )
2
Parameters used:
g  9.8 m/s 2 , l  2m
Pendulum: Equation of Motion and Solution
Lagrange function and equation motion:
L
m 2 2
d L L
g
l   mgl cos  

 0   (t )  sin  (t )   0
2
dt  
l
Initial conditions (radian):
 (0)  0.3,  (0)  0
Parameters:
g  9.8 m/s 2 , l  2m
 (0)  1.8,  (0)  0
 (0)  3.12,  (0)  0
Bouncing Ball: Lagrange Function Derivation
General form
L  T ( y)  U ( y)
Kinetic energy
m 2
T y
2
Potential energy
y
U  mgy   e  y
Purple dashed line: the first term (gravitational
interaction between the ball and the Earth).
Blue dotted line: the second term (repulsion between
the ball and the bouncing surface).
Red solid line: the total potential.
β is a relative scaling factor: the ball-surface repulsive
potential growths in eβ times for a unit length
ball/surface penetration (from y = 0 to y = – 1m).
The total Lagrangian:
m 2
L( y, y )  y  mgy   e   y
2
Parameters used:
g  9.8 m/s 2 , m  1kg,
  1J,   4m 1
Bouncing Ball: Equation of Motion and Solution
Lagrange function and equation motion:
m 2
d L L
   y (t )
 y
L  y  mgy   e


 0  y (t )  g 
e
0
2
dt y y
m
Parameters:
g  9.8m/s 2 , m  1kg,   1J,   4m 1
Initial conditions:
y(0)  5m
y(0)  0
y(0)  10m
y(0)  0
y(0)  15m
y(0)  0
Particle in a Circular Cavity: Lagrange Function Derivation
General form
L  T ( x, y )  U ( x, y )
Kinetic energy
m
T  ( x2  y 2 )
2
Potential energy
U  e
   R r 
 e
y
r
x

  R  x2  y 2
R

The potential energy grows quickly and becomes larger than the typical kinetic energy, when the
distance r between the particle and the center of the cavity approaches value R.
R is the effective radius of the cavity. At r < R, U does not alter the trajectory.
β is a relative scaling factor: the potential energy growths in eβ times between r = R and r = R+1.
The total Lagrangian
   R
m 2
2
L( x, y, x, y)  ( x  y )   e
2
x2  y 2

Particle in a Circular Cavity: Equation of Motion and Solution
Lagrangian function and equations:
   R
m
L  ( x2  y 2 )   e
2
x2  y 2

 m x(t )   e  ( R r (t )) x(t )/ r (t )
 
  ( R  r ( t ))
y(t )/ r (t )
m y(t )   e
  1027 J, r (t )  x 2 (t )  y 2 (t )
Parameters:   4nm1 , R  10nm, m  109 kg
Initial conditions:
Potential barrier:
Summary of the Lagrangian Method
1.
The choice of s generalized coordinates (s – number of degrees of freedom).
2.
Derivation of the kinetic and potential energy in terms of the generalized
coordinates
3.
The difference between the kinetic and potential energies gives the Lagrangian
function.
4.
Substitution of the Lagrange function into the Lagrangian equation of motion and
derivation of a system of s second-order differential equations to be solved.
5.
Solution of the equations of motion, using a numerical time-integration
algorithm.
6.
Post-processing and visualization.
Hamiltonian Mechanics
Description of mechanical systems in terms of generalized coordinates and velocities is not
unique.
Alternative formulation formulation in terms of generalized coordinates and momenta can be
utilized. This formulation is used in statistical mechanics.
 q, q 
Legendre’s transformation
(the passage from one set of independent
variables to another)
Differential of the Lagrange function:
Generalized momentum (definition):
L
p 
q

 q, p 
L  L ( q, q ) 
L
L
dL  
dq  
dq ,   1, 2,..., s
 q
 q
 dL   p dq   p dq


A mass point in Cartesian coordinates:
px  mx, p y  m y, pz  mz
Hamiltonian Equations of Motion
dL   p dq   p dq




dL   p dq  d   p q    q dp

 
 


d   p q  L    q dp   p dq

 
 
Hamiltonian of the system:
H ( p, q, t )   p q  L

dH   q dp   p dq


Equations of motion:
H
q 
,
p
H
p  
q
Hamiltonian Equations of Motion: Pendulum
Lagrange function and the generalized momentum:
L
m 2 2
L
l   mgl cos   p 
 ml 2
2

Hamiltonian of the system:
p2
p2
T
, U  mgl cos   H 
 mgl cos 
2
2
2ml
2ml
Equations of motion:

H
H
, p
p



p
, p  mgl sin 
2
ml
Conversion to the Lagrangian form (elimination of p):
p  ml 2  mgl sin    
g
sin   0
l
l

Hamiltonian Equations of Motion: Bouncing Ball
Lagrange function and the generalized momentum:
m 2
L
 y
L  y  mgy  e
 p
 my
2

Hamiltonian of the system:
p2
p2
 y
T
, U  mgy  e
 H
 mgy  e   y
2m
2m
Equations of motion:
y
H
H
, p
p
y

y
p
, p   mg   e   y
m
Conversion to the Lagrangian form (elimination of p):
p  my  mg   e   y  y  g 

m
e  y  0
y
Summary of the Hamiltonian Method
1.
The choice of s generalized coordinates (s – number of degrees of freedom).
2.
Derivation of the kinetic and potential energy in terms of the generalized
coordinates.
3.
Derivation of the generalized momenta.
4.
Expression of the kinetic energy in terms of the generalized momenta.
5.
The sum the kinetic and potential energies gives the Hamiltonian function.
6.
Substitution of the Hamiltonian function into the Lagrangian equation of motion
and derivation of a system of 2s first-order differential equations to be solved.
7.
Solution of the equations of motion, using a numerical time-integration
algorithm.
8.
Post-processing and visualization.
Predictable and Chaotic Systems
Divergence of trajectories
in predictable systems
q  q  t , q(0)  : q t , q(0)    q t , q(0)   C t
Variance of the trajectory depends linearly on time
Divergence of trajectories
in chaotic (unpredictable) systems
q  q  t , q(0)  : q  t , q(0)    q t , q(0)    e t
Variance of the trajectory depends exponentially on time (J.M. Haile,
Molecular Dynamics Simulation, 2002)
Trajectories in MD systems are unpredictable/unstable; they are
characterized by a random dependence on initial conditions.
Example Quasiperiodic System: Particle in a Circular Cavity
Initial conditions
x(0)  2.522 nm, x(0)  0
y (0)  0, y (0)  10 nm/s
x(0)  2.500 nm, x(0)  0
y (0)  0, y (0)  10 nm/s
Dependence of solutions on initial conditions is smooth and predictable in stable systems.
Example Unstable System: Discontinuous Boundary
Initial conditions
x(0)  2.522 nm, x(0)  0
y (0)  0, y (0)  10 nm/s
x(0)  2.500 nm, x(0)  0
y (0)  0, y (0)  10 nm/s
No predictable dependence on initial conditions exists unstable systems.
Comparison for Longer-Time Simulation
Stable quasiperiodic trajectory
Unstable chaotic trajectory
Divergence of Trajectories
Slow and predictable divergence
in stable systems
Relative variance of x(0):
 (0)  0.00 to 0.06
Quick and random divergence
in unstable systems
 (0)  0.0000 to 0.0003
Transition from Stable to Chaotic Motion: Example
Initial conditions
l1  l2  0.6m
1 (0)  1.9rad, 2 (0)  1.55rad
1 (0)  2 (0)  0
1 (0)  1.9rad, 2 (0)  1.7rad
1 (0)  2 (0)  0
Periodic motion
Chaotic motion
A transition from a stable to chaotic motion can be observed in this system
The Phase Space Trajectory
The 2s generalized coordinates represent a phase vector (q1, q2, …, qs, p1, p2, …, ps)
in an abstract 2s-dimensional space, the so-called phase space.
A particular realization of the phase vector at a given time provides a phase point in the
phase space. Each phase point uniquely represents the state of the system at this time.
In the course of time, the phase point moves in phase space, generating a phase trajectory,
which represents dynamics of the 2s degrees of freedom.
p (t )
q(t2 ), p(t2 )
q(t1 ), p(t1 )
q (t )
The Phase Space Trajectories: Examples
Examples of phase space trajectories:
(Projection to the plane x, px)
2.2 Molecular Forces
• Newtonian dynamics
• Interatomic potentials
• Molecular dynamics simulations
• Boundary conditions
• Post-processing and visualization
Newtonian Dynamics
Molecular forces and positions change with time.
In principle, an MD simulation is a solution of a system of Newtonian
equations of motion.
The Newtonian equation (the Newton’s second law) is a special case of
the Lagrangian equation of motion for mass points in a Cartesian
system.
For such a system, the Lagrangian function is given by
Z
mi 2
L   ( xi  yi2  zi2 )  U (r )
2
i
rij
ri
r   r1 , r2 ,... , ri   xi , yi , zi 
rj
X
Y
Newtonian Dynamics
By utilizing the Lagrangian equation of motion at qα= xi, qα+1= yi, qα+2= zi,
d L L

 0,
dt q q
d L
L

 0,
dt q 1 q 1
d L
L

0
dt q  2 q  2



mi xi  Fx ,i ,
mi yi  Fy ,i ,
mi zi  f z ,i
This can be rewritten in the vector form:
Fi  mi ri ,
U (r )
Fi   Fx ,i , Fy ,i , Fz ,i   
ri
Interatomic Potential
The most general form of the potential is given by the series
U (r1 , r2 ,..., rN )  W1 (ri )   W2 (ri , r j ) 
i
i , j i

i , j i , k  j
W3 (ri , r j , rk )  ...
The one-body potential W1 describes external force fields (e.g.
gravitational filed), and external constraining fields (e.g. the “wall
function” for particles in a circular chamber)
z
j
The two-body potential describes dependence of the
potential energy on the distances between pairs of
atoms in the system:
W2 (ri , r j )  W2 (r ),
r  | rij |  | ri  r j |
The three-body and higher order potentials (often ignored)
provide dependence on the geometry of atomic
arrangement/bonding.
For instance, a dependence on the angle between
three mass points is given by
r ji  r jk
W3 (ri , r j , rk )  W3 (cos ijk ), cos ijk 
|r ji | |r jk |
i
rij
ri
rj
y
x
z
rji
θijk
j
rj
x
i
rjk
rk
k
y
Pair-Wise Potentials: Lennard-Jones
  12   6 
WLJ (r )  4      
 r  
 r 
WLJ (r )

FLJ (r )  
 24
r

   13   7 
2      
  r   r  
Here, ε is the depth of the potential energy well and σ is
the value of r where that potential energy becomes zero;
the equilibrium distance ρ is given by
FLJ (  )  0    6 2
• The first term represents repulsive interaction
•
At small distances atoms repel due to quantum
effects (to be discussed in a later lecture)
• The second term represents attractive interaction
•
This term represents electrostatic attraction at
large distances
Pair-Wise Potentials: Morse
Another pair-wise potential model, effective for
modeling crystalline solids, is the Morse potential,
WM (r )   e 2  (   r )  2e  (   r ) 
FM (r )  2  e 2  (   r )  e  (   r ) 
Here, ε is the depth of the potential energy well and β
is a scaling factor; the equilibrium distance is given
by ρ
FM (  )  0
Similarly to the earlier example,
• The first term represents repulsive interaction
• The second term represents attractive interaction
Truncated Potential
• In system of N atoms
1
N ( N  1) unique pair interactions
2
•
accumulates
•
if all pair interactions are sampled, the number increases with square of the
number of atoms
• Saving computer time
•
Neglect pair interactions beyond some distance
•
Example: Lennard-Jones potential used in simulations
  
4 
WLJ (r )    r


rc
12
6
   
     r  rc
  r  
0
r  rc
Instability of Trajectories
Due to essential non-linearity of the molecular interaction, classical
equations of motion yield non-stable (non-predictable) trajectories.
Therefore, the results of MD simulations are intriguing.
Fluctuations in MD Simulation
Macroscopic properties are determined by behavior of individual molecules, in
particular:
•
Any measurable property can be translated into a function that depends on the
positions of phase points in phase space
•
The measured value of a property is generated from finite duration experiments
As a phase point travels on a hypersurface of constant energy, most quantities are not
constant, but fluctuating (discussed in more detail Week 3 lecture).
While kinetic energy Ek and potential energy U fluctuate, the value of total energy E
is preserved,
E  Ek  U  const
Fluctuations in MD Simulation…
6
averaged observable
6
observable
5
4
3
2
1
0
2
4
time
6
8
5
4
3
2
1
0
2
4
time
6
Averaging of a fluctuating value F(t) over period of time t1 to t2 is accomplished according to
t
1 2
F (t ) 
F (t )dt
t2  t1 t1
8
Periodic Boundary Conditions
Simulated system encompasses boundary conditions
Periodic Boundary conditions for particles in a box
•
Box replicated to infinity in all three Cartesian directions;
Motivation for periodic boundary conditions: domain
reduction and analysis of a representative substructure only.
Particles in all boxes move simultaneously, though only one
modeled explicitly, i.e. represented in the computer code
•
•
•
Each particle interacts with other particles in the box and with images in
nearby boxes
Interactions occur also through the boundaries
No surface effects take place
Periodic Boundary Conditions
Original box
Translation
image
boxes
Periodic Boundary Conditions: Example
Example simulation of an
atomic cluster with periodic
boundary conditions.
A particles, going through a
boundaries returns to the box
from the opposite side:
This model is equivalent to a
larger system, comprised of the
translation image boxes:
Dynamic Molecular Modeling
Develop Model
Model Molecular
Interactions
Molecular Dynamics Simulation
Develop Equations
of Motion
Initialize
Parameters
Initialization
Initialize Atoms
Analysis of Trajectories
Generate Trajectories
Static and dynamic
(macro) properties
Re-initialization
Predictions
Initialization
• Decisions concerning preliminaries
•
Systems of units in which calculations will be carried out
•
Numerical algorithm to be used
•
Assignment of values to all free parameters
• Initialization of atoms
•
Assignment of initial positions
•
Assignment of initial velocities
Post-Processing
After simulation is completed, macroscopic properties of the system are
evaluated based on (microscopic) atomic positions and velocities:
1. Macroscopic thermodynamic parameters
- temperature
- internal energy
- pressure
- entropy (will be discussed in week 3 lecture)
2. Thermodynamic response functions, e.g. heat capacity
3. Other properties (e.g. viscosity, crack propagation speed, etc.)
Temperature and Potential Energy
Absolute temperature is proportional to the average kinetic energy
2
T
E kin ,
fk
E
kin
m( xi2  yi2  zi2 )
1
 
N i
2
k Boltzman constant
N total number of atoms
f number of DOF per atom
The time-averaged potential energy
U int
t

1 
 U   W1  ri ( )   W2  rij ( )   ... d
t 0 i
i j i

W1 one-body potential, W2 two-body potential,
t total time
Pressure and Mean Square Force
• Pressure is defined by the atomic velocities
mN 2
mN 2
2 N kin
Px 
vx ,
Px 
vx 
Ex
V
V
V
V - volume, m - mass, N - number of particles
• Mean square force is defined by the time-averaged derivative of the
potential function
F12 
 W (r1 j ) 
j 1
2
, W 
u (r )
r
Entropy
For an adiabatically isolated system, the entropy is related to the
phase volume integral:
S  k ln 
k Boltzmann constant
 phase volume integral (over all allowed states)

1
  E  E kin (p)  U (r ) dr1...drN dp1...dp N
3N

h N!
0, x  0
 ( x)  
1, x  0
h
Planck's constant
In greater detail, the properties and calculation of entropy for various systems will be
discussed on weeks 3-5 lectures.
Thermodynamic Response Functions
The TD response functions reveal how simple thermodynamic quantities respond to
changes in measurables, usually either pressure or temperature. Thus, they are derivative
quantities (coefficients):
• Heat capacity (how the system internal energy responds to an isometric change
in temperature):
 E 
Cv  


T

V
• Thermal pressure (how pressure responds to an isometric change in
temperature):
 P 
v   
 T V
• Adiabatic compressibility (how the system volume responds to an isentropic,
change in pressure):
1  V 
s   

V  P S
Equilibration
In order to evaluate the averaged macroscopic parameters, the simulated system must
achieve a thermodynamic equilibrium. Indeed, the thermodynamic parameters, such as
temperature, internal energy, etc., are defined for equilibrium systems.
In the equilibrium:
1) the macroscopic parameters fluctuate around their statistically averaged values;
2) the property averages are stable to small perturbations;
3) different parts of the system yield the same averaged values of the macroscopic
parameters.
Equilibration of a macroscopic parameter is achieved in distinct ways in closed adiabatic
and isothermal systems (surrounded by a heat bath).
Closed systems: value of the macroscopic parameter fluctuates about the averaged value
with a decaying fluctuation amplitude.
Isothermal systems: value of the macroscopic parameter both fluctuates, and also
asymptotically approaches the averaged statistical value (examples are to follow).
Adiabatic Example: Interactive Particles in a Circular Chamber
Repulsive interaction between the particles and the
wall is described by the “wall function”, a one-body
potential that depends on ri – distance between the
particle i and the chamber’s center):
Wwl (ri )   e
  ( R ri )
 e

  R  xi2  yi2

Interaction between particles is modeled with the
two-body Lennard-Jones potential (rij – distance
between particles i and j):
  12  6 
WLJ (rij )  4  12  6 
r
rij 
 ij
The total potential:
U   Wwl (ri )   WLJ (rij ),
i
i
j i
ri  xi2  yi2
rij  ( xi  x j ) 2  ( yi  y j ) 2
y
ri
rij
R
x
rj
Three Particles: Equation of Motion and Solution
The total potential:
U  Wwl (r1 )  Wwl (r2 )  Wwl (r3 )
 WLJ (r12 )  WLJ (r13 )  WLJ (r23 ),
ri  xi2  yi2
rij  ( xi  x j ) 2  ( yi  y j ) 2
Equations of motion:
mi xi  
U
U
, mi yi  
, i  1, 2,3
xi
yi
Parameters:
  4nm1 , R  10nm, m  109 kg
Initial conditions (nm, m/s):
x1 (0)  0, x1 (0)  25, y1 (0)  3.0, y1 (0)  0
x2 (0)  0, x2 (0)  30, y2 (0)  0.5, y2 (0)  0
x3 (0)  0, x3 (0)  20, y3 (0)  2.5, y3 (0)  0
Five Particles: Equations of Motion and Solution
The total potential:
U  Wwl (r1 )  Wwl (r2 )  Wwl (r3 )  Wwl (r4 )  Wwl (r5 )
 WLJ (r12 )  WLJ (r13 )  WLJ (r14 )  WLJ (r15 )
 WLJ (r23 )  WLJ (r24 )  WLJ (r25 )
 WLJ (r35 )  WLJ (r35 )
 WLJ (r45 )
Equations of motion:
mi xi  
U
U
, mi yi  
, i  1, 2,3, 4,5
xi
yi
Parameters:
  4nm1 , R  10nm, m  109 kg
Initial conditions (nm, nm/s):
x1 (0)  0, x1 (0)  25, y1 (0)  5.6, y1 (0)  0
x2 (0)  0, x2 (0)  30, y2 (0)  3.0, y2 (0)  0
x3 (0)  0, x3 (0)  20, y3 (0)  0.5, y3 (0)  0
x4 (0)  0, x4 (0)  24, y4 (0)  2.5, y4 (0)  0
x5 (0)  0, x5 (0)  22, y5 (0)  4.9, y5 (0)  0
Post-Processing: Kinetic Energy, Temperature and Pressure
Averaged kinetic energy vs. time (three particles)
Time averaged kinetic energy of
particles is approaching the value
E kin  3.057 1025 J
which corresponds to temperature
1 kin 3.057 1025
T E 
k
1.38 1023
0.02K
Note: a low temperature system was chosen in
order to observe the real-time atomic motion.
Pressure in the system is due to the
radial components of velocities:
P
mN 2
vrad , vrad  vx cos   v y sin   P  3.23 109 Pa
V
Kinetic energy, and therefore temperature and pressure are due to motion of the particles.
Post-Processing: Potential Energy
Averaged potential energy vs. time (three particles)
Time averaged potential
energy of the system is
approaching the value
U  0.480 1025 J
that gives the potential energy
of the system.
Potential energy is due to interaction of particles with each other and with
external constraining fields.
Post-Processing: Total Energy
Kinetic, potential and total energy vs. time (three particles)
E E
tot
E
kin
kin
E
pot
m 3 2 2
  ( xi  yi )
2 i 1
E pot  U
The total energy (solid red line): 9.65x10-25 J. This value does not vary in time
Isothermal Example: 1D Lattice with a “Cold” Region
“Cold” region
y
Total number of particles is large.
Interaction between particles is modeled with the two-body harmonic potential:
Wh (rij ) 
1 2
k rij
2
here, rij = yi – yj is the relative distance between particles i and j in the vertical
(y-axis) direction, k – linear interaction coefficient (similar to spring stiffness).
Several atoms in the middle of the chain (between the yellow dashed lines)
represent the initially “cold” region. Initial velocities and displacements for
these atoms are zero. The remaining atoms have randomly distributed initial
velocities and displacements.
1D Lattice: Equations of Motion and Solution
Number of particles simulated: 50.
Boundary conditions are periodic, so that the coupling between the 50th and 1st particles is
established. For a more symmetric view, the simulation shows the 1 st particle at both ends of
the lattice.
The total potential (the last term is due to periodic boundary conditions):):
49
U   Wh ( yi  yi 1 )  Wh ( y1  y50 ), Wh ( yi  y j ) 
i 1
Equations of motion:
m yi  
1
k ( yi  y j ) 2
2
U
, i  1, 2,...50
yi
Parameters:
m  1021 kg, k  1021 N/m
Initial conditions:
middle atoms 22..27 : yi (0)  yi (0)  0,
other atoms : yi (0), yi (0)  random within [  1,1]
Post-Processing: Kinetic Energy and Temperature
Averaged kinetic energy per particle vs. time
(for the initially “cold” subsystem)
Time averaged kinetic
energy of particles in the
“cold” subsystem is
approaching the value
 0.286 10
21
J
which corresponds to
temperature
2
T  E kin
k
2  0.286 1021

1.38 1023
Kinetic energy, 10
E
kin
21
J
0.4
0.3
0.2
Time averaged
0.1
Equilibrium value
41.4K
0
100
200
300
400
Time, s
In contrast to the adiabatic system example, the kinetic energy for the open isothermal
subsystem both fluctuates and approaches asymptotically the statistical average.
500
Post-Processing: Potential Energy
U  0.480 1025 J
Potential energy, 10
Time averaged potential
energy of the “cold”
subsystem is approaching
the value
21J
Averaged potential energy vs. time (“cold” subsystem)
2
1.5
1
Time averaged
0.5
Internal energy asymptote
0
100
200
300
400
Time, s
In contrast to the adiabatic system example, the internal energy for the open isothermal
subsystem both fluctuates and approaches asymptotically the statistical average.
500
Limitations of MD
• Realism of forces
•
Simulation imitates the behavior of a real system only to the extend that
interatomic forces are alike to those that real nuclei would experience when
arranged in same configuration
•
In simulation forces are obtained as a gradient of a potential energy function,
which depends on the positions of the particles
•
Realism depends on the ability of potential chosen to replicate the conduct of
the material under the circumstance at which the simulation is governed
Limitations of MD
• Time Limitations
•
Simulation is “safe” when duration of the simulation is much greater than
relaxation time
•
Systems have a propensity to become slow and sluggish near phase transitions
•
Relaxation times order of magnitude larger than times reachable by simulation
Limitations of MD
• System Size Limitations
•
Correlation lengths can increase or deviate near phase transitions when
comparing the size of the MD cell with one of the spatial correlation functions
•
Partially assuaged by Finite Size Scaling
•
Compute physical property, A, using many different size boxes, L,
then relating the results to:
A( L)  A0 
and using
•
A0 , c, n
c
Ln
as fitting parameters
where A0  lim L A( L) , and should be taken as the best
estimate for true physical quantity
2.3 Molecular Dynamics Applications
Tensile failure of a gold nanowire
Carbon nanotube immersed into
monoatomic helium gas
Dislocation dynamics of nanoindentation
Deposition of an amorphous carbon film
Molecular Dynamics Applications (II)
MD fracture simulations
Shear-dominant crack propagation