Computational Quantum chemistry

Download Report

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