Ch121a Atomic Level Simulations of Materials and Molecules

1. Room BI 115 Hours: Monday, Wednesday, Friday 2-3pm Lecture 4, April 9, 2014 FF2: standard FF Ch121a Atomic Level Simulations of Materials and Molecules Presented by Caitlin Scott William A. Goddard III, wag@wag.caltech.edu Charles and Mary Ferkel Professor of Chemistry, Materials Science, and Applied Physics, California Institute of Technology TA’s Caitlin Scott and Andrea Kirkpatrick

2. 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 5weeks is to turn in a one page report on progress with the project The final is a research report describing the calculations and conclusions

3. Dreiding Dreiding is a generic or rule-based force field. Parameters were based on general principles, not fitted to specific molecules Bond distance = RA + RB -0.01A where RA are bond radii based on the A-CH3 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)

4. Dreiding Force constants” Kbond = BO*700 kcal/mol/Å2 where BO bond order (1,2,3) Kangle = 100 kcal/mol/rad2. Inversion barrier for planar molecules = 40 kcal/mol/rad2 torsion barrier = 2.0 kcal/mol for single bonds 45 kcal/mol for double bonds Dreiding is a generic or rule-based force field. Parameters were based on general principles, not fitted to specific molecules Bond distance = RA + RB -0.01A where RA are bond radii based on the A-CH3 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)

5. Dreiding Atom Types C_1, C_2, C_3 indicate sp,sp2, and sp3 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 sp2 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

6. Example Atom Types tetrapeptide Lys-Met-Phe-Pro

7. Geometric Valence parameters for Dreiding

8. The van der Waals Parameters for Dreiding R0 is total bond distance, vdw radius is ½ this size R0=3.195 for H was a bad choice. Should have been 2.6-2.7. it is compensated by a small D0 Parameters in the same row are similar. General trend to larger Re and De as go down the columns

9. Dreiding parameters for dihedrals E(φ) = ½ B[1+(-1)CM cos(Pφ)] a Two sp3 atoms: B=2, P=3, 0° is max b One sp3-one sp2: B=1, P=6, 0° is min c two sp2(BO=2): B=45, P=2, 0° is min d two sp2R (BO=1.5): B=25, P=2, 0° is min e two sp2(BO=1): B=5, P=2, 0° is min f two sp2R (BO=1.0): B=10, P=2, 0° is min g one or two sp1, monovalent (F..), metals (Fe..) all have B=0 h sp3-sp3 in O column: B=2, P=2, 0° is max i sp3 O column-sp3 other column: B=2 P=2, 0° is min j exception to b: if sp2 is bonded to one H: then B=2, P=3, 0° to H is max (eg propene) k barriers should decrease by factor of 3 as go down column, but this trend was ignored in Dreiding B=barrier (kcal/mol) P = periodicity CM=+1 if 0° is min = 0 if 0° is max

10. Validation: experimental barriers (kcal/mol)

11. Compare conformations to experiment

12. Dreiding validate using 1st 76 molecules in Cambridge Crystallographic Data Base (in 1988)

13. Dreiding accuracy over 76 molecules Note that the test minimized the molecule in a vacuum. The experiment was in a crystal

14. Complex case Linear CF3-(CF2)n-CF3 L K 150º Get F—F steric interactions for F on C1 and C3 164º 180º Not all trans as for CH3-(CH2)n-CH3 Get 163º torsion not 180º 14

15. Energy ComponentsMolecular origin of helicity in Teflon FF electrostatic QM vdW

16. High Pressure forms of C20F42 The figure shows predicted stable helical conformations for C20F42. 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.

17. DREIDING DREIDING Internal Energies Bond, Angle, Torsion, Inversion Electrostatics Point charges Dielectric=1.0 VDW Lennard-Jones 12-6, exp-6, Morse… Hydrogen Bonds 3-body term, Lennard-Jones 12-10 with a cos4 term

18. 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 QHQAcoulomb 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

19. DREIDING Hydrogen Bond Term General Form of a 3-Body DREIDING Hydrogen Bond Term: cos 4th power Lennard-Jones 12-10 Dhb: Hydrogen Bond Well Depth Rhb: Equilibrium Hydrogen Bond Distance RDA: 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

20. Dreiding Hydrogen Bond H AHD A The Dreiding HB is factored as EHB(RAD, AHD) = Ed(RAD)Ea(AHD) RAD distance term: D LJ12-10 (Dreiding II): where Morse: Here Re is the equilibrium distance between acceptor and donor (A-D), and De is the energy well depth. Angle term: =0 if AHD < 90° Dreiding II: =0 if AHD < 90° Dreiding III: Ea(AHD) = (cos AHD)2

21. Caitlin Scott stopped here Apr. 9

22. Equilibrium Water Dimer Structure • X3LYP/aug-cc-pvtz(-f) minimized structure • Binding Energy (BE): • Ebind = Edimer - 2*Ewater • For water dimer: • 5.00 kcal/mol Donor D A H Acceptor

23. Charge Assignment Mulliken charges: assign partial atomic charges

24. Angular Dependence for water dimer BE (cos 4th power) BE (cos 3rd power) BE (cos 2nd power) VDW QM_BE Coulomb BE (cos 1st power) • QM constrain O-O distances and rotating donor hydrogen bond water • Plot cos(θAHD), cos2(θAHD), cos3(θAHD), cos4(θAHD) • cos2(θ) chosen: • Better fit compared to cos4(θ) • With vanishing derivative at 90 degrees Tod redo this plot from 0 to 90

25. Radial Dependence of HB for Water dimer Constrain O-H…O angle = 180° Use 6-31G** charges for DREIDING Morse-Potential (best fit) is shown γ=9.70, R0=3.10, D0=1.75 Modify off-diagonal VDW terms VDW HB Coulomb QM Total (DREIDING) water dimer O-O (A)

26. Final Form of New Hydrogen Bond Radial and angular considerations leads to the following updated DREIDING hydrogen bond term: Updated DREIDING Hydrogen Bond Term: where Fitted with γ=9.70, Rhb and Dhb

27. Dreiding III HB Parameters for Amino Acids fitted to QM O_3W O_3 Ser,Thr • 30 pairs of HB donor-acceptor parameters • 7 atom types • 30 pairs of model compounds QM data • 30 pairs of Dhb, Rhb • Parameters fitted to within 0.01A and 0.1 kcal/mol of QM values • Mulliken charges from B3LYP with 6-31G** O_2 N_A N_R N_A His, Trp Asn, Gln, Amide S_3 O_R Cys,Met Tyr Model compounds for amino acids neutral at standard pH Red: Atom types involved in HB Donor/Acceptors

28. Examples Parameters for Dreiding III HB

29. HB parameters in Dreiding II and III *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 II 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. TYPE = 1: LJ 12-10 TYPE = 2: Morse

30. 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 SimulationsJ 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 CompoundsJ Am Chem Soc 114: 10046-10053 (1992)

31. Universal Force Field (UFF) : Bond stretch RIJ = equilibrium distance = Force Constant = XI = electronegativity element I XI = electronegativity element I XI = electronegativity element I Harmonic Morse n = BO, l = 0.1332 DIJ = bond energy (Morse) =BO*70 kcal/mol

32. 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()=(KIJK/4)[1+cos(2)] O: For all cases, K uses the bond force constant between the 1-3 neighbors, based on ZI* andZK*. No new force constant for angles!

33. UFF Z=H-Ne

34. Universal Force Field (UFF) : Dihedral sp3-sp3 Same cases as for Dreiding, but different barriers sp3-sp3 barriers based on the hydride, writing Vsp3 =Sqrt(VJVK) validation

35. Universal Force Field (UFF) : Dihedral -other Uj C row 2 Si row 1.25 Ge row 0.7 Sn row 0.2 Pb row 0.1 sp2-sp2 Both in O column Period 2, minimum at 90, VJ = 2.0 kcal/mol O, VJ = 6.8 kcal/mol S, Se, Te One atom not main group: VJ=0 One atom sp1: VJ=0

36. 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: C0 = 1,C1 = -1, C2 = 0; K=6 kcal/mol, except K=50 kcal/mol if J, K, or L =O_2 For N column C0 = 1,C1 = 0, C2 = 1; K=0 for N and K=22 for other column 15 (to fit inversion barrier of hydride) All other atoms set K=0

37. UFF Na-Ca

38. UFF Sc-Tc

39. UFF Ru-Eu

40. UFF Gd-Tl

41. UFF Pb-Lr Note: LwLr Element 103

42. Compare UFF to experiment

43. Compare UFF to experiment

44. Compare UFF to experiment

48. 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

49. 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)

50. Optimized Potentials for Liquid SimulationsOPLS-aa Jorgensen, Yale geometric combination rules Intramolecular non-bonded interactions (Eab) counted for atoms three or more bonds apart 1, 4 interactions are scaled down by fij = 0.5; otherwise, fij = 1.0. Jorgensen WL, Tirado-Rives J (1988). JACS 110: 1657–1666. Jorgensen WL, Maxwell DS, Tirado-Rives J (1996). JACS.118: 11225–11236.