C7790 Introduction to Molecular Modelling -1C7790
Introduction to Molecular
Modelling
TSM Modelling Molecular Structures
Petr Kulhánek
kulhanek@chemi.muni.cz
National Centre for Biomolecular Research, Faculty of Science
Masaryk University, Kamenice 5, CZ-62500 Brno
PS/2021 Present Form of Teaching: Rev1
Lesson 25
Molecular Dynamics I
C7790 Introduction to Molecular Modelling -2-
Context
microworldmacroworld
equilibrium (equilibrium constant)
kinetics (rate constant)
free energy
(Gibbs/Helmholtz)
partition function
phenomenological thermodynamics
statistical thermodynamics
microstates
(mechanical properties, E)
states
(thermodynamic properties, G, T,…)
microstate ≠ microworld
Description levels (model chemistry):
• quantum mechanics
• semiempirical methods
• ab initio methods
• post-HF methods
• DFT methods
• molecular mechanics
• coarse-grained mechanics
Structure
EnergyFunction
Simulations:
• molecular dynamics
• Monte Carlo simulations
• docking
• …
C7790 Introduction to Molecular Modelling -3System
Evolution in Time
t
)(tM
How to simulate time evolution of the system?
𝑀 =
1
𝑡𝑡𝑜𝑡
න
𝑜
𝑡 𝑡𝑜𝑡
𝑀(𝑡)𝑑𝑡
snapshots of the system are a microstates
Mechanical Description (classical physics)*:
➢ Newton's laws of motion
➢ First law states that an object at rest will stay at rest, and an object in motion
will stay in motion unless acted on by a net external force.
➢ Second law states that the acceleration of a body over time is directly
proportional to the force applied and occurs in the same direction as the
applied force.
➢ Third law states that all forces between two objects exist in equal magnitude
and opposite direction.
𝑭𝑖 = 𝑚𝑖 𝒂𝑖
the system is composed
of N atoms
* time evolution can also be described by QM but at cost of theoretical and computational
complexity
C7790 Introduction to Molecular Modelling -4-
Forces
Origin of interatomic forces:
R
E(R)
R
E(R)➢ the most stable local
configuration
➢ the system at rest
➢ the system tends to
reach more energy
favorable configuration
𝜕𝐸(𝑹)
𝜕𝑹
= 0
𝜕𝐸(𝑹)
𝜕𝑹
≠ 0
Interatomic forces:
𝑭𝑖 = −
𝜕𝐸(𝑹)
𝜕𝒓𝑖
negative value of potential energy gradient is force
total potential energy
force acting on atom i position of atom i
The only forces that can act on atoms in the system are from interatomic interactions.
C7790 Introduction to Molecular Modelling -5Equation
of Motions
𝑭𝑖 = 𝑚𝑖 𝒂𝑖 𝑭𝑖 = −
𝜕𝐸(𝑹)
𝜕𝒓𝑖
Second Newton's Law Forces in molecular systems
Final equations of motions (EM):
𝑚𝑖 𝒂𝑖 = −
𝜕𝐸 𝑹
𝜕𝒓𝑖
𝑚𝑖
𝑑2
𝒓𝑖
𝑑𝑡2
= −
𝜕𝐸 𝑹
𝜕𝒓𝑖
To describe evolution of the system in time, it is necessary to solve system of N (number of
atoms) second order differential equations or motions.
𝑹 𝒕 = 𝒓1(𝑡), 𝒓2(𝑡), … , 𝒓 𝑁(𝑡)
Result: position of atoms in time (trajectory)
C7790 Introduction to Molecular Modelling -6Numerical
Integration
The solution of EM can be obtained by integration of differential equations. Unfortunately,
the analytical solution is not feasible even for small systems (three and more atoms).
Numerical integrations
➢ Finite difference methods
➢ leap-frog algorithm (a variant of Verlet algorithm)
➢ velocity Verlet algorithm
➢ Gear corrector-predictor methods
➢ Runge-Kutta methods
most often used algorithms in MD simulations
of (bio)chemical systems
C7790 Introduction to Molecular Modelling -7Leap-frog
algorithm
𝒓 𝑡 + 𝑑𝑡 = 𝒓 𝑡 + 𝒗(𝑡 + 𝑑𝑡/2) ∙ 𝑑𝑡
1) Initial conditions:
2) Molecular dynamics (MD loop)
1) Calculation of forces and accelerations
2) Update velocities
1) Update positions
𝒗 𝑡 + 𝑑𝑡/2 = 𝒗 𝑡 − 𝑑𝑡/2 + 𝒂(𝑡) ∙ 𝑑𝑡
𝒓 𝑡 ; 𝒗 𝑡 − 𝑑𝑡/2
𝒂 𝑡 =
1
𝑚
𝜕𝐸 𝑹(𝑡)
𝜕𝒓
velocities
positions
accelerations
𝒂 𝑡 =
𝒗 𝑡 + 𝑑𝑡/2 − 𝒗 𝑡 − 𝑑𝑡/2
𝑑𝑡
𝒗 𝑡 + 𝑑𝑡/2 =
𝒓 𝑡 + 𝑑𝑡 − 𝒓 𝑡
𝑑𝑡
𝑡
𝑑𝑡
time
time step (integration step)
C7790 Introduction to Molecular Modelling -8Time
Step
𝑑𝑡
numerical accuracycomputational efficiency
➢ The time step size is usually taken as 1/10 of the fastest motions.
➢ The fastest motions are X-H vibrations (higher PES curvature, light atom (hydrogen)).
➢ Then, the typical size of the integration step is 1 fs (10-15 s)
Strategies how to increase the integration time step:
➢ remove the fastest motions by the constraining X-H distances, which allows a 2-fs
step size
➢ SHAKE, RATTLE, SETTLE, LINCS algorithms
➢ in addition, constrain valence angles (mathematically too complex, not use)
➢ hydrogen mass repartitioning (up to 4 fs)
➢ multiple time-step integrators
➢ computationally cheap short-range forces (shorter integration time step)
➢ computationally expensive long-range forces (longer integration time step)
C7790 Introduction to Molecular Modelling -9HMR
- Hydrogen Mass Repartitioning
Hopkins, C. W.; Le Grand, S.; Walker, R. C.; Roitberg, A. E. Long-Time-Step Molecular
Dynamics through Hydrogen Mass Repartitioning. J. Chem. Theory Comput. 2015,
11 (4), 1864–1874. https://doi.org/10.1021/ct5010406.
Since the molecular dynamics uses the classical physics, each degree of freedom is
thermalized to
1
2
𝑘 𝐵 𝑇, which is independent to atom masses.
C7790 Introduction to Molecular Modelling -10Initial
Conditions, Equilibration
For integration of EM, we need initial geometry and velocities:
𝒓 𝑡 ; 𝒗 𝑡 − 𝑑𝑡/2
initial geometry (structure) of model
velocities can be generated randomly to satisfy MaxwellBoltzmann
distribution for given temperature
The initial geometry of models derived from experimental structures (X-RAY, NMR,
CryoEM, etc.) is usually of low quality.
Equilibration:
➢ The aim of the equilibration is to bring model to desired thermodynamic state
(temperature, pressure, density, etc.).
Initial
model
Equilibrated
model
temperature adjustment
density adjustment
production
dynamics
analysis
C7790 Introduction to Molecular Modelling -11Thermostats
and Barostats
TB (target temperature)
PB (target pressure)
simulated system
thermal bath
pressure bath
(usually virtual)
T (temperature)
P (pressure)
At equilibrium:
𝑇 = 𝑇𝐵
P = 𝑃𝐵
The capital P (pressure) is employed to
avoid ambiguities with the small p
(momentum).
C7790 Introduction to Molecular Modelling -12-
Thermostats
Equipartition principle:
3𝑁 − 𝑐
2
𝑘 𝐵 𝑇 = 𝐸 𝑘 = ෍
𝑖=1
𝑁
1
2
𝑚𝑖 𝒗𝑖
2
degrees of freedom (DOF)
(c - constrained DOF) temperature
mean kinetic energy
The temperature can be controlled by a thermostat:
➢ weak coupling thermostat, Berendsen thermostat
➢ simple, incorrect ensemble
➢ temperature is controlled by velocity scaling
➢ dangerous to use for simulations in vacuum
➢ susceptible to various artefacts (flying ice cube, etc.)
➢ Langevin thermostat (stochastic, correct ensemble)
➢ it thermalizes each degree of freedom by random collisions
➢ Nosé-Hoover barostat (correct ensemble)
C7790 Introduction to Molecular Modelling -13-
Barostats
Virial theorem (Clausius 1870):
2 𝐸 𝑘 = − ෍
𝑖=𝑁
𝑁
𝑭𝑖 𝒓𝑖
time average of
kinetic energy
the virial (it reflects
potential energy)
NpT ensemble (ideal gas model):
2 𝐸 𝑘 = −3𝑃𝑉 − ෍
𝑖=𝑁
𝑁
𝑭𝑖 𝒓𝑖 𝑃 =
1
𝑉
𝑁𝑘 𝐵 𝑇 −
1
3
෍
𝑖=𝑁
𝑁
𝑭𝑖 𝒓𝑖
pressure
The pressure can be controlled by a barostat:
➢ weak coupling barostat (simple, incorrect ensemble)
➢ Monte-Carlo barostat (stochastic, correct ensemble)
➢ Nosé-Hoover barostat (correct ensemble)
the pressure change is achieved by changing the size of the simulation box.
from equipartition principle
C7790 Introduction to Molecular Modelling -14Output
from MD
NSTEP = 60000 TIME(PS) = 3970719.913 TEMP(K) = 300.74 PRESS = -209.4
Etot = -68239.6682 EKtot = 14401.7354 EPtot = -82641.4035
BOND = 225.2198 ANGLE = 498.4282 DIHED = 712.6121
1-4 NB = 259.0139 1-4 EEL = -4287.9441 VDWAALS = 10154.5260
EELEC = -90203.2595 EHBOND = 0.0000 RESTRAINT = 0.0000
EKCMT = 6899.5998 VIRIAL = 7977.8787 VOLUME = 238531.9588
Density = 1.0214
actual temperature (K) and pressure (atm)
(they DO NOT represent thermodynamical properties)
the same property, which is the actual kinetic energy,
expressed in different units
C7790 Introduction to Molecular Modelling -15Output
from MD, cont.
A V E R A G E S O V E R
NSTEP = 5000000 TIME(PS) = 3990599.911 TEMP(K) = 299.90 PRESS = 2.0
Etot = -68232.2253 EKtot = 14361.2461 EPtot = -82593.4714
BOND = 248.6326 ANGLE = 517.2225 DIHED = 724.2102
1-4 NB = 253.7846 1-4 EEL = -4299.0145 VDWAALS = 10259.0679
EELEC = -90297.3769 EHBOND = 0.0000 RESTRAINT = 0.0023
EAMBER (non-restraint) = -82593.4736
EKCMT = 6859.7668 VIRIAL = 6849.2815 VOLUME = 238422.2094
Density = 1.0219
R M S F L U C T U A T I O N S
NSTEP = 5000000 TIME(PS) = 3990599.911 TEMP(K) = 1.53 PRESS = 153.2
Etot = 17.0060 EKtot = 73.4906 EPtot = 75.4230
BOND = 12.9980 ANGLE = 16.9666 DIHED = 10.9737
1-4 NB = 5.9275 1-4 EEL = 19.2799 VDWAALS = 126.9605
EELEC = 159.5253 EHBOND = 0.0000 RESTRAINT = 0.0359
EAMBER (non-restraint) = 75.3871
EKCMT = 57.0484 VIRIAL = 788.8336 VOLUME = 174.8435
Density = 0.0007
thermodynamic temperature (K) and pressure (atm)
Thermostat: T = 300 K (weak coupling)
Barostat: p = 1 atm (weak coupling)