1 / 49

Ch121a Atomic Level Simulations of Materials and Molecules

Room BI 115 Lecture: Monday, Friday 2-3pm Lab Session: Wednesday 2:30-3:30pm Extra help: Monday 3-3:30pm and Wednesday 2-2:30pm. Ch121a Atomic Level Simulations of Materials and Molecules. Lecture 3 , April 8, 2013 FF1: valence, QEq , vdw. William A. Goddard III, wag@wag.caltech.edu

yardan
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. Room BI 115 Lecture: Monday, Friday 2-3pm Lab Session: Wednesday 2:30-3:30pm Extra help: Monday 3-3:30pm and Wednesday 2-2:30pm Ch121a Atomic Level Simulations of Materials and Molecules Lecture 3, April 8, 2013 FF1: 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

  2. Room BI 115 Hours: Monday, Wednesday, Friday 2-3pm Lecture 3, April 8, 2013 FF1: valence, QEq, vdw Ch121a Atomic Level Simulations of Materials and Molecules 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

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

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

  5. The Bending potential surface for CH2 C  H H 1B1 1Dg 1A1 3B1 3Sg- 9.3 kcal/mol linear Note that E(p-) = E(p+) Symmetric about  = p=180°

  6. 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 parameters 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

  7. 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 QM 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  ∞

  8. Bond Stretch Terms: Morse oscillator form R Re De 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

  9. Morse Potential – get analytic solution for quantum vibrational levels where X = exp[-a(R-Re)] = exp[-(g/2)(R/Re -1)] EMorse(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 D0 De

  10. Angle Bend Terms: cosine harmonic J  K I Expand E() as Fourier series. Note only cos(n) since must be symmetric about =0 and p E() = C0 + C1cos() + C2cos(2) + C3cos(3) + …. E() ~ ½ Ch [cos  - cos e]2 E=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 k = E”(e) = Ch(sin e)2 Using cos(2) = 2 cos2 - 1

  11. The Bending potential surface for CH2 C  H H 1B1 1Dg 1A1 3B1 3Sg- 9.3 kcal/mol linear Note that E(p-) = E(p+) Symmetric about  = p=180°

  12. Barriers for angle term The energy is a maximum at  = p. Thus energy barrier to “linearize” the molecule becomes By symmetry the angular energy satisfies 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, ..

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

  14. 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 = CP2where P=4 is the periodicity, so we can write C=K/P2

  15. 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°

  16. Dihedral barriers 7.6 kcal/mol Cis barrier HOOH Part of such 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° φe ~ 111° 1.19 kcal/mol Trans barrier HSSH: φe ~ 92° 5.02 kcal/mol trans barrier 7.54 kcal/mol cis barrier

  17. 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 kn and dnfor each n.

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

  19. 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 3 kcal/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.

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

  21. Twisting potential surface for ethene The ground state (N) of ethene prefers φ=0º to obtain the highest overlap of the pp orbitals on each C The rotational barrier is 65 kcal/mol. 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 The triplet excited state (T) prefers φ =90º to obtain the lowest overlap. N φ = 90° φ = 0° φ = -90°

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

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

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

  25. Summary: Valence Force Field Terms example Description Typical Expressions Harmonic Stretch Bond Stretch Harmonic-cosine Angle bend dihedral Torsion E(φ) = ½ B [1 – cos(3 φ)], Umbrella inversion Inversion

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

  27. 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 not the idealized charges of +1 and -1 We need a method to estimate such charges in order to calculate properties of materials. A side note 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-10esu A Where m(D) = 2.5418 m(au)

  28. QEq Electrostatics energy 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

  29. 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 The 2nd order expansion in Q is ok Get minimum at Q=-0.887 Emin = -3.676 = 8.291 = 9.352

  30. QEq: the optimum charges lead to Equilibrationthat is, equal chemical potential Expand the energy of an atom in a power series of 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)

  31. QEq parameters

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

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

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

  35. 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 shielded unshielded The plot shows the shielding for C atoms using various Slater orbitals Using RC=0.759a0 And l = 0.5

  36. QEq results for alkali halides

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

  38. QEq for deoxy adenosine Amber charges in parentheses

  39. QEq for polymers Nylon 66 PEEK

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

  41. Noble gas dimers 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 Ar2 s Re De Most popular form: Lennard-Jones 12-6 E=A/R12 –B/R6 = De[r-12 – 2r-6] where r= R/Re Note that E is a min at Re, note the 2 • A 2nd form for LJ 12-6 is • = 4 De[t-12 – t-6] where t=R/s Note that E=0 at R = s, note the 1 Here s = Re(1/2)1/6 =0.89 Re

  42. van der Waals interaction • What do vdW interactions account for? • Short range: Pauli Repulsion of overlapping electron pairs Pauli repulsion: Born repulsion: • Long range attraction  London dispersion due to Instantaneous fluctuations in the QM charges Dipole-dipole Dipole-quadrupole Quadruole-quadrupole We usually neglect higher order (>R6) terms.

  43. 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 Here Dvdw=well depth, ndRvdw= Equilibrium distance for vdwdimer Alternative form E(R) = 4Dvdw{[svdw/R)12]- [svdw/R)6]} Where s = 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.

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

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

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

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

  48. 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 so that we do not want to include them again in the nonbond list. Thus normally we exclude from the vdw and coulomb sums the contributions from bonds (1-2) and next nearest neighbor (1-3) interactions

  49. 1-4 Nonbond exclusions There is disagreement about 1-4 interactions. In fact the origin 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 as assumed in the vdw. 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