Module 5:

Molecular Dynamics Simulations

Systems
 Why Not Quantum Mechanics?
• Modeling the motion of a complex molecule by solving the wave functions of
the various subatomic particles would be accurate…

− 2 2
  + U ( x, y, z ) ( x, y, z ) = E ( x, y, z )
2m

• But it would also be very hard to program and take more computing power
than anyone has!
 Classical or Molecular Mechanics

• Instead of using Quantum mechanics, we can use classical Newtonian mechanics

to model our system.

• Calculate the structure and energy of molecules based on nuclear motions.

• electrons are not considered explicitly

• Born-Oppenheimer approximation of Schrödinger equation

➢ nuclei are heavier and move slower than electrons

➢ nuclear motions (vibrations, rotations) can be studied separately

➢ electrons move fast enough to adjust to any nuclei movement

 Molecular Mechanics

Molecular mechanics treats a molecule as a collection of

weights connected with springs, where the weights represent the
nuclei and the springs represent the bonds.

• Collection of atom types (to define the atoms in

a molecule), parameters (for bond lengths, bond
angles, etc.) and equations (to calculate the
energy of a molecule)
 MM: Atom Types

In MM, a given element may have several atom types.

• sp3-Hybridized carbons have a tetrahedral bonding geometry
• aromatic carbons have a trigonal bonding geometry.
• C-C bond in the ethyl group differs from a C-C bond in the phenyl ring
• C-C bond between the phenyl ring and the ethyl group differs from all other C-C bonds in ethylbenzene.
 Force
V ( R) = Ebonded + Enon−bonded Bonded Terms

Bond Stretch

Bond Angle

1 1 2
E str = k str ,ij (r ij − r 0)2 Ebend = k b , ijk ( θijk − θ0)
2 2

Dihedral

Improper
Dihedral

1
E d = k d, ijkl (1+cos (nΦ− Φs)) Eid = k id ,ijkl (Φijkl − Φ0)2
2
Non-bonded Terms

van der Waals Interaction Electrostatic Interaction

q1q2
F=
4 0 r 2

[( ) ( )]
12 6
E LJ = 4 ϵ σ − σ
r ij r ij
 MM: Force Field
Force Field used to calculate the energy and geometry of a molecule.

• Total energy of a molecule is divided into several parts called force potentials, or potential energy
equations.
• Force potentials are calculated independently, and summed to give the total energy of the molecule.
• These equations define the potential energy surface of a molecule.
 Equation of Motion

• Newton’s second law or the equation of

motion, F=ma, where F is the force exerted on
the particle, m is its mass and a is its
acceleration.

• Knowing the force on each atom, one can

determine the acceleration of each.

• Integration of the equations of motion then

yields a trajectory that describes the positions,
velocities and accelerations of the particles as
they vary with time.

• From this trajectory, the average values of

properties can be determined.
Cars cruise down an expressway at 25 m/s. Engineers are designing an off-ramp in an
interchange with a deceleration of −2.0 m/s2 that lasts 3.0 s.
• What velocity will cars have at the end of the off-ramp?
• What minimum length should the ramp have?

v = v0 + at
v = (25 m/s) + (−2.0 m/s2)(3.0 s)
v = 19 m/s

Δs = v0t + ½at2
Δs = (25 m/s)(3.0 s) + ½(−2.0 m/s2)(3.0 s)2
Δs = 66 m
 Integration Algorithms

• The potential energy is a function of the atomic positions (3N) of all the atoms in
the system.

• Due to the complicated nature of this function, there is no analytical solution to the
equations of motion; they must be solved numerically.

• Numerous numerical algorithms have been developed for integrating the equations
of motion, such as
Verlet algorithm
Velocity Verlet
Leap-frog algorithm
Three Atoms
Progress in time
Atoms move forward in time
 The idea
Molecular Dynamics
Simulations
c what atoms do in real life, assuming a
n potential
• energy
Build realistic function
atomistic model of the system
e energy
• function
Simulate allowsofusyour
the behavior to system
calculate
over the
time force
using specific conditions
perienced by any pressure,
(temperature, atom given
volume,theetc.)
positions of the
er atoms
•
wton’s MD
lawssimulations modelthose
tell us how the motion of some
forces group ofthe
will affect particles by solving the classical
equations
tions of of motion.
the atoms
Energy (U)

Position Position 4
The Simulation Cell
Periodic Boundary Conditions
 Periodic Boundary Conditions
⚫ Simulation box replicated to form an infinite
lattice.

⚫ Periodic images of a particle in the

simulation box move with the same
orientation.

⚫ Nullify the surface effects and study the bulk

properties of the system with small number
of particles.

⚫ Linear momentum of the system under PBC

is conserved.

⚫ Minimum image convention is adopted in

order to prevent particles interacting with
 Periodic boundary conditions

Consider box of size L, repeat

box infinitely many times in all
directions

• Each particle interacts (in principle) with all particles in all boxes → problems for long-
range interactions (infinite re-summation necessary)

• short-range interactions: minimum image convention: at most the closest of all images of
a particle j can interact with a given particle i → great simplification: pick the closest
image and use this to calculate V(rij) 22
Minimum Image Convention
 forward_arrow_colour

24 of 6 © Boardworks Ltd 2012

 Molecular Dynamics

MD is a computational microscope. Understanding how things work on a

detailed molecular level ranging from electronic structures to long-time phase
behavior of molecules
 MD algorithm

Set Initial conditions (positions, velocities)

Calculate forces

Solve equations of motion

Move atoms

Increase time

Repeat as long as you need

Calculate the desired physical quantities

 Initial Positions
• X-ray, NMR, simulation or analytical calculation

– Protein Databank ([Link]

 Velocity
• The initial velocities are assigned taking them from a Maxwell
distribution at a certain temperature

[ ]
3 /2 2 2 2

( )
m m(v x +v y +v z )
P(v )= exp −
2 πk B T 2kBT
 Time Step t

• Time step t should be as large as possible to still get accurate trajectories (on the time
scale needed) and conserve of energy.

• In general, t should be ≈ 0.01 x the fastest behavior of your system

Atoms oscillate about once every 10-12 s in a solid 

MD time steps are ≈ 10-14 s in simulations of solids.
 Solvation
• Biological activity is the result of interactions between molecules and occurs at the
interfaces between molecules (protein-protein, protein-DNA, protein-solvent, DNA-
solvent, etc).

• Many biological processes occur in aqueous solution

• Solvation effects play a crucial role in determining molecular conformation, electronic

properties, binding energies, etc

• How to model solvation?

• explicit treatment: solvent molecules are added to the molecular system

• implicit treatment: solvent is modeled as a continuum dielectric
 Topology Files
Topolgy files contain:
• atom types are assigned to identify different elements and different molecular orbital environments
• charges are assigned to each atom
HN
• connectivity between atoms are established N

HB1

CB
CA HB2

HA
C HB3

O
 MASS HS 1.0080 ! thiol hydrogen
MASS C 12.0110 ! carbonyl C, peptide backbone
MASS CA 12.0110 ! aromatic C
........ (missing data here)
Example of
!-----------------------------------------------------------
RESIDUE ALA Topology File
GROUP
ATOM N TYPE=NH1 CHARGE= -.4700 END ! |
ATOM HN TYPE=H CHARGE= .3100 END ! N--HN
ATOM CA TYPE=CT1 CHARGE= .0700 END ! | HB1
ATOM HA TYPE=HB CHARGE= .0900 END ! | /
GROUP ! HA-CA--CB-HB2
HN
ATOM CB TYPE=CT3 CHARGE= -.2700 END ! | \ N
ATOM HB1 TYPE=HA CHARGE= .0900 END ! | HB3
ATOM HB2 TYPE=HA CHARGE= .0900 END ! O=C HB1
ATOM HB3 TYPE=HA CHARGE= .0900 END ! |
GROUP !
ATOM C TYPE=C CHARGE= .5100 END
ATOM O TYPE=O CHARGE= -.5100 END CB
!END GROUP CA HB2
BOND CB CA
BOND N HN
BOND N CA
BOND O C
HA
BOND C CA C HB3
BOND CA HA
BOND CB HB1
BOND CB HB2
BOND CB HB3
DONOR HN N O
ACCEPTOR O C
END {ALA }
 Parameter Files

Parameter files contain:

• force constants necessary to describe the bond energy, angle energy, torsion energy,
nonbonded interactions (van der Waals and electrostatics)
• suggested parameters for setting up the energy calculations
!BOND PARAMETERS: Force Constant, Equilibrium Radius
BOND C C 600.000 {SD=.022} 1.335 ! ALLOW ARO HEM Example of
BOND CA CA 305.000 {SD=.031} 1.375 ! ALLOW ARO

!ANGLE PARAMETERS: Force Constant, Equilibrium Angle,

Parameter File
Urie-Bradley Force Const., U.-B. equilibrium (if any)
ANGLE CA CA CA 40.00 {SD=.086} 120.0000 UB 35.000 2.416
ANGLE CP1 N C 60.00 {SD=.070} 117.0000 ! ALLOW PRO

!DIHEDRAL PARAMETERS: Energy Constant, Periodicity, Phase Shift, Multiplicity

DIHEDRAL C CT2 NH1 C 1.60 {SD=.430} 1 180.0000 ! ALLOW PEP
DIHEDRAL C N CP1 C .80 {SD=.608} 3 .0000 ! ALLOW PRO PEP

!IMPROPER PARAMETERS: Energy Constant, Periodicity(0), Phase Shift(0)

! Improper angles are introduced for PLANARITY maintaining
IMPROPER HA C C HA 20.00 {SD=.122} 0 .0000 ! ALLOW PEP POL ARO
IMPROPER HA HA C C 20.00 {SD=.122} 0 180.0000 ! ALLOW PEP POL ARO

! -----NONBONDED-LIST-OPTIONS-------------------------------
-----ELECTROSTATIC OPTIONS--------------------------------
EPS= 1.000 E14FAC= 1.000 CDIELECTRIC SHIFT
NBXMOD= 5
! ------------
! EPS SIGMA EPS(1:4) SIGMA(1:4)

NONBONDED C .1100 4.0090 .1100 4.0090 ! ALLOW PEP POL ARO

NONBONDED CA .0700 3.5501 .0700 3.5501 ! ALLOW ARO
 Simulation Procedures
Minimization
• Minimize the energy of the system in order to reach the most favorable
configuration
Equilibration
• The point of the equilibration phase is to run the simulation until the
structure, pressure, temperature and energy become stable with respect to
time.
Dynamics
• Normal/Periodic boundary condition
• Single/Multiple time stepping
• Integrators
• Electrostatics
Analyses
 MD – Shortcomings
• Quality of the force-field - systems that have been parameterized

• Size and Time – atomistic simulations can be performed only for systems of a few tenths
of angstroms on the length scale and for a few nanoseconds on the time scale

• Conformational freedom of the molecule – the number of possible conformations a

molecule can adopt is enormous, growing exponentially with the number or rotatable
bonds.

• Connectivity of atoms cannot change during dynamics – no chemical reactions.

 MD

What you have up to now:

• pdb file,
• topology file
• parameter file

What you need next: a program capable of reading and manipulating this information

Programs:

AMBER
CHARMm
VASP (DFT)
XMD
CPMD(DFT) (Car-Parrinello MD)
NAMD
 Small Molecule Parameters and Topology

Web-servers:
• [Link]

• [Link]

• [Link]

• [Link]

• [Link]
 Web-servers for running MD Simulations

• [Link]
• [Link]
• [Link]
• [Link]
 forward_arrow_colour

40 of 6 © Boardworks Ltd 2012