Force Fields and Molecular Dynamics
1
Mathieu Linares – PDC, KTH�linares@kth.se
2
I) Introduction to Molecular Mechanics and Molecular Dynamics
3
Non relativistic, time independent Schrödinger equation
Quantum chemistry
Born Oppenheimer approximation
Electronic Equation
Equation for nuclear motion on the potential energy surface
Force-Field Theory
HΨ (r,R) = EΨ (r,R)
[Tn + U]Φ(R) = EΦ(R)
Heψe (r; R) = Ee(R) ψe (r; R)
H = Tn + Te + Vnn + Ven + Vee
Semiempirical fit to the PES cad U(R):
Force field (U)
with U = Ee + Vnn
4
Conformational changes = movements on the Potential Energy Surface (PES)
Simulation of the PES
Morphology Study: Theory
Exploration of the PES
Local Algorithms
Global Algorithms
Mechanics
Dynamics
Monte Carlo
V(r)
r
V(r)
r
5
Truncated Taylor’s expansion:
E
E0
How is built a force field ?
∙ No coupling between internal coordinates
⇒ only diagonal terms
r1
r2
θ
⇒ Separation of the energy contributions
Valence Terms Non-bonded Terms Other
∙ Zero of energy arbitrary ⇒ E0 = 0
⇒ Relative Energies
∙ E0 Equilibrium ⇒ first derivatives = 0
∙ Small deformations close to the minimum
⇒ only the quadratic term (n = 2)
6
1) Bond Stretching
∙ Small deformations around the minimum
⇒ only the quadratic form is considered (n = 2)
∙ For bond lengths unusually long
⇒ also terms beyond the quadratic one (n = 3,4)
Rem: at large deformations cubic functions invert. Disaster during an energy minimization from a poor starting geometry.
Morse Fonction :
with
+ : Finite energy of bonding dissociation D0
Description of a large range of situations
- : Three parameters for each bond, calculation time
Vibrational frequency
Reduced mass
or
7
2) Bending
often, n = 2 ⇒ harmonic potential
From a truncated Taylor expansion :
-: Large deformation not well reproduced
Here, the energy must go to infinity
Problem at this point
⇒ Fourier Expansions in cosine:
+: singularity-free derivatives
reproduce large angular deformations
8
3) Torsion Angles
Fourier Expansion in cosine
(generally up to 6 terms):
A[1+ cos(nϕ - ϕ0]
ϕ0
ϕ
ϕ0
ϕ0
½ barrier height
Periodicity
Phase
Note that Erotational barrier = Etorsion + ......
Torsions
Soft degree of freedom
Bond-stretching and angle-bending
Hard degrees of freedom
Responsible of conformational changes; important to reproduce them correctly
9
4) Inversion
Reproduce the tendency of I to stay in the JKL plane
∙ Fourier expansion in cosine:
ω = 0° ⇒ minimum
ω = 180° ⇒ maximum
ω = 0° ⇒ maximum
ω = 90° ⇒ minimum
With ω the angle between the IL bond and the IJK plane
IJ ILK
IK IJL
CH2=CH2
PH3
∙ Harmonic potential (for planar systems):
10
5) Valence Cross-Terms
Coupling between internal coordinates
Off-diagonal terms, generally functions of two coordinates with atoms in common.
(Δθ)2
ΔθΔr2
ΔθΔr1
Δθ
Δr2Δθ
(Δr2)2
Δr2Δr1
Δr2
Δr1Δθ
Δr1 Δr2
(Δr1)2
Δr1
Δθ
Δr2
Δr1
r1
r2
θ
+ : Description of the vibrational frequencies, conformational changes
- : Calculation time
11
6) van der Waals interactions
Lennard-Jones Potential in 6-12:
-: Repulsive term too “hard”
Lennard-Jones in 6-9:
+: Decrease the hardness of the repulsive term
-: Too attractive at large distance
Induced dipole – induced dipole attraction
Electron repulsion
due to exchange forces
with D0 the well depth
Exponential-6:
-: more calculation time than a power of r
12
7) Electrostatic Interactions
A) Coulombic interactions between pairs of fractional point charges
⇒ Better description of the electrostatic properties
Different methods to attribute partial charges
Ex: bond increments
δij = charge transfer between the i and j atoms
⇒ Partial charge in i:
B) Bond dipole representation
μi
χ
μj
αi
αj
Rij
Unit conversion
Which value to use in the simulations (microscopic level) ?
Sometimes, a relative dielectric constant proportional to r is used
13
8) Description of hydrogen bonds
LJ potentials and electrostatic interactions have a radial symmetry
of the hydrogen bonds
⇒ Difficult analysis of the results.
Explicit energy term with an angular
dependency (CHARMM and Dreiding)
With static charges, no ability of polarization
AMBER uses a LJ10-12 potential, going to zero much more quickly than a LJ6-12
D
A
Donor directionality
ϕ
C=O
H
Acceptor directionality
14
∙ Valence terms proportional to N (number of atoms)
How to manage the non-bonding terms?
∙ Non-bonding terms proportional to N2 (i.e. n(n-1)/2 each)
⇒ Calculation time quickly with the dimension of the system
⇒ Cutoff distance. Beyond, no interactions
System with 5000 atoms.
Number of non-bonding terms in
function of the cutoff distance
Slow the decrease of the potential to avoid discontinuities in potential and forces
Avoid periodical discontinuities when the list of interactions is updated
Lennard-Jones, fast decrease with the distance: 1/r6
Electrostatic interactions, slow decrease: 1/r
⇒ Large cutoff distance (minimum 10(+4) Å)
∙ Generally excluding 1,2 and 1,3 interactions
15
Sources of parameterization data
Molecular structures
X-ray diffraction on crystals
- standard deviation σ ≥ 0.01 Ǻ
- determine the position of electron densities
⇒ X-H bonds underestimated (not an issue for other elements)
Neutron diffraction on crystals
Electron diffraction of gas-phase
- mesures the average distance between nuclear positions
⇒ vibrational effects (bond elongation due to anharmonicity)
differences in a few 0.001 Ǻ
Microwave spectroscopy of gas-phase
- determines molecular moments of inertia from which distances
are determined (+ vibrational corrections)
NMR spectroscopy on solution
Force constants
Vibrational spectroscopy
Molecule | w(cm-1) | Kij (kcal/mol.Ǻ2) |
H2 | 4401 | 821 |
HI | 2309 | 449 |
I2 | 215 | 248 |
16
Sources of parameterization data
Energy differences and barriers
NMR coalescence data
Vibrational spectroscopy
Temperature dependence of thermodynamic data
Charges
Electric moments (dipole, quadrupole, octupole,…) to fit atomic partial charges
-: the dipole moment being a Cartesian vector, cannot be used when > 2 atoms
problem with lone pairs
vdW
Attraction: atomic polarizability
Repulsion: crystal data
Theoretical methods can also provide all these information
+
Detailed shape of the PES
For 2 atoms:
17
Ex: Determination of the torsion parameters
3) The relative torsion energy for a given angle ϕ is thus:
The energy profile for
is fit by a Fourier cosine series:
Erotational barrier = Etorsion + ......
Responsible of the conformational changes
(overcoming the rotational barriers)
⇒ Pay extreme attention to the torsion energy profiles
A method
18
Parameterization approaches
A) Simple assignment
Assume that force constant are indeed constant and independent of atom type
Admit that data and time are lacking
Ex: Dreiding: k = 700 kcal/mol.Ǻ2 to all stretching terms
k = 100 kcal/mol.rad2 to all bending terms
Assume that stretching and bending are less important than torsions and non-bonding interactions in determining the conformations
B) Fitting experimental and theoretical data
Adjust parameters to reproduce experimental observables or data obtained by other theoretical methods (quantum chemistry) such as molecular structure, vibrational frequencies, internal rotation and inversion barriers, conformational energy difference
-: availability and quality of data
19
How many parameters ?
The equilibrium values for each bond length, angle, … in a molecule are different (unless they are related by symmetry).
However, for practical purposes, a force-field will have a minimum number of transferable parameters, and bond length and angles are given standard transferable values.
One standard bond length between primary, secondary, and tertiary carbon atoms
One single force constant for all deformations of a given type
The equilibrium geometry is reached by allowing the molecule to relax; it is the result of the contribution of different energy terms
Harmonic stretching potential
Nonbonded and other energy terms
Stretching potential in strained hydrocarbon
20
Classification of the Force Fields
Generic Force Fields
1) Parameters limited to simple rules (derived from electronegativity, atomic size)
2) Simple potential energy functions
3) Not practical for specific task of high precision.
Classical Force Fields of first generation
1) Highly parameterized (many atom types)
Better prediction for specific applications, mainly in biochemistry Limited to a given combination of atoms, not easy to fit a few new parameters to a new system (each parameter derived with some dependence on other parameters)
2) Simple functional forms similar to those of the first category
⇒ AMBER, CHARMM
Force Fields of second generation
21
Ex: UFF
Atom types are classifications based on element and bonding environment
+
+
+
+
+
Expression of the Force-Field
Minimal number of terms, no cross terms, few parameters
22
Minimal number of terms,
No cross terms,
Many parameters
(3) Dihedrals + maintain the tetrahedral nature of sp3 centers with united-atoms
(6) augments the electrostatic description of HB, by adding about 0.5 kcal mol-1 to the HB energy (option)
Ex: AMBER
The atom types in AMBER are quite specific to amino acids and DNA bases.
23
United atoms: Atoms removed (H) by adding their mass to the atoms they are connected to
⇒ Remove hydrogen vibrations (increase time step in MD)
24
Ex: CFF family
CFF91 is useful for hydrocarbons, proteins, protein-ligand interactions. Parameterization based on quantum mechanics calculations and molecular simulations
Many terms,
Cross terms,
Many parameters
25
Limits and advantages of the force-fields
∙ Electronic Transitions
(photon absorption)
∙ Charge transport properties
∙ Acid-base reactions
∙ Reactions breaking or
building bonds
Limitations:
Advantages:
∙ Conceptual simplicity
∙ Computational quickness
∙ Study of large systems (> 100000 atoms)
∙ Decompose the potential energy into its individual components
Presence probability
Advantages to have different force-fields
Large accuracy computational resources available
∙ Allow studying a large number of systems
∙ Allow determining the degree of independence of the result with respect to the
force-field
∙ Weight the benefits of a highly parameterized FF to the utility of a more generic,
thus more extendable FF
∙ Different functional forms increase the simulation flexibility
26
Conformational changes ≡ movements on the Potential Energy Surface (PES)
Force Fields : semi-empirical fit with equations derived from classical mechanics
Simulation of the PES
Morphology Study: Theory
Exploration of the PES
Local Algorithms
Global Algorithms
Mechanics
Dynamics
V(r)
r
27
Molecular Mechanics – Principle of Energy Minimization
B) Conformational Changes
Exemple:
A) Energy Evaluation
Initial coordinates: a (x0, y0)
Energy Expression: E(x,y) = x2 + 5 y2
⇒ initial direction: gradient in a: ∇E = (2x, 10y)
(one-dimensional minimization)
⇒ Adapt α to minimize E (x',y') (iterative process)
At the minimum (c), the search line is tangent
to the energy contour
⇒ the next derivative vector is orthogonal the previous one
28
Line search direction along the local gradient -∇E(xi, yi)
∙ Direction change when the one-dimensional minimum is reached
Molecular Mechanics – Minimization Algorithms
A) Steepest descent
∙ Direction change when a point of lower energy is reached on the line search
12 iterations
8 (typically 3-10) evaluations of E/iteration
12 iterations
1.3 evaluations of E/iteration
+: Very robust
- : Slow convergence close to the minimum: gradient close to 0; oscillating directions, when each segment of the path is orthogonal to the previous one, thus tending to revert the progression obtained in a previous iteration
The truncated line search reduces a lot the number of function evaluation per iteration (gain of calculation time)
29
Base: To accelerate the convergence: at each iteration, the energy is minimized on a hyperspace that contains all the previous search directions.
Direction of the line search: local gradient + previous search line:
hi+1
hi
gi+1
gi+1, is orthogonal to g0, g1, g2, ..., gi,
hi+1, is conjugated to h0, h1, h2, ..., hi
or
+: In theory, the solution is found in N iterations (N freedom degrees)
Molecular Mechanics – Minimization Algorithms
B) Conjugate gradient
with
30
Molecular Mechanics – Minimization Algorithms
C) Newton-Raphson
Methods using first and second derivatives
Local minimum when
Minimization = iterative determination of x for which f’(x)= 0
Taylor expansion …
= 0
Line approximation
With 3n – 6 coordinates:
i = 1 to 3n-6
-: storage of the (3n)2 second derivative matrix elements
determination of the (3n)2 second derivatives
obtaining the inverse of the matrix
Time consuming
In cartesian coordinates
Steepest descent
31
Molecular Mechanics – Minimization Algorithms
Comparison of the methods: minimization of lactic acid
CPU time/step
SD: 0.08 s
CG: 0.22 s
NR: 0.99 s
Not yet converged after 500 steps
CPU time to reach
E = 14 kcal/mol
SD: 4.11 s
CG: 3.07 s
NR: 4.95 s
SD or CG
To finish off the minimization: Newton-Raphson
In one shot: CG a good compromise
32
∙ Derivative of the function
- Root Mean Square (RMS) (+: better take into account the large derivatives)
- Average
- Highest value
∙ Energy changes between each iteration
∙ Coordinate changes between each iterations
Other termination criterion: Maximum number of iterations
To reach the minimum, the number of iterations depends on the algorithm, energy expression, dimension of the system, and convergence criteria.
⇒ Minimization stopped when convergence criteria met
Molecular Mechanics – Convergence criteria
How close to the absolute convergence is close enough?
Depends on the aim:
Relaxation of overlapping atoms before starting a dynamics: 1.0 kcal mol-1Å-1
Analysis of normal modes: < 10-5 kcal mol-1Å-1
with n atoms
33
Conformational changes ≡ movements on the Potential Energy Surface (PES)
Force Fields : semi-empirical fit with equations derived from classical mechanics
Simulation of the PES
Morphology Study: Theory
Exploration of the PES
Local Algorithms
Global Algorithms
Mechanics
Dynamics
V(r)
r
300 K
Solve the Newton’s equations of motion
V(r)
r
Generate Statistical Ensemble (NVT, NPT,...)
Study Molecular Motions
Conformational search
34
Verlet “Leap Frog” Algorithm
v(t-1/2Δt)
r(t)
a(t)
v(t+1/2Δt)
r(t+Δt)
or
Integration ⇒ r, v, a of the atoms with time: Trajectory
⇒ Solved numerically (method of finite difference)
1) r0.
2) v0 generated randomly, according to the target temperature.
3) r, v, a determined at a time t+Δt from their value at time t.
Newton’s equation:
Molecular Dynamics – Integration Algorithms
Criteria of a good integrator in MD:
• Quickness, ideally only one energy evaluation per Δt.
• Low memory requisites.
• Using a large Δt.
• Good energy conservation.
-: Calculation of the positions and velocities out of synchrony
35
Molecular Dynamics – Selection of the Time Step
Large trajectory ⇒ Large Δt
Verlet’s Supposition: v and a constant within Δt.
⇒ Decompose the atomic motion ⇒ small Δt.
Vibrations C-H: 10-14 s.
⇒ Δt of 10-15 s.
⇒ Many iterations.
⇒ Time consuming, depends on the dimension of the system and the time scale of the dynamic process investigated
Local Motions: de 0.01 to 5 Å, 10-15 to 10-1 s.
Atomic Fluctuations, motions of small molecular segments.
Medium Scale Motions: from 1 to 10 Å, 10-9 to 1s (106 and 1015 time steps).
Molecular parts locally organized,
Rigid bodies: helices and domains in proteins.
Large Scale Motions: > 5 Å, 10-7 to 104 s (108 and 1019 time steps).
Folding and Unfolding of proteins,
Dissociation/Association of Molecular Complexes.
36
II) Conformational Search
37
Conformational Search – Meaning of the energy minimum
The calculated energy is relative
Eethylamine
Ebiphenyl
E1biphenyl
E2biphenyl
The calculated energy is the classic enthalpy at 0K
(ignoring the vibrational motions at 0K)
Boltzmann equation
Population of each conformation, Pi:
Or the relative populations of two conformations:
38
Conformational Search – Selection of a procedure
In practice:
Small molecules: we can obtain the global minimum.
Large Molecules: A few conformers for understanding the conformational space
Many minima ?? ⇒ many minimizations to explore the whole potential energy surface and obtain the global minimum
39
Conformational Search - Systematic Exploration
Modify a few selected torsion angles according to a grid of equally spaced values.
N = 32
Ex. :
Vary ϕ1, ϕ 2; each one being given the following values: 60°, 180° and 300°
ϕ1 | 30° | 30° | 30° | 60° | 60° | 60° | 90° | 90° | 90° |
ϕ2 | 30° | 60° | 90° | 30° | 60° | 90° | 30° | 60° | 90° |
⇒ Limit the number of torsions to vary
and the number of values given to the torsions.
+: Efficient to localize the global minimum.
-: Becomes quickly time consuming with the number of angles to vary.
60°
180°
300°
number of torsion values
number of torsions to vary
40
Relative Energy (kcal/mol)
Torsion angle (degrees)
Exemple 1: Conformational search of a polyphenylene dendrimer
monomer
In biphenyl, a 360°rotation generates 4 minima (~40°, ~140°, ~220°, ~320°) and 4 transition states (0°, 90°, 180° y 270°):
But for peripheral phenyl groups, two orientations are redundant
212x43, i.e. 262,144 structures are generated!!
45°
225°
135°
315°
Conformational Search - Systematic Exploration
41
Relative energy (kcal/mol)
Torsion angle (degrees)
monomer
How to reduce the number of conformations generated?
By eliminating the variations of the peripheral phenyl groups ⇒ 43, i.e. 64 structures.
Due to symmetry, one monomer is given one orientation ⇒ 42, i.e. 16 structures.
Exemple 1: Conformational search of a polyphenylene dendrimer
Conformational Search - Systematic Exploration
When optimizing the structure as built, one conformer is generated, showing that the orientations of the phenyl groups are correlated
42
Exemple 2: poly(propylene imine) dendrimer
R = NH2
60°
180°
300°
There is no key to determine the redundant combinations of angles
⇒ Use another method
Conformational Search - Systematic Exploration
117 torsion angles (around the CC and CN bonds).
Bonds between sp3 atoms ⇒ 3 values separated by 120°.
⇒ obtain 3117, i.e. 6.5 1050 structures.
How to reduce the number of conformations generated?
Conclusion about the systematic explorations:
Very efficient for systems with - a few atoms/freedom degrees
- steric constraints
Not practical for studying large flexible systems
43
Heating at high temperature (800K)
for overcoming the potential energy barriers.
Cooling to low temperature (5 K)
relaxation of the system, to avoid being trapped
in a local minimum of high energy.
Energy Minimization
leads to a local minimum.
MD
MM
Conformational Search – Simulated Annealing
Combination of the MD and MM methods
44
Example: Conformational search of poly(benzylether dendrimers)
Conformational Search – Simulated Annealing
Ideal for large flexible molecules:
⇒ Generating conformers among the most stable ones.
Not practical for rigid molecules:
⇒ Limited conformational changes.
conformers
Relative Energy (kcal/mol)
In both cases: Extraction of information during the exploration of the PES