Molecular Dynamics

Molecular Dynamics Simulation
A Short Introduction
Michel Cuendet
Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008
Plan
Introduction
The classical force field
Setting up a simulation
Connection to
statistical mechanics
Usage of MD simulation
Why we do simulation
In some cases, experiment is :
1.
impossible
Inside of stars
Weather forecast
2.
too dangerous
Flight simulation
Explosion simulation
3.
expensive
High pressure simulation

Windchannel simulation
4.
blind
Some properties cannot be

observed on very short time-scales
and very small space-scales
Simulation is a useful complement, because it can :
replace experiment
provoke experiment
explain experiment
aid in establishing intellectual property
Molecular modeling
What is Molecular Modeling?
Molecular Modeling is concerned with the description of the atomic
and molecular interactions that govern microscopic and macroscopic
behaviors of physical systems
What is it good for?

The essence of molecular modeling resides in the connection between
the macroscopic world and the microscopic world provided by the
theory of statistical mechanics
Macroscopic
observable
(Solvation energy,
affinity between two
proteins, H-H distance,
conformation, ... )
Average of observable
over selected microscopic
states
Computational tools
Quantum Mechanics (QM)
Electronic structure (Schrdinger)
ACCURATE
EXPENSIVE small system
10-100 atoms
10-100 ps
Classical Molecular Mechanics (MM)

Empirical forces (Newton)
LESS ACCURATE
FAST
104-105 atoms
10-100 ns
Mixed Quantum/Classical (QM/MM)
MM
QM
104-105 atoms
10-100 ps
5
Types of phenomena
Goal : simulate/predict processes such as
1.
2.
3.
4.
polypeptide folding
biomolecular association
partitioning between solvents
membrane/micelle formation
thermodynamic equilibria
governed by weak
(non-bonded) forces
5.
6.
7.
chemical reactions, enzyme catalysis

enzyme catalysis
photochemical ractions, electron transfer
chemical transformations
governed by strong forces
characteristics (1-4):
degrees of freedom:
equations of motion:
governing theory:
characteristics (5-7):
degrees of freedom:
equations of motion:
governing theory:
atomic (solute + solvent)

classical dynamics
electronic, nuclear
quantum dynamics
quantum statistical mechanics
classical MD
quantum MD
Processes: Thermodynamic Equilibria

Folding
folded/native
Micelle Formation
denatured
Complexation
bound
unbound
micelle
mixture
Partitioning
in membrane
in water in mixtures
Plan
Introduction
Connection to
Definition of a model for molecular simulation

Every molecule consists of atoms that are very strongly bound to each other
Degrees of freedom:
atoms as elementary
particles + topology
Forces or
interactions
between atoms
Boundary conditions
MOLECULAR
MODEL
Force field =
physico-chemical
knowledge
Methods to generate
configurations of
atoms: Newton
system
temperature
pressure
walls
external forces
Classical force fields

Goals of classical (semi-empirical) force fields
- Definition of empirical potential energy functions V(r) to model the
molecular interactions
- These functions need to be differentiable in order to compute the
forces acting on each atom: F=-V(r)
Implementation of calssical potential energy functions

1. Theoretical functional forms are derived for the potential energy V(r).
2. Definition of atom types that differ by their atomic number and
chemical environment, e.g. the carbons in C=O or C-C are of
different types.
3. Parameters are determined so as to reproduce the interactions between
the various atom types by fitting procedures
- experimental enthalpies (CHARMM)
- experimental free energies (GROMOS, AMBER)
Parametrization available for proteins, lipids, sugars, ADN, ...
10
Covalent bonds and angles

Type: CT1-CT1
Bonds
r
E bond = K b (r " r0 )
r0
Type: CT1-CT1-CT3
Angles
E angle = K" (" # " 0 )

11
Dihedrals and Improper torsions

Type: X-CT1-CT1-X
Dihedral angles
E dihedral = K" [1+ cos(n" # $ )]

Improper angles
!
Type: OC-OC-CT1-CC
O
C
H2N
E improper = K" (# $ # 0 )
R
2
12
Van der Waals interactions

Lennard -Jones potential :
: collision parameter
: well depth
16
rm : distance at min rm = 2 "
+$ # '12 $ # ' 6 .
+$ r '12 $ r ' 6 .
E VdW = 4" -& ) * & ) 0 = " -& m ) * 2& m ) 0
%r(/
,% r ( % r ( /
,% r (
Combination rule for two different atoms i, j :
rm = rm,i + rm, j
Type: CT3-CT3
" = "i" j
!
Repulsive :
Pauli
! exclusion principle
1
" 12
r
EVdW
Attractive:
induced dipole / induced dipole
!
rm
1
"# 6
r
13
Electrostatic interactions
Coulomb law
qi q j
E elec =
4 "#0# rij
where is the dielectric constant :

1
for vacuum,
4-20
for protein core,
80
for water
The Coulomb energy decreases only as 1/r

Despite dielectric shielding effects, it is a long range interaction
Special techniques to deal with this :
- PME : for stystems with periodic boundary conditions
- Reaction Field : suppose homogeneous dielectric outside cutoff
14
Derived Interactions
Some interactions are often referred to as particular interactions, but they
result from the two interactions previously described, i.e. the electrostatic
and the van der Waals interactions.
1) Hydrogen bonds (Hb)

- Interaction of the type D-H A
- The origin of this interaction is a
dipole-dipole attraction
- Typical ranges for distance and angle:
2.4 < d < 4.5 and 180 < < 90
H
d
2) Hydrophobic effect
- Collective effect resulting from the energetically
unfavorable surface of contact between the water
and an apolar medium (loss of water-water Hb)
- The apolar medium reorganizes to minimize the
water exposed surface (compaction, association... )
Water
Oil
15
The total potential energy function

E=
# K b (r " r0 ) +
bonds
# K$ ($ " $0 ) +
angles
# K% [1+ cos(n% " &)] +

dihedrals
# K' (' " '0 )
impropers
/) r ,12 ) r , 6 2
qiq j
m
m
+ #(1+ . " 2+ . 4 + #
* * r - 3 i> j 4 5(0( r
i> j 0 r
For a system with 1500 atoms

~ 106 pairs of interacting atoms
Introduction of cutoff
16
Cutoff for non-bonded interactions

For an atom A, only non-bonded interactions
with atoms within are calculated
Non-bonded neighbour lists
Generally, = 8 to 14
Three cutoff schemes: strict, shift, switch

Shift and switch:
cutofnb
E'(r) = E(r) " S(r)

S(r) differentiable
cutonnb cutnb
!
17
Effect of cutoff
No cutoff
8 cutoff
Elec
VdW
Total
Elec
VdW
Total
18
Force field parametrization

E=
# K (r " r ) + #K ($ " $ )
b
bonds
angles
#K [1+ cos(n% " &)] + # K (' " ' )

%
'
impropers
dihedrals
0* ) -12 * ) -6 3
qiq j
+ # 4(2, / " , / 5 + #
1+ r . + r . 4 i> j 46(0( r
i> j
Phase
Type of properties
Force field
parameters
structural data
(exp.)
Type of
system
small
molecules
crystalline solid
phase
r0, 0, 0
spectroscopic data
(exp.)
small
molecules
gas phase
molecular geometry:
bond lengths, angles
intra-molecular
vibrations: force
constants
small
molecules
gas phase
torsional-angle
rotational profiles
K, , n
small
molecules
gas phase
atom charges
charges qi (initial)
thermodynamic
data
(exp.)
molecules in
solution,
mixtures
condensed phase
heat of vaporisation,
density, free energy of
solvation
v. d. Waals : i, i
charges qi (final)
dielectric data
(exp.)
small
molecules
condensed phase
dielectric permittivity,
relaxation
charges
qi
transport data
(exp.)
small
molecules
condensed phase
transport coefficients:
diffusion, viscosity
v. d. Waals : i, i
charges qi
Type of data
quantum-chemical
calculations :
energy profiles
(theor.)
quantum-chemical
calculations :
electron densities
(theor.)
K b , K , K
19
Solvation
Fundamental influence on the structure, dynamics and
thermodynamics of biological molecules
Effect through:
Solvation of charge
Screening of charge - charge interactions
qi q j
E elec =
4 "#0#(r) r
=1
-
=80
Hydrogen bonds between water molecules and

polar functions of the solute
Taken into account via:

Explicit solvation. Water molecules are included.
Stochastic boudary conditions
Periodic boundary conditions
Implicit solvation. Water effect is modeled.
Screening constant
Implicit solvation models (Poisson Boltzmann, Generalized Born)
20
Explicit solvation schemes

Stochastic boundary conditions
The region of interest is solvated in a water
sphere at 1atm. The water molecules are
submitted to an additional force field that
restrain them in the sphere while maintaining
a strong semblance to bulk water.
Periodic boundary conditions

The fully solvated central cell is simulated, in
the environment produced by the repetition of
this cell in all directions.
21
Implicit solvation
Screening constant
" = Nr
N=4,8. The dielectric constant is a function of atom distance.

Mimic screening effect of solvent. Simple, unphysical but efficient.
In CHARMM:
NBOND RDIE EPS 4.0
Implicit solvent models

Poisson Boltzmann (PB) equation.
" {#(r)"$ (r)} % & 'sinh[$ (r)] = %4 '( (r)
f(r) : electrostatic potential,

r(r) : charge density
Equation solved numerically. Very time consuming.

In CHARMM: PBEQ module.
!
Generalized born (GB) equation.
=80
Gelec = $
i> j
=1
qiq j
qi q j
1 & 1)
% (1% +$ $
4 "#0 r 2 ' # * i j r 2 + ai a j exp(% r 2 4ai a j )
ai : Born radius
solvation energy
Others: EEF1, SASA, etc...

!
Explicit hydrogen bonds with water molecules are lost!

22
Introduction to molecular surfaces

Hydrophobic effect :
Simply modelled by a non-polar solvation (free) energy term,
proportional to the solvent accessible surface area (SASA) :
"Gnp, solv = #SASA + b
= 0.00542 kcal mol 1 2

b = 0.92 kcal mol-1
Rotating probe: radius 1.4

Connolly Surface
Solvent accessible
Surface (SASA)
Van der Waals
Surface
Definitions:
- Van der Waals:
- Connolly:
- Solvent:
ensemble of van der Waals sphere centered at each atom

ensemble of contact points between probe and vdW spheres
ensemble of probe sphere centers
23
Examples of molecular surfaces

Van der Waals
Connolly
(Contact)
Solvent accessible
24
Limitations of classical MD
Problems
Solutions
1) Fixed set of atom types

2) No electronic polarization:
- fixed partial charges allow
for conformational polarization
but not electronic polarization
3) Quadratic form of potentials:

- problematic far from equilibrium
values
- no bond creation or deletion
Fluctuating charges treated as

dynamical parameters
Charges on springs representing
e- clouds
QM-MM
Full QM simulations
QM-MM
Full QM simulations
25
Coarse grain models

All-atom model
16 (CH2 or CH3) atoms
Map
to all-atom
configurations
Centre of mass
A1 A4
Coarse-grained model
4 atoms
Centre of mass
B1 B4
Centre of mass
C1 C4
Centre of mass
D1 D4
Centre of mass
W1 W4
B
C
D
26
Molecular dynamics software

Package name
CHARMM
www.charmm.org
Amber
supported force fields

CHARMM (E / I; AA / UA), Amber
Amber (E / I ; AA)
amber.scripps.edu
GROMOS
www.igc,ethz.ch/GROMOS
Gromacs
www.gromacs.org
NAMD
www.ks.uiuc.edu/Research/namd
E = explicit solvent
I = implicit solvent
Gromos (E / vacuum ; UA)
Amber, Gromos, OPLS - (all E)
CHARMM, Amber, Gromos, ...

AA = all atom
UA = united atom (apolar H omitted)
27
Plan
Introduction
Connection to
28
Minimal input for MM

1) Topological properties:
description of the covalent connectivity of the molecules to be modeled
2) Structural properties:
the starting conformation of the molecule, provided by an X-ray structure,
NMR data or a theoretical model
3) Energetical properties:
a force field describing the force acting on each atom of the molecules
4) Thermodynamical properties:
T0
NV
a sampling algorithm that generates the thermodynamical ensemble that

matchese experimental conditions for the system, e.g. N,V,T , N,P,T, ...
29
MD flowchart
Initial coordinates (X-ray, NMR)

Structure minimization (release strain)
Solvation (if explicit solvent)
Initial velocities assignment
Heating dynamics (Temp. to 300K)
Equilibration dynamics (control of Temp. and structure)
Rescale velocities
Temp. Ok?
Structure Ok?
Production dynamics (NVE, NVT, NPT)

Analysis of trajectory
Calculation of macroscopic values
http://www.ch.embnet.org/MD_tutorial/
30
Coordinate files
Cartesian coordinates
(x,y,z)
ATOM
ATOM
ATOM
ATOM
ATOM
ATOM
ATOM
1
2
3
4
5
6
7
N
HT1
HT2
HT3
CA
HA
CB
VAL
VAL
VAL
VAL
VAL
VAL
VAL
1
1
1
1
1
1
1
-0.008
-0.326
-0.450
-0.172
1.477
1.777
2.038
-0.022
0.545
-0.956
0.566
-0.077
-0.598
-0.740
-0.030
0.778
-0.084
-0.876
0.073
0.971
-1.193
Internal coordinates (IC)
1.00
1.00
1.00
1.00
1.00
1.00
1.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
PEP
PEP
PEP
PEP
PEP
PEP
PEP
Given 4 consecutive atoms A-B-C-D, the IC are:
D
B
RAB, ABC, ABCD, BCD,

1
2
3
4
5
6
7
8
1
1
1
1
1
1
1
1
N
N
N
CG1
CG1
CA
HG11
HG11
1
1
1
1
1
1
1
1
C
C
CA
CA
CA
CB
CB
CB
1
1
1
1
1
1
1
1
*CA
*CA
CB
*CB
*CB
CG1
*CG1
*CG1
1
1
1
1
1
1
1
1
CB
HA
CG1
CG2
HB
HG11
HG12
HG13
1.4896
1.4896
1.4896
1.5421
1.5421
1.5353
1.1099
1.1099
105.11 117.88
105.11 -118.25
108.86 176.05
110.92 121.67
110.92 -118.15
110.92
56.96
111.11 -119.80
111.11 120.81
111.68
108.30
110.92
110.44
109.36
111.11
110.60
110.60
RCD
1.5353
1.0807
1.5421
1.5454
1.1177
1.1099
1.1134
1.1103
It is possible to calculate missing cartesian coordinates

from the existing ones and the IC
31
Definition of atom types

The residue topology file (RTF) contains the atom types
and the standard topology of residues.
Example: the CHARMM force field, version 22
file : top_all22_prot.inp
Atom types section :

MASS
MASS
MASS
MASS
MASS
MASS
[...]
MASS
MASS
MASS
MASS
MASS
MASS
1
2
3
4
5
6
20
21
22
23
24
25
No.
H
HC
HA
HT
HP
HB
C
CA
CT1
CT2
CT3
CPH1
atom
type
1.00800
1.00800
1.00800
1.00800
1.00800
1.00800
H
H
H
H
H
H
!
!
!
!
!
!
polar H
N-ter H
nonpolar H
TIPS3P WATER HYDROGEN
aromatic H
backbone H
12.01100
12.01100
12.01100
12.01100
12.01100
12.01100
C
C
C
C
C
C
!
!
!
!
!
!
carbonyl C, peptide backbone

aromatic C
aliphatic sp3 C for CH
aliphatic sp3 C for CH2
aliphatic sp3 C for CH3
his CG and CD2 carbons
mass
90 atom types
atom
32
Decomposition into residues

In CHARMM, molecules are decomposed into residues.
A molecule may be composed of one to several hundreds or thousands of
residues.
Residues correspond to amino acids for proteins or to nucleotides for DNA.
Topologies for individual residues are pre-defined in CHARMM
Easy construction of the protein topology from the sequence.

Only informations about 20 amino acids needed to construct the
topology of all proteins.
O
NH
NH
O
R
NH
NH
R
NH
R
Residue definition in
CHARMM
33
Residue topology definition
file : top_all22_prot.inp
Total charge
RESI ALA
0.00
GROUP
ATOM N
NH1
-0.47 !
|
ATOM HN
H
0.31 ! HN-N
ATOM CA
CT1
0.07 !
|
HB1
ATOM HA
HB
0.09 !
|
/
GROUP
! HA-CA--CB-HB2
ATOM CB
CT3
-0.27 !
|
\
ATOM HB1 HA
0.09 !
|
HB3
ATOM HB2 HA
0.09 !
O=C
ATOM HB3 HA
0.09 !
|
GROUP
!
ATOM C
C
0.51
ATOM O
O
-0.51
BOND CB CA N HN N CA
BOND C CA C +N CA HA CB HB1 CB HB2 CB HB3
DOUBLE O C
IMPR N -C CA HN C CA +N O
DONOR HN N
ACCEPTOR O C
IC -C
CA
*N
HN
1.3551 126.4900 180.0000
IC -C
N
CA
C
1.3551 126.4900 180.0000
IC N
CA
C
+N
1.4592 114.4400 180.0000
IC +N
CA
*C
O
1.3558 116.8400 180.0000
IC CA
C
+N
+CA
1.5390 116.8400 180.0000
IC N
C
*CA CB
1.4592 114.4400 123.2300
IC N
C
*CA HA
1.4592 114.4400 -120.4500
IC C
CA
CB
HB1
1.5390 111.0900 177.2500
IC HB1 CA
*CB HB2
1.1109 109.6000 119.1300
IC HB1 CA
*CB HB3
1.1109 109.6000 -119.5800
Atom name
Atom type
Atom charge
Bond
Improper
angle
Previous residue
Next residue
115.4200
114.4400
116.8400
122.5200
126.7700
111.0900
106.3900
109.6000
111.0500
111.6100
0.9996
1.5390
1.3558
1.2297
1.4613
1.5461
1.0840
1.1109
1.1119
1.1114
IC
Angles and dihedrals can be generated automatically from this.

34
Patching residues
Treat protein special features.
To make N- and C-termini:
PRES NTER
1.00 ! standard N-terminus
GROUP
! use in generate statement
ATOM N
NH3
-0.30 !
ATOM HT1 HC
0.33 !
HT1
ATOM HT2 HC
0.33 !
(+)/
ATOM HT3 HC
0.33 ! --CA--N--HT2
ATOM CA
CT1
0.21 !
|
\
ATOM HA
HB
0.10 !
HA
HT3
DELETE ATOM HN
BOND HT1 N HT2 N HT3 N
DONOR HT1 N
DONOR HT2 N
DONOR HT3 N
IC HT1 N
CA
C
0.0000 0.0000 180.0000 0.0000
IC HT2 CA
*N
HT1
0.0000 0.0000 120.0000 0.0000
IC HT3 CA
*N
HT2
0.0000 0.0000 120.0000 0.0000
0.0000
0.0000
0.0000
To make disulfide bridges:

PRES DISU
GROUP
ATOM 1CB CT2
ATOM 1SG SM
GROUP
ATOM 2SG SM
ATOM 2CB CT2
DELETE ATOM 1HG1
DELETE ATOM 2HG1
BOND 1SG 2SG
IC 1CA 1CB 1SG
IC 1CB 1SG 2SG
IC 1SG 2SG 2CB
-0.36 ! patch for disulfides. Patch must be 1-CYS and 2-CYS.

! use in a patch statement
! follow with AUTOgenerate ANGLes DIHEdrals command
-0.10 !
-0.08 !
2SG--2CB-!
/
-0.08 ! -1CB--1SG
-0.10 !
2SG
2CB
2CA
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
180.0000
90.0000
180.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
Etc...
35
The parameter file

Contains force field parameters. file : par_all22_prot.inp
Bonds:
BONDS
!
!V(bond) = Kb(b - b0)**2
!
!Kb: kcal/mole/A**2
!b0: A
!
!atom type Kb
b0
!
!Carbon Dioxide
CST OST
937.96
1.1600 ! JES
!Heme to Sulfate (PSUL) link
SS
FE
250.0
2.3200 !force constant a guess
!equilbrium bond length optimized to reproduce
!CSD survey values of
!2.341pm0.01 (mean, standard error)
!adm jr., 7/01
C
C
600.000
1.3350 ! ALLOW ARO HEM
! Heme vinyl substituent (KK, from propene (JCS))
CA
CA
305.000
1.3750 ! ALLOW
ARO
! benzene, JES 8/25/89
CE1 CE1
440.000
1.3400
!
! for butene; from propene, yin/adm jr., 12/95
CE1 CE2
500.000
1.3420
!
! for propene, yin/adm jr., 12/95
CE1 CT2
365.000
1.5020
!
! for butene; from propene, yin/adm jr., 12/95
CE1 CT3
383.000
1.5040
!
! for butene, yin/adm jr., 12/95
CE2 CE2
510.000
1.3300
!
Kb
r0
Idem for angles, dihedrals, impropers, ...

36
Setting up a system in CHARMM

$
charmm < input_file.inp > output_file.out
*
!------ Standard Topology and Parameters
OPEN UNIT 1 CARD READ NAME top_all22_prot.inp
READ RTF CARD UNIT 1
CLOSE UNIT 1
OPEN UNIT 1 CARD READ NAME par_all22_prot.inp
READ PARA CARD UNIT 1
CLOSE UNIT 1
! Generate actual topology
OPEN UNIT 1 READ CARD NAME 1aho-xray.pdb
READ SEQUENCE PDB UNIT 1
GENE 1S SETUP FIRST NTER LAST CTER
REWIND UNIT 1
READ COOR PDB UNIT 1
CLOSE UNIT 1
! Fill IC table and build missing coordinates

IC FILL
IC PARA ALL
IC BUILD
! Build better coordinates for hydrogens
HBUILD SELE TYPE H* END
! Write coordinates in CHARMM format
OPEN UNIT 1 WRITE CARD NAME 1aho.crd
WRITE COOR CARD UNIT 1
CLOSE UNIT 1
! Write coordinates in PDB format
OPEN UNIT 1 WRITE CARD NAME 1aho.pdb
WRITE COOR PDB UNIT 1
CLOSE UNIT 1
! Write Protein Structure File for subsequent use
OPEN UNIT 1 WRITE CARD NAME 1aho.psf
WRITE PSF CARD UNIT 1
CLOSE UNIT 1
! Make disulfide bridges

PATCH DISU 1S 12 1S 63
AUTOgenerate ANGLes DIHEdrals
PSFSUM> Summary of the structure file counters :

Number of segments
=
1
Number of residues
=
Number of atoms
=
962
Number of groups
=
Number of bonds
=
980
Number of angles
=
Number of dihedrals
=
2591
Number of impropers =
Number of cross-terms
=
0
Number of HB acceptors =
98
Number of HB donors =
Number of NB exclusions =
0
Total charge =
1.00000
64
298
1753
169
118
37
The protein structure file (PSF)

The PSF is generated by CHARMM from the sequence of the proteins,
the ligands, the water molecules, etc...,
using the information present in the residue topology file (RTF)
It contains all the information needed for future simulations :
Residues and segments. How the system is divided into
residues and segments.
Atom information. Names, types, charges, masses.
Bond, angle, dihedral and improper dihedral lists
Electrostatic groupings. How some numbers of atoms are grouped
for the purpose of calculating long range electrostatic
A segment is a group of molecules, for example:
one single protein
a collection of water molecules
a collection of ions
a ligand
38
Performing a calculation in CHARMM

Second step : perform calculation, e.g. energy evaluation
input
output
*
CHARMM>

CLOSE UNIT 1
CLOSE UNIT 1
!------ Actual topology
OPEN UNIT 1 READ CARD NAME 1aho.psf
READ PSF CARD UNIT 1
CLOSE UNIT 1
!------ Coordinates
OPEN UNIT 1 READ CARD NAME 1aho.pdb
CLOSE UNIT 1
!------ Energy calculation
ENERGY
!------ End of input file
STOP
ENERGY
NONBOND OPTION FLAGS:

ELEC
VDW
ATOMs
CDIElec SHIFt
VATOm
VSWItch
BYGRoup NOEXtnd NOEWald
CUTNB = 14.000 CTEXNB =999.000 CTONNB = 10.000 CTOFNB = 12.000
WMIN
= 1.500 WRNMXD = 0.500 E14FAC = 1.000 EPS
= 1.000
NBXMOD =
5
There are
0 atom pairs and
0 atom exclusions.
There are
0 group pairs and
0 group exclusions.
<MAKINB> with mode
5 found
2733 exclusions and
2534 interactions(1-4)
<MAKGRP> found
886 group exclusions.
Generating nonbond list with Exclusion mode = 5
== PRIMARY == SPACE FOR
276432 ATOM PAIRS AND
0 GROUP PAIRS
General atom nonbond list generation found:
224678 ATOM PAIRS WERE FOUND FOR ATOM LIST
9439 GROUP PAIRS REQUIRED ATOM SEARCHES
ENER ENR: Eval#
ENERgy
ENER INTERN:
BONDs
ENER EXTERN:
VDWaals
-----------------ENER>
0 -1222.13834
ENER INTERN>
29.79474
ENER EXTERN>
-302.00504
------------------
Delta-E
ANGLes
ELEC
--------0.00000
88.24015
-1285.42224
---------
GRMS
UREY-b
HBONds
--------4.26768
5.92868
0.00000
---------
DIHEdrals
ASP
---------
IMPRopers
USER
---------
239.13668
0.00000
---------
2.18868
0.00000
---------
39
Energy landscape
Landscape for / plane of dialanine
Complex landscape for a protein

E
3N Spatial coordinates
40
Minimization
Landscape for / plane of dialanine
Finding minimum energy conformations

given a potential energy function:
"E
=0
"x i
2
"
E
and
2 >0
"x i
Used to:
relieve strain in experimental conformations
! (energetically)
find
stable states
!
Huge number of degrees of freedom in macromolecular systems.

Huge amount of local minima
Impossible to find true global minimum
Different minimization methods available:
Steepest descent (SD) relieve strain, find close local minimum
Conjugated gradient (CONJ)
Adopted Basis Newton Raphson (ABNR)
Find lower
energy mimima
41
CHARMM input for minimization

*
CLOSE UNIT 1
CLOSE UNIT 1
OPEN UNIT 1 READ CARD NAME val.psf
CLOSE UNIT 1
!------ Coordinates
OPEN UNIT 1 READ CARD NAME val.pdb
CLOSE UNIT 1
!------ ABNR minimization
MINI ABNR NSTEP 200
STOP
ABNER> An energy minimization has been requested.

EIGRNG =
0.0005000
MINDIM =
5
NPRINT =
50
NSTEP =
200
PSTRCT =
0.0000000
SDSTP =
0.0200000
STPLIM =
1.0000000
STRICT =
0.1000000
TOLFUN =
0.0000000
TOLGRD =
0.0000000
TOLITR =
100
TOLSTP =
0.0000000
FMEM
=
0.0000000
MINI MIN: Cycle
ENERgy
Delta-E
GRMS
MINI INTERN:
BONDs
ANGLes
UREY-b
MINI EXTERN:
VDWaals
ELEC
HBONds
---------------------------------MINI>
0
-18.69945
0.00000
3.66846
MINI INTERN>
0.49721
2.67684
0.28823
MINI EXTERN>
7.68221
-36.34723
0.00000
---------------------------------MINI>
50
-26.40712
7.70767
0.60709
MINI INTERN>
0.65246
3.04259
0.34085
MINI EXTERN>
4.90369
-36.38436
0.00000
---------------------------------MINI>
100
-27.90707
1.49995
0.38850
MINI INTERN>
0.84536
3.91951
0.53225
MINI EXTERN>
6.11013
-41.28270
0.00000
---------------------------------UPDECI: Nonbond update at step
103
Generating nonbond list with Exclusion mode = 5
== PRIMARY == SPACE FOR
172 ATOM PAIRS AND
Step-size
DIHEdrals
ASP
--------0.00000
6.34896
0.00000
--------0.00545
0.97825
0.00000
--------0.00279
1.87179
0.00000
---------
ABNER> Minimization exiting with number of steps limit (
ENERgy
BONDs
VDWaals
---------28.09805
0.92928
6.19000
---------
Delta-E
ANGLes
ELEC
--------0.00062
3.80437
-41.86221
---------
GRMS
UREY-b
HBONds
--------0.00826
0.58415
0.00000
---------
0.15433
0.00000
--------0.05940
0.00000
--------0.09659
0.00000
---------
0 GROUP PAIRS
General atom nonbond list generation found:

120 ATOM PAIRS WERE FOUND FOR ATOM LIST
0 GROUP PAIRS REQUIRED ATOM SEARCHES
MINI>
150
-28.09742
0.19035
0.04299
MINI INTERN>
0.93508
3.79489
0.58387
MINI EXTERN>
6.18507
-41.85983
0.00000
---------------------------------MINI>
200
-28.09805
0.00062
0.00826
MINI INTERN>
0.92928
3.80437
0.58415
MINI EXTERN>
6.19000
-41.86221
0.00000
----------------------------------
ABNR MIN: Cycle

ABNR INTERN:
ABNR EXTERN:
---------ABNR>
200
ABNR INTERN>
ABNR EXTERN>
----------
IMPRopers
USER
---------
0.00027
2.19377
0.00000
--------0.00005
2.18679
0.00000
---------
0.06973
0.00000
--------0.06957
0.00000
---------
200) exceeded.
Step-size
DIHEdrals
ASP
--------0.00005
2.18679
0.00000
---------
IMPRopers
USER
--------0.06957
0.00000
---------
42
Simple Molecular dynamics

Newtons law of motion:
Fi = mia i = "
dE i
dri
In discete time : integration algorithm.

!
Example:
Verlet algrorithm
1
r(t + "t) = r(t) + v(t) # "t + a(t) #!"t 2
2
1
r(t " #t) = r(t) " v(t) $ #t + a(t) $ #t 2
2
!
MD algorithm
F(t)
r(t + "t) = 2r(t) # r(t # "t) +
$ "t 2
m
calculate
forces F(t)
t = t + "t
!
update
r(t + "t)
t ~ 1 fs = 10-15 s
Propagation of time: position at time t+dt is a determined by position at time t and t-dt, and by
the acceleration at time t (i.e., the forces at time t)
The equations of motion are deterministic, e.g., the positions and the velocities at time zero
!
determine
the positions and velocities at all other times, t.
43
CHARMM input file for MD simulation

*
CLOSE UNIT 1
CLOSE UNIT 1
OPEN UNIT 1 READ CARD NAME val.psf
CLOSE UNIT 1
Control propagation of time.

Timestep in ps (i.e. 1fs here)
!------ Coordinates
OPEN UNIT 1 READ CARD NAME val.pdb
CLOSE UNIT 1
!------ MD
OPEN WRITE
OPEN WRITE
OPEN WRITE
DYNA
simulation
UNIT 31 CARD NAME traj/dyna.rst
UNIT 32 FILE NAME traj/dyna.dcd
UNIT 33 CARD NAME traj/dyna.ene
VERLET
IUNWRI
IPRFRQ
FIRSTT
IEQFRQ
INBFRQ
Generated files:
Restart file (rst). To restart after
crash or to continue MD sim.
Collection of instantaneous
coordinates along trajectory (dcd)
! Restart file
! Coordinates file
! Energy file
START NSTEP 1000 TIMESTEP 0.001 31 IUNCRD 32 KUNIT 33 100 NPRINT 100 NSAVC 100 NSAVV 100 ISVFRQ 2000 300.0 FINALT 300.0 ICHEW 1 10 IASORS 0 ISCVEL 0 IASVEL 1 ISEED 8364127 10

STOP
Control output
Control temperature
Control non bonded list update
frequency
44
Output of MD simulation
DYNA DYN: Step
DYNA PROP:
DYNA INTERN:
DYNA EXTERN:
DYNA PRESS:
---------DYNA>
0
DYNA PROP>
DYNA INTERN>
DYNA EXTERN>
DYNA PRESS>
----------
Time
GRMS
BONDs
VDWaals
VIRE
--------0.00000
1.79744
0.49721
4.45183
0.00000
---------
TOTEner
HFCTote
ANGLes
ELEC
VIRI
---------6.35935
-6.35862
2.67684
-37.44124
2.21273
---------
TOTKe
HFCKe
UREY-b
HBONds
PRESSE
--------19.31479
19.31700
0.28823
0.00000
0.00000
---------
ENERgy
EHFCor
DIHEdrals
ASP
PRESSI
---------25.67414
0.00074
3.69866
0.00000
0.00000
---------
TEMPerature
VIRKe
IMPRopers
USER
VOLUme
--------381.16244
-3.31909
0.15433
0.00000
0.00000
---------
TOTKe
HFCKe
UREY-b
HBONds
PRESSE
--------16.21830
16.53681
3.70007
0.00000
0.00000
---------
ENERgy
EHFCor
DIHEdrals
ASP
PRESSI
---------8.40109
0.10617
3.31010
0.00000
0.00000
---------
TEMPerature
VIRKe
IMPRopers
USER
VOLUme
--------320.05568
-73.15251
0.21480
0.00000
0.00000
---------
[...]
DYNA DYN: Step
Time
TOTEner
DYNA PROP:
GRMS
HFCTote
DYNA INTERN:
BONDs
ANGLes
DYNA EXTERN:
VDWaals
ELEC
DYNA PRESS:
VIRE
VIRI
-------------------------DYNA>
10
0.01000
7.81721
DYNA PROP>
14.53782
7.92338
DYNA INTERN>
3.28801
14.20105
DYNA EXTERN>
9.54825
-42.66337
DYNA PRESS>
0.00000
48.76834
-------------------------UPDECI: Nonbond update at step
10
[...]
45
Newtonian dynamics in practice

In theory, Newtonian dynamics conserves the total energy (isolated system) :
p
m
p = F(r)
r =
pi2
H (r ,p) = "
+ V(r ) =
i 2mi
cste
In practice, constant energy dynamics is not often used for two reasons :
1) Inaccuracies of the MD algorithm tend to heat up the system

! We can couple the system to a heat reservoir to absorb the excess heat
Integrator
System :
reservoir
cutoffs
constraints
barostat
2) The constant energy dynamics (NVE) does rearely represents the

experimental conditions for the system simulated.
46
Plan
Introduction
Connection to
47
Thermodynamic ensembles
A macroscopic state is described by :
- number of particles :
- volume :
- energy :
N
V
E
- chemical potential :
- pressure :
- temperature :
P
T
Definition: a thermodynamical ensemble is a collection of microscopic

states that all realize an identical macroscopic state
A microscopic state of the system is given by a point (r, p) of the phase
space of the system, where r= (r1, ..., rN ) and p= (p1, ..., pN ) are the positions
and the momenta of the N atoms of the system.
Examples of thermodynamical ensembles:
- Microcanonical:
- Canonical:
- Constant P-T:
- Grand Canonical:
fixed N, V, E
fixed N, V, T
fixed N, P, T
fixed , P, T
often used in MD
often used in MD
48
Bolzmann (canonical) distribution

Boltzmann showed that the canonical probability of the microstate i is given by
1 " #E i
Pi = e
Z
= 1/KBT
KB = Boltzmann constant
Z is the partition function,
Z =$ e" #Ej
such that :
"P =1
i
The partition function is a very complex function to compute, because it represents

a measure of the whole space accessible to the system.
!
Illustration:
If a system can have two unique states, state (1) and state (2),
then the ratio of systems in state (1) and (2) is
P 1 e__E
= _ _E = e_ __E _ E _= e_ __ E
P2 e
1
at 300K, a E of 1.3 kcal/mol

results in a P1/P2 of 1 log10.
(1)
(2)
Cave: if state (1) and (2) are composed of several microscopic states, E G
49
The partition function

The determination of the macroscopic behavior of a system from a
thermodynamical point of view is tantamount to computing the
partition function, Z, from which all the properties can be derived.
Z = $ e" #E i
i
O =
1
# $E i
O
e
" i
Z i
...
Expectation Value
!
"
ln( Z ) = U
E =
"#
Internal Energy
# "ln( Z ) &
p = kT%
(
$ "V ' N,T
Helmoltz free energy
Pressure
A = "kT ln(Z)
50
The ensemble average
Microscopic
E1, P1 ~ e-E1
Macroscopic
Expectation value
O =
Where
!
1
# $E
Oie i
"
Z
Z = $e
E2, P2 ~ e-E2
E3, P3 ~ e-E3
" #E i
is the partition function
E4, P4 ~ e-E4
E5, P5 ~ e-E5
51
Ergodic Hypothesis
The ergodic hypothesis is that the ensemble averages used to compute expectation
values can be replaced by time averages over the simulation.
O
1
Z
Ergodicity
ensemble
1
" #E(r, p )
drdp =
$ O(r, p)e
%
time
$ O(t)dt
0
The microstates sampled by molecular dynamics are usually a small subset

of the entire thermodynamical ensemble.
The validity
! of this hypothesis depends on the quality of the sampling produced by
the molecular modelling technique. The sampling should reach all important
minima and explore them with the correct probability,
- NVE simulations
- NVT simulations
- NPT simulations
Microcanonic ensemble
Canonical ensemble
Isothermic-isobaric ensemble
P = cst.
P(E) = e-E
P(E) = e-(E+PV)
" #E
Note that the Bolzmann weight e

is not present in the time average because
it is assumed that conformations are sampled from the right probability.
52
Ergodic hypothesis : intuitive view

MD Trajectory
E
ensemble
NVE simulation in a
local minimum
3N Spatial coordinates
1
=
Z
$ O(r, p)e
" #E (r, p )
drdp
?=
1
%
$ O(t)dt =
time
Two main requirements for MD simulation :
1) Accurate energy function E(r,p)

2) Appropriate algorithm, which
- generates the right ensemble

- samples efficiently
53
Sampling of the various ensembles

1)
2)
Microcanonical ensemble (constant N,V,E)

sampling is obtained by simple integration
of the Newtonian dynamics:
- Verlet, Leap-Frog, Velocity Verlet, Gear
Canonical ensemble (constant N,V,T)
sampling is obtained using thermostats :
1)
Berendsen: scaling of velocities to
obtain an exponential relaxation of the
temperature to T
2)
Nose-Hoover: additional degree of
freedom coupled to the physical system
acts as heat bath.
NVE
NV
T
(infinite E
reservoir)
fixed P
3)
Isothermic-isobaric ensemble (constant N,P,T)

In addition to the thermostat, the volume of the
system is allowed to fluctuate, and is regulated
by barostat algorithms.
T
(infinite E
reservoir)
54
The Nos-Hoover thermostat

Phase space extended by two extra variables :
physical variables
Newton
friction term
temperature regulation
One can demonstrate that
the
canonical distribution is reproduced
for the physical variables
Conserved quantity :
Non-Hamiltonian dynamics...
55
Other sampling methods I

Langevin Dynamics (LD)
In Langevin Dynamics, two additional forces are added to the standard force field:
- a friction force:
-i pi
whose direction is opposed to the velocity of atom i
- a stochastic (random) force:
(t)
such that <(t)> = 0.
This leads to the following equation for the motion of atom i:
ri =
pi
mi
p i = Fi (r) + "pi + # (t)
This equation can for example simulate the friction and stochastic effect of the
solvent in implicit solvent simulations. The temperature is adjusted via and ,
using the dissipation-fluctuation
! theorem.
The stochastic
term can improve barrier crossing and hence sampling.
!
LD does not reproduce dynamical properties

56
Other sampling methods II

Monte Carlo Simulations and the Metropolis criterion
In this sampling method, instead of computing the forces on each atom to
solve its time evolution, random movements are assigned to the system and
the potential energy of the resulting conformer is evaluated.
To insure Boltzmann sampling, additional criteria need to be applied on the
new conformer. Let C be the initial conformer and C' the randomly modified:
- if V(C') < V(C), the new conformer is kept and C' becomes C for next step
- if V(C') > V(C), a random number, , in the [0,1] interval is generated and
if e-(V'-V) > , the new conformer is kept and C' becomes C for next step
Using this algorithm, one insures Boltzmann statistics,
P(C ")
= e# $ (V "#V )
P(C)
57
Plan
Introduction
Connection to
58
The importance of entropy

Mechanics:
A state is characterised by one minimum energy structure

(global minimum)
Statistical mechanics:
A state is characterised by an ensemble of structures
Very small energy differences between states (~kBT = 2.5 kJ/mol)

resulting from summation over very many contributions
Entropic effects :
Not only energy minima are of importance but whole range of

x-values with energies ~kBT
energy
U(x)
The free energy (F)

governs the system
kBT
S(x1)
S(x2)
F = U - TS
U(x2)
x2 may have higher energy but lower
free energy than x1
U(x1)
x1
x2
coordinate x
Energy (U) entropy (S)

compensation at finite
temperature T
59
Dynamical behavior of proteins

Biological molecules exhibit a wide range of time scales over which
specific processes occur; for example:
Local Motions (0.01 to 5 , 10-15 to 10-1 s)
Atomic fluctuations
Sidechain Motions
Loop Motions
Rigid Body Motions (1 to 10 , 10-9 to 1 s)
Helix Motions
Domain Motions
Subunit motions
Large-Scale Motions (> 5 , 10-7 to 104 s)

Helix coil transitions
Dissociation/Association
Folding and Unfolding
60
Types of problems
Molecular dynamics simulations permit the study of complex, dynamic
processes that occur in biological systems. These include, for example:
Protein stability
Conformational changes
Protein folding
Molecular recognition: proteins, DNA, membranes, complexes
Ion transport in biological systems
and provide the mean to carry out the following studies,

Drug Design
Structure determination: X-ray and NMR
61
Historical perspective
Theoretical milestones:
Newton (1643-1727):
Boltzmann(1844-1906):
Schrdinger (1887-1961):
Classical equations of motion: F(t)=m a(t)

Foundations of statistical mechanics
Quantum mechanical eq. of motion: -ih t (t)=H(t) (t)
Molecular mechanics milestones:

Metropolis (1953):
Rahman (1964):
First MD simulation of proteins

The CHARMM general purpose FF & MD program
The AMBER general purpose FF & MD program
First full QM simulations
First QM-MM simulations
Proteins
Karplus (1977) &

McCammon (1977)
Karplus (1983):
Kollman(1984):
Car-Parrinello(1985):
Kollmann(1986):
Liquids
Wood (1957):
Alder (1957):
First Monte Carlo (MC) simulation of a liquid

(hard spheres)
First MC simulation with Lennard-Jones potential
First Molecular Dynamics (MD) simulation of
a liquid (hard spheres)
First MD simulation with Lennard-Jones potential
62
System sizes
Simulations in explicit solvent
BPTI (VAC), bovine pancreatic
trypsin inhibitor without solvent
BPTI, bovine pancreatic
trypsin inhibitor with solvent
RHOD, photosynthetic reaction
center of Rhodopseudomonas
viridis
HIV-1, HIV-1 protease
ES, estrogenDNA
STR, streptavidin
DOPC, DOPC lipid bilayer
RIBO, ribosome
Solid curve, Moores law
doubling every 28.2 months.
Dashed curve, Moores law
doubling every 39.6 months.
K.Y. Sanbonmatsu, C.-S. Tung, Journal of Structural Biology 157(3) : 470, 2007
63
Largest system 2006

satellite tobacco mosaic virus
1 million atoms
Simulation time : 50 ns
system size : 220 A
P. L. Freddolino, A. S. Arkhipov, S. B. Larson, A. McPherson, and K. Schulten, Molecular dynamics

simulations of the complete satellite tobacco mosaic virus, Structure 14 (2006), 437.
64
The tradeoff we can afford
65
Exemple : Aquaporin
Selective translocation of water across a membrane
S. Hub and Bert L. de Groot. Mechanism of selectivity in aquaporins and aquaglyceroporins PNAS. 105:1198-1203 (2008)
66
Exemple : Aquaporin
Selective translocation of water across a membrane
Propose a model for
selectivity at the atomic
level
S. Hub and Bert L. de Groot. Mechanism of selectivity in aquaporins and aquaglyceroporins, PNAS 105:1198-1203 (2008)
67

Molecular Dynamics

Uploaded by

Document Informationclick to expand document information

Copyright:

Available Formats

Molecular Dynamics

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Molecular Dynamics

Uploaded by

Copyright:

Available Formats

Molecular Dynamics Simulation

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

High pressure simulation

Some properties cannot be

Simulation is a useful complement, because it can :

What is it good for?

Classical Molecular Mechanics (MM)

Mixed Quantum/Classical (QM/MM)

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

chemical reactions, enzyme catalysis

atomic (solute + solvent)

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

Processes: Thermodynamic Equilibria

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

Definition of a model for molecular simulation

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

Classical force fields

Implementation of calssical potential energy functions

Covalent bonds and angles

E angle = K" (" # " 0 )

Dihedrals and Improper torsions

E dihedral = K" [1+ cos(n" # $ )]

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

Van der Waals interactions

where is the dielectric constant :

The Coulomb energy decreases only as 1/r

1) Hydrogen bonds (Hb)

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

The total potential energy function

# K% [1+ cos(n% " &)] +

# K' (' " '0 )

For a system with 1500 atoms

Cutoff for non-bonded interactions

Three cutoff schemes: strict, shift, switch

E'(r) = E(r) " S(r)

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

Force field parametrization

#K [1+ cos(n% " &)] + # K (' " ' )

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

Hydrogen bonds between water molecules and

Taken into account via:

Explicit solvation schemes

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008

Periodic boundary conditions

N=4,8. The dielectric constant is a function of atom distance.

NBOND RDIE EPS 4.0

Implicit solvent models

f(r) : electrostatic potential,

Equation solved numerically. Very time consuming.

Others: EEF1, SASA, etc...

Explicit hydrogen bonds with water molecules are lost!

Introduction to molecular surfaces

"Gnp, solv = #SASA + b

= 0.00542 kcal mol 1 2

Rotating probe: radius 1.4

ensemble of van der Waals sphere centered at each atom

Molecular Dynamics Simulation - Michel Cuendet - EMBL 2008