The Force Field Toolkit ( fg TK) Christopher G. Mayne, Emad - - PowerPoint PPT Presentation
Parametrizing Small Molecules Using: The Force Field Toolkit ( fg TK) Christopher G. Mayne, Emad Tajkhorshid Beckman Institute for Advanced Science and Technology University of Illinois, Urbana-Champaign Klaus Schulten University of Illinois,
Beckman Institute for Advanced Science and Technology University of Illinois, Urbana-Champaign
Georgia Institute of Technology
University of Illinois, Urbana-Champaign
Computational Biophysics Workshop | DICP | July 11 2018
bonds bond
angles ki
angle
dihedralski
dihedral
i Σ j≠i 4∈ij （
12
6
i Σ j≠i
bonds bond
angles ki
angle
dihedralski
dihedral
i Σ j≠i 4∈ij （
12
6
i Σ j≠i
O O P O O O O O O P O O O O O N N N N O NH2 N N NH O NH2 HN O O Me P O O O O N N N N H2N
Peptides and Proteins Nucleic Acids repetitive backbone unit limited set of isolated building blocks
H N N H H N N H H N N H Me O O O O O O NH2 O Me Me OH O O S Me
N O HN N HN Me Me N N N N O Me Me O O OH S S
Imatinib (Gleevec) Tiotropium (Spiriva)
H2N COOH O
N H COOH HO
Coenzyme A
O N O OH P O O O N N N NH2 O P O O O P O O O N H N H Me Me OH O O HS
H2N O S NH2 O S S H2N O S NH2 O Fe
Water Interaction En. (QM) Charge Optimization Hessian Calculation (QM) Bond & Angle Optimization Torsion Scan (QM) Torsion Optimization Find Missing Parameters Geometry Optimization (QM) update PSF update PAR update PAR build init PAR update PDB
C.G. Mayne et al., J. Comput. Chem. 2013, 34, 2757-2270.
update PSF update PAR update PAR build init PAR update PDB
Hessian Calculation (QM) Bond & Angle Optimization Torsion Scan (QM) Torsion Optimization Find Missing Parameters Geometry Optimization (QM)
Water Interaction En. (QM) Charge Optimization
Load QM optimized geometry | Auto-detect interaction sites | Generat
VMD main window fgTK GUI
mized geometry | Auto-detect interaction sites | Generate Gaussian Input Files | Run Q
VMD main window fgTK GUI
raction sites | Generate Gaussian Input Files | Run QM | Inspect water optimization
fgTK GUI
iles | Run QM | Inspect water optimization
fgTK GUI
Load QM Target Data Prepare Optimization Optimizer: Assign Charges Compute UMM, dMM, µMM Compute Objective Function Return Optimized Charges Analyze Performance Write Charges to PSF Setup Optimization
0.2 0.4
dMM-QM (Å)
2 4 6 8 10
UMM-QM (kcal/mol)
Initial Charges Intermediate Charges Final Optimized Charges Literature Charges
0.2 0.4
1 2 3 4 5 6
0.2 0.4 -0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4 500 1000 1500 2000 2500 3000 3500 4000 4500
N H N H N H H H N H H H
2
2
2
2
In practice, it is impossible to fit all of these perfectly! Often we decrease wd and wµ to improve the fit to the energies
50 100 150 200 250 300
Optimization Iteration
1 2 3 4 5 6 7 8
ΔE from Target (kcal/mol)
export data
https://studynights.blogspot.com/2015/03/the- single-point-energy-of-mnh2o6-and.html
An alternative to water interactions for charges, commonly used in Amber The QM electrostatic potential is calculated and then fit by optimizing the MM charges Has problems with buried atoms, which may not noticeably affect the ESP RESP fitting is supported by FFTK, requires downloading the resp program as part of AmberTools (free)
Bond and angle are fit by creating a small distortion of the bond/angle and calculating the QM energy and the MM energy, then choosing the force constants to match
bonds bond
angles ki
angle
periodicity (1-6 possible) phase (always 0 or 180 deg.)
Dihedrals are scanned in QM in 10-15 deg. increments Energies of each conformation are fit in MM according to the pre-determined dihedral terms included
dihedralski
dihedral
periodicity phase (always 0 or 180 deg.)
Energies above a threshold (e.g., 8-10 kcal/mol) are ignored
QM PES: black initial: red First fit: blue After 7 rounds of simulated annealing, the fit is much better
One multiplicity and locked phase Several multiplicities and free phase Pro: very good fit of QM PES Cons: possible incorrect behavior multiple sets of force constants Pros: limits incorrect behavior, sets of force constants Cons: fit to QM PES not always possible
imidazole-pyrridine moiety of antibiotic telithromycin
50 100 150 200
Configuration 1 2 3 4 5 6 7 8 9 10
E (kcal/mol)
QM PES MM Fit
Dihedral Fit
3 multiplicities, free phases for each dihedrals
k1[1 + cos(φ + δ1)] + k2[1 + cos(2φ + δ2)] + k3[1 + cos(3φ + δ3)]
multiple dihedrals scanned; plotted simultaneously to reference a common global minimum
Pavlova, Gumbart. Parametrization of macrolide antibiotics using the force field toolkit. J. Comp. Chem. 2015, 36, 2052−2063.
imidazole-pyrridine moiety of antibiotic telithromycin
Pavlova, Gumbart. Parametrization of macrolide antibiotics using the force field toolkit. J. Comp. Chem. 2015, 36, 2052−2063.
imidazole-pyrridine moiety of antibiotic telithromycin
planar dihedrals have multiplicity 2 and phase 180 deg.
50 100 150 200
Configuration 1 2 3 4 5 6 7 8 9 10
E (kcal/mol)
QM PES MM Fit
Dihedral Fitting
The fit looks (surprisingly) better despite using fewer terms (why?)
k2[1 + cos(2φ + π)]
Fitting a lot of parameters simultaneously cannot always find the best fit!
Pavlova, Gumbart. Parametrization of macrolide antibiotics using the force field toolkit. J. Comp. Chem. 2015, 36, 2052−2063.
imidazole-pyrridine moiety of antibiotic telithromycin
Restraining phase of CH dihedrals to 0 prevents eclipsed conformations
Pavlova, Gumbart. Parametrization of macrolide antibiotics using the force field toolkit. J. Comp. Chem. 2015, 36, 2052−2063.
planar dihedrals have multiplicity 2 and phase 180 deg.
k2[1 + cos(2φ + π)]
cobalamin is also known as vitamin B12, is a large, cobalt-containing compound; inability to absorb vitamin B12 causes pernicious anemia its large size and metal center make it particularly challenging for simulation
Pavlova, Parks, Gumbart. Development of CHARMM-Compatible Force-Field Parameters for Cobalamin and Related Cofactors from Quantum Mechanical Calculations. J. Chem. Theory Comput. 2018, 14, 784−798.
N N N N NH O Co O N N P O O O R HO OH O NH2 O NH2 O NH2 O NH2 O NH2 O H2N H
α β
A
R = 5’-deoxyadenosyl, Me, CN, OH
C1R C4 C5 N21 C6 N22 C3 C2 C1 C19 N24 C7 C8 C9 C10 C11 N23 C18 C17 C16 C15 C14 C12 C13 C55 C57 N59 O58 Co1 C1P C2P O3 C5B C4B C9B C8B C7B C6B N3B C2B N1B C6M C5M P O5 O2 O4 R C54 C60 C26 C37 C46 C48 C47 C41 C36 C30 C25 C20 C35 C38 C27 C61 C42 C43 C49 C50 C31 C32 C3R C4R O6R C2R C5R O8R O7R C53 C56 O33 N34 O39 N40 O44 N45 O51 N52 O62 N63 O28 N29
B
C3P
The first (tedious) step is to assign unique names to all the atoms! (and make a nice ASCII schematic if you are so inclined!)
Pavlova, Parks, Gumbart. JCTC. 2018, 14, 784−798.
A
N N N N Co N N R O NH2 O NH2 O NH2 O NH2 H
R = 5’-deoxyribose, Me, CN, OH
A B C D
Often the molecule of interest is too large (> 50 atoms) and/or flexible for direct application of QM optimization We create one or more molecule fragments for independent parametrization, combining them all at the end (possibly needing to create fragments for linker regions) Black is the corrin ring; w/blue and red side chains were separately used as well
Pavlova, Parks, Gumbart. JCTC. 2018, 14, 784−798.
He Co
Optimizing van der Waals parameters is especially challenging, but rarely necessary as existing atom types are almost always appropriate
Used a Helium probe (no charge!) to fit interaction energies in QM and MM
*Yin, D.; MacKerell, A. D., Jr. Combined ab initio/empirical approach for optimization of Lennard-Jones parameters. J. Comput. Chem. 1998, 19, 334−348. *Chen, I. J.; Yin, D.; MacKerell, A. D., Jr. Combined ab initio/empirical approach for optimization of Lennard-Jones parameters for polar-neutral compounds.
1 2 3 4
E (kcal/mol)
Cbl(I)
BP86 PBE B3LYP
QM EDiff
1 2 3 4
MM Start QM MM Fit
MM Fit ε = 0.133 rmin=4.0
A
B C 2 3 4 5
2 3 4 5
He-Co Distance (Å)
2 3 4 5 6-1
Pavlova, Parks, Gumbart. JCTC. 2018, 14, 784−798.
tried three different DFT methods (typically use MP2 or HF for CHARMM)
CHARMM focuses on interaction with waters - but where to place them for buried atoms? RESP approach was tried, but gave unphysical charges due to problems with buried atoms as well Instead, a hybrid approach was used: Natural Population Analysis (NPA) for buried Co and N atoms; RESP for all others Cobalt atom Recommended approaches can and will fail! Don’t be afraid to experiment!
Pavlova, Parks, Gumbart. JCTC. 2018, 14, 784−798.
A
N N N N Co N N R O NH2 O NH2 O NH2 O NH2 H
R = 5’-deoxyribose, Me, CN, OH
A B C D
Only n = 2 and n = 4 dihedrals for the imidazole were fit
30 60 90 120
Configuration 2 4 6 8 10 E(kcal/mol)
QM MM No Fit MM Fit
Ring Dihedrals
10 20 30 40
Configuration 2 4 6 8 10
Imidazole Dihedrals
Many (but not all!) dihedrals within the corrin ring set to 0 (avoid overfitting)
Pavlova, Parks, Gumbart. JCTC. 2018, 14, 784−798.
QM MM (no dih.) MM fit
fit to low energies is more important than fit to high (> 8 kcal/mol) energies
Multiple rounds of fitting and tested were required!
Minimization in NAMD using final parameters produced excellent agreement with QM minimized geometry (BP86/Def2-SVP)
A B C
QM MM
MeCbl(III) Cbl(II) Cbl(I)
2 4 6 8 Time (ns)
1 2 3 4 5 6
Distance 2 (Å)
2 4 6 8 Time (ns) 1 2 3 4 5 6
Pavlova, Parks, Gumbart. JCTC. 2018, 14, 784−798.
Also ran simulations of Cbl bound to proteins, monitoring various interactions over time (each run twice)
Water Interaction En. (QM) Charge Optimization Hessian Calculation (QM) Bond & Angle Optimization Torsion Scan (QM) Torsion Optimization Find Missing Parameters Geometry Optimization (QM)
fgTK is available as a VMD Plugin (1.9.1 or newer) Mayne et al.; J. Comp. Chem. 2013, 34, pp. 2757-2770 http://www.ks.uiuc.edu/Research/vmd/plugins/fftk
Volume 34 | Issues 31–32 | 2013 Included in this print edition: Issue 31 (December 5, 2013) Issue 32 (December 15, 2013)Research in Systems Neuroscience
www.c-chem.org
Journal of