Classical Molecular Dynamics

1
Classical Molecular Dynamics
Computational Chemistry 5510 Spring 2006 Hai Lin
2
Motions of Atoms
The motions of atoms are determined by classical
mechanics (Newtons Laws).
F ma
http//ca.geocities.com/spacephysicsisu/newton.htm
l
3
A Deterministic Scheme
The state of the system at any future time can be
predicted from its current state. mi(d2ri/dt2)
- dV/dri ri(t 0) ri0 dri/dt (t 0)
vi0 The positions and momenta of the atoms at t
t1 are completely determined by the positions
and momenta at t 0.
Initial conditions
4
Generate the Trajectory
Solve the differential equations numerically A
small finite time step Dt. Conserve the
energy Kinetic and potentail energies can
fluctuate, but the total energy should be
conserved ( 0.01 variation). Computational
effort Generate 104 to 106 or more points in
the phase space at a reasonably short time.
5
Verlet Algorithm
Taylor expansion of postion ri 1 ri vi Dt
½ ai (Dt)2 (1/6) bi (Dt)3 ri - 1 ri - vi Dt
½ ai (Dt)2 - (1/6) bi (Dt)3 ? Position ri
1 2ri - ri 1 ai (Dt)2 Velocity vi
ri1 ri-1 / 2Dt Acceleration ai Fi / mi
To start, one needs r-1 ! Can be estimated by
r-1 r0 v0Dt There are other algorithms to
propagate the trajectory, but the Verlet
algorithm is quite popular.
6
Let us Do the Simulation
Describe the atoms, functional groups, and
interactions.
System Representation
Select the ensembles, boundary treatments,
initial configurations, time steps, and
constraints.
Simulation Setup
Heat the system up to the desired temperature and
let it reach equilibration ( 104 to 105 steps).
Heat up Equilibration
Let the system run ( 105 to 106 steps),
producing trajectories where representative
configurations are extracted.
Production Phase
Results and Error Estimation
Calculate thermodynamic quantities, estimate the
statistical errors.
7
Represent the System
Atomic level Accurate, but expensive United
Atoms Represent a functional group as a whole

• Interactions
• Greatly simplified Hard ball
• Realistic Van der Waals electrostatic
• Advanced Polarization included

8
Boundary Effect
Experiments measure macroscopic or bulk
properties. (Avogadro constant NA 6.0221023
mol1)
A liter of water one in 1.5 million water
molecule is affected by the interaction with the
walls of the container. A drop of water having
1000 molecules the majority of the waters are
within the influence of the boundary.
Surface region
Bulk region
Can we minimize the boundary effect and make our
simulations more relevant to experimental
measurements?
9
Periodic Boundary

• The cell is duplicated in x, y, and z directions.
• A molecule leaving the cell will enter the cell
  from the other side.

• Cut-off for non-bonded interactions should not
  let an atom see its image --- they are indeed the
  same atoms!
• Results are quasi-periodic. You might not want it
  sometimes.

B'
B
A'
A
10
Non-Periodic Boundary
Active region normal MD simulations Buffer
region atoms fixed or restrained to certain
positions
Not as safe as the periodic boundary treatment,
but it reduces the number of atoms in simulation.
Active region
Buffer region
11
Non-bonded Interactions

• Non-bonded interactions are expensive to
  calculate (scaled as N 2).
• Approaches to obtain a better scaling behavior
• Use of Cut-off
• Evaluate exactly for near-field contributions but
  approximate the far-field contributions

• Fast Multipole Moment (FMM) method
• Divide the system into boxes, each box containing
  a number of atoms
• Approximate the far-field interactions by
  interactions between multi-poles centered at the
  boxes.
• The farther away the atoms are, the larger box
  can be used
• Scaled closed to linearly ( N)!

12
Time Step size
The time step size is determined by the fastest
process in the system. Too small expensive,
subject to round-off errors Too large
inaccurate, producing artefacts
Too small Too large Just fine
13
SHAKE
Stretches of X-H bonds are fast motions, but
usually do not have significance in our
simulations. It is convenient to freeze the X-H
stretches, i.e., to constrain their bond lengths
to equilibrium bond lengths during the
simulations, so that one can use a larger time
step.
i
j
sij (ri rj)2 dij2 0
d
The atoms are first allowed to move under the
influence of the forces, then their positions are
adjusted iteratively to satisfy the constraints
to an acceptable accuracy. Sometimes also apply
to bends, but rarely to torsions.
14
Heat up and Equilibration

• Example
• A protein solvated in water
• Start from the X-ray structure.
• Surround the protein by water.
• Rise the temperature up gradually (e.g., 300 K).
• Let the system run freely for a while (e.g., 30
  ps).

15
Monitoring Equilibration
How can I know if the system is in equilibration?
Unfortunately, you can never really know. But
you can show that the system is likely in
equilibration.

• Monitoring several quantities over
• time
• Energy
• Temperature
• Pressure
• Volume
• Density
• Mean squared displacement
• Dr2(t) (1/N) Si1,,N ri(t) - ri(0)2

t ps
16
Production Run Result Analysis

• Production Run
• After the system has reached equilibration,
  continue the MD rum for 105 to 106 steps.
• Result Analysis
• View the trajectories.
• Calculate thermodynamics quantities
• Energy, Heat capacity, density, Redial
  distribution function, Vibrational spectra,
• Estimate the errors and reliability.
• Compare with experimental results.

17
Thermodynamics Quantities

• Kinetic Energy
• Internal Energy
• Temperature
• Pressure
• Heat Capacity
• Radial Distribution Function
• Correlation Function

18
Kinetic Energy
For a given point on the trajectory in the phase
space, the kinetic energy is Ek ½ Si1,,N mi
vi2 where the sum is over N atoms.
The ensemble averaged kinetic energy is
? Ek ?M (1/M) Sa 1,,M Ek,a
where average is over M points on the MD
trajectory (or M states generated by the MD
simulation).
19
Internal Energy
For a given point on the trajectory in the phase
space, the internal energy is U Siltj Eij where
the sum is over all pairwise interactions.
The ensemble averaged internal energy is
? U ?M (1/M) Sa 1,,M Ua
where average is over M points on the MD
trajectory (or M states generated by the MD
simulation).
20
Temperature
Temperature is related to the ensemble averaged
kinetic energy by ½ (3N Nc) kB T ? Ek
?M where N is the number of atoms and Nc is the
number of constraints. Typically we require the
total linear momentum of the system is
constrained to zero (the center of mass of the
system does not move), and Nc is 3.
21
Pressure
Pressure is related to the product of the
positions and forces (for pairwise
interactions) PV NkBT (1/3) ? Siltj rij fij
?M ideal gas contribution where N is the
number of atoms, rij is the distance between a
pair of interacting atoms, fij is the
corresponding force, and the sum is over all
pairwise interactions.
22
Heat Capacity

• Calculate internal energies at different
  temperatures and take the partial derivatives
• CV (?U / ?T)V
• (U2 U1) / (T2 T1) at constant V
• Calculate the fluctuation of internal energy
  around its mean value
• kBT2 CV ? (U ?U?M)2 ? ? U2 ?M ?U?M2
• It requires a longer simulation time for one
  simulation at one temperature. (A trade-off!)

23
Radial Distribution Function
A radial distribution function measures the
probability of finding a particle as a function
of distance from a given particle. g(r, Dr) ? ?
N(r, Dr) ?M / 4p r2Dr
Number of particles between r and r Dr from the
given particle
The volume of a spherical shell with thickness Dr
r
24
Correlation Function

• A correlation function measures the relationship
  between two variables
• Cxy ? x y ?M / (? x2 ?M ? y2 ?M) ½
• If x (or y) fluctuate about a non-zero mean
  value, replace x (or y) in the above equation by
  x ? x ?M (or y ? y ?M).
• If x y, Cxx is called an auto-correlation
  function.

25
Time Correlation Function
Correlation functions are often time-dependent,
and they can be normalized by initial values at t
t0 Cxy (t) ? x(t) y(t0) ?N / ? x(t0) y(t0)
?N

• The average is over the collection of N particles
  at a given time t. The value often drops from 1
  at t t0 to 0 at t ?.
• Fourier transform of the auto-correlation
  function is sometimes related to experimentally
  measured quantities.

Cxy
1
0
t
26
Estimate the Error

• Two kinds of errors
• Statistical Error
• Err ? 1 / N ½
• A 1 error requires N 10,000
• A 0.1 error requires N 1,000,000
• Systematic Error
• Inaccurate model
• Finite precision in calculations
• Equilibration not reached

27
Summary

• Generate the Trajectory
• (Verlet Algorithm)
• Boundary Effect
• (Periodic and non-periodic boundary treatments)
• Non-bonded Interactions
• (Cut-off and approximated far-field
  contributions)
• Time Step Size
• (SHAKE)
• Heat up, Equilibration, and Production
• Thermodynamic Quantities
• Error Estimations

28
Your Homework

• Read the slides
• Read textbook (take notes when you read!)
• 16.1 (excluding 16.1.1, 16.1.2, and 16.1.3)
• 16.2 (up to page 387)
• Questions
• Point out two differences between MC and MD
  simulations.
• What is a periodic boundary?
• Supposed that I want to study protein folding.
  Can I use the SHAKE algorithm? Why?
• Could you tell me how I should estimate the error
  bar in my MD simulations?

29
Practise with Tinker

• An anion in a water box
• Make your own subdirectory (e.g., hailin), and
  copy two files anion.xyz and anion.key from the
  example subdirectory to your own subdirectory.
• Use Tinker to load the anion.xyz, the anion.key
  will be loaded automatically.
• Set the color mode to partial charge. Set the
  display mode for the ion to Spacefill and for the
  water molecules to Wireframe.
• Go to the Keyword Editor section
• Set Output Control ARCHIVE and VERBOSE.
  (Important!)
• Check the Active Keywords (No change is needed at
  the moment).
• Go to the Modeling Commands section, set
• Number of Dynamic Steps 1000
• Time Step Length in Femtoseconds 1.0
• Time between Dumps in Picoseconds 0.1
• Ensemble NVT
• Temperature in Kelvin 300

30
Practise with Tinker (2)

• An anion in a water box (2)
• Launch the dynamics run
• At the bottom the Graphic window, monitor the
  time elapsed and energy.
• When it is done, read the log file.
• Check how the energies, temperature, and pressure
  change over time.
• Check the statistical errors for the
  thermodynamic quantities and their variations
  over time.
• Check the system linear velocity and
  translational kinetic energy. (Does the center of
  mass move?)
• Go to the Graphic window, set the display mode
  for anion.xyz to invisible.
• Load anion.arc. Set the color mode to partial
  charge. Set the display mode for the ion to
  Spacefill and for the water molecules to
  Wireframe.
• Go to the Trajectory menu, and click on the Play
  button.

31
Practise with Tinker (3)

• An anion in a water box (3)
• Stop the trajectory play, and close the anion.arc
  file.
• Go to your own subdirectory, rename anion.arc to
  anion_NVT.arc.
• Now redo the dynamics again with the NVE option,
  using the procedure described as steps 2 to 12.
• Number of Dynamic Steps 1000
• Time Step Length in Femtoseconds 1.0
• Time between Dumps in Picoseconds 0.1
• Ensemble NVE
• Temperature in Kelvin 300
• When it is done, read the log file. Compare the
  results (e.g., the statistical errors) with the
  NVT simulation you previously done.
• Load anion.arc. Set the color mode to partial
  charge. Set the display mode for the ion to Ball
  Stick and for the water molecules to Wireframe.
• Go to the Trajectory menu, and click on the Play
  button.

32
Your Homework

• If you have not yet finished the practice about
  an anion in a water box (3), please finish it.
• Practise 50 butane molecules
• Copy the butanes.xyz and butanes.key files from
  the test subdirectory to your own directory.
• Go to Keyword section,
• Set Output Control ARCHIVE and VERBOSE.
• Check the Active Keywords (No change is needed at
  the moment).
• Go to Modeling Commands section, redo the
  dynamics again using the NVT or NVE options.
• When it is done, read the log file. Use Excel or
  what ever programs that you have to plot the
  energies and temperature as functions of time.
• View the Trajectory.
• Change the length of the cubic cell (A-AXIS) from
  20 Å to 100 Å, and redo the simulations. What do
  you found? Why? Is it reasonable? Is it an
  artefact? Or does the program has a bug?