molecular dynamics:
integrate Newton's equations of motion for all atoms in the system
Applications of molecular dynamics
-conformational searches
-generate statistical ensembles
to calculate energetic, thermodynamic, structural, and dynamic (time-dependent) properties
-study motions of molecules (i.e. time evolution)
Newton's equations of motion:
: force on particle i
: mass of particle i
: acceleration of particle i
V: potential energy
trajectory: coordinates and velocities for a molecular dynamics run
Properties of Newton's equations of motion
- Energy conservation: V+K is constant
- Equations of motion are time reversible: if we change the signs of all velocities,
the atoms will retrace their trajectories
- Numerical methods are required to solve these equations for systems larger
than 2 independent particles
A. Finite difference methods
-numerically solve Newton's equations of motion
-positions and velocities at time t 
positions and velocities at a later time t+dt
Desirable properties of an integrator
-fast
-requires little computer memory
-permits the use of a relatively long time step
-duplicates classical trajectory as closely as possible
-conserves energy and momentum
-time-reversible
-simple in form and easy to program
No method will duplicate the classical trajectory for a long time:
-any 2 trajectories that are initially very close will eventually diverge from one another exponentially with time
-a small perturbation will cause a computer-generated trajectory to diverge from the true classical trajectory
Molecular dynamics has 2 roles
1. need essentially exact solutions of equations of motion for times comparable
with correlation times of interest to accurately calculate time correlation functions
2. use the method to generate states sampled from the microcanonical ensemble;
do not need exact classical trajectories to do this, but energy conservation is
crucial
A shorter time step should be used for
-higher temperatures (i.e. larger velocities of particles)
-more rapidly varying potential surfaces
-lighter mass of particles
appropriate time step energy conservation, time-reversibility
1. Verlet algorithm
Taylor expansion about 
:
Add these 2 equations up:
Note: velocities do not appear at all. They can be calculated from:
The errors in Eq. 21 are of order 
,
and the velocities (Eq. 22) are subject to errors of order
.
Algorithm:
Start with and 
Repeat these steps:
1. calculate 
(from the force
)
2. use Eq. 21 to calculate 
3. calculate 
if desired
4. replace 
with
and with 
Advantages:
-the advancement of positions is all in one step (using Eq. 21)
-compact and simple to program
-time-reversible
-conserves energy well even with relatively long time steps
Disadvantages:
-the velocities at t can be calculated only after are known
-to start a trajectory, need to know and
rather than and
2. Velocity Verlet algorithm
The Verlet algorithm can be recovered by eliminating the velocities.
Algorithm:
Start with and and calculate
Repeat the following steps:
1. calculate using Eq. 23
2. calculate velocities at mid-step using
3. calculate 
4. complete the velocity move using
Advantages:
-kinetic energy at time 
is available
-can start with positions and velocities at time t (i.e.
and )
-numerically stable
-simple to program
-time-reversible
-conserves energy well even with relatively long time steps
Disadvantages:
-the advancement of velocities takes 2 steps
Storage requirements of Verlet and Velocity Verlet identical (9N).
Other methods: Gear predictor-corrector, Beeman
3. Time step
Larger time step decrease computation time
Too large a time step instability and inaccuracy in numerical integration
Velocity Verlet and Verlet algorithms assume that velocities and accelerations are constant over a given time step
the time step is limited by the highest-frequency motion that must be simulated since fast vibrations imply rapidly changing velocities and accelerations
A vibrational period should be split into at least 
segments to satisfy
the assumption that velocities and accelerations are constant over
a given time step
Example: if the fastest vibration is C-H stretch, which
has a period of about 10-14 seconds, then the time step should be about
10-15 seconds (1 fs)
Test the stability of numerical integration
-energy conservation
-time reversibility
B. Constraint dynamics
Constrain some intramolecular bond lengths and/or angles to fixed values
Reasons for use:
-intramolecular bond vibrations are typically the highest frequencies in the system and therefore determine the largest time step that can be used if bonds are constrained, a larger time step can be used, which speeds up the computation
-the quantum mechanical nature of bond vibration makes the classical description questionable
1. Fundamental concepts of constraint dynamics
where 
is a constraint force
Use a set of undetermined multipliers to represent the magnitudes of forces
directed along the bonds to keep the bond lengths constant. Solve the
equations of motion for one time step in the absence of the constraint
forces, and then determine the magnitudes of the constraint forces and
correct the atomic positions (and/or velocities)
Iterative method:
Go through constraints one by one, cyclically, adjusting coordinates to satisfy
each constraint in turn until all are satisfied to within a specified
tolerance.
SHAKE (Verlet algorithm): only atomic positions are constrained
RATTLE (velocity Verlet): both positions and velocities are constrained
2. RATTLE algorithm
H.C. Andersen, J. Comp. Chem. 52, 24 (1983).
Bond constraints:
where 
is the fixed distance between atoms i and j (whose coordinates are
and 
).
The time derivatives of the constraint equations give the velocity constraints:
The equations for constrained dynamics are
where
is the force due to unconstrained interactions and
is the force on atom i due to the constraints:
where the prime indicates a sum over only atoms j that are
connected to atom i by a constraint, 
indicates the gradient with respect to ,
and 
are Lagrange multipliers associated with the constraints.
The RATTLE equations of motion are analogous to the velocity Verlet equations, with
the addition of the constraint forces:
where 
.
The Lagrange multipliers 
are chosen so that the
constraint equations 25 are satisfied at time .
The Lagrange multipliers 
are chosen so that the time derivatives of the constraint equations
26 are satisfied at
time .
Algorithm (slightly simplified for clarity)
Start with and and calculate
Repeat the following steps:
1. calculate in the absence of the constraints (as in Eq. 21 )
2. calculate velocities at mid-step using
3. start the iterative loop for calculating the Lagrange multipliers
(change so that all bond constraints (Eq. 25) are satisfied within tolerance)
a. pick a constraint involving atoms i and j
b. if the current bond length differs from the fixed bond length
by more than the prescribed tolerance, then determine the
value of necessary to satisfy the constraint more closely
and add this term to 
,
and 
,
according to Eqs.29 and 30, respectively;
otherwise go back to step 3a and pick a new constraint
Continue this iterative procedure until all bond constraints are satisfied
within the prescribed tolerance
4. calculate
5. complete the velocity move in the absence of constraints using
6. start the iterative loop for calculating the Lagrange multipliers
(change so that all velocity constraints (Eq. 26) are satisfied within tolerance)
a. pick a constraint involving atoms i and j
b. if the dot product of 
and 
differs from zero by more than the
prescribed tolerance then determine the value of
necessary to satisfy
the velocity constraint more closely
and add this term to 
,
according to Eq. 30;
otherwise go back to step 6a and pick a new constraint
Continue this iterative procedure until all velocity constraints are satisfied
within the prescribed tolerance
Constraining bond lengths
-time step can be increased
-overhead for iterative procedure is not significant
speeds up computation
Constraining angles
-time step can be increased
-overhead for iterative procedure is significant (requires longer time to converge)
typically does not speed up computation (and sometimes slows it down)
Thus, constraining bond lengths is worthwhile, but bond angles should be free to evolve
C. Temperature: Maxwell-Boltzmann distribution of velocities
Temperature:
-proportional to the kinetic energy of the system, which can be expressed
in terms of the atomic velocities
-a measure of an average distribution
-meaningful (thermodynamically) only at equilibrium
Boltzmann:
For a dilute gas,
: probability of finding a molecule with energy 
: degeneracy,
(i.e. the number of different states of a molecule
corresponding to energy )
k: Boltzmann's constant
This applies to any type of energy - translational kinetic energy, vibrational, rotational, electronic, etc.
Maxwell:
derived special case of Boltzmann equation for translational kinetic energies ( 
)
P(v) dv = the probability of finding a molecule with speed between v and dv
F(v) dv = the fraction of molecules with speed between v and dv
(where molecule has mass m and temperature T)
Note: 
If we do the integral in the denominator (the normalization factor), then we get the
Maxwell-Boltzmann equation:
Increasing the temperature (or decreasing the mass)
-a broader distribution of speeds
-a shift to higher speeds
The x, y, and z components of the velocities have a Gaussian distribution:
Calculate average speed ( 
indicates the average of quantity x):
For a system of N atoms and no internal constraints:
Kinetic energy: 
Average kinetic energy: 
Equipartition principle:
there is an average energy of kT/2
associated with each degree of freedom that appears as a squared term in the Hamiltonian
(In this case with no constraints, there are 3N velocity degrees of freedom, and
each one appears as 
in the kinetic energy K of the Hamiltonian)
If there are N atoms and 
internal constraints, then the number of degrees of freedom is 
.
In this case,
instantaneous kinetic temperature:
thermodynamic temperature T:
average of the instantaneous kinetic temperature over many molecular dynamics
time steps ( 
)
C. Initialization and equilibration
Initialization
initial coordinates:
can be obtained from a previous simulation
Otherwise:
1. liquid
-start out in lattice fcc configuration with appropriate density and random
orientation of dipole moments
-equilibrate to melt lattice
2. proteins
-start with X-ray structure
-minimize energy to eliminate any large forces due to repulsive van der Waals contacts or poor geometries
-place constraints on all atoms (tethering) to prevent large shifts from starting
constraints due to very large forces
initial velocities:
Gaussian distribution of velocities 
, 
, and 
appropriate to a given temperature
(i.e. Maxwell-Boltzmann distribution)
Procedure to generate a Maxwell-Boltzmann distribution of velocities:
1. use a random number generator to generate a normal distribution with zero mean and unit variance (i.e. a set of 3N numbers corresponding to such a normal
distribution)
2. multiply by 
to get proper distribution of velocities
3. correct so there is no overall momentum
There is no unique set of initial velocities - this depends on the random number seed
An ensemble of such initial velocity sets satisfies
Equilibration
preliminary molecular dynamics simulation to equilibrate the system at the desired T
Even if initial velocities satisfy a Maxwell-Boltzmann distribution at the desired
temperature, the distribution will change during simulation, especially
if initial configuration is far from equilibrium
Procedure for equilibration:
1. Generate initial velocities for specified temperature
2. Perform molecular dynamics simulation, scaling velocities occasionally
to reach the desired temperature
Test if system is equilibrated:
- Plot thermodynamic quantities such as potential energy, temperature, etc. as a function of time. If equilibrated,
1. quantities should fluctuate around average values
2. average should remain constant over time (no systematic drift)
-Different initial conditions (spatial coordinates and velocities) should produce the same thermodynamic quantities after equilibration
Can first run at higher temperature than desired to ensure full equilibration
E. Temperature control
Uses:
-During equilibration: to attain the desired temperature
-During data collection: to maintain the correct ensemble (i.e. canonical NVT)
1. Velocity scaling
velocities of all atoms scaled to match the target temperature
: target temperature
: either the instantaneous temperature (i.e. the temperature with 
) or an average temperature
Useful during equilibration to attain the desired temperature but does
not create a canonical ensemble
2. Andersen method
H.C. Andersen, J. Chem. Phys. 72, 2384 (1980)
Periodically replace all velocities with velocities selected from a Maxwell- Boltzmann distribution for the desired temperature.
-Useful for studying equilibrium properties
(generates a canonical ensemble)
-Problematic for studying dynamical quantities
sudden changes in the velocities discontinuous trajectories
3. Nosé-Hoover dynamics
S. Nosé, J. Chem. Phys. 81, 511 (1984);
W.G. Hoover, Phys. Rev. A 31, 1695 (1985).
-produces a true canonical ensemble both in coordinate and momentum space
-smooth, deterministic, and time-reversible trajectories
Method:
-add an additional (fictitious) degree of freedom to the real physical system
-solve the equations of motion for the extended system (real plus fictitious)
-the constant-energy dynamics (microcanonical ensemble NVE) of the extended system produces the canonical ensemble (NVT) of the real physical system
The extended Hamiltonian that will accomplish this is
where the first two terms make up the Hamiltonian of the real system (with potential
energy 
) and the last two terms are the kinetic and potential energies, respectively,
of the fictitious degree of freedom.
: coordinates of real system
: momenta of real system
S: the fictitious coordinate
Q: the fictitious ``mass"
: the fictitious momentum
g: the number of degrees of freedom in the real system
: 
After some algebraic manipulations (which eliminate the fictitious coordinate S), it can
be shown that the equations of motion for this extended Hamiltonian are:
Q is an arbitrary parameter that should be chosen carefully
too small a Q
the frequency of the fictitious degree of freedom is too high
a smaller timestep is required for the numerical integration
too large a Q
there is not enough energy exchange between the real system and the heat bath (the fictitious degree of freedom)
thermalization is inefficient
F. Ensembles
microcanonical (NVE)
-no temperature or pressure adjustment
-cannot use to equilibrate to desired temperature or pressure
-can use during data collection
-energy will drift due to roundoff and truncation errors in integration
canonical (NVT)
-control temperature using velocity scaling during equilibration and Nos dynamics (or temperature-bath coupling) during data collection
-can use to equilibrate to desired temperature
isothermal-isobaric (NPT)
-control both temperature and pressure
-can use only with periodic boundary conditions
-can use to equilibrate to desired temperature and pressure
Ensemble averages should be identical for all ensembles
 |
 |
 |
 |
 |
 |
 |
 |
 |
