atomistische simulationen vorlesungen zur molekularen dynamikmeuwly/pdfs/nanoiii.pdf ·...
Post on 20-Feb-2020
3 Views
Preview:
TRANSCRIPT
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
top related