1 / 46

Ch121a Atomic Level Simulations of Materials and Molecules

Ch121a Atomic Level Simulations of Materials and Molecules. BI 115 Hours: 2:30-3:30 Monday and Wednesday Lecture or Lab: Friday 2-3pm (+3-4pm). Lecture 3, April 6, 2011 Force Fields – 1 valence, QEq, vdw. William A. Goddard III, wag@wag.caltech.edu

bart
Download Presentation

Ch121a Atomic Level Simulations of Materials and Molecules

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Ch121a Atomic Level Simulations of Materials and Molecules BI 115 Hours: 2:30-3:30 Monday and Wednesday Lecture or Lab: Friday 2-3pm (+3-4pm) Lecture 3, April 6, 2011 Force Fields – 1 valence, QEq, vdw William A. Goddard III, wag@wag.caltech.edu Charles and Mary Ferkel Professor of Chemistry, Materials Science, and Applied Physics, California Institute of Technology Teaching Assistants Wei-Guang Liu, Fan Lu, Jose Mendozq, Andrea Kirkpatrick

  2. The full Quantum Mechanics Hamiltonian: Htotal(1..Nel,1..Mnuc) Hel(1..Nel) = N M M,N M N nuclei electron-nucleus electrons nucleus-nucleus electron-electron Born-Oppenheimer approximation, fix nuclei (Rk=1..M), solve for Ψel(1..N) Hel(1..Nel) Ψel(1..Nel)= EelΨel(1..Nel) Ψel(1..Nel)and Eel(1..Mnuc) are functions of the nuclear coordinates. Next solve the nuclear QM problem Hnuc(1..Mnuc) Ψnuc(1..Mnuc)= EnucΨnuc(1..Mnuc) Hnuc(1..Mnuc) = +Eel(1..Mnuc) FF is analytic fit to this EPE(1..Mnuc)

  3. Force Field EPE(1..Mnuc) = Eel(1..Mnuc) The Force Field (FF) is an analytic fit to such expressions but formulated to be transferable so can use the same terms for various molecules and solids General form: Ecross-terms Involves covalent bonds and the coupling between them (2-body, 3-body, 4-body, cross terms) E nonbond = Evdw + Eelectrostatic Involves action at a distance, long range interactions

  4. Bond Stretch Terms: Harmonic oscillator form R Taylor series expansion about the equilibrium, Re, d = R-Re E(d)=E(Re)+(dE/dR)e[d]+½(d2E/dR2)e[d]2+(1/6)(d3E/dR3)e[d]3 + O(d4) ignore ignore zero ignore Kb Re = 0.9 – 2.2 Å Ke = 700 (Kcal/mol)/Å 2 E(R) = (ke/2)(R-Re)2, Simple Harmonic Units: multiply by 143.88 to convert mdynes/ Å to (Kcal/mol)/Å 2 eigenstates, En = n + ½ , where n=0,1,2,… Ground state wavefunction (Gaussian) 0(d) = (mw/pЋ) exp[- (mw/2Ћ)d2] Re Some force fields, CHARRM, Amber use k=2Ke to avoid multiplying by 2 Problem: cannot break the bond, E  ∞

  5. Bond Stretch Terms: Morse oscillator form R R0 D0 Want the E to go to a constant as break bond, popular form is Morse function where X = exp[-a(R-Re)] = exp[-(g/2)(R/Re -1)] Evdw (Rij) = De[X2-2X] Note that E(Re) = 0, the usual convention for bond terms De is the bond energy in kcal/mol Re is the equilibrium bond distance a is the Morse scaling parameter calculate from Ke and De • = sqrt(ke/2De), where ke = d2E/dR2 at R=Re (the force constant) E = 0

  6. Morse Potential – get analytic solution for quantum vibrational levels where X = exp[-a(R-Re)] = exp[-(g/2)(R/Re -1)] Evdw (Rij) = De[X2-2X] a = sqrt(ke/2De), where ke = d2E/dR2 at R=Re (the force constant) Ev = hv0(n + 1/2 ) – [hv0(n + ½)]2/4De n0 = (a/2p)sqrt(2De/m) Level separations decrease linearly with level En+1 – En = hv0 – (n+1)(hv0)2/2De Write En/hc = we(n + ½ ) - xewe(n + ½ )2

  7. The Bending potential surface for CH2 1B1 1Dg 1A1 3B1 3Sg- 9.3 kcal/mol

  8. Angle Bend Terms: cosine harmonic J  K I E() = C0 + C1 cos() + C2 cos(2) + C3 cos(3) + …. ~ ½ Ch [cos  - cos e]2 = 0 at  = e E”() = dE()/d = - Ch [cos  - cos e] sin  = 0 at  = e E”() = d2E()/d2 = - Ch [cos  - cos e] cos  + Ch (sin )2 = Ch (sin e)2 at  = e Thus k = Ch (sin e)2

  9. Barriers for angle term The energy barrier to “linearize” the molecule becomes By symmetry the angular energy satisfy This is always satisfied for the cosine expansion but A second more popular form is the Harmonic theta expansion E() = ½ K  + [ - e]2 However except for linear molecules this does NOT satisfy the symmetry relations, leading to undefined energies and forces for  = 180° and 0°. This is used by CHARMM, Amber, ..

  10. Angle Bend Terms for linear molecule, K J  I If e = 180° then we write E() = K [1 + cos()] since for  = 180 – d this becomes E(d) = K [1 + (-1 + ½ d2)] ~ 1/2 Kd2

  11. Angle Bend Terms Simple Periodic Angle Bend Terms For an Octahedral complex The angle function should have the form E() = C [1 - cos4] which has minima for  = 0, 90, 180, and 270° With a barrier of C for  = 45, 135, 225, and 315° Here the force constant is K = 16C= CP2 where P=4 is the periodicity, so we can write C=K/P2

  12. Angle Bend Terms Simple Periodic Angle Bend Terms • For an Octahedral complex • The angle function should have the form • E() = C [1 - cos4] which has minima for  = 0, 90, 180, and 270° With a barrier of C for  = 45, 135, 225, and 315° Here the force constant is K = 16C= CP2 where P=4 is the periodicity, so we can write C=K/P2 However if we wanted the minima to be at • = 45, 135, 225, and 315° with maxima at  = 0, 90, 180, and 270° then we want to use E() = ½ C [1 + cos4] Thus the general form is E() = (K/P2)[1 – (-1)Bcos4] Where B=1 for the case with a minimum at 0° and B=-1 for a maximum at 0°

  13. Rotational barriers 7.6 kcal/mol Cis barrier HOOH Part of these barriers can be explained as due to vdw and electrostatic interactions between the H’s. But part of it arises from covalent terms (the pp lone pairs on each S or O This part has the form E(φ) = ½ B [1+cos(2φ)] Which is E=0 for φ=90,270° and E=B for φ=0,180° 1.19 kcal/mol Trans barrier HSSH: 5.02 kcal/mol trans barrier 7.54 kcal/mol cis barrier

  14. dihedral or torsion terms J L φ L I J K K I Given any two bonds IJ and KL attached to a common bond JK, the dihedral angle φ is the angle of the JKL plane from the IJK plane, with cis being 0 E(φ) = C0 + C1 cos(φ) + C2 cos(2 φ) + C3 cos(3 φ) + ….which we write as Where each kn is in kcal/mol, n = 1, is the periodicity of the potential and d=+1 the cis conformation is a minimum d=-1 the cis conformation is thea maximum Input data is K and d for each n.

  15. Consider the central CC bond in Butane L I J K There are nine possible dihedrals: 4 HCCH, 2 HCCC, 2 CCCH, and 1 CCCC. Each of these 9 could be written as E(φ) = ½ B [1 – cos(3 φ)], For which E=0 for φ = 60, 180 and 300° and E=B the rotation barrier for φ = 0, 120 and 240°.

  16. Consider the central CC bond in Butane L I J K There are nine possible dihedrals: 4 HCCH, 2 HCCC, 2 CCCH, and 1 CCCC. Each of these 9 could be written as E(φ) = ½ B [1 – cos(3 φ)], For which E=0 for φ = 60, 180 and 300° and E=B the rotation barrier for φ = 0, 120 and 240°. For ethane get 9 HCCH terms. The total barrier in ethane is 3kcal/mol but standard vdw and charges of +0.15 on each H will account for ~1 kcal/mol. Thus only need explicit dihedral barrier of 2 kcal/mol.

  17. Consider the central CC bond in Butane L I J K There are nine possible dihedrals: 4 HCCH, 2 HCCC, 2 CCCH, and 1 CCCC. Each of these 9 could be written as E(φ) = ½ B [1 – cos(3 φ)], For which E=0 for φ = 60, 180 and 300° and E=B the rotation barrier for φ = 0, 120 and 240°. For ethane get 9 HCCH terms. The total barrier in ethane is 3kcal/mol but standard vdw and charges of +0.15 on each H will account for ~1 kcal/mol. Thus only need explicit dihedral barrier of 2 kcal/mol. Can do two ways. Incremental: treat each HCCH as having a barrier of 2/9 and add the terms from each of the 9 to get a total of 2 (Amber, CHARRM) Or use a barrier of 2 kcal/mol for each of the 9, but normalize by the total number of 9 to get a net of 2 (Dreiding)

  18. Twisting potential surface for ethene The ground state (N) of ethene prefers φ=0º to obtain the highest overlap with a rotational barrier of 65 kcal/mol for prefers φ =90º to obtain the lowest overlap. We write this as E(φ) = ½ B [1-cos2 φ)]. Since there are 4 HCCH terms, we calculate each using the full B=65, but divide by 4. T N φ = 90° φ = 0° φ = -90°

  19. Inversions L J I K L J I K When an atom I has exactly 3 distinct bonds IJ, IK, and IL, it is often necessary to include an exlicit term in the force field to adjust the energy for “planarizing” the center atom I. Where the force constant in kcal/mol is Umbrella inversion: ω is the angle between The IL axis and the JIK plane

  20. AMBER Improper Torsion JILK J L I K Amber describes inversion using an improper torsion. N=2 for planar angles (0=180°) and n=3 for the tetrahedral Angles (0 =120°). Improper Torsion:  is the angle between the JIL plane and the KIL plane

  21. CHARMM Improper Torsion IJKL L J I K In the CHARMM force field, inversion is defined as if it were a torsion. For a tetrahedral carbon atom with equal bonds this angles is φ=35.264°. Improper Torsion: φ is the angle between the IJK plane and the LJK plane

  22. Summary: Valence Force Field Terms Description Illustration Typical Expressions Harmonic Stretch Bond Stretch Harmonic-cosine Angle bend Dreiding dihedral Torsion Umbrella inversion Inversion

  23. Coulomb (Electrostatic) Interactions Most force fields use fixed partial charges on each nucleus, leading to electrostatic or Coulomb interactions between each pair of charged particles. The electrostatic energy between point charges Qi and Qk is described by Coulomb's Law as EQ(Rik) = C0 Qi Qk /(e Rik) where Qi andQk are atomic partial charges in electron units C0 converts units: if E is in eV and R in A, then C0 = 14.403 If E is in kcal/mol and R in A, then C0 = 332.0637 e = 1 in a vacuum (dielectric constant) TA check numbers

  24. How estimate charges? Even for a material as ionic as NaCl diatomic, the dipole moment  a net charge of +0.8 e on the Na and -0.8 e on the Cl. We need a method to estimate such charges in order to calculate properties of materials. First a bit more about units. In QM calculations the unit of charge is the magnitude of the charge on an electron and the unit of length is the bohr (a0) Thus QM calculations of dipole moment are in units of ea0 which we refer to as au. However the international standard for quoting dipole moment is the Debye = 10-10 esu A Where m(D) = 2.5418 m(au)

  25. ReaxFF Electrostatics energy – Critical element I atomic interactions Keeping: Jij Hardness (IP-EA) Electronegativity (IP+EA)/2 1/rij rij ri0 +rj0 • Self-consistent Charge Equilibration (QEq) • Describe charges as distributed (Gaussian) • Thus charges on adjacent atoms shielded • (interactions  constant as R0) and include interactions over ALL atoms, even if bonded (no exclusions) • Allow charge transfer (QEq method) 2 Three universal parameters for each element: 1991: used experimental IP, EA, Ri; ReaxFF get all parameters from fitting QM

  26. Charge Equilibration First consider how the energy of an atom depends on the net charge on the atom, E(Q) Including terms through 2nd order leads to • Charge Equilibration for Molecular Dynamics Simulations; • K. Rappé and W. A. Goddard III; J. Phys. Chem. 95, 3358 (1991) (2) (3)

  27. Charge dependence of the energy (eV) of an atom E=12.967 E=0 E=-3.615 Cl+ Cl Cl- Q=+1 Q=0 Q=-1 Harmonic fit Get minimum at Q=-0.887 Emin = -3.676 = 8.291 = 9.352

  28. QEq parameters

  29. Interpretation of J, the hardness Define an atomic radius as RA0 Re(A2) Bond distance of homonuclear diatomic H 0.84 0.74 C 1.42 1.23 N 1.22 1.10 O 1.08 1.21 Si 2.20 2.35 S 1.60 1.63 Li 3.01 3.08 Thus J is related to the coulomb energy of a charge the size of the atom

  30. The total energy of a molecular complex Consider now a distribution of charges over the atoms of a complex: QA, QB, etc Letting JAB(R) = the Coulomb potential of unit charges on the atoms, we can write Taking the derivative with respect to charge leads to the chemical potential, which is a function of the charges or The definition of equilibrium is for all chemical potentials to be equal. This leads to

  31. The QEq equations Adding to the N-1 conditions The condition that the total charged is fixed (say at 0) leads to the condition Leads to a set of N linear equations for the N variables QA. AQ=X, where the NxN matrix A and the N dimensional vector A are known. This is solved for the N unknowns, Q. We place some conditions on this. The harmonic fit of charge to the energy of an atom is assumed to be valid only for filling the valence shell. Thus we restrict Q(Cl) to lie between +7 and -1 and Q(C) to be between +4 and -4 Similarly Q(H) is between +1 and -1

  32. The QEq Coulomb potential law We need now to choose a form for JAB(R) A plausible form is JAB(R) = 14.4/R, which is valid when the charge distributions for atom A and B do not overlap Clearly this form as the problem that JAB(R)  ∞ as R 0 In fact the overlap of the orbitals leads to shielding The plot shows the shielding for C atoms using various Slater orbitals Using RC=0.759a0 And l = 0.5

  33. QEq results for alkali halides

  34. QEq for Ala-His-Ala Amber charges in parentheses

  35. QEq for deoxy adenosine Amber charges in parentheses

  36. QEq for polymers Nylon 66 PEEK

  37. Polarizable QEq Allow each atom to have two charges: A fixed core charge(+4 for Si)with a Gaussian shape A variable shell charge with a Gaussian shape but subject to displacement and charge transfer Electrostatic interactions between all charges, including the core and shell on same atom Allow Shell to move with respect to core, to describe atomic polarizability Self-consistent charge equilibration(QEq) Four universal parameters for each element: Get from QM

  38. Noble gas dimers Ar2 s Re De All nonbonded atoms and molecules exhibit a very repulsive interaction at short distances due to overlap of electron pairs (Pauli Repulsion) and a weak attractive interaction scaling like 1/R6 at long R (London Dispersion. Together these are called van der Waals (vdW) interactions • Most popular form: LJ 12-6 • E=A/R12 –B/R6 • = De[r-12 – 2r-6] where r= R/Re • = 4 De[t-12 – t-6] where t=R/s Here s = Re(1/2)1/6 =0.89 Re

  39. van der Waals interaction • What vdW interactions account for? • Short range repulsion  Overlap of orbitals and the Pauli Principle Pauli repulsion: Born repulsion: • Long range attraction  Instantaneous multipolar interaction (London dispersion) Dipole-dipole Dipole-quadrupole Quadruole-quadrupole We often neglect higher order (>R6) terms.

  40. Popular vdW Nonbond Terms: LJ12-6 NonBond R • Lennard-Jones 12-6 • E(R)=A/R12 – B/R6 • =Dvdw{1/r12 – 2/r6} where r = R/Rvdw • = 4Dvdw{[svdw/R)12]- [svdw/R)6]} • Dvdw=well depth, • Rvdw = Equilibrium distance for vdw dimer • = point at which E=0 (inner wall, s ~ 0.89 Rvdw) Dimensionless force constant k = {[d2E/dr2]r=1}/Dvdw = 72 The choice of 1/R6 is due to London Dispersion (vdw Attraction) However there is no special reason for the 1/R12 short range form (other that it saved computation time for 1950’s computers) LJ 9-6 would lead to a more accurate inner wall.

  41. Popular vdW Nonbond Terms- Exponential-6 NonBond R Buckingham or exponential-6 E(R)= A e-CR - B/R6 which we write as whereR0 = Equilibrium bond distance; r = R/R0 is the scaled distance Do=well depth z = dimensionless parameter related to force constant at R0 We define a dimensionless force constant as k= {(d2E/dr2)r=1}/D0 k = 72 for LJ12-6 z = 12 leads to –D0/r6 at long R, just as for LJ12-6 z = 13.772 leads to k = 72 just as for LJ12-6 A problem with exp-6 is that E- ∞ as R 0. To avoid this BioGraf, LinGraf, CeriusII calculate the inner maximum and reflect E about this point for smaller R so that E and E’ are continuous. When z>10 this point is well up the inner wall and not important

  42. Converting between X-6 and LJ12-6 Usual convention: use the same R0 and D0 and set z = 12. I recall that this leads to small systematic errors. Second choice: require that the long range form of exp-6 be the same as for LJ12-6 (ie -2D0/r6) and require that the inner crossing point be the same. This leads to This was published in the Dreiding paper but I do not know if it has ever been used

  43. Popular vdW Nonbond Terms: Morse NonBond R E(R) = D0 {exp[-b(R-R0)] - 2 exp[(-(b/2)(R-R0)]} At R=R0, E(R0) = -D0 At R=∞, E(R0) = 0 D0 is the bond energy in Kcal/mol R0 is the equilibrium bond distance b is the Morse scaling parameter I prefer to write this as E(R) = D0 [X2 - 2X] where X = exp[(a/2)(1-r)] At R=R0, r = 1 X = 1  E(R0) = -D0 At R=∞, X = 0  E(R0) = 0 k= {(d2E/dr2)r=1}/D0 = a2/2 Thus a = 12  k= 72 just as for LJ12-6 D0 R0 Few theorists believe that Morse makes sense for vdw parameters since it does not behave as 1/R6 as R∞ However, the vdw curve matches 1/R6 only for R>6A, and for our systems there will be other atoms in between. So Morse is ok.

  44. vdW combination rules Generally the vdw parameters are provided for HH, CC, NN etc (diagonal cases) and the off-diagonal terms are obtain using combination rules D0IJ = Sqrt(D0II D0JJ) R0IJ = Sqrt(R0II R0JJ) zIJ = (zIJ + zIJ)/2 Sometimes an arithmetic combination rule is used for vdw radii R0IJ = (R0II +R0JJ)/2 but this complicates vdw calculations (the amber paper claims to do this but the code uses geometric combinations of the radii)

  45. 1-2 and 1-3 Nonbond exclusions For valence force fields, it is assumed that the bond and angle quantities include already the Pauli Repulsion and electrostatic contributions included in the nonbond interactions. Thus normally we exclude from the vdw and coulomb sums the contributions from bonds (1-2) and next nearest neighbor (1-3) interactions

  46. 1-4 Nonbond exclusions There is disagreement about 1-4 interactions. In fact the original of the rotational barrier in ethane is probably all due to Pauli repulsion orthogonality effects. Thus a proper description of the vdW interactions between 1-4 atoms should account for the barrier. In fact it accounts only for 1/3 of the barrier. I believe that this is because the CH bond pairs are centered at the bond midpoint not on the H atoms. Thus using atom centered vdw accounts for only part of the barrier necessitating an explicit dihedral term. Thus I believe that the1-4 vdw and electrostatic terms should be included. However some FF, such as Amber, CHARRM reduce the 1-4 interactions by a factor of 2. I do not know of a justification for this.

More Related