Transcript Computational Quantum chemistry
Ch121a Atomic Level Simulations of Materials and Molecules
Room BI 115 Hours: Monday, Wednesday, Friday 2-3pm
Lecture 4, April 9, 2014 FF2: standard FF
Presented by Caitlin Scott
William A. Goddard III, [email protected]
Charles and Mary Ferkel Professor of Chemistry, Materials Science, and Applied Physics, California Institute of Technology TA’s Caitlin Scott and Andrea Kirkpatrick 1 L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved
Homework and Research Project
First 5 weeks: The homework each week uses generally available computer software implementing the basic methods on applications aimed at exposing the students to understanding how to use atomistic simulations to solve problems. Each calculation requires making decisions on the specific approaches and parameters relevant and how to analyze the results. Midterm: each student submits proposal for a project using the methods of Ch121a to solve a research problem that can be completed in the final 5 weeks.
The homework for the last 5 weeks is to turn in a one page report on progress with the project The final is a research report describing the calculations and conclusions L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 2
Dreiding
Dreiding is a generic or rule-based force field. Parameters were based on general principles, not fitted to specific molecules
Bond distance
= R A + R B -0.01A where R A are bond radii based on the A-CH 3 bond distance.
Bond angles
based on hydride: HOH=104.5
°, HSH=92.2°, H--C_3--H = 109.5
°, H--C_2--H = 120°, H--C_1--H = 180°, Steve Mayo, Barry Olafson, WAG, “DREIDING - A Generic Force-field for Molecular Simulations,” J Phys Chem 94 8897 (1990) L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 3
Dreiding
Dreiding is a generic or rule-based force field. Parameters were based on general principles, not fitted to specific molecules Bond distance = R A + R B -0.01A where R A are bond radii based on the A-CH 3 bond distance.
Bond angles based on hydride: HOH=104.5
°, HSH=92.2°, H--C_3--H = 109.5
°, H--C_2--H = 120°, H--C_1--H = 180°, Force constants” K bond = BO*700 kcal/mol/Å 2 where BO bond order (1,2,3) K angle = 100 kcal/mol/rad 2 .
Inversion barrier for planar molecules = 40 kcal/mol/rad 2 torsion barrier = 2.0 kcal/mol for single bonds 45 kcal/mol for double bonds Steve Mayo, Barry Olafson, WAG, “DREIDING - A Generic Force-field for Molecular Simulations,” J Phys Chem 94 8897 (1990) L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 4
Dreiding Atom Types
C_1, C_2, C_3 indicate sp,sp 2 , and sp 3 hybridized carbon atoms That is C_3 is a tetrahedral C with 4 bonds C_2 is a trigonal C with 3 bonded atoms (one is a double bond) C_1 has two bonds, one of which is a triple bond C_R is sp 2 but in an aromatic ring These may have different FF parameters but vdw and Q depend only on the element • Atom Type rules make possible the correct assignment of force field parameters throughout the molecule • The rules are easy for “chemists” to understand and easy to code • Also they allow automatic autotyping L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 5
Example Atom Types
tetrapeptide Lys-Met-Phe-Pro L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 6
Geometric Valence parameters for Dreiding
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 7
The van der Waals Parameters for Dreiding
R 0 is total bond distance, vdw radius is ½ this size R 0 =3.195 for H was a bad choice. Should have been 2.6-2.7. it is compensated by a small D 0 Parameters in the same row are similar. General trend to larger R e and D e as go down the © copyright 2013 William A. Goddard III, all rights reserved 8
a Two sp 3 c two sp 2
Dreiding parameters for dihedrals
E( φ) = ½ B[1+(-1) CM atoms: B=2, P=3, 0 ° is max b One sp 3 -one sp 2 : B=1, P=6, 0 ° is min (BO=2): B=45, P=2, 0 ° is min d two sp 2 R (BO=1.5): B=25, P=2, 0 ° is min cos(P φ)] B=barrier (kcal/mol) P = periodicity CM=+1 if 0 ° is min = 0 if 0 ° is max e two sp 2 (BO=1): B=5, P=2, 0 ° is min f two sp 2 R (BO=1.0): B=10, P=2, 0 ° is min g one or two sp 1 , monovalent (F..), metals (Fe..) all have B=0 h sp 3 -sp 3 in O column: B=2, P=2, 0 ° is max i sp 3 O column-sp 3 other column: B=2 P=2, 0 ° is min j exception to b: if sp 2 is max (eg propene) is bonded to one H: then B=2, P=3, 0 ° to H k barriers should decrease by factor of 3 as go down column, but this trend was ignored in Dreiding L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 9
Validation: experimental barriers (kcal/mol)
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 10
Compare conformations to experiment
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 11
Dreiding validate using 1
st
76 molecules in Cambridge Crystallographic Data Base (in 1988)
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 12
Dreiding accuracy over 76 molecules
Note that the test minimized the molecule in a vacuum. L04-Ch121a-Goddard 13
Complex case Linear CF
3
-(CF
2
)
n
-CF
3 C6F14 Dihedral QM energy (B3LYP, 631G*)
3,5 3 2,5 2 1,5 1 0,5 0 145 -0,5 150 155 160 165 170 175 180 185
Dihedral angle
Get F —F steric interactions for F on C1 and C3 Not all trans as for CH 3 -(CH 2 ) n -CH 3 Get 163 º torsion not 180º L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 150 º L K 164 º 180 º
Energy Components Molecular origin of helicity in Teflon
3 2,5 2 1,5 1 0,5 0 150 -0,5 -1 155 160 165 170 f torsion vdw electrostatic Force Field QM B3LYP O
(degrees)
FF electrostatic 175 180 vdW QM L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 15
High Pressure forms of C
20
F
42
The figure shows predicted stable helical conformations for C 20 F 42 . From left to right t+, t-, g+, g-, h+, and h enantiomeric pair conformations. The atoms are colored to facilitate the viewing of their helical nature. The tighter the dihedral angle (from 164 to 60) the shorter the molecule gets. Fluorine atoms of each color would be located on the same side if the molecule were prepared in the all-trans conformation. L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 16
DREIDING
E
E Internal
E VDW
E HB
E Coulomb
DREIDING – Internal Energies • Bond, Angle, Torsion, Inversion – Electrostatics • Point charges • Dielectric=1.0
– VDW
E VDW
( – Hydrogen Bonds
R
)
D VDW
[(
R VDW
/
E Coulomb R
) 12 • Lennard-Jones 12-6, exp-6, Morse… 2 (
R c
0
VDW
• 3-body term, Lennard-Jones 12-10 with a cos 4 /
q
1
r
12
R q
) 6 2 ] term
E hb
D hb
[ 5 (
R hb
/
R DA
) 12 6 (
R hb
/
R DA
) 10 ] cos 4 (
DHA
) © copyright 2013 William A. Goddard III, all rights reserved L04-Ch121a-Goddard 17
Hydrogen Bond terms
Why do we need “special” HB terms? When a hydrogen atom is bonded to very electronegative atoms HD (donor: F, O, Cl, N, S), the charge is moved toward the bond midpoint so that there is much less charge remaining at the center of the H. This leads to a strong coulomb 2-body interaction with other electonegative atoms (A for acceptors) which we include in the Q H Q A coulomb terms. However the H-A vdw interaction should be reduced since most of the charge on H has moved toward D. Thus we need to modify the H--A vdW term. The standard terms would push the H away from A.
Various FF have different strategies to handle this problem L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 18
DREIDING Hydrogen Bond Term
General Form of a 3-Body DREIDING Hydrogen Bond Term:
E hb
(
R DA
,
AHD
) ( ) (
AHD
)
E hb
D hb
[ 5 (
R hb
/
R DA
) 12 6 (
R hb
/
R DA
) 10 ] cos 4 (
DHA
) Lennard-Jones 12-10 cos 4 th power
D hb : Hydrogen Bond Well Depth R hb : Equilibrium Hydrogen Bond Distance R DA : Donor-Acceptor Distance θ DHA : Donor-Hydrogen-Acceptor Angle
• With constraints set forth by VDW and point charges, difficult to accurately describe polar interactions without a HB Term • Examine water dimer structure to determine best radial and angular functional form L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 19
Dreiding Hydrogen Bond
H The Dreiding HB is factored as EHB(R AD , AHD ) = E d (R AD )E a ( AHD ) A AHD R AD distance term: LJ12-10 (Dreiding II):
b
AR
12
BR
10
D e
[5(
R e R
) 12
b R e
10 ) ]
D e
[ 2
R
D Morse:
b
D e
[ 2 where
e
(
e
)
e
2 (
R R e
1) Here R e where is the equilibrium distance between acceptor and donor (
R
(A-D), and D e
e
e
)
e
2
R e
1) Angle term: Dreiding II:
E a
4 =0 if AHD < 90 ° Dreiding III: L04-Ch121a-Goddard E a ( AHD ) = (cos AHD ) 2 =0 if AHD < 90 ° 20
Caitlin Scott stopped here Apr. 9
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 21
Equilibrium Water Dimer Structure
• X3LYP/aug-cc-pvtz(-f) minimized structure • Binding Energy (BE): – E bind = E dimer - 2*E water – For water dimer: • 5.00 kcal/mol Donor
D
L04-Ch121a-Goddard
H A
Acceptor © copyright 2013 William A. Goddard III, all rights reserved 22
Charge Assignment
Basis Set DFT-B3LYP 6-311G** 6-31G** Aug-cc-pvtz(-f) 6-31G/HF Experiment Mulliken Charge on Monomer Oxygen (e) -0.48
Dipole Moment with Mulliken Charges (Debye) 1.35
-0.61
1.72
-0.67
-0.87
n/a 1.91
2.43
1.8
Mulliken charges: assign partial atomic charges L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 23
Angular Dependence for water dimer
Angle Dependence
4 3 2 -4 -5 -6 1 0 -1 0 -2 -3 BE (cos 4 th BE (cos 3 rd BE (cos 2 nd power) power) power) 50 QM_BE BE (cos 1 st 100 power) 150 200 250 300 VDW 350 400 Coulomb
Rotation Angle
• QM constrain O-O distances and rotating donor hydrogen bond water – Plot cos(θ AHD ), cos 2 ( θ AHD ), cos 3 ( θ AHD ), cos 4 ( θ AHD ) • cos 2 ( θ) chosen: – Better fit compared to cos 4 ( θ) Tod redo this plot from 0 to 90 – With vanishing derivative at 90 degrees L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 24
2 1
Radial Dependence of HB for Water dimer
Constrain O H…O angle = 180 ° Use 6-31G** charges for DREIDING
Water Dimer Binding
Morse-Potential (best fit) is shown γ=9.70, R 0 =3.10, D 0 =1.75
Modify off-diagonal VDW terms -3 -4 0 2 -1 -2 2,2 -5 L04-Ch121a-Goddard 2,4 2,6 2,8 HB 3 Coulomb 3,2 3,4 3,6 QM Total (DREIDING)
O-O Dist (in Angstroms)
water dimer O-O (A) © copyright 2013 William A. Goddard III, all rights reserved 4 25
Final Form of New Hydrogen Bond
Radial and angular considerations leads to the following updated DREIDING hydrogen bond term:
Updated DREIDING Hydrogen Bond Term:
E hb
D hb
[ 2 2 ] cos 2 (
DHA
)
where
e
2
R R HB
1 Fitted with γ=9.70, R hb and D hb L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 26
Dreiding III HB Parameters for Amino Acids fitted to QM
H O H
O_3W
Water O H
O_3
Methanol Ser,Thr O
O_2
H N H
N_R
Amide Asn, Gln, Amide H N
N_A
Imidazole
S_3
S H H O
O_R
N
N_A
His, Trp Cys,Met • 30 pairs of HB donor acceptor parameters – 7 atom types – 30 pairs of model compounds QM data – 30 pairs of D Phenol Thiol Tyr Model compounds for amino acids neutral at standard pH Red: Atom types involved in HB Donor/Acceptors L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved hb , R • Parameters fitted to within 0.01A and 0.1 B3LYP with 6-31G** hb kcal/mol of QM values • Mulliken charges from 27
Examples Parameters for Dreiding III HB
Donor Atom Type O_3 Acceptor Atom O_3W O_3 O_2 O_R N_A S_3 Model (Donor Acceptor) MeOH - H 2 O MeOH - MeOH QM BE/R kcal/mol/Ǻ D HB BE/R kcal/mol/Ǻ 4.77 / 2.93 1.5 / 2.925
5.16 / 2.90 0.8 / 2.85
CH 3 OH – Amide 6.35 / 2.83 1.3 / 2.75
MeOH – C 6 H 5 OH MeOH – Me-Im 3.434 / 3.00
0.40 / 3.09
6.33 / 2.89 2.70 / 2.79
MeOH – MeSH 3.79 / 3.40 2.50 / 3.25
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 28
HB parameters in Dreiding II and III
Dreiding II *HBOND TYPE -DE HB RE HB *X -X 1 -9.0000 2.7500 ! no charges obsolete *X -X 1 -7.0000 2.7500 ! Gasteiger charges obsolete *X -X 1 -4.0000 2.7500 ! "experimental" charges Dreiding III MPSIM_HB (A-H-D) TYPE -DE HB RE HB O_3 -H___A-O_3 1 -4.8000 2.7500
O_3 -H___A-O_3M 1 -4.8000 2.7500
O_3M -H___A-O_3 1 -4.8000 2.7500
O_3M -H___A-O_3M 1 -4.8000 2.7500
O_2 -H___A-O_3 1 -4.8000 2.7500
O_2 -H___A-O_3M 1 -4.8000 2.7500
* Note: (D-H-A) is correct order.
L04-Ch121a-Goddard 29
Universal Force Field (UFF)
Generic force field for all elements from H (Z=1) to Lr (Z=103) Bond, angle, dihedral, inversion, vdw, electrostatics 6 constants per atom describes all interactions A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard III, and W. M. Skiff;
UFF, A Full Periodic-table Force-field For Molecular Mechanics And Molecular-dynamics Simulations
J Am Chem Soc 114 10024-10035 (1992) Casewit C J, Colwell K S, Rappe A K,
Application Of A Universal Force-field To Main Group Compounds
J Am Chem Soc 114: 10046-10053 (1992) L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 30
Universal Force Field (UFF) : Bond stretch
Harmonic Morse R IJ = equilibrium distance =
X
I
=
electronegativity element I
n
= BO, l = 0.1332
Force Constant = D IJ = bond energy (Morse) =BO*70 kcal/mol L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 31
Universal Force Field (UFF) : Angle Bend
Based on general form Linear (n=1), trigonal-planar (n=3), square planar (n=4), octahedral (n=4) S, Se, Te: E( )=(K IJK /4)[1+cos(2 )] O: For all cases, K uses the bond force constant between the 1-3 neighbors, based on Z I * and Z K * . No new force constant for angles!
( ) L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 32
UFF Z=H-Ne
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 33
Universal Force Field (UFF) : Dihedral
sp 3 -sp 3 Same cases as for Dreiding, but different barriers sp 3 -sp 3 barriers based on the hydride, writing V sp3 =Sqrt(V J V K ) validation L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 34
Universal Force Field (UFF) : Dihedral -other
sp 2 -sp 2 C row Uj 2 Si row 1.25
Ge row 0.7
Sn row 0.2
Pb row 0.1
Both in O column Period 2, minimum at 90, V J = 2.0 kcal/mol O, VJ = 6.8 kcal/mol S, Se, Te One atom not main group: One atom sp 1 : V J =0 L04-Ch121a-Goddard V J =0 © copyright 2013 William A. Goddard III, all rights reserved 35
Inversion Terms
Atom I bonded to J, K, L angle ω is angle of IL bond from JIK plane or IJ bond from LIK plane or of IK bond from JIL plane. We do all 3 cases and average.
C_2 and C_R: C 0 = 1,C 1 = -1, C 2 = 0; K=6 kcal/mol, except K=50 kcal/mol if J, K, or L =O_2 For N column C 0 = 1,C 1 = 0, C 2 = 1; K=0 for N and K=22 for other column 15 (to fit inversion barrier of hydride) All other atoms set K=0 L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 36
UFF Na-Ca
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 37
UFF Sc-Tc
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 38
UFF Ru-Eu
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 39
UFF Gd-Tl
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 40
UFF Pb-Lr
Note: Lw Lr © copyright 2013 William A. Goddard III, all rights reserved 41
Compare UFF to experiment
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 42
Compare UFF to experiment
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 43
Compare UFF to experiment
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 44
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 45
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 46
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 47
Valence Force Fields
AMBER, Assisted Model Building and Energy Refinement AMBER/OPLS CHARMM, Chemistry at HARvard Macromolecular Mechanics DISCOVER, force fields of the Insight/Discover DREIDING, force fields of POLYGRAF/BIOGRAF GROMOS, GROningen MOlecular Simulation package MM2/MM3/MM4 , Allinger molecular mechanics FF MMFF94, the Merck Molecular Force Field Tripos, force field of the Sybyl molecular modeling program UFF: Universal Force Field Used in Maestro, Cerius, BioGraf L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 48
Some Existing Force Fields
AMBER (Assisted Model Building with Energy Refinement)
Proteins and nucleic acids
CHARMM (Chemistry at Harvard Macromolecular Mechanics)
Proteins and nucleic acids
CFF (Consistent Force Field)
Conformations, vibrational spectra, strain energy, and vibrational enthalpy of proteins. (Variations: UBCFF, CVFF, Lynghy CFF)
CHEAT (Carbohydrate Hydroxyls represented by External AToms)
Specifically designed for modeling Carbohydrates
DREIDING
All-purpose organic or bio-organic FF
GROMOS (Gronigen Molecular Simulation)
Predicting the dynamical motion of molecules, bulk liquids, and bio-molecules.
MMx
MM1,MM2(MMx, MM+),MM3, and MM4 are general-purpose organic FF.
MOMEC
FF for describing transition metal coordination compounds.
OPLS (Optimized Potentials for Liquids Simulations)
Bulk liquids, and bio-molecules.
Tripos (Tripos Inc.)
Organic and bio-organic molecules.
UFF
Full periodic table. Most widely used for systems containing inorganic elements.
YETI
Accurate representation of nonbonded interactions. Modeling interactions between biomolecules and small substrate molecules (not for geometry optimization but docking) 49
Optimized Potentials for Liquid Simulations OPLS-aa Jorgensen, Yale
geometric combination rules Intramolecular non-bonded interactions (
E
ab) counted for atoms three or more bonds apart 1, 4 interactions are scaled down by
f
ij = 0.5; otherwise,
f
ij = 1.0. Jorgensen WL, Tirado-Rives J (1988). JACS 110 : 1657 –1666. Jorgensen WL, Maxwell DS, Tirado-Rives J (1996). JACS. 118 : 11225 –11236. © copyright 2013 William A. Goddard III, all rights reserved 50
OPLS
A distinctive feature of the OPLS parameters is that they were optimized in Monte Carlo calculations of the liquid at 330K to fit experimental density and heat of vaporization of the liquid, in addition to fitting gas-phase torsional profiles. The problem is that the parameters are not transferable to other molecules.
Cannot describe new compounds. OPLS simulations in aqueous solution typically use the TIP4P or TIP3P water model.
Later discussed in class. I recommend Levitt’s F3C (Flexible 3 charge model originally from Ferguson rigid) with Q(H)=0.39697, H-bond (Ro=2.5, Do=3.2 LJ12-10) (density, heat of vaporization, Surface tension, dielectric constant, self-diffusion coefficient, radial Distribution functions) L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 51
AMBER (Assisted Model Building with Energy Refinement) Force field
The functional form of the AMBER force field
^
Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM Jr, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA "A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules".
J.
Am. Chem. Soc. 117 : 5179 –5197 (1995) L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 52
AMBER Force Field
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 53
Amber HB terms
AMBER 1984: The charges assigned to the H and A are used to account for the bulk of the attractive HB interactions. This is supplemented by a weak 12-10 potential (dependent upon the H-A distance) designed to help adjust the resulting H-A distance. The normal 12-6 interaction between H and A is ignored.
Tod I believe that later Amber changed this to a 12-6 potential Please update L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 54
CHARMM force field (Chemistry at HARvard Macromolecular Mechanics)
Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M (1983). "
CHARMM
: A program for macromolecular energy, minimization, and dynamics calculations". J Comp Chem 4 : 187 –217. MacKerell, A.D., Jr.; Brooks, B.; Brooks, C. L., III; Nilsson, L.; Roux, B.; Won, Y.; Karplus, M. (1998). "
CHARMM
: The Energy Function and Its Parameterization with an Overview of the Program". in Schleyer, P.v.R.; et al..
The Encyclopedia of Computational Chemistry
.
1
. Chichester: John Wiley & Sons. pp. 271 –277. MacKerell, Jr. AD,
et al.
(1998). "All-atom empirical potential for molecular modeling and dynamics studies of proteins".
J Phys
Chem B 102 : 3586 –3616.
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 55
CHARMm Force Field
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 56
CHARMM
HB terms
CHARMM 1984 The normal van der Waals and electrostatic interactions of the H with all other atoms is ignored, and replaced with a special HB potential involving the D-A distance and the D-H-A angle (described later).
I believe that later CHARMM changed this to a 12-6 potential Please update L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 57
Consistent Valence Force Field
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 58
Add in discussion of water FF
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 59
Practical issues Cutoffs
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 60
Need for a Non-bond cutoff
Calculation of the energy and forces constitute the main computational time in structure optimization, Molecular Dynamics and Monte Carlo calculations. The valence terms scale linearly with the number of atoms since they involve near neighbor atoms (1-2,1-3,1-4 interactions), but there are N*(N-1)/2 non-bonded interactions (both Coulomb and vdW). Thus for hemoglobin with 6000 atoms, there are 18,000,000 NB interactions to evaluate.
In order to reduce the computational cost, we would like to consider only the closer interactions by assuming that the non bonded interactions becomes negligible beyond cutoff distance.
R ij < R cut 61 L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved
Rcut
The first question is how large to take Rcut. Thus two unit charges separated by 10 A have a Coulomb interaction of 33.2 kcal/mol, while charges separated by 100 A still have 3.32 kcal/mol Coulomb energy. However common practice is a cut off of 12A. This is because the number of positive and negative interactions tend to cancel at larger distances General practice is to generate a list of atom pairs within the cutoff distance (NB List). This NB list must be updated periodically (generally every 50 to 100 steps) during the of dynamics. For small systems (< 500), it is practical and more accurate to include
all
nonbonded interactions in the calculations. Since it is the forces that are important in the MD, one cannot just truncate the forces at a finite distance because when atoms pass from beyond Rcut to below or vice versa, the forces would 62 L04-Ch121a-Goddard
Switching function: Cubic Spline
has the form
on
,
on
,
on
,
x off x off x off
) ) )
S x x on on on
, ,
x x x off off off
) )
x off
x
(
x
2
on on x
,
x off off
)
off x off x
) 2 )
on
)
x
3 (
x off
3
x on x
( ) 2 ) ( if 2
x x on
x off
x on
f )
x x
3
x on
x x on x
)
x x x on
f
x off
where 1 if
x x on
2 (
x off
x
( 2 ) (
x off x off
x on
2 ) 3
x
3
x on
) x = R ij if
x on x off
x = R x = R 2 on = R 0 i x off 2 off off ij 2 on 2 off
x on
x
Generally R off is taken to be ~ 1 to 2 A smaller than R cut x = R 2 atom pairs that are initially greater than R off so that apart do not move closer than R off during the time that the NB List is not updated. on on Common choice is R off cut 2 off = 12 A, R on = 10 A, R off = 8 A.
L04-Ch121a-Goddard © copyright 2013 William A. Goddard III, all rights reserved 63
x off
Cubic Spline
1.2
S r 1.0
0.8
0.6
0.4
0.2
0.0
0 2 Cubic Spline 4 S(r)/r 6 S(r) 1/r 8 R off 1.2
S r 1.0
0.8
0.6
0.4
0.2
10 r A 0.0
0 1/r 6 S(r)/r 2 6 4 Cubic Spline 6 S(r) 8 R off 10 r A L04-Ch121a-Goddard R on = 0 A, R off = 9 A © copyright 2013 William A. Goddard III, all rights reserved 64
S r 1.2
1.0
0.8
0.6
0.4
0.2
0.0
0
Higher order switching function: Taper function
Taper function S r 1.2
1.0
0.8
0.6
0.4
0.2
0.0
0 1/r 6 Tap(r) 1/r 2 Tap(r)/r 4 Taper function 6 Tap(r) 8 R off 10 r A Tap(r)/r 6 2 4 6 8 R on = 0 A, R off = 9 A R off 10 r A 7 th order polynomial:
Tap
n
7 0
Tap R n ij n
where
Tap
0 1
Tap Tap
4 1
Tap
2
Tap
3 35 /
R
4
cut Tap Tap Tap
5 6 7 84 / 5
R cut
70 / 6
R cut
20 /
R
7
cut
0 65
Effect of NB cutoff (216 water molecules)
Eele EvdW 690 -2600 688 -2700 686 -2800 684 -2900 682 -3000 680 -3100 678 -3200 5 10 15 30 35 5 10 15 20 Cutoff (A) 25 30 20 Cutoff(A) 25 Etot -3200 -3300 -3400 -3500 220000 200000 180000 160000 140000 -3600 120000 100000 -3700 80000 -3800 5 10 15 30 35 60000 5 10 15 20 Cutoff (A) 25 30 20 Cutoff (A) 25 Electrostatic interaction is most sensitive to NB cutoff (R off ) since it falls off more slowly than vdW. Energy converges at R off L04-Ch121a-Goddard = ~15 A (R on = 0 A, R cut =R off +1A). © copyright 2013 William A. Goddard III, all rights reserved 35 Nnb 35 66