Transcript QM/MM - Australian National University

A coupled (hybrid) potential
QM/MM/MD simulations in Amber
Dr. Vladislav Vasilyev
Supercomputer Facility
The Australian National University
Hybrids are popular from the
ancient time
Why Do We Need a Hybrid QM/MM
Approach?
Quantum chemical methods are generally applicable and
allow the calculation of ground and excited state properties
(molecular energies and structures, energies and structures
of transition states, atomic charges, reaction pathways etc.)
Molecular Mechanical methods are restricted to the
classes of molecule it have been designed for and their
success strongly depends on the careful calibration of a large
number of parameters.
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:
Molecular Mechanical methods
~ N2
Semiempirical Quantum Chemical methods
~ N2
Ab initio Quantum Chemical methods
~ N4
A Hybrid QM/MM Approach
The development of hybrid QM/MM approaches is guided by the general idea that
large chemical systems may be partitioned into an electronically important region
which requires a quantum chemical treatment and a remainder which only acts in a
perturbative fashion and thus admits a classical description.
The Simplest Hybrid QM/MM Model
Hamiltonian for the 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

“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
O
P O
N
O
O
O
N
O
1800
RHF/6-31G*
36228
1
PM3
1
F
F
QM Methods in Amber 9
• Available semi- empirical Hamiltonians
are PM3, AM1, MNDO, PDDG/PM3,
PDDG/MNDO.
• They can be used for gas phase,
generalized Born and PME periodic
simulations.
QM Methods in Amber 9
•
Support is also available, on a functionally
limited basis:
1. The Density Functional Theory-based-tightbinding (DFTB) Hamiltonian
2. The Self-Consistent-Charge version, SCCDFTB
The DFTB/SCC-DFTB implementation does not currently support generalized
Born, PME or Ewald calculations,
The elements supported by QM
methods in Amber 9
•
•
•
•
•
•
•
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
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
6-31G(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
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
Example of Application of the
QM/MM Method to the Enzyme
Catalysis
Tetrahedral Intermediate Formation in the
Acylation Step of Acetylcholinesterases
• Acetylcholinesterases (AChE) are the serine protease
enzymes which hydrolyze the neurotransmitter
acetylcholine (Ach)
CH3COOCH2CH2N+(CH3)3 + AchE  CH3CO-AchE +
HOCH2CH2N+(CH3)3  CH3COO- + H+ + AChE
and which function at a rate approaching that of a diffusioncontrolled reaction.
• Remark: Human Cathepsin G is also a serine protease
enzyme
V.V. Vasilyev, J.Mol.Struct. (Theochem), 304(1994), 129.
Steps along the reaction pathway of serine
protease catalyzed bond cleavage
Asp-102
His-57
C
 “Catalytic triad”
Ser-195
O
R H
H
O
O
H
R3 C N C C N
O
N
N
H
C
-
R1
R2 O
Asp-102
His-57
C
Ser-195
H
O
R
N
+
N
C
H
O
O
H
H
C
R3 C N C
O
R2
N
R1
O Tetrahedral Intermediate
Asp-102
His-57
C
Ser-195
N
R H O
R3 C N C C
O
R2 O
Acyl-Enzyme
O
H
H
N
C
H
O
O
 Rate limiting step of the
reaction
Partition into the QM and MM Parts
Asp-102
His-57
C
Ser-195
O
R H
H
O
O
H
R3 C N C C N
O
N
N
H
C
-
R1
R2 O
Asp-102
 QM Region (PM3) – CH3OH
(as Ser-200), methylimidazole (as
His-440), CH3CH2COOH (as
Glu-327), and CH3COOCH3 (as
a substrate)
His-57
C
Ser-195
H
O
R
N
+
N
C
H
O
O
H
H
C
R3 C N C
O
R2
N
R1
O Tetrahedral Intermediate
Asp-102
His-57
C
Ser-195
N
R H O
R3 C N C C
O
R2 O
Acyl-Enzyme
O
H
H
N
C
H
O
O
 MM Region – 5161 enzyme
atoms
Activation of the Serine Residue
O
H
N
N
O
-
H
N
+
Critical Step of the Reaction is activation of the Serine residue
In Gas Phase:
PA(CH3O-) = ~ -350 kcal/mole
PA(Imidazole) = ~ -220 kcal/mole
where PA is a Proton Affinity
N
Activation of the Serine Residue
Energetics of the proton transfer
from Ser to His in the enzyme and in the gas phase
The energy of the first point, Ser(0)-His(0)-Glu(-), is taken as zero
A proton transfer from the Serine
to Histidine residue is very
unfavorable in the gas phase
(~34.6 kcal/mol).
40
34.6
Relative Energy, kcal/mol
30
In Gas Phase
21.5
In Enzyme
20
17.2
10
0
0
Ser(0)
Ser(-)
A proton transfer is less
unfavorable in the enzyme
(energy barrier is about 21.5
kcal/mol).
Activation of the Serine Residue
Energetics of the proton transfer
from Ser to His in the enzyme and in the gas phase
The energy of the first point, Ser(0)-His(0)-Glu(-), is taken as zero
40
O
34.6
N
H
-0.8
Relative Energy, kcal/mol
30
In Gas Phase
21.5
In Enzyme
20
-28.4
17.2
10
-54.3
Electrostatic
Potential (kcal/mol)
0
0
Ser(0)
Ser(-)
Enzyme environment creates electrostatic potential which favours the proton
transfer from the Serine to Histidine residue.
N
Tetrahedral Intermediate Formation in the
Acylation Step of Acetylcholinesterases
 Incorporation of enzyme
environment in simulation
changes drastically the flow of
the reaction
40
Relative energy, kcal/mol
20
 A decrease in the activation
energy of TI formation in the
enzyme environments versus
the uncatalysed gas-phase
reaction is about 27 kcal/mol
Reaction path
0
In Enzyme
TI
In Gas Phase
-20
-40
1.00
2.00
3.00
Serine-Substrate Distance, Angstrom
4.00
(experimental estimation of the
reduction in activation energy for
the whole enzyme reaction
versus uncatalysed neutral,
basic, and acidic hydrolysis of
AchE is 14, 11, and 18 kcal/mol,
respectively.)
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
semi-empirical or DFTB quantum
description.
Hints for running QM/MM calculations
Choosing the QM region
Hints for running QM/MM calculations
Parallel Simulations
• At present all parts of the QM simulation are parallel
except the density matrix build and the matrix
diagonalisation.
• For small QM systems these two operations do not take a
large percentage of time and so acceptable scaling can be
seen to around 8 cpus.
• However, for large QM systems the matrix diagonalization
time will dominate and so the scaling will not be as good.
Amber 9 QM/MM Example
Resume
• Amber 9 features new and significantly improved
QM/MM support
• The QM/MM facility supports gas phase, implicit
solvent (GB) and periodic boundary (PME)
simulations
• Compared to earlier versions, the QM/MM
implementation offers improved accuracy, energy
conservation, and performance.