Atomistische

Simulationen

Vorlesungen zur Molekularen Dynamik

Prof. M. Meuwly

Departement Chemie

Universitaet Basel

1. Einführung

2. Elektronenstruktur

3. Kraftfelder

4. Optimierungsstrategien

5. Molekulare Dynamik

6. Statistische Mechanik von Proteinen

• Fundamentales Verständnis atomistischer Prozesse

• Interpretation von Experimenten komplexer Systeme kann sehr schwierig sein.

• Gefährliche Experimente

• Aufwendige Experimente

• Stellen weiterführender Fragen

Warum Simulationen?

Leitmotiv The underlying physical laws necessary for the mathematical theory of a

large part of physics and the whole of chemistry are thus completely known,

and the difficulty is only that the exact application of these laws leads to equations

much too complicated to be soluble. It therefore becomes desirable that

approximate practical methods of applying quantum mechanics should be

developed, which can lead to an explanation of the main features of complex

atomic systems without too much computation.

Chemie-Nobelpreise 1998 (Kohn und Pople) und 2013 (Karplus, Levitt und Warshel)

P. A. M. Dirac (1929)

Water Transport in Channels

Water Transport in Channels

Blue curve: water-water interaction; green curve: protein-water interaction (H-bonding);

black curve: sum of green and blue.

Chaperones

Proteine sind “liquid-like”

F. Schotte et al., Science 300, 1944 (2003)

IACS, Kolkata February 2014

Charge Exchange in N2+ + N2 → N2 + N2

+

Tong et al., Chem. Phys. Lett. (2012)

Statistics of several 1000 reactive trajectories from exciting v9=5 and v9=6

Vibrationally Induced Dissociation of Sulfuric Acid

Reyes and Meuwly JPCA (2011)

Empirical Force Fields

and their optimization

Ziel: Beschreibung der Wechselwirkungen in

einem komplexen System.

Problem: Anzahl Freiheitsgrade (3N mit N104)

Ansatz: Ausgangspunkt ist Struktur des Systems

Bindungen, Valenzwinkel, Diederwinkel

Parametrisierung der Wechselwirkung

PEF/PES: Potential Energie Funktion

Potential Energy Surface

Topologie (“Landkarte”) der Wechselwirkung

Niedrig-dimensionale Darstellung der WW

Projektion!

Komplizierte “Landschaft”

Potential energy function

(mathematical equations)

Empirical force field

(equations and parameters that relate chemical structure

and conformation to energy)

Class I

CHARMM

AMBER

OPLS/AMBER/Schrödinger

ECEPP (free energy force field)

GROMOS

Class II

CFF95 (Biosym/Accelrys)

MM3

MMFF94 (CHARMM, Macromodel, elsewhere)

UFF, DREIDING

Common empirical force fields

Class I Potential Energy Function

Intramolecular (internal, bonded terms)

2,3,13,1

2

22)cos(1

o

BradleyUrey

UBo

impropers

torsions

o

angles

o

bonds

b

rrKK

nKKbbK

Intermolecular (external, nonbonded terms)

qiq j

4Drij ij

nonbonded

Rmin,ij

rij

12

2Rmin,ij

rij

6

Class II force fields (e.g. MM3, MMFF, UFF, CFF)

Kb, 2 b bo 2 Kb, 3 b bo

3 Kb, 4 b bo

4 bonds

K , 2 o 2 K , 3 o

3 K , 4 o

4 angles

K ,1 1 cos K , 2 1 cos 2 K , 3 1 cos 3 dihedrals

K

impropers

2

Kbb'bonds'

bonds

b bo b 'bo ' K '

angles'

angles

o ' o '

Kbangles

bonds

b bo o

b bo K , b1 cos K , b2 cos 2 K , b3 cos 3 dihedrals

bonds

b 'bo ' K , b '1 cos K , b' 2 cos 2 K , b ' 3 cos 3 dihedrals

bonds'

o K , 1 cos K , 2 cos 2 K , 3 cos 3 dihedrals

angles

o ' o ' cosdihedrals

angles'

angles

Equilibrium terms

bo: bonds

o: angles

n: dihedral multiplicity

o: dihedral phase

o: impropers

r1,3o: Urey-Bradley

Force constants

Kb: bonds

K: angles

K: dihedral

K: impropers

KUB: Urey-Bradley

Intramolecular parameters

Kbbonds

b bo 2 K

angles

o 2 K

torsions

1 cos(n )

K

impropers

o 2 KUBUreyBradley

r1,3 r1,3,o 2

Vbond Kb bbo 2

1 2

3 4

Vangle K o 2

VdihedralK (1 (cosn ))

Intermolecular interactions between bonded atoms

1,2 interactions: 0

1,3 interactions: 0

1,4 interactions: 1 or scaled

> 1,4 interactions: 1

H

C

HH

H

VimproperK o 2

VUreyBradley KUB r1,3 r1,3o 2

Bond Energy versus Bond length

0

100

200

300

400

0.5 1 1.5 2 2.5

Bond length, Å

Pote

nti

al En

ergy, kcal/

mol

Chemical type Kbond bo

C-C 100 kcal/mole/Å2 1.5 Å

C=C 200 kcal/mole/Å2 1.3 Å

C=-C 400 kcal/mole/Å2 1.2 Å

Vbond Kb bbo 2

Dihedral energy versus dihedral angle

0

5

10

15

20

0 60 120 180 240 300 360

Dihedral Angle, degrees

Pote

nti

al En

ergy, kcal/

mol

VdihedralK (1 (cosn ))

= 0˚

qi: partial atomic charge

D: dielectric constant

: Lennard-Jones (LJ, vdW) well-depth

Rmin: LJ radius (Rmin/2 in CHARMM)

Combining rules (CHARMM, Amber)

Rmin,ij = Rmin,i + Rmin,j

i,j = SQRT(i * j )

Intermolecular parameters

qiq j

4Drij ij

nonbonded

Rmin,ij

rij

12

2Rmin,ij

rij

6

Electrostatic Energy versus Distance

-100

-80

-60

-40

-20

0

20

40

60

80

100

0 1 2 3 4 5 6 7 8

Distance, Å

In

tera

cti

on

en

erg

y,

kca

l/m

ol

q1=1, q2=1

q1=-1, q2=1

Lennard-Jones Energy

-0 .5

-0 .3

-0 .1

0 .1

0 .3

0 .5

0 .7

0 .9

1 2 3 4 5

Distance, Å

In

tera

cti

on

En

erg

y,

kca

l/m

ol

Series 1

Rmin,ij

i,j

ijRmin,ij

rij

12

2Rmin,ij

rij

6

Lennard-Jones Energy versus Distance

-0.5

-0.3

-0.1

0.1

0.3

0.5

0.7

0.9

1 2 3 4 5 6 7 8

Distance, Å

In

tera

cti

on

En

erg

y,

kca

l/m

ol

e=0.2,Rmin=2.5

Rmin,i,j

eps,i,j

Extent of Parameter Optimization

Minimal optimization

by analogy (i.e. direct transfer of known parameters)

Maximal optimization

time-consuming

requires appropriate target data (expt, calculations)

Choice based on goal of the calculations Minimal

database screening

NMR/X-ray structure determination

Maximal

free energy calculations (perturbations, potential of mean force)

mechanistic studies

subtle environmental effects

lead optimization

Intermolecular Optimization

Partial Atomic Charges

VDW Parameters

Intramolecular Optimization

Bonds

Angles

Torsions

Impropers, Urey-Bradley

Parameter Optimization Complete

if intermolecular and intramolecular changes < convergence criteria

Initial Geometry

if intermolecular change > conv.crit.

if intermolecular microscopic and macroscopic change < convergence criteria

if intramolecular change > conv.crit.

if intramolecular and intermolecular change > conv.crit.

1) Identify previously parameterized

compounds

2) Access topology information

i) Assign atom types

ii) Connectivity (bonds)

iii) Charges

CHARMM topology (parameter files)

top_all22_model.inp (par_all22_prot.inp)

top_all22_prot.inp (par_all22_prot.inp)

top_all22_sugar.inp (par_all22_sugar.inp)

top_all27_lipid.rtf (par_all27_lipid.prm)

top_all27_na.rtf (par_all27_na.prm)

top_all27_na_lipid.rtf (par_all27_na_lipid.prm)

top_all27_prot_lipid.rtf (par_all27_prot_lipid.prm)

top_all27_prot_na.rtf (par_all27_prot_na.prm)

toph19.inp (param19.inp)

From top_all22_model.inp

RESI PHEN 0.00 ! phenol, adm jr.

GROUP

ATOM CG CA -0.115 !

ATOM HG HP 0.115 ! HD1 HE1

GROUP ! | |

ATOM CD1 CA -0.115 ! CD1--CE1

ATOM HD1 HP 0.115 ! // \\

GROUP ! HG--CG CZ--OH

ATOM CD2 CA -0.115 ! \ / \

ATOM HD2 HP 0.115 ! CD2==CE2 HH

GROUP ! | |

ATOM CE1 CA -0.115 ! HD2 HE2

ATOM HE1 HP 0.115

GROUP

ATOM CE2 CA -0.115

ATOM HE2 HP 0.115

GROUP

ATOM CZ CA 0.11

ATOM OH OH1 -0.54

ATOM HH H 0.43

BOND CD2 CG CE1 CD1 CZ CE2 CG HG CD1 HD1

BOND CD2 HD2 CE1 HE1 CE2 HE2 CZ OH OH HH

DOUBLE CD1 CG CE2 CD2 CZ CE1

HG will ultimately be deleted.

Therefore, move HG (hydrogen) charge

into CG, such that the CG charge

becomes 0.00 in the final compound.

Use remaining charges/atom types

without any changes.

Do the same with indole

Top_all22_model.inp contains all protein

model compounds. Lipid, nucleic acid and

carbohydate model compounds are in the full

topology files.

Intermolecular Optimization Target Data

Local/Small Molecule

Experimental

Interaction enthalpies (MassSpec)

Interaction geometries (microwave, crystal)

Dipole moments

Quantum mechanical

Mulliken Population Analysis

Electrostatic potential (ESP) based

CHELPG (g98: POP=(CHELPG,DIPOLE)

Restricted ESP (AMBER)

Dimer Interaction Energies and Geometries (OPLS, CHARMM)

Global/condensed phase (all experimental)

Pure solvents (heats of vaporization, density, heat capacity, compressibility)

Aqueous solution (heats/free energies of solution, partial molar volumes)

Crystals (heats of sublimation, lattice parameters, interaction geometries)

Additive Models: account for lack of explicit

inclusion of

polarizability via “overcharging” of atoms.

RESP: HF/6-31G overestimates dipole moments (AMBER)

Interaction based optimization (CHARMM, OPLS)

local polarization included

scale target interaction energies (CHARMM)

1.16 for polar neutral compounds

1.0 for charged compounds

Partial Atomic Charge Determination

For a particular force field do NOT change the QM level of theory. This

is necessary to maintain consistency with the remainder of the force

field.

Comparison of analogy and optimized charges

N

NHO

Name Type Analogy Optimized

C1 CT3 -0.27 -0.27

H11 HA3 0.09 0.09

H12 HA3 0.09 0.09

H13 HA3 0.09 0.09

C2 C 0.51 0.58

O2 O -0.51 -0.50

N3 NH1 -0.47 -0.32

H3 H 0.31 0.33

N4 NR1 0.16 -0.31

C5 CEL1 -0.15 -0.25

H51 HEL1 0.15 0.29

C6 CT3 -0.27 -0.09

H61 HA 0.09 0.09

H62 HA 0.09 0.09

H63 HA 0.09 0.09

Geometries (equilibrium bond, angle, dihedral, UB and improper terms)

microwave, electron diffraction, ab initio

small molecule x-ray crystallography (CSD)

crystal surveys of geometries

Vibrational spectra (force constants)

infrared, raman, ab initio

Conformational energies (force constants)

microwave, ab initio

Intramolecular optimization target data

Note that the potential energy surface about a given torsion is the sum of the

contributions from ALL terms in the potential energy function, not just the dihedral

term

Addition of simple functional groups is generally

straightforward once the full compound parameters

have been optimized.

C

H

HH

N

H

H

N

H

HH

C

O

O

O

H

C

O

NH2

O

CH3F

Lead Optimization

Summary: Force Fields

1) Junk in, junk out: Parameter optimization effort

based on application requirements.

2) Follow standard protocol for the force field of

interest (higher level QM is not necessarily

better).

3) Careful parameter optimization of lead molecules

4) Simple substitutions often require minimal or no

optimization.

Differences between Force Fields

Differences between Force Fields

Differences between Force Fields

Molecular Dynamics

Simulations

MD: Verlet Method

Newton’s equation represents a set of N second order

differential equations which are solved numerically at

discrete time steps to determine the trajectory of each

atom.

Advantage of the Verlet Method:

requires only one force evaluation per timestep

Energy function:

used to determine the force on each atom:

MD: Velocity Verlet Method

Advantage of Velocity Verlet:

No differences

Velocities available directly

tttattvttv

ttatvttv

ttattvtrttr

2

121

2

121

2

1 2

Molecular Dynamics

Ensembles

Constant number of particles, energy, volume (NVE)

(microcanonical)

Constant number of particles, temperature, volume (NVT)

(canonical)

Constant number of particles, temperature, pressure (NPT)

(isothermal-isobaric)

Constant temperature, volume, chem.Potential (mVT)

(Grand canonical)

Transformation between different ensembles via Legendre Transformation

Steps in

Molecular Dynamics Simulations

1) Build realistic atomistic model of the system

2) Simulate the behavior of your system over time using

specific conditions (temperature, pressure, volume,

etc)

3) Analyze the results obtained from MD and relate to

macroscopic level properties

Example: KscA channel

solvent

solvent

KcsA channel protein

(in blue) embedded in a

(3:1) POPE/POPG

lipid bilayer. Water

molecules inside the

channel are shown

in vdW representation.

Summary of simulations:

• protein/membrane system contains 38,112 atoms,

including 5117 water molecules, 100 POPE and 34 POPG

lipids, plus K+ counterions

• CHARMM26 forcefield

• periodic boundary conditions, PME electrostatics

• 1 ns equilibration at 310K, NpT

• 2 ns dynamics, NpT

Program: NAMD2

Platform: Cray T3E (Pittsburgh Supercomputer Center)

Simulating the system:

Free MD

MD Results

RMS deviations for the KcsA protein and its selectivity filer indicate that the protein is stable during the

simulation with the selectivity filter the most stable part of the system.

Temperature factors for individual residues in the four monomers of the KcsA channel protein indicate

that the most flexible parts of the protein are the N and C terminal ends, residues 52-60 and residues 84-

90. Residues 74-80 in the selectivity filter have low temperature factors and are very stable during the

simulation.

Simulation Procedure Overview

• Setup

1. PDB file

– Protein Databank (http://www.rcsb.org/pdb/)

2. PSF file

– Generated specifically for the molecule

– Contains the detailed composition and

connectivity of the molecule(s) of interest

Simulation Procedures

• Setup

1. Topology file

– information for putting molecules together,

such as what atoms are to be used, which of

these atoms are bonded to each other, and

the sets of atoms that form bond angles

2. Parameter file

– physical parameters (force constants, van der

Waals forces, bonds, angles, etc.)

Simulation Procedures

• Solvation

– Create water box or shell to enclose the

molecule

• Minimization

– Minimize the energy of the system in

order to reach the most favorable

configuration

Simulation Procedures

• Heating

– Initial velocities are assigned at a low

temperature. Periodically, new velocities

are assigned at a slightly higher

temperature and the simulation is allowed

to continue. This is repeated until the

desired temperature is reached.

Simulation Procedures

• Equilibration

– The point of the equilibration phase is to

run the simulation until the structure,

pressure, temperature and energy become

stable with respect to time.

Simulation Procedures

• Dynamics

– Normal/Periodic boundary condition

– Single/Multiple time stepping

– Integrators

– Electrostatics

Simulation Procedures

0obs )(

1lim dttAAA

time

Computer Simulations generate information at the

atomistic/microscopic level (positions, velocities, etc.).

Conversion of this data to observable/macroscopic

quantities (pressure, internal energy, infrared spectra,

etc.) is domain of statistical mechanics. In MD the

system evolves in time. Connection between

computer simulation and experimental observable is

provided by Gibbs’ theorem:

It states that the ensemble average approaches the

time average for infinite simulation time (complete,

ergodic sampling)

Converting simulations to information

– Mean energy

– RMS difference between two structures

Converting simulations to information

NVTB

V ETk

C 2

2

1

Temperature

Converting simulations to information

Specific heat

N

i i

i

BB mNkNk

KT

1

2

3

1

3

2 p

Diffusion coefficient

0

03

1tdtD ii vv

Infrared spectrum ttidtI μμ 0exp

Converting simulations to information

Visualization

VMD

MolMol

Shortcomings of MD

• Quality of the forcefield

• Size and Time: atomistic simulations can be performed only

for systems of a few tens of angstroms (length scale) and

for a few nanoseconds (time scale).

• Conformational freedom of the molecule: the number of

possible conformations a molecule can adopt is enormous,

growing exponentially with the number or rotatable bonds.

• Only applicable to systems that have been parameterized

• Connectivity of atoms: can not change during dynamics, i.e.

no chemical reactions