Transcript Amber 8

Introduction to Amber
The theory and practice of biomolecular
simulations using the Amber suite of
programs
Dr. Vladislav Vassiliev
NCI National Facility, The Australian National University,
ACT 0200, Canberra, Australia
February 2011
Presentation Outline
Introduction to Amber 12
Hands-on
• Setting up a standard Amber MD Run
• Building non-standard Residues
• QM/MM: Using Amber-Gaussian
Interface
• QM/MM: Using Amber inbuilt QM
methods
What is AMBER?
AMBER
Assisted
Model
Building with
Energy
Refinement
What is Amber?
“Amber” refers to two things:
1) a set of molecular mechanical force fields for the
simulation of biomolecules
2) a package of molecular simulation programs (about 50 )
which includes source code and demos
The current version of the code is Amber version 12, which is
distributed by UCSF (University of California, San
Francisco) subject to a licensing agreement
Amber Home Page: http://ambermd.org/
What is Amber
Amber is distributed in two parts:
AmberTools12 and Amber 12:
•
AmberTools12 could be used without Amber12, but not vice
versa
• AmberTools12 currently consists of several independently
developed packages that work well by themselves, and with
Amber itself
• Amber 12 centered around the sander and pmemd simulation
programs and continues to be licensed as before, under a more
restrictive license (Academic/non-profit/government: $400.
Industrial (for-profit): $20,000 for new licensees, $15,000 for
licensees of Amber 10).
AmberTools
NAB
build molecules; run MD or distance geometry,
etc.
antechamber & MCPB
Create force fields for general organic molecules
ptraj & cpptraj
Analyze trajectories from Amber or CHARMM
tleap and xleap
Basic preparation program for Amber simulations
3D-RISM
Solves integral equation models for solvation
sqm
semiempirical and DFTB quantum chemistry
pbsa
Performs numerical solutions to PoissonBoltzmann models
Mdgx
Code for explicit solvent molecular dynamics
simulations
MMPBSA.py & amberlite
Energy-based analyses of MD trajectories
AmberTools
• AmberTools is released under the GNU General Public License
(GPL)
• A few components are included that are in the public domain or
which have other, open-source, licenses.
• AmberTools is distributed in source code format, and must be
compiled in order to be used. One needs C, C++, and fortran
compilers to compile the AmberTools programs.
• The source code of AmberTools could be obtained here:
http://ambermd.org/AmberTools-get.html
Versions of Amber
Version
Released
12
2012
11
2010
10
2008
9
2006
8
2004
AMBER Home
• Have a look at the Amber Home Page: http://ambermd.org/
Amber Main References
A general overview of the Amber codes:
D.A. Case, T.E. Cheatham, III, T. Darden, H. Gohlke, R. Luo, K.M.
Merz, Jr., A. Onufriev, C. Simmerling, B. Wang and R. Woods. The
Amber biomolecular simulation programs. J. Comput. Chem. 26,
1668-1688 (2005)
An overview of the Amber protein force fields, and how they were developed:
W. Ponder and D.A. Case. Force fields for protein simulations. Adv. Prot.
Chem. 66, 27-85 (2003).
E. Cheatham, III and M.A. Young. Molecular dynamics simulation of
nucleic acids: Successes, limitations and promise. Biopolymers 56,
232-256 (2001).
What is Amber?
“Amber” is a software package for
modelling of Large Molecular Systems
Why Do We Need A Special Treatment
for Large Molecular Systems?
Or, Why Do We Need Amber?
Quantum Chemistry Methods Provide a
Rigorous Description of Molecular Systems
They solve Schrödinger equation
And they are generally applicable:
But… they are very time consuming…
To Treat Large Molecular Systems We Need
to Reduce the Complexity of the System
Molecular Mechanics is a non-quantum
mechanical technique for treating Large
Molecular Systems
As a result Molecular Mechanics methods are
thousands times faster than Quantum
Chemistry methods
Molecular Mechanics vs
Quantum Mechanics
Quantum Mechanics
Molecular Mechanics
Considers atoms as collections of
electrons and nuclei
Considers atoms as soft or hard spheres.
Covalent bonds are treated as springs
Solves quantum Schrödinger equation
Uses classical potential energy
equations
Force Fields
The potential energy equations to calculate the energy in
Molecular Mechanics methods and the parameters/constants
used in the equations are known as a Force Field
There are many force fields designed for different purposes.
QVBMM
SIBFA
UFF
MM2, MM3, MM4
COSMOS-NMR
AMBER
DRF90
PIPF
OPLS
MMFF
ECEPP/2
ENZYMIX
GROMACS
X-Pol
CVFF
CHARMM
CFF
CHARMm
QCFF/PI
GROMOS
AMOEBA
Amber Force Field
The total Energy in Amber force field consists of
1. bonded terms relating to atoms linked by
covalent bonds and
2. nonbonded terms describing the long-range
electrostatic and van der Waals interactions:
Etotal = Ebonded + Enonbonded
Amber Force Field: Bonded Terms
The bonded Energy in Amber force field consists of
bond stretching, angle bending, and torsion terms:
Ebonded = Estretch + Ebend + Etorsion
Amber Force Field: Bond Stretching
Amber force field treats covalent bonds between atoms
as springs (Hooke's Law, F = -kx)

K
(
r

r)

2
bond
stretching
r
eq
bonds
E
where Kr is the empirical stretching force constant, r is the
actual bond length and req is the “natural” (empirical) bond
length
Amber Force Field: Bond Stretching
Amber Force Field: Angle Bending
Amber force field treats angles that are bonded to the
same central atom as springs (Hooke's Law, F = -kx)

K
(


)

2
angle
bending

eq
angles
E
where Kθ is the empirical bending force constant, θ is the actual
bond angle and θeq is the “natural” (empirical) bond angle
Amber Force Field: Angle Bending
Amber Force Field: Torsion Energy
Torsion Energy: torsional (dihedral) angle rotation between
atoms that are vicinal (bonded to adjacent atoms) to each other

V
n


E

(
1

cos
n

)

dihedral
2
dihedrals
where Vn is the barrier to free rotation for the “natural” bond,
n is the periodicity of the rotation (number of cycles in 360°),
φ is the torsion angle and γ is the angle where the potential
passes through its minimum value
Amber Force Field: Torsion Energy
Amber Force Field: Nonbonded Terms
The nonbonded Energy terms in Amber force field
describe the long-range electrostatic and van der
Waals interactions:
Enonbonded = Eelectrostatic + EvdW
Amber Force Field: Electrostatic Energy
The Electrostatic Energy in the Amber force field
represents the pair-wise sum of the electrostatic
energies of all possible interacting non-bonded atoms i
and j:
q
q
atoms
i j
electrosta
tic
i
j
ij
E


R
where qi and qj are the point
charges on atoms, Rij is the
interatomic distance and ε is the
dielectric constant
Amber Force Field: van der Waals Energy
The van der Waals Energy in the Amber force field
represents the pair-wise sum of the van der Waals
energies of all possible interacting non-bonded atoms i
and j:
B A

ij
  
E

6
R R
ij

atoms
ij
vdW
12
i
j
ij
where the Aij and Bij parameters
control the depth and position
(interatomic distance) of the
potential energy well for a given
pair of non-bonded interacting
atoms and Rij is the interatomic
distance
Amber Force Field
- Empirical Parameters
Where Do Empirical
Parameters Come From?
Parameter Derivation: Partial Charges
In AMBER:
Connolly
Connolly
1) Partial atomic charges are
static
2) Quantum chemical methods
(B3LYP/ccpVTZ//HF/631G**) are used to generate
an electrostatic potential
(ESP) around a molecule on
the spheric grid
3) RESP (Restrained Electrostatic
Potential) Method is used to derive
the partial charges
QM = ab initio, DFT, semi-empirical
Parameter Derivation:
Van der Waals Parameters
It is the most difficult part…
1) Optimizing van der Waals parameters to reproduce the
experimental or high-level Quantum Chemical data
Could be computationally expensive
2) Optimizing van der Waals parameters through the
Monte Carlo or MD simulations to reproduce the
experimental properties of bulk solvent (density, etc.).
For example, OPLS van der Waals parameters
Could be computationally expensive
3) Reusing existing van der Waals parameters for similar atom types from the
same or other force field
The simplest approach
Parameter Derivation:
Bond and Angle Interactions
req and θeq come either from experimental data (X-ray, neutron diffraction) or
Quantum Chemical calculations (geometry optimization)
Kr and Kθ force constants are usually optimized to reproduce the vibration
frequencies calculated using high-level Quantum Chemical methods.
Or (the simplest approach)
Kr and Kθ force constants could be derived from the existing bond/angle
parameters for similar bond/angle types from the same or other force field
Parameter Derivation:
Dihedral Angle Interactions
Vn, n, and γ are derived to reproduce
the rotational profile from the highlevel Quantum Chemical calculations.
Or (the simplest approach)
Vn, n, and γ could be derived from
the existing dihedral angle
parameters for similar dihedral angle
types from the same or other force
field
J.Wang et al., Development and testing of a
general amber force field, Journal of
Computational Chemistry, 25 (2004), 1157
Force Fields in Amber 12
ff99SB
ff10
ff12SB
Proteins
ff99 +
backbone torsion
Modifications
no change from
ff99SB
ff99SB + new
backbone and
sidechain torsions
DNA
ff99
ff99 +
“Barcelona”
backbone torsion
modifications
no change from
ff10
RNA
ff99
ff99 +
no change from
“Barcelona”
ff10
backbone
changes + “OL3”
J.Wang et al., Development and testing of a
changes
for c
general amber force field, Journal of
Computational Chemistry, 25 (2004), 1157
Force Fields in Amber 12
Lipid11: A modular lipid force field
A new modular force field for the simulation of
phospholipids and cholesterol designed to be compatible
with the other pairwise additive Amber force field
J.Wang et al., Development and testing of a
general amber force field, Journal of
Computational Chemistry, 25 (2004), 1157
Other Force Fields in Amber:
Inclusion of Polarization
Non-additive" force fields based on atom-centered dipole
polarizabilities can also be used.
These add a "polarization" term to what was given above
where μi is an induced atomic dipole.
In addition, charges that are not centered on atoms, but are off-center (as for
lone-pairs or "extra points") can be included in the force field.
Other Force Fields in Amber: AMOEBA
(Atomic Multipole Optimized Energetics for Biomolecular Applications)
• Atomic Multipoles: The model uses a polarizable atomic multipole description of
electrostatic interactions. Multipoles through the quadrupole are assigned to each atomic
center based on a distributed multipole analysis (DMA) derived from large basis set
molecular orbital calculations at the MP2/aug-cc-pVTZ level and the experimental geometry
of the gas-phase monomer.
• Polarization is treated via self-consistent induced atomic dipoles. Atomic dipole
polarizabilities can be derived from an empirical fit to experimentally known molecular
polarizabilities.
The induced dipole at each atomic site is computed as
where αi is the atomic polarizability and Ei,α is the sum of the fields generated by both
permanent multipoles and induced dipoles
Other Force Fields: AMOEBA
(Atomic Multipole Optimized Energetics for Biomolecular Applications)
• The functional forms for bond stretching and angle bending were taken from
the MM3 force field:
• A Urey-Bradley functional form was chosen for the stretch-bend term:
Other Force Fields: AMOEBA
(Atomic Multipole Optimized Energetics for Biomolecular Applications)
• Repulsion-Dispersion. The buffered 14-7 potential has been applied to model
pairwise additive vdW interactions
where εij is the potential well depth, ρij= Rij/R0ij with Rij as the i-j separation and R0ij the
minimum energy distance. n = 14, m = 7, δ = 0.07, γ=0.12. The combining rules are:
• The buffered 14-7 function yields a repulsive region softer than the Lennard-Jones
6-12 function but steeper than typical Buckingham exp-6 formulations.
• The buffered 14-7 form was found to outperform Lennard-Jones and Buckingham
potentials in simultaneously reproducing gas phase ab initio results and liquid
thermodynamic properties of noble gases and a series of diatomic species.
What we can do with Amber?
Molecular Dynamics Simulations
The Molecular Dynamics simulation method is based
on Newton’s second law or the equation of motion,
F=ma,
where F is the force exerted on the particle, m is its
mass and a is its acceleration
Integration of the equations of motion then yields a
trajectory that describes the positions, velocities and
accelerations of the particles as they vary with time.
From this trajectory, the average values of properties
can be determined.
MD: Melting of Ice
Human carboxyl
esterasecomplexed with morphine
MD: Translocation of DNA
This movie shows the electrophoretically-driven translocation of a 58-nucleotid DNA strand through the
transmembrane pore of alpha-hemolysin
Molecular Dynamics:
Amber MD Workhorses
SANDER - Simulated Annealing with NMR-Derived Energy
Restraints
PMEMD - Particle Mesh Ewald Molecular Dynamics
SANDER
PMEMD is up to 55% faster
than SANDER
GPU
PMEMD
Molecular Dynamics:
What Are Current Simulation
Capabilities?
Time scales of biological processes
Femtosecond (fs)
Picosecond (ps)
Nanosecond (ns)
Microsecond (μs)
= 10-15 second
= 10-12 second
= 10-9 second
= 10-6 second
Molecular Dynamics Trajectory
Molecular Dynamics Simulation
Time
Snapshots
Snapshots of
Representative
Structures
Molecular dynamics trajectory is a file
containing snapshots of the simulated system
For each Snapshot Amber Saves
Structure and Energy Decomposition
Energy Decomposition for each snapshot is written in the form:
NSTEP =
Etot
BOND
1-4 NB
EELEC
EKCMT
100
TIME(PS) =
-73238.8859 EKtot
654.2822 ANGLE
817.8912 1-4 EEL
-109230.1680 EHBOND
7732.1695 VIRIAL
10.200 TEMP(K) =
297.23 PRESS = -1257.4
=
=
18131.6434 EPtot
=
-91370.5294
=
=
1929.4666 DIHED
=
775.1417
=
=
4242.6763 VDWAALS
=
9440.1805
=
=
0.0000 RESTRAINT =
0.0000
=
=
16916.4636 VOLUME
=
338304.5955
Density
=
0.9012
Ewald error estimate:
0.1550E-03
------------------------------------------------------------------------------
Snapshots of
Representative
Structures
Plotting Molecular Dynamics
Properties
Equilibration step
MD Phase
Equilibration step
allows atoms and
molecules to find more
natural positions with
respect to one another
During the MD Phase
molecular properties
(structures, energies,
etc.) are accumulated
for future analysis
Not all system properties
reach equilibrium at the
same time
Advanced MD Analysis in Amber
Conformational clustering tools is available in ptraj
ptraj uses several different algorithms for clustering trajectory frames into groups based on
pairwise similarity
Clustering
Trajectory snapshots after MD
Statistical Ensembles
•
The microscopic state of a system is defined by the atomic positions, q, and momenta,
p; these can also be considered as coordinates in a multidimensional space called
phase space
•
A single point in phase space, denoted by G, describes the state of the system
•
An ensemble is a collection of points in phase space satisfying the conditions of a
particular thermodynamic state.
•
A Molecular Dynamics simulations generates a sequence of points in phase space as a
function of time;
•
These points belong to the same ensemble, and they correspond to the different
conformations of the system and their respective momenta
Statistical Ensembles Supported in
Amber
• Microcanonical ensemble (NVE) : The thermodynamic state
characterized by a fixed number of atoms, N, a fixed volume, V,
and a fixed energy, E. This corresponds to an isolated system.
• Canonical Ensemble (NVT): This is a collection of all systems
whose thermodynamic state is characterized by a fixed number
of atoms, N, a fixed volume, V, and a fixed temperature, T.
• Isobaric-Isothermal Ensemble (NPT): This ensemble is
characterized by a fixed number of atoms, N, a fixed pressure, P,
and a fixed temperature, T.
Advanced MD Techniques in Amber
Adaptively Biased Molecular Dynamics (ABMD) method
ABMD is a method for the computation of the free energy surface of a reaction
coordinate using non-equilibrium dynamics.
• Chemical reactions,
conformational transitions, etc,
occur when the system migrates
from one local equilibrium
minimum to another, overcoming
the usually large energy barriers
that separate reagents from
products.
• The probability of such an event
occurring spontaneously depends
exponentially on the energy barrier
and easily exceeds the
computational time regime that
present-day computer technology
can afford.
Advanced Molecular Dynamics
Techniques in Amber:
Path integral molecular dynamics
Path integral molecular dynamics simulations can be used to sample equilibrium canonical
distributions using quantum dynamics rather than Newton's equations for nuclear
motion. Both equilibrium and kinetic isotope effects can be estimated via
thermodynamic integration over mass.
•
Centroid Molecular Dynamics (CMD) is an approximate method for calculating realtime quantum correlation functions.
•
Ring Polymer Molecular Dynamics (RPMD).
•
Both CMD and RPMD simulations provide an efficient route for the calculation of
approximate correlation functions, which can then be related to the true quantum
correlation functions.
How to Treat Bulk System?
Bulk (“infinite”) solvent
We run computer simulation to predict and
study the properties of a system in bulk
(very big or “infinite” system)
How to Treat Bulk System?
Treatable system
But…
We can simulate only a relatively small
number of particles in order not to slow
down the computation.
Artificial surface effect
Problem: we are not interested
in surface effects 
Possible Solutions?
1) The system size should
be extremely large to
ensure that the surface has
only a small influence on
the bulk properties
But…
Such system is too
big to simulate… 
Surface effect has small
influence on the bulk
properties
Possible Solutions?
2) Surface effects can be ignored for all system sizes if we use periodic boundary
conditions.
The cubical simulation box (Central
Box) is replicated throughout space to
form an infinite lattice
All other boxes are identical to the
Central Box (its copies)
Central Box
No surface
Periodic Boundary Conditions
2) Then one of its images will
enter through the opposite face
1) If a molecule leaves the
Central Box
Periodic Boundary Conditions
Periodic Boundary Conditions in
AMBER
Rectangular parallelepiped
1)
Truncated octahedron
2)
Truncated octahedron has the advantage of being more
nearly spherical than most other MD cells. This can be
very useful when simulating a large molecule in
solution, where fewer solvent molecules are required for
a given simulation cell width.
Estimation of Binding Energies in
Non-Covalent complexes
In general, non-covalent bonding refers to attractive
intermolecular forces that are not covalent in nature.
Non-covalent interactions may include ionic bonds,
hydrophobic interactions, hydrogen bonds and van der
Waals forces.
Estimation of Binding Energies in
Non-Covalent complexes
Protein-ligand complex
Evaluating Free Energies of
Binding using Amber
Evaluating Free Energies of
Binding: MM-PBSA
• The acronym MM-PBSA stands for
Molecular Mechanics- Poisson Bolzmann
Surface Area
• The MM-PBSA approach represents the
postprocessing method to evaluate free
energies of binding or to calculate absolute
free energies of molecules in solution.
Acc. Chem. Res. 2000, 33, 889-897
Evaluating Free Energies of
Binding: MM-PBSA
1. One carries out a molecular dynamics simulation,
typically in a periodic box with water and counterions
(“regular” MD simulation), and correct
representation of long-range electrostatic effects such
as PME, saving a set of representative structures.
2. After MD Simulation any solvent and counterion
molecules are removed, and the free energy, G, is
calculated according to the following equation:
Evaluating Free Energies of
Binding: MM-PBSA
any solvent and
counterion molecules
are removed
Evaluating Free Energies of
Binding: MM_PBSA
where G is the calculated average free energy, and EMM is the average molecular
mechanical energy:
where these correspond to the bond, angle, torsion, van der Waals, and
electrostatic terms in the molecular mechanical force field, evaluated with no
nonbonded cutoff.
Evaluating Free Energies of
Binding: MM_PBSA
GPBSA is the solvation free energy calculated with a numerical solution of the
Poisson-Bolzmann equation and an estimate of the nonpolar free energy with a
simple surface area term.
-TSMM is the solute entropy, which can be estimated by quasi harmonic analysis
of the trajectory or, in selected cases, by using normal-mode analysis. This final
term is likely to be much smaller than the other two in many applications of
estimating relative free energies.
Evaluating Free Energies of
Binding: Thermodynamic integration
The thermodynamic integration (TI) technique allows to calculate the free
energy difference between two systems, A and B, by slowly interconverting the
Hamiltonian HA (representing system A) into the Hamiltonian HB (representing
system B), during the course of the simulation.
This process could involve the annihilation or creation of atoms
(“Computational alchemy” ).
Examples:
Atom → nothing
Group of Atoms (or Molecule) → nothing
Charge on Atom → No charge on Atom
Charge on Group of Atoms (or Molecule) →
→ No charge on Group of Atoms (Molecule)
“Computational alchemy”
• One common application of this model is pKa calculations,
where the charges are mutated from the protonated to the
deprotonated form
O
C
O
H
O
C
O
Disappears during simulation
Evaluating Free Energies of
Binding: Thermodynamic integration
The free energy difference is then given by

H
 d

A

0   
1
The subscript λ at the pointed angles indicates that the average should be
taken over an ensemble with Hamiltonian Hλ .
In MD simulations the integral is often replaced by a sum over a discrete set of
values of λ:

H
 

A




 i
i
where Δλ is chosen such that the result is statistically accurate while using a
minimum of computer time.
Inclusion of Solvation Effects in
Amber
Practically all important biological
processes take place in solvent
Solvation methods can be devided into two
main categories: explicit
(supermolecule) and implicit solvation
methods
Explicit solvent model
Molecular solvent models
employ hundreds or
thousands of discrete
solvent molecules
Pros: Many of the
properties of solutions
and solutes can be
reproduced
Cons: Such calculations converge only slowly to precise answers
because of the large number of particles and states involved;
expensive computationally.
Implicit Solvation Methods in Amber
Implicit solvation schemes speed up the
calculations by orders of magnitude and
are assumed to compromise little on
essential features of the solvation
phenomenon.
Continuum solvation models
Continuum model treat the
solvent as a continuous
medium having the average
properties of the real
solvent and surrounding the
solute beginning at or near
its van der Waals or
Solvent-accessible surface.
Pros: Faster than
molecular solvation
models
Cons: Obtaining accurate numerical solutions for a large
system such as a protein still has a significant
computational cost
Implicit Solvation Methods in Amber:
The Generalized Born/Surface Area Model
To estimate the total solvation free energy of a molecule, ΔGsolv ,
one typically assumes that it can be decomposed into the
"electrostatic" and "non-electrostatic“ parts:

G


G


G
solv
el
nonel
where ΔGnonel is the free energy of solvating a molecule from which all charges
have been removed (i.e. partial charges of every atom are set to zero), and
ΔGel is the free energy of first removing all charges in the vacuum, and then
adding them back in the presence of a continuum solvent environment.
ΔGnonel comes from the combined effect of two types of interaction: the
favorable van der Waals attraction between the solute and solvent molecules,
and the unfavorable cost of breaking the structure of the solvent (water)
around the solute.
Implicit Solvation Methods in Amber:
The Generalized Born/Surface Area
Model
Calculating ΔGnonel:
In the Amber code ΔGnonel is taken to be
proportional to the total solvent accessible
surface area (SASA) of the molecule, with a
proportionality constant derived from
experimental solvation energies of small
non-polar molecules, and uses a fast Linear
Combinations of Pairwise Overlaps (LCPO)
algorithm [J. Comput. Chem. 20, 217-230
(1999)] to compute an analytical
approximation to the surface accessible area
of the molecule.

G



Area

b
nonel
Implicit Solvation Methods in Amber:
The Generalized Born/Surface Area
Model
Calculating ΔGel:
Within Amber GB models, each atom in a molecule is represented as a sphere of
radius ρi with a charge qi at its center; the interior of the atom is assumed to be
filled uniformly with a material of dielectric constant of 1. The molecule is
surrounded by a solvent of a high dielectric εw (80 for water at 300 K)
q
q
1
i
j

G


G


el
gb 


2
f
r
,
R
,
R
ij
gb
ij
i
j
where rij is the distance between atoms i and j, the Ri are the so-called effective
Born radii of atoms i and j, and fgb is a certain smooth function of its arguments.
A common choice of fgb is
1
/
2





r
2


f

r

R
R


gb
ij
i
jexp
 4

R
R
i
j




2
ij
Implicit Solvation Methods in Amber:
ALPB (Analytical Linearized PoissonBoltzmann)
Based on an approximate analytical solution of the linearized Poisson-Bolzmann
equation for a sphere (Kirkwood, 1934).
The basic ALPB equation that approximates the electrostatic part of the solvation
free energy is:




1
1
1
1 
1





G


G



q
q


el
alpb
i
j




2
1

f
A
ij 
in
ex
gb



where β = εin /εex is the ratio of the internal and external dielectrics, α = 0. 571 412,
and A is the so-called effective electrostatic size of the molecule. fgb is the same
smooth function as in the GB model.
The GB approximation is then just the special case of ALPB when the solvent
dielectric is infinite; however, for finite values of solvent dielectric the ALPB tends to
be more accurate.
Grigori Sigalov, Andrew Fenley, and Alexey Onufriev, J. Chem. Phys. 124, 124902 (2006)
Grigori Sigalov, Peter Scheffel, and Alexey Onufriev, J. Chem. Phys. 122, 094511 (2005)
Implicit Solvation Methods in Amber:
Poisson-Boltzmann solver
An efficient finite-difference numerical solver is implemented for various
applications of the Poisson-Boltzmann (PB) method.
The electrostatic potential φj at atomic charge site is computed by solving the PB
equation:
 












r

r


4
r

4
z
exp

z
r
/
k
Y

i
i
B
i
where ε(r) is the dielectric constant, φ(r) is the electrostatic potential, ρ(r) is the solute
charge, zi is the charge of ion type i, ci is the number density of ion type i far from the
solute, kB is the Boltzmann constant, and T is temperature; the summation is over all
different ion types.
• This is the most rigorous method for treatment of implicit solvent in Amber
• It can be used for both static (single point) and dynamic applications.
• However, it is much slower than GB and ALPB and memory intensive for
macromolecules.
Inclusion of Solvation Effects in
Amber: RISM
RISM - Reference Interaction Site Model
RISM is an approximate solution to the OrnsteinZernike (OZ) equation:
where r12 is the separation between particles 1 and 2 while Ω1 and Ω2 are
their orientations relative to the vector r12. The two functions in this
relation are h, the total correlation function, and c, the direct correlation
function.
RISM: Practical Considerations
• Calculating a 3D-RISM solution for a single
solute conformation typically requires about
100 times more computer time than the same
calculation with explicit solvent or PB.
• Memory: anywhere from a few megabytes
for the smallest solutes to gigabytes for large
complexes
Exploring Conformational Space
of Biomolecules
Conformational Space
of Biomolecules Can Be Very Complex
Exploring Conformational Space
of Biomolecules
• Due to this property of the free energy landscape,
efficient computational approaches for searching
for low-energy minima in these complex systems
present a great challenge.
Exploring Conformational Space:
Simulating Annealing
Temperature
Time
Energy Profile
Local Minima
Exploring Conformational Space:
REMD
REMD stands for the Replica Exchange Method Dynamics
In REMD several noninteracting copies (replicas) are independently and
simultaneously simulated at different temperatures.
Replica 1, T1
Replica 2, T2
Replica N, TN
At intervals during the otherwise standard simulations, conformations of the
system being sampled at different temperatures are exchanged based on a
Metropolis-type criterion
Exploring Conformational Space:
REMD
Replica 1, T1
Replica 2, T2
Replica N, TN
As a result, the low temperature simulations (replicas) have the potential to
escape kinetic traps by jumping to minima that are being sampled by the
higher-temperature replicas where kinetic trapping is less prevalent.
Treating Long-Range
Electrostatic Interactions
Eelectrostatic 
atoms atoms

i
j i
qi q j
rij
Treating Long-Range
Electrostatic Interactions
Cut-off
Treating Long-Range
Electrostatic Interactions
Treating Long-Range
Electrostatic Interactions in
Amber
• The particle-mesh Ewald (PME) procedure
(or, optionally, a "true" Ewald sum) is used
to handle long-range electrostatic
interactions.
Doing Semi-Empirical Quantum
Chemistry with Amber
• Amber 12 is packaged with sqm - a linear
scaling semi-empirical program for
calculation of energies, charges and
geometries of systems up to ˜20,000 atoms.
Doing Semi-Empirical Quantum
Chemistry with Amber
sqm’s Available features include:
• Linear scaling Divide and Conquer (D&C) calculations.
• Single point AM1, PM3, MNDO, MNDO/d or PDDG-PM3
calculations.
• Geometry Optimization (steepest decent, conjugate
gradient, BFGS, and LBFGS available)
• Mulliken, CM1 and CM2 charge analysis
• Nuclear Magnetic Resonance prediction and simulation
• Mixed quantum mechanics/molecular mechanics
(QM/MM) linear scaling Semi-Empirical calculations.
Doing Semi-Empirical Quantum
Chemistry with Amber
• Amber 12 is packaged with SQM semiempirical program.
MNDO: H, Li, Be, B, C, N, O, F, Al, Si, P, S, Cl, Zn, Ge, Br, Sn, I,
Hg, Pb
AM1: H, C, N, O, F, Al, Si, P, S, Cl, Zn, Ge, Br, I, Hg
PM3: H, Be, C, N, O, F, Mg, Al, Si, P, S, Cl, Zn, Ga, Ge, As, Se, Br,
Cd, In, Sn, Sb,
PDDG/PM3: H, C, N, O, F, Si, P, S, Cl, Br, I
PDDG/MNDO: H, C, N, O, F, Cl, Br, I
RM1: H, C, N, O, P, S, F, Cl, Br, I
PM3CARB1: H, C, O
PM6: H, He, Li, Be, B, C, N, O, F, Ne, Na, Mg, Ar, K, Ca, Zn, Ga,
Ge,Kr, Rb,
DFTB/SCC-DFTB: (Any atom set available from the
www.dftb.org website)
A Hybrid Quantum
Mechanical/Molecular Mechanical
(QM/MM) Approach
Why Do We Need a Hybrid QM/MM
Approach?
Quantum Mechanics
Molecular Mechanics
generally applicable
restricted to the classes of molecule
it have been designed for
allow the calculation of ground and
excited state properties: molecular
energies and structures, energies and
structures of transition states, atomic
charges, reaction pathways etc.
allow the calculation of ground state
properties: relative molecular energies
and structures
CPU and memory hungry.
Computationally efficient
Suitable for small and medium size
systems
Suitable for large molecular
systems
Why Do We Need a Hybrid QM/MM
Approach?
The main bottleneck of quantum chemical methods is that they are CPU and memory
hungry.
For example, for small peptide of 126 atoms one energy evaluation requires:
CPU Time
Method
Memory
Seconds
Time units
KB
Memory units
Quantum
chemical*
273.00
1820
4889
85
Molecular
Mechanical
0.15
1
58
1
*Semi-empirical PM3 method
In general, CPU and memory requirements (N – number of atoms):
Molecular Mechanical methods
~ N2
Semiempirical Quantum Chemical methods
~ N2
Ab initio Quantum Chemical methods
~ N4
A Hybrid QM/MM Approach
The general idea of a hybrid QM/MM approach is that large chemical systems may be
partitioned into 1) an electronically important region (QM region) which requires a
quantum chemical treatment and 2) a remainder which only acts in a perturbative fashion
and thus admits a classical description (MM region).
The Simplest Hybrid QM/MM Model
Hamiltonian for molecular system in the Born-Oppenheimer approximation:
1 electrons 2 electronsnuclei Z j electronselectrons 1 nuclei
H      
  
 
2 i
R
r
i
j
i
j i
i
ij
ij
nuclei

Zi Z j
j i
Rij
1 electrons 2 electronsnuclei Z j electronselectrons 1 nuclei
H      
  
 
2 i
R
r
i
j
i
j i
i
ij
ij
nuclei
Zi Z j
j i
Rij

The MM region is viewed in the QM
calculations as a set of point charges
“Standard” QM hamiltonian

electrons ch arg es
 
i
k
Qk nuclei ch arg es Z iQk
  
Rik
Rik
i
k


Effectof External Ch arg es
The main drawbacks of this simple QM/MM model are:
 it is impossible to optimize the position of the QM part relative to the external
charges because QM nuclei will collapse on the negatively charged external charges.
 some MM atoms possess no charge and so would be invisible to the QM atoms
 the van der Waals terms on the MM atoms often provide the only difference in
the interactions of one atom type versus another, i.e. chloride and bromide ions both
have unit negative charge and only differ in their van der Waals terms.
A Hybrid QM/MM Model
So, it is quite reasonable to attribute the van der Waals parameters (as it is in the MM
method) to every QM atom and the Hamiltonian describing the interaction between the
QM and MM atoms can have a form:
Hˆ QM / MM  
electrons MM atoms
 
i
j
Qj
rij

nuclei MM atoms
 
i
j
Z iQ j
Rij
nuclei MM atoms

i

j

 Aij Bij 


 12
6
R
R

ij 
 ij

The van der Waals term models also electronic repulsion and dispersion interactions,
which do not exist between QM and MM atoms because MM atoms possess no
explicit electrons.
A. Warshel, M. Levitt // Theoretical Studies of Enzymic Reactions: Dielectric,
Electrostatic and steric stabilization of the carbonium ion in the reaction of
lysozyme. // J.Mol.Biol. 103(1976), 227-49
The Hybrid QM/MM Model
Now we can construct a “real” hybrid QM/MM Hamiltonian:
Hˆ  Hˆ QM  Hˆ QM / MM  Hˆ MM
Hˆ QM / MM  
electrons MM atoms
 
i
j
Qj
rij

nuclei MM atoms
 
i
j
Z iQ j
Rij
nuclei MM atoms

i

j

 Aij Bij 

 12  6 
Rij 

 Rij

A “standard” MM force field can be used to determine the MM energy. For example,
AMBER-like force field has a form:
ETotal 
 Aij Bij qi q j 
 Cij Dij qi q j 
2






K
(
R

R
)
  K (   0 ) 
 12






b
0
6
12
10
Rij
Rij  H bonds  Rij
Rij
Rij  bonds
nonbonded 
angles
 Rij
V
(1  cos(n ))

dihedrals 2
Choice of QM method
... is a compromise between computational efficiency and practicality and the desired
chemical accuracy.
The main advantage of semi-empirical QM methods is that their computational
efficiency is orders of magnitude greater than either the density functional or ab initio
methods
Calculation times (in time units)
F
O
O
O
Ab initio method
O
P O
N
O
O
O
N
F
F
O
1800
RHF/6-31G*
36228
1
PM3
1
Semi-empirical method
Calibration of the QM/MM potential
Crucial aspect is how the interaction between
QM and MM parts is determined.
In choosing the appropriate form, it is required
that the balance between attractive and
repulsive forces must be preserved and the
QM/MM interactions must be of the correct
magnitude with respect to the separate QM
and MM contributions
Calibration of the QM/MM potential:
Parameterizations
Hˆ QM / MM  
electrons MM atoms
 
i
j
Qj
rij

nuclei MM atoms
 
i
1)
j
Z iQ j
Rij
nuclei MM atoms

i

j

 Aij Bij 


 12
6
R
R

ij 
 ij

2)
1) Modification of the one-electron terms arising from
interaction of the electron cloud of the QM fragment
with the point charge of an MM atom.
2) By varying the radii in the van der Waals terms.
3) By varying 1)+2)
Calibration of the QM/MM potential
Hˆ QM / MM  
electrons MM atoms
 
i
j
Qj
rij

nuclei MM atoms
 
i
j
Z iQ j
Rij
nuclei MM atoms

i

j

 Aij Bij 

 12  6 
Rij 

 Rij

1) By hand, to find the optimum values of the parameters by
calculating interaction curves for charge/ion systems and
comparing them with the MP2/6-311++G** ab initio results.
M.J. Field, P.A. Bash, M. Karplus, J.Comp.Chem., 11(1990), 700-733.
2) Fitting calculated H-bond energies to experimental data
on ion-molecular complexes in the gas phase.
V.V. Vasilyev, A.A. Bliznyuk, A.A. Voityuk, Int.J.Quant.Chem. 44(1992), 897930.
Calibration of the QM/MM potential
Hˆ QM / MM  
electrons MM atoms
 
i
j
Qj
rij

nuclei MM atoms
 
i
j
Z iQ j
Rij
nuclei MM atoms

i

j

 Aij Bij 

 12  6 
Rij 

 Rij

3) Optimizing van der Waals parameters on QM atoms to reproduce the 631G(d) interaction energies for H-bonded complexes in the gas phase.
P.A. Bash, L. Lawrence, A.D. MacKerell, Jr., D. Levine, P. Hallstrom, PNAS USA, 93(1996),
3698-703.
4) Optimizing van der Waals parameters on QM atoms to reproduce the
MP2/6-31G(dp) interaction energies for H-bonded complexes in the gas
phase.
J. Gao // Toward a molecular Orbital Derived Empirical Potential for Liquid Simulations //
J.Phys.Chem. B 101(1997), 657-63
5) By varying the radii in the van der Waals terms to reproduce experimental
free energies of solvation using MD simulations.
P.L. Cummins, J.E. Gready, J.Comp.Chem., 18(1997), 1496-512.
Dividing Covalent Bonds across the
QM and MM Regions
In many simulations it is necessary
to have the QM/MM boundary cut
covalent bonds, and a number of
additional approximations have to
be made.
Dividing Covalent Bonds across the
QM and MM Regions
Using a hybrid orbital on the frontier MM atom
MM Region
O
QM Region
O
Frontier MM
Atom
Frontier QM
Atom
Frontier MM
Atom
A. Warshel, M. Levitt // Theoretical Studies of
Enzymic Reactions: Dielectric, Electrostatic and
steric stabilization of the carbonium ion in the
reaction of lysozyme. // J.Mol.Biol.
103 (1976), 227-249
V. Thery, D. Rinaldi, J.-L. Rivail, B. Maigret,
G.G. Ferenczy, J.Comp.Chem. 15 (1995), 269
Dividing Covalent Bonds across the
QM and MM Regions
Using “link” atoms
“Link” atoms are used to gracefully cap
the electron density.
This approach is used in Amber
Implementation of “link” Atom
Approach in Amber 9 & 10
The link atom is placed along the bond
vector joining the QM and MM atom
The default link atom type is hydrogen
It interacts with MM region only
electrostatically (no VDW term).
WdV interaction between QM and MM
atoms which form 1-2 and 1-3 “bonded”
pairs is not calculated.
Bond stretching, angle bending, and
torsion interactions between QM and MM
regions are calculated as those in MM if
1-2, 1-2-3, or 1-2-3-4 terms contain at
least one MM atom
Reviews on QM/MM
•H. Hu and W. Yang, Free energies of chemical reactions in solution and in
enzymes with ab initio quantum mechanics/molecular mechanics methods, Annu
Rev Phys Chem. 2008;59:573-601
•C. Bo and F. Maseras, QM/MM methods in inorganic chemistry, Dalton Trans.,
2008, 2911–2919
• H.M. Senn and W. Thiel, QM/MM studies of enzymes, Current Opinion in
Chemical Biology, 2007(11), 182-187
• R.A. Friesner and V. Guallar, Ab initio Quantum Chemical and Mixed
Quantum Mechanics/Molecular Mechanics (QM/MM) Methods for Studying
Enzymatic Catalysis, Annual Review of Physical Chemistry, 2005 (56), 389-427
• G. Monard, X. Prat-Resina, A. González-Lafont, J.M. Lluch, Determination of
enzymatic reaction pathways using QM/MM methods, Int. J Quant Chem, 2003,
93 Issue 3, Pages 229 - 244
Hints for running QM/MM calculations
Choosing the QM region
There are no good universal rules here
One might want to have as large a QM region
as possible
However, having more than 80-100 atoms in
the QM region will lead to simulations that
are very expensive.
Hints for running QM/MM calculations
Choosing the QM region
For many features of conformational analysis, a
good MM force field may be better than a semiempirical or DFTB quantum description.
Hints for running QM/MM calculations
Choosing the QM region
QM Methods in Amber 12
Available semi- empirical Hamiltonians are
MNDO, AM1, PM3, RM1, PDDG/PM3,
PDDG/MNDO, and PM3CARB1, PM3MAIS, MNDO/d, AM1/d (Mg from AM1/d
and H, O, and P from AM1/d-PhoT) and
PM6
They can be used for gas phase, generalized
Born and PME periodic simulations.
QM Methods in Amber 12
Support is also available the DFT methods:
1. The Density Functional Theory-based-tightbinding (DFTB) Hamiltonian
2. The Self-Consistent-Charge version, SCC-DFTB
In Amber 9 the DFTB/SCC-DFTB implementation does not support generalized
Born, PME or Ewald calculations,
The elements supported by QM
methods in Amber 12
MNDO: H, Li, Be, B, C, N, O, F, Al, Si, P, S, Cl, Zn,
Ge, Br, Sn, I, Hg, Pb
AM1: H, C, N, O, F, Al, Si, P, S, Cl, Zn, Ge, Br, I,
Hg
PM3: H, Be, C, N, O, F, Mg, Al, Si, P, S, Cl, Zn,
Ga, Ge, As, Se, Br, Cd, In, Sn, Sb, Te, I, Hg,
Tl, Pb, Bi
PDDG/PM3: H, C, N, O, F, Si, P, S, Cl, Br, I
PDDG/MNDO: H, C, N, O, F, Cl, Br, I
PM3CARB1: H, C, O
DFTB/SCC-DFTB: H, C, N, O, S, Zn
QM/MM calculations: ab initio
and DFT methods
Amber can support QM/MM simulations via an
interface to external QM software packages:
ADF
Gaussian
(Amsterdam Density Functional)
GAMESS-US
Orca
NWChem
TeraChem
QM/MM calculations: ab initio
and DFT methods
Mechanical and electrostatic embedding:
• Gaussian
• Orca
• TeraChem
Mechanical embedding:
• ADF
• GAMESS-US
• NWChem
Importance of Visualization
One quick look at the structure can help to detect
errors and save days or weeks of your time
Freeware Visualization Programs:
RasMol
http://www.openrasmol.org/
Freeware Visualization Programs:
VMD (Visual Molecular Dynamics)
VMD is a molecular visualization program for displaying, animating, and analyzing
http://www.ks.uiuc.edu/Research/vmd/
large biomolecular systems using 3-D graphics and built-in scripting.
Freeware Visualization Programs:
gOpenMol
http://www.csc.fi/gopenmol/
http://sirius.sdsc.edu/
Freeware Visualization Programs:
Chimera
http://www.cgl.ucsf.edu/chimera/
Freeware Visualization Programs:
MD Display
A Multi-platform 3D Stereo Molecular Dynamics Trajectory Visualization Package
http://www.cgl.ucsf.edu/chimera/
Commercial Programs
… they represent an expert molecular modeling
environment which provides construction,
editing, and visualization tools for both large
and small molecules
Tripos (www.tripos.com)
Accelrys (http://www.accelrys.com)
and others…
Learning Amber
Amber Basic Tutorials
• http://ambermd.org/tutorials/
Simulating a small
fragment of DNA
Basic introduction to LEaP, sander, and ptraj, to build,
solvate, run MD and analyze trajectories.
Using VMD with AMBER
Brief introduction to using VMD for visualising AMBER
inpcrd, restrt and trajectory files
Folding TRP Cage
Vreating structures using XLeap followed by running heating
and long MD simulations to conduct protein folding
experiments. Advanced analysis: RMSd fitting, mdcrd to
binpos conversion, average structure calculation, hydrogen
bond analysis and dihedral angle tracking using ptraj
Demo of Ptraj Commands
How to use AMBER's ptraj analysis program to analyse a
peptide simulation and gather a range of statistics from the
trajectory.
Visualizing Amber
Trajectories with Sirius
how to use Sirius visualization software to display and analyze
AMBER MD trajectory files
Amber Advanced Tutorials
• http://ambermd.org/tutorials/
Setting up an Advanced
System (Including Charge
Derivation)
Preparing a system, for simulation with sander, that contains
several non-standard residues
A simple coupled potential
QM/MM/MD simulation.
How to set up a simple QM/MM/MD simulation of NMA in
solution using AMBER 9
MM-PBSA
Step by step explanation of using the mm_pbsa script in
AMBER 9 to calculate the binding energy of the RAS-RAF
protein complex
Nudged Elastic Band
(NEB) method
How use the NEB method to predict a pathway for a
conformational change in alanine dipeptide.
pKa Calculations using
Thermodynamic
Integration
How to calculate the pKa value of the ASP residue in the
protein thioredoxin
… and other
Resume
• Amber package represents an expert molecular
modelling environment with a reach
functionality and good computer performance.
Hands-on
Web: sf.anu.edu.au/~vvv900/monash
Tutorial files (AMBER_INTRO_COURSE/):
• standard-setup.tar – “Standard” setup (long)
• nonstandard.tar – Handling non-standard
residues (long)
• amber-gaussian.tar – QM/MM using AmberGaussian interface (short)
• qm-mm.tar – QM/MM using Amber inbuilt
semiempirical methods