Transcript Essential Bioinformatics and Biocomputing (LSM2104

SMA5233
Particle Methods and Molecular Dynamics
Lecture 2: Physical Principles and Design Issues of MD
A/P Chen Yu Zong
Tel: 6874-6877
Email: [email protected]
http://bidd.nus.edu.sg
Room 08-14, level 8, S16
National University of Singapore
When is classical mechanics a reasonable approximation?
In Newtonian physics, any particle may possess any one of a
continuum of energy values. In quantum physics, the energy is
quantized, not continuous. That is, the system can accomodate
only certain discrete levels of energy, separated by gaps.
At very low temperatures these gaps are much larger than
thermal energy, and the system is confined to one or just a few
of the low-energy states. Here, we expect the `discreteness' of
the quantum energy landscape to be evident in the system's
behavior.
As the temperature is increased, more and more states become
thermally accessible, the `discreteness' becomes less and less
important, and the system approaches classical behavior.
2
When is classical mechanics a reasonable approximation?
For a harmonic oscillator, the quantized energies are
separated by DE=hf, where h is Planck's constant and f
is the frequency of harmonic vibration.
Classical behavior is approached at temperatures for
which KBT>>hf, where KB is the Boltzmann constant and
KBT = 0.596 kcal/mol at 300 K. Setting hf = 0.596
kcal/mol yields f = 6.25/ps, or 209cm-1. So a classical
treatment will suffice for motions with characteristic times
of a ps or longer at room temperature.
3
Classical Mechanics - in a Nutshell
Building on the work of Galileo and others, Newton unveiled his
laws of motion in 1686. According to Newton:
1.
A body remains at rest or in uniform motion (constant velocity both speed and direction) unless acted on by a net external
force.
2.
In response to a net external force, F, a body of mass m
accelerates with acceleration a=F/m.
3.
If body i pushes on body j with a force Fij, then body j pushes
on body i with a force Fji=-Fij.
4
Classical Mechanics - in a Nutshell
For energy-conserving forces, the net force Fi on particle i is the
negative gradient (slope in three dimensions) of the potential energy
with respect to particle i's position: Fi= - ▼iV(R) where V(R)
represents the potential energy of the system as a function of the
positions of all N particles, R.
In three dimensions, ri is the vector of length 3 specifying the position
of the atom, and R is the vector of length specifying all coordinates.
In the context of simulation, the forces are calculated for energy
minimizations and molecular dynamics simulations but are not needed
in Monte Carlo simulations.
Classical mechanics is completely deterministic: Given the exact
positions and velocities of all particles at a given time, along with the
function , one can calculate the future (and past) positions and
velocities of all particles at any other time. The evolution of the
system's positions and momenta through time is often referred to as a
trajectory.
5
Statistical Mechanics - Calculating Equilibrium Averages
According to statistical mechanics, the probability that a given state
with energy E is occupied in equilibrium at constant particle number N,
volume V, and temperature T (constant , the `canonical' ensemble) is
proportional to exp{-E/ KBT}, the `Boltzmann factor.'
Probability ~ exp{-E/ KBT}
The equilibrium value of any observable O is therefore obtained by
averaging over all states accessible to the system, weighting each
state by this factor.
6
Statistical Mechanics - Calculating Equilibrium Averages
Classically, a microstate is specified by the positions and velocities
(momenta) of all particles, each of which can take on any value. The
averaging over states in the classical limit is done by integrating over
these continuous variables:
where the integrals are over all phase space (positions and
momenta _) for the N particles in 3 dimensions.
When all forces (the potential energy V) and the observable O are
velocity-independent, the momentum integrals can be factored and
canceled:
where
is the total kinetic energy,
and E=K+V. As a result, Monte Carlo simulations
compare V's, not E's.
7
What is molecular mechanics?
Molecular Mechanics (MM) finds the geometry that corresponds
to a minimum energy for the system - a process known as
energy minimization.
A molecular system will generally exhibit numerous minima,
each corresponding to a feasible conformation. Each minimum
will have a characteristic energy, which can be computed. The
lowest energy, or global minimum, will correspond to the most
likely conformation.
8
Energy Minimization
U
Torsion angles
Are 4-body
Non-bonded
pair
1
2


K
b

b

b
0
2
all bonds
1
2
 
K   0 
all angles 2

 K 1  cosn 
all torsions
Angles
Are 3-body
Bonds
Are 2-body
6
 R 12
 Rij  
ij
   ij    2  
r  
 rij 
i , j nonbonded
 ij  



i , j nonbonded
qi q j
4 0rij
U is a function of the conformation C of the protein. The problem of
“minimizing U” can be stated as finding C such that U(C) is minimum.
9
Units in force fields
Most force fields use the AKMA (Angstrom – Kcal – Mol – Atomic Mass Unit) unit
system:
Quantity
AKMA unit
Equivalent SI
Energy
1 Kcal/Mol
4184 Joules
Length
1 Angstrom
10-10 meter
Mass
1 amu
(H=1amu)
1.6605655 10-27 Kg
Charge
1e
1.6021892 10-19 C
Time
1 unit
4.88882 10-14 second
Frequency
1 cm-1
18.836 1010 rd/s
10
Energy Minimization
Local minima:
A conformation X is a local minimum if there exists a
domain D in the neighborhood of X such that for all Y≠X in D:
U(X) <U(Y)
Global minimum:
A conformation X is a global minimum if
U(X) <U(Y)
for all conformations Y ≠X
11
Molecular Modeling:
Energy minimization:
p2
H
atoms
Force guided approach:
vn
[1  cos(n   )]

2

bond  rotation
 [V (1  e
Initialize:
Change atom position: xi -> xi+dxi
1
1
k r (r  req ) 2 
k (   eq ) 2 

2
bond  stretch
bond  angle bending 2
 2m  
H bond
0
 [V (1  e
S bond
 a ( r  r0' ) 2
)  V0 ] 

non bonded
 a ( r  r0' ) 2
)  V0 ] 
0
[
Aij
rij12

Bij
rij6

qi q j
 ij rij
]
Compute potential energy change:
V -> V +dV
Determine next movement:
Force: Fxi=-dV/dxi; Fyi=-dV/dyi; Fzi=-dV/dzi
Atom displacement: dxi=C*Fxi
New position: xi=xi+dxi
12
Notes on computing the energy
Bonded interactions are local, and therefore their computation has a linear
computational complexity (O(N), where N is the number of atoms in the
molecule considered.
The direct computation of the non bonded interactions involve all pairs of
atoms and has a quadratic complexity (O(N2)).
This is usually prohibitive for large molecules.
U NB
6
 R 12



R
qi q j
ij
ij
   ij    2    




r
r

i , j nonbonded
i , j nonbonded 4 0rij
ij 
ij  




Reducing the computing time: use of cutoff in UNB
13
Cutoff schemes for faster energy computation
U NB
6
 R 12



R
qi q j
ij
ij
  ij S rij  ij    2     ij S rij 
 r   i, j
4 0rij
 rij 
i, j
ij 



ij : weights (0< ij <1). Can be used to exclude bonded terms,
or to scale some interactions (usually 1-4)
S(r) : cutoff function.
Three types:
1) Truncation:
1 r  b
S (r)  
0 r  b
b
14
Cutoff schemes for faster energy computation
2. Switching
ra
 1

S ( r )  1  y ( r ) 2 2 y ( r )  3 a  r  b
 0
rb

a
r2  a2
y(r)  2
b  a2
with
b
3. Shifting
  r 2 
S1 ( r )  1    
  b  
2
rb
or
 r
S2 (r)  1  
 b
2
rb
b
15
The minimizers
Some definitions:
Gradient:
The gradient of a smooth function f with continuous first and second
derivatives is defined as:
 f
f
f  X   
...
xi
 x1
f 
xN 
...
Hessian
The n x n symmetric matrix of second derivatives, H(x), is called
the Hessian:
 2 f
2 f
2 f 

2
 x1
 ...
 2 f
H ( x)  
 xi x1
 ...
 2 f

 xN x1
...
...
...
...
...
x1x j
...
2 f
xi x j
...
2 f
xN x j
...
...
...
...
...
x1xN
...
2 f
xi xN
...
2 f
xN2










16
The minimizers
Minimization of a multi-variable function is usually an iterative process, in which
updates of the state variable x are computed using the gradient, and in some
(favorable) cases the Hessian.
Iterations are stopped either when the maximum number of steps (user’s input)
is reached, or when the gradient norm is below a given threshold.
Steepest descent (SD):
The simplest iteration scheme consists of following the “steepest
descent” direction:
( sets the minimum
along the line defined
k 1
k
k
by the gradient)
x
 x  f ( x )
Usually, SD methods leads to improvement quickly, but then exhibit
slow progress toward a solution.
They are commonly recommended for initial minimization iterations,
when the starting function and gradient-norm values are very large.
17
The minimizers
Conjugate gradients (CG):
In each step of conjugate gradient methods, a search vector pk is
defined by a recursive formula:
pk 1  f xk   bk 1 pk
The corresponding new position is found by line minimization along
pk:
xk 1  xk  k pk
the CG methods differ in their definition of b:
- Fletcher-Reeves:
- Polak-Ribiere
- Hestenes-Stiefel
f ( xk 1 )  f ( xk 1 )
b 
f ( xk )  f ( xk )
f ( xk 1 )  f ( xk 1 )  f ( xk )
b kPR

1
f ( xk )  f ( xk )
FR
k 1
b kHS1 
f ( xk 1 )  f ( xk 1 )  f ( xk )
pk  f ( xk 1 )  f ( xk )
18
Newton’s methods:
The minimizers
Newton’s method is a popular iterative method for finding the 0 of a
one-dimensional function:
g xk 
xk 1  xk 
g ' xk 
x3 x2
x1
x0
It can be adapted to the minimization of a one –dimensional function, in which
case the iteration formula is:
g ' x 
xk 1  xk 
k
g ' ' xk 
The equivalent iterative scheme for multivariate functions is based on:
xk 1  xk  H
1
k
xk f xk 
Several implementations of Newton’s method exist, that avoid computing
the full Hessian matrix: quasi-Newton, truncated Newton, “adopted-basis
Newton-Raphson” (ABNR),…
19
What is molecular dynamics simulation?
Molecular Dynamics (MD) applies the laws of
classical mechanics to compute the motion of the
particles in a molecular system.
Simulation that shows how the atoms in the system
move with time
Typically on the nanosecond timescale
Atoms are treated like hard balls, and their motions
are described by Newton’s laws.
20
What is molecular dynamics simulation?
Molecular dynamics simulations generate information at the microscopic
level, including atomic positions and velocities.
The conversion of this microscopic information to macroscopic observables
such as pressure, energy, heat capacities, etc., requires statistical
mechanics.
MD Simulation
Microscopic information
Statistical Mechanics
Macroscopic properties
21
Characteristic protein motions
Type of motion
Timescale
Amplitude
Local:
bond stretching
angle bending
methyl rotation
0.01 ps
0.1 ps
1 ps
<1Å
Periodic (harmonic)
Medium scale
loop motions
SSE formation
ns – ms
Global
protein tumbling
(water tumbling)
protein folding
20 ns
(20 ps)
ms – hrs
1-5 Å
>5Å
Random (stochastic)
22
Molecular Dynamics Simulation
Molecule: (classical) N-particle system
Newtonian equations of motion:
d2  
mi 2 ri  Fi ( r )
dt

Fi (r )  iV (r )
with


r  ( r1 ,..., rN )
Integrate numerically via the „leapfrog“ scheme:
with
Δt  1fs!
(equivalent to the Verlet algorithm)
23
How do you run a MD simulation?
Get the initial configuration
From x-ray crystallography or NMR spectroscopy (PDB)
Assign initial velocities
At thermal equilibrium, the expected value of the kinetic energy of the
system at temperature T is:
Ekin
1 3N
1
2
  mi vi  (3N )kBT
2 i 1
2
This can be obtained by assigning the velocity components vi from a
random Gaussian distribution
with mean 0 and standard deviation (kBT/mi):
2
i
v
kBT

mi
24
How do you run a MD simulation?
For each time step:
Compute the force on each atom:
E
F ( X )   E ( X )  
X
X: cartesian vector
of the system
Solve Newton’s 2nd law of motion for each atom, to get new
coordinates and velocities
M diagonal mass matrix

M X  F(X )
.. means second order
differentiation with
respect to time
Store coordinates
Stop
Newton’s equation cannot be solved analytically:
Use stepwise numerical integration
25
What the integration algorithm does
Advance the system by a small time step Dt during which
forces are considered constant
Recalculate forces and velocities
Repeat the process
If Dt is small enough, solution is a reasonable
approximation
26
Choosing a time step Dt
Not too short so that conformations are efficiently sampled
Not too long to prevent wild fluctuations or system ‘blowup’
An order of magnitude less than the fastest motion is ideal
– Usually bond stretching is the fastest motion:
C-H is ~10fs so use time step of 1fs
– Not interested in these motions?
Constrain these bonds and double the time step
27
Molecular Dynamics Ensembles
The method discussed above is appropriate for the
micro-canonical ensemble: constant N (# of particles) V
(volume) and ET (total energy = E + Ekin)
In some cases, it might be more appropriate to simulate
under constant Temperature (T) or constant Pressure
(P):
– Canonical ensemble: NVT
– Isothermal-isobaric: NPT
– Constant pressure and enthalpy: NPH
How do you simulate at constant temperature and pressure?
28
Common Ensembles
Name
Probability distribution
Microcanonical
(EVN)
 i  1
Canonical
(TVN)
 ( Ei )  Q1 e b Ei
Schematic
Isothermal-isobaric  ( Ei ,Vi )  D1 e b ( Ei  PVi )
(TPN)
Grand-canonical
(TVm)
 ( Ei , Ni )  1 e b ( Ei  m Ni )
Note: b  1/ kT
29
MD in different thermodynamic ensemble
MD naturally samples the NVE (microcanonical
ensemble): this is not very useful for most
applications.
• More useful ensembles include:
– Canonical (NVT): requires the particles to interact
with a thermostat
– Isobaric-Isothermal (NPT): requires the particles to
interact with a thermostat and a baristat
30
Temperature
Thermodynamics definition
1  S 
 
T  E V , N
1


ln ( E ,V , N )
kT E
Measuring T
1 N pi2
2
kT 

K

Nd i 1 m Nd
31
Canonical (NVT)
Several possible approaches
– Berendsen Method (velocities are rescaled
deterministically after each step so that the system
is forced towards the desired temperature)
– Stochastic collisions (Anderson: a random
particle’s velocity is reassigned by random
selection from the Maxwell-Boltzmann distribution
at set intervals)
– Extended Lagrangians (Nose and Nose-Hoover:
the thermal reservoir is considered an integral part
of the system and it is represented by an
additional degree of freedom: set of deterministic
equations that give sampling in NVT)
32
Simulating at constant T:
the Berendsen scheme
system
Heat bath
Bath supplies or removes heat from the system as
appropriate
Exponentially scale the velocities at each time step
by the factor :
Dt  Tbath 
  1  1 

 
T 
T: “kinetic” temperature
where  determines how strong the bath influences
the system
Berendsen et al. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684 (1984)
33
Simulating at constant P:
Berendsen scheme
system
pressure bath
Couple the system to a pressure bath
Exponentially scale the volume of the simulation box
at each time step by a factor :
  1
Dt
P
P  Pbath 
where
N
2

P   Ekin   xi  Fi 
3v 
i 1

where  : isothermal compressibility
P : coupling constant
u : volume
xi: position of particle i
Fi : force on particle i
Berendsen et al. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684 (1984)
34
Nose-Hoover Thermostat
Introduce a new coordinate, s, representing reservoir in the Lagrangian
N
LNose
2
mi ( sri )
Q 2
N

 U (r )  s  gkT ln s
2
2
i 1
With momenta:
L
2
pi 
 mi s ri
ri
L
ps 
 Qs
s
Q=effective mass associated
with s
g is a parameter we will fix
later
35
Nose-Hoover Thermostat
Extended-system Hamiltonian:
N
H Nose  
2
pi
 U (r ) 
N
2
2
ps
 gkT ln s
2Q
i 1 2mi s
it can be shown that the molecular positions and
momenta follow a canonical (NVT) distribution if g =
3N+1

 p( t )

p 
A  lim dtA
, r ( t )   A , r 
s 
 s( t )

0
 
 A p' , r 
N VT
Nose
With p’=p/s
36
Nose-Hoover Thermostat r '  r
Real and virtual variables are related by:
p
p' 
s
s'  s
Dt
Dt ' 
s
s can be interpreted as a scaling factor of the time step
– treal = tsim/s
– since s varies during the simulation, each “true” time step
is of varying length
37
Nose-Hoover Thermostat
From the Hamiltonian, can derive the equations of motion
for the virtual variables p, r and t.
Scaled-variables
equations of motion
– constant simulation Dt
– fluctuating real Dt
Real-variables equation
of motion (here, I have
Dropped the prime)
H
p
ri 
 i2
pi mi s
ri 
H
pi  
 Fi
ri
sps
pi  Fi 
pi
Q
s sps

s Q
H ps
s

ps Q

H 1  N pi
ps  
 
 gkT 
2
s s  i 1 mi s

pi
mi

 ( sps / Q ) 1  N pi
    gkT 
38
t
Q  i 1 mi

IMPLEMENTATION
Simplify expression (Nose-Hoover) by introducing a
thermodynamic
Friction coefficient:
s' ps'

Q
So equations of motions (dropping the prime):
ri 
pi
mi
pi  Fi   pi
s

s

1  N pi
     gkT 
Q  i 1 mi

(g=3N)
39
An experiment is usually made on a macroscopic sample that
contains an extremely large number of atoms or molecules
sampling an enormous number of conformations. In statistical
mechanics, averages corresponding to experimental
observables are defined in terms of ensemble averages;
40
AVERAGES
Thermodynamic (Ensemble) average: average over all points in
phase space at a Single time
Dynamic average: average over a single point in phase space at
All times
Ergodic hypothesis: for infinitely long trajectory, thermo average=
Dynamic average
A time  A ensemble
A large number of observations made on a single system at N
arbitrary instants in time have the same statistical properties
as observing N arbitrarily chosen systems at the same time.
41
To get a thermo average, need to know the probability of
finding the system at each point (state) in phase space.

 r ,p
N
N

 H (r N , p N ) 
1
 exp 

Q
k
T
B


Probability density
 H (r N , p N ) 
Q   dp dr exp 

k
T
B


N
Partition function
N
Ensemble (thermodynamic) average of observable A(rN,pN):
A
ens


  dp N dr N A r N , p N  r N , p N

42
In an MD simulation, we determine a time average of A,
which is expressed as
where  is the simulation time, M is the number of time
steps in the simulation and A(pN,rN) is the instantaneous
value of A.
43
THERMODYNAMIC PROPERTIES FROM MD SIMULATIONS
Thermodynamic (equilibrium) averages can be calculated
via time averaging:
1
A 
N
N
 A(t )
i 1
i
Examples:
1 N
Average potential energy
U   U (t i )
N i 1
N M m
1
j
K  
v j t i v j t i  Average kinetic energy
N i 1 j 1 2
44
Other quantities obtained from MD:
Heat Capacity
Expressed in terms of fluctuations in the energy
2


E

( b A) 
 
2
Cv   
 k b 

2
 T V , N
 b 
V ,N
 1
N
N
b E
 k b
dr
dp
Ee
b Q( b ) 
2
Cv  k b
2

E
2
 E
2

45
Other Quantities
Root Mean Square Deviation (for instance,
over C atoms)
RMSD 
N
1
N
B
A 2
|
r

r
 i i |
i 1
Radius of Gyration
Rg 
1
N
N
2
|
r

r
|
 i cm
i 1
46
From PDB orXtal
Relieve local stresses,
Relax bond length/angle distortions
velocities (vx, vy,vz) are randomly
assigned according to a Maxwellian
distribution,
T(t) = (1/kB(3N-n)) Si (mivi2)
Scale velocities by (T'/T(t))1/2 to get new
temperatureT
P(vi)dv = (m/2pikBT)1/2 exp[ -mv2/2kBT ]dv
to overcome “hot spots” due to
Initial high velocities
47
Types of MD Systems
Empirical Potentials
Empirical potentials either ignore quantum mechanical effects, or
attempt to capture quantum effects in a limited way through entirely
empirical equations. Parameters in the potential are fitted against
known physical properties of the atoms being simulated, such as elastic
constants and lattice parameters.
Empirical potentials can further be subcategorized into pair potentials,
and many-body potentials
– Pair potentials: the total potential energy of a system can be calculated
from the sum of energy contributions from pairs of atoms.
– Many-body potentials: the potential energy cannot be found by a sum over
pairs of atoms.
Tersoff Potential involves a sum over groups of three atoms, with the angles
between the atoms being an important factor in the potential.
Tight-Binding Second Moment Approximation (TBSMA), the electron density of
states in the region of an atom is calculated from a sum of contributions from
surrounding atoms, and the potential energy contribution is then a function of
this sum.
48
Types of MD Systems
Semi-Empirical Potentials
Semi-empirical potentials make use of the matrix representation from
quantum mechanics. However, the values of the matrix elements are
found through empirical formulae that estimate the degree of overlap of
specific atomic orbitals. The matrix is then diagonalized to determine
the occupancy of the different atomic orbitals, and empirical formulae
are used once again to determine the energy contributions of the
orbitals.
There are a wide variety of semi-empirical potentials, known as tightbinding potentials, which vary according to the atoms being modeled.
49
Types of MD Systems
Ab-initio (First Principles) Methods
Ab-initio (first principles) methods use full quantum-mechanical formula
to calculate the potential energy of a system of atoms or molecules.
This calculation is usually made "locally", i.e., for nuclei in the close
neighborhood of the reaction coordinate. Although various
approximations may be used, these are based on theoretical
considerations, not on empirical fitting. Ab-Initio produce a large
amount of information that is not available from the empirical methods,
such as density of states information.
Perhaps the most famous ab-initio package is the Car-Parrinello
Molecular Dynamics (CPMD) package based on the density functional
theory..
50
Types of MD Systems
Coarse-Graining: "Reduced Representations"
The earliest MD simulations on proteins, lipids and nucleic acids
generally used "united atoms", which represents the lowest level of
coarse graining. For example, instead of treating all four atoms of a
methyl group explicitly (or all three atoms of a methylene group), one
represents the whole group with a single pseudoatom. This
pseudoatom must, of course, be properly parameterized so that its van
der Waals interactions with other groups have the proper distancedependence. Similar considerations apply to the bonds, angles, and
torsions in which the pseudoatom participates. In this kind of united
atom representation, one typically eliminates all explicit hydrogen
atoms except those that have the capability to participate in hydrogen
bonds ("polar hydrogens").
51
Types of MD Systems
Coarse-Graining: "Reduced Representations"
When coarse-graining is done at higher levels, the accuracy of the
dynamic description may be less reliable. But very coarse-grained
models have been used successfully to examine a wide range of
questions in structural biology. Following are just a few examples.
– protein folding studies are often carried out using a single (or a few)
pseudoatoms per amino acid;
– DNA supercoiling has been investigated using 1-3 pseudoatoms per
basepair, and at even lower resolution;
– Packaging of double-helical DNA into bacteriophage has been investigated
with models where one pseudoatom represents one turn (about 10
basepairs) of the double helix;
– RNA structure in the ribosome and other large systems has been modeled
with one pseudoatom per nucleotide.
– The parameterization of these very coarse-grained models must be done
empirically, by matching the behavior of the model to appropriate
experimental data.
52
Electrostatics and the `Generalized Born' Solvent Model
All problems in electrostatics boil down to the solution of a single
equation, Poisson's equation:
where is the Laplacian operator,  is the electrostatic potential,
 is the charge density (total charge per unit volume including all
`free' and `polarization' charges), and 0 is the permittivity of
free space. In Cartesian (xyz) coordinates,
The electrostatic potential at a given point in space is the
potential energy per unit charge that a test charge would have if
positioned at that point in the electric field E specified by  :
53
Solvent Model: Structure of water molecule
54
Solvent Model: Estimate the average
distance between water molecules
3
Vmole
Mass / Mole 18*10


3
dengsity
10
V1 _ molecule
Vmole
 29
3


3
.
0
*
10
m
6,022.1023
29
L  3.0 *10
L  3.1A
3
 1.8 *105 m 3
m
3
55
Solvent Model: Layout of water molecule
Divide the cubic evenly
into meshes, according to
the distances between
water molecules
Place oxygen atom of
each water molecule on
the node of the mesh,
generate random
direction for hydrogen
atoms.
56
Add-in a Molecule (Nano-tube)
57
MD Simulation Result
58
MD-Experiments with Argon Gas
59
Role of environment - solvent
explicit
or
implicit?
box
or
droplet?
60
Surface (tension) effects?
periodic boundary conditions
and the minimum image convention
61
Proteins jump between many, hierarchically ordered
„conformational substates“
62
H. Frauenfelder et al., Science 229 (1985) 337
Molecular Dynamics Simulations
• External coupling:
Temperature (potential truncation, integration errors)
Pressure (density equilibration)
System translation/rotation
• Analysis
Energies (individual terms, pressure, temperature)
Coordinates (numerical analysis, visual inspection!)
 Mechanisms
63
Molecular Dynamics Simulations
Design Constraints:
•
Design of a molecular dynamics simulation can often encounter limits of
computational power. Simulation sizes and time duration must be selected
so that the calculation can finish within a useful time period.
•
The simulation's time duration is dependent on the time length of each
timestep, between which forces are recalculated. The timestep must be
chosen small enough to avoid discretization errors, and the number of
timesteps, and thus simulation time, must be chosen large enough to
capture the effect being modeled without taking an extraordinary period of
time.
•
The simulation size must be large enough to express the effect without the
boundary conditions disrupting the behavior. Boundary constraints are often
treated by choosing fixed conditions at the boundary, or by choosing periodic
boundary conditions in which one side of the simulation loops back to the
opposite side.
64
Molecular Dynamics Simulations
Design Constraints:
•
The scalability of the
simulation with respect to the
number of molecules is
usually a significant factor in
the range of simulation sizes
which can be simulated in a
reasonable time period.
•
In Big O notation, common
molecular dynamics
simulations usually scale by
either O(nlog(n)), or with good
use of neighbour tables, O(n),
with n as the number of
molecules
65