Unformatted text preview: Molecular Mechanics Overview Molecular Mechanics Molecular mechanics can be used to study molecules that are too large for quantum mechanical models. Molecules are treated as simple mechanical systems governed by Newtonian mechanics. The potential energy of the molecule is specified as a function of the nuclear coordinates. This function depends on parameters determined using data or ab initio calculations (empirical force fields). The distribution of electrons is usually not explicitly modeled. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 1 / 39 Molecular Mechanics Overview Applications Molecular mechanics focuses on static aspects of molecular structure: Find the most stable/lowest energy structure. Perform Monte Carlo simulations to sample structures according to the Boltzman distribution: P(X) = e V (X)/kB T . e V (Y)/kB T dY However, the MM potential energy function can also be used to construct the force fields that are used in molecular dynamics simulations: m d 2X = F (X) = V (X). dt 2 Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 2 / 39 Molecular Mechanics Overview Assumptions To be useful, the potential energy function should satisfy the following conditions: Additivity: The total energy can be written as a sum of potentials arising from interactions between small numbers of atoms or functional groups.
Allatom fields treat all atoms explicitly. Unified atom fields depend only on the relative positions of functional groups (e.g., bases, amino acids). Transferability: Each term in this sum can be parameterized using data obtained from representative structures. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 3 / 39 Molecular Mechanics Overview Potential Energy The potential energy is usually expressed as a sum of bonded and nonbonded contributions V (R) = rb Vbond (rb ) + + where i=j Vangle () + VCoul (rij ) VvdW (rij ) + i=j Vtor ( ) rb , b,b and are bond lengths, bond angles and dihedral angles; rij is the distance between atoms i and j. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 4 / 39 Molecular Mechanics Bonded interactions Bonded Interactions There are usually three or four contributions to the local or bonded energy: bond length; bond angle (between bonds from a common atom); torsion/dihedral angle (along a series of four bonded atoms); improper torsion angles/outofplane bending. The number of terms in the local potential energy of a molecule containing N atoms is O(N). Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 5 / 39 Molecular Mechanics Bonded interactions Bond Length Potential Small changes in bond length r are usually modeled by a harmonic potential 1 Vbond (r ) = kH (r  r0 )2 2 where r0 is the equilibrium bond length and kH is a constant that depends on the bond type and the atoms. This is adequate only for deviations of up to 0.1 . A Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 6 / 39 Molecular Mechanics Bonded interactions Hooke's Law The harmonic potential is equivalent to Hooke's law d 2r = F (r ) = V (r ) = kH (r  r0 ) dt 2 where k is the spring constant and m is the reduced mass of the two atoms. This solution is just the simple harmonic oscillator m r (t) = r0 + A sin(t + ) where A, and = kH /m are the amplitude, the phase, and the frequency of the oscillations. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 7 / 39 Molecular Mechanics Bonded interactions Parametrization The spring constant kH can be estimated from the bond stretching frequencies measured using IR spectroscopy. Molecules absorb photons at frequencies that match internal vibrations. Spectroscopic frequencies are usually reported as wavenumbers (), which can be converted to Hz using: =c Bond HO HN HC CN CC Frequency (cm1 ) 36003700 34003500 29003000 1250 1000 Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 8 / 39 Molecular Mechanics Bonded interactions Morse Potential The Morse potential can be used to model larger changes in bond length: 2 VMorse (r ) = D 1  e kM (r r0 ) 7 4 2 2 3 3 4 D kM (r  r0 )  kM (r  r0 ) + kM (r  r0 ) 12 D and kM are adjustable parameters. Spectroscopic data can be used to estimate kM by setting: kM = kH /D. Taylor series approximations can be used to speed calculations. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 9 / 39 Molecular Mechanics Bonded interactions Morse Potential The Morse potential reflects the fact the potential energy should level off at large distances as the bond breaks. This asymptotic value is called the dissociation energy and can be measured experimentally. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 10 / 39 Molecular Mechanics Bonded interactions Bond Angle Potential Bond angle potentials are often taken to be harmonic functions of either the bond angle or of the bond angle cosine: Vangle () = kh (  0 )2 2 ~ Vangle () = kt cos()  cos(0 ) where kh , kt , 0 are constants that depend on the trio of bonded atoms. The simple harmonic potential is cheaper to compute but is unbounded. The trigonometric potential can be expanded as a power series to speed calculations. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 11 / 39 Molecular Mechanics Bonded interactions Bond Angle Parameters The stiffness coefficients in the bond angle potential can be estimated from the bending vibrational frequencies measured with IR spectroscopy. Frequency (cm1 ) 1600 1600 1500 500 500 Equilibrium angles can be estimated using Solved structures. Valence shell electron pair repulsion (VSEPR) theory. Bonds HOH HNH HCH OC=O CC=O Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 12 / 39 Molecular Mechanics Bonded interactions Torsional Potential Recall that the dihedral angle of a chain i  j  k  l of atoms is the angle between the planes i  j  k and j  k  l. ijkl ijk n = cos1 ( ijk jkl ) n n jk ij =  jk ij Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 13 / 39 Molecular Mechanics Bonded interactions Torsional Potential Because the torsion potential may have several local minima and maxima, it is often expanded as a trigonometric series Vtor ( ) =
n k=1 Vk 1 + cos(k ) where n is usually 2 or 3 and the Vk are adjustable parameters. n determines the number of minima/maxima. Multiple terms are required to allow for asymmetric maxima/minima. The Vk 's can be estimated using energy differences between maxima and minima. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 14 / 39 Molecular Mechanics Bonded interactions Example: Torsional Potential of Butane syn: = 0 gauche: = 60 eclipsed: = 120 anti: = 180 Here, is the dihedral angle along the four carbons. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 15 / 39 Molecular Mechanics Bonded interactions Force field parameters are not independent. Suppose that Vtor = Then syn/anti anti/gauche V3 V2 1 + cos(2 ) + 1 + cos(3 ) . 2 2 = Vtot ( = 0 )  Vtot ( = 180 ) = V3 + [Vtot  Vtor ]( = 0 )  [Vtot  Vtor ]( = 180 ) = Vtot ( = 180 )  Vtot ( = 60 ) 3 = V2 + [Vtot  Vtor ]( = 180 )  [Vtot  Vtor ]( = 60 ) 4 from which V2 and V3 can be calculated. Notice, however, that these parameters will depend on the remainder Vtot  Vtor .
Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 16 / 39 Molecular Mechanics Bonded interactions Improper Torsion Outofplane potentials are included for some configurations. If atoms i, l, k are bonded to j, then the improper Wilson angle is the angle between bond j  l and the plane i  j  k. The potential energy is usually modeled as a harmonic function 1 Vimp () = kimp (  0 )2 , 2 where taking 0 = 0 favors planarity of the groups attached to j. This can be used to maintain planarity about the C atom in an amino acid; the glycosyl nitrogen in nucleic acid bases. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 17 / 39 Molecular Mechanics Bonded interactions Cross Terms Crossterms are sometimes added to couple the distortions in different internal coordinates. Three examples include V (r , r ) = krr (r  r0 )(r  r0 ) (stretchstretch); (bendbend); (stretchbend). V (, ) = k (  0 )(  0 ) V (r , ) = kr (r  r0 )(  0 ) These can improve the accuracy of the potential energy function, but are used sparingly to limit computational costs. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 18 / 39 Molecular Mechanics Nonbonded interactions Nonbonded Interactions There are usually two contributions to the nonlocal energy: pairwise van der Waals forces; pairwise Coulombic forces. The number of terms in the nonlocal potential energy of a molecule containing N atoms is O(N 2 ). Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 19 / 39 Molecular Mechanics Nonbonded interactions van der Waals Forces van der Waals forces arise from quantum mechanical interactions between the electron distributions about nonbonded atoms. vdW forces are repulsive at shortrange because of the Pauli exclusion principle (exchange repulsion). vdW forces are attractive at longrange because fluctuations in electron density near any one atom induce dipoles in neighboring atoms (London disperson forces). While van der Waals forces are usually weaker than electrostatic forces between charged particles, they are important in nonpolar molecules and increase in importance with molecular size. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 20 / 39 Molecular Mechanics Nonbonded interactions The LennardJones Potential van der Waals forces are often modeled by the LennardJones potential: VLJ (rij ) =  Aij Bij + 12 6 rij rij where the constants Aij and Bij determine the strength of the attractive and repulsive forces, respectively.
6 The rij term can be derived using quantum mechanics. 12 The rij term has been chosen for computational convenience. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 21 / 39 Molecular Mechanics Nonbonded interactions Parametrization of the LJ Potential
0 The parameters Aij and Bij can be calculated from the energy Vij and 0 radius rij at which the LJ potential has its minimum: Aij Bij
0 rij 0 Vij 0 0 = 2(rij )6 Vij 0 0 = (rij )12 Vij 1 0 0 = r11 + r22 2 0 0 1/2 = Vii Vjj 0 0 The quantities Vii and rii pertain to a monoatomic substance and are estimated using Xray crystallography or QM calculations. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 22 / 39 Molecular Mechanics Nonbonded interactions Alternatives to the LJ potential Because the LJ potential overestimates the repulsion at short distances, the vdW potential is sometimes calculated using other functions: Aij + Bij e rij 6 rij VBuck (rij ) =  V7/14 (rij ) =  (Buckingham potential) Aij Bij + 0 0 7 0 (rij + 1 rij )7 (rij + 1 rij )7 [rij + 2 (rij )7 ] (buffered 7/14 attraction/repulsion) Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 23 / 39 Molecular Mechanics Nonbonded interactions Electrostatic Interactions Electrostatic interactions between nonbonded atoms are modeled using the Coulomb potential: Vcoul (rij ) = qi qj rij where qi is the partial atomic charge of atom i and is the dielectric constant for a vacuum. The Coulomb potential falls off slowly with distance. This approach assumes that the charge distribution in the molecule can be modeled as a sum of point charges centered on the nuclei. The partial atomic charges can be estimated using QM methods. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 24 / 39 Molecular Mechanics Nonbonded interactions Screened Coulomb Potentials If the molecule is not immersed in a vacuum, then the Coulomb potential should be modified to take into account screening of charges: VCoul (rij ) = qi qj . (rij )rij Here (r ) is an increasing function of distance. Two standard choices include: (r ) = D0 e r Ds + D 0 (r ) =  D0 . 1 + ke (Ds +D0 )r Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 25 / 39 Molecular Mechanics Nonbonded interactions Some Established Force Fields for MM and MD CHARMM (Chemistry at Harvard Macromolecular Mechanics): developed by Martin Karplus and colleagues. AMBER (Assisted Model Building with Energy Refinement): developed by Peter Kollman and colleagues. GROMOS (Groningen Molecular Simulation System) GROMACS (Groningen Machine for Chemical Simulations): freeware, used in [email protected] Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 26 / 39 Molecular Mechanics Optimization Minimization of the Potential Energy One of the objectives of molecular mechanics is to identify the global minimum of the potential energy function, which should correspond to the most stable conformation of the molecule. Here we are effectively assuming the Thermodynamic hypothesis (C. Anfinsen): Biological macromolecules have a unique kinetically accessible stable structure that minimizes the free energy. This is clearly false in general, as is shown by prions and knotted or linked molecules, but appears to be a good working hypothesis for many proteins and nucleic acids. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 27 / 39 Molecular Mechanics Optimization Multivariate Optimization Several algorithms can be used to search for a local minimum of the potential energy function: Steepest descent algorithms Trust region algorithms Newton's method QuasiNewton methods Conjugate gradient methods TruncatedNewton methods Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 28 / 39 Molecular Mechanics Optimization Steepest Descent Algorithms Algorithm: Choose x0 and then
1 2 3 4 Calculate the descent direction: pk = V (xk ). Set xk+1 = xk + k pk . Determine the step length k using a line search strategy. Test xk+1 for convergence and then return to (1). Often the goal of the line search strategy is to identify a k that gives a local minimum for the onedimensional function: F () = V (xk + pk ). This can be done using Newton's method or by minimizing a polynomial interpolant for F .
Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 29 / 39 Molecular Mechanics Optimization Trust Region Algorithms Algorithm: Choose x0 and R(x0 ) and then
1 Calculate the gradient gk = V (xk ) and the Hessian matrix Hk = ( 2 V (xk )/xi xj ). Find the minimum sk R(xk ) of the quadratic form 1 qk (s) = V (xk ) + gT s + sT Hk s. k 2 2 3 4 Set xk+1 = xk + sk . Test xk+1 for convergence and then return to (1). R(xk ) is the trust region, which is chosen so that the quadratic form q() is a good approximation to V in a neighborhood of xk . Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 30 / 39 Molecular Mechanics Optimization Newton's Method Algorithm: Choose x0 and then
1 Calculate the gradient gk = V (xk ) and the Hessian matrix Hk = ( 2 V (xk )/xi xj ). Calculate pk = H1 gk . k Set xk+1 = xk + pk . Test xk+1 for convergence and then return to (1). 2 3 4 Newton's method converges quadratically in a neighborhood of the stationary point provided H is locally invertible. pk can be found by solving the system Hk pk = gk . Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 31 / 39 Molecular Mechanics Optimization QuasiNewton Methods QuasiNewton methods eliminate the cost of calculating the Hessian matrix by using an approximate Hessian Bk which satisfies the secant equation: V (xk + dx)  V (xk ) = Bk dx. The step vector is then taken to be pk = B1 gk . Bk itself is usually k chosen so that Bk  Bk1  is small. Bk  Bk1 has low rank. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 32 / 39 Molecular Mechanics Optimization Conjugate Gradient Methods The CG method was originally developed to solve linear systems Ax = b, where A is a symmetric, positivedefinite n n matrix. The method is motivated by the following observations: Solving Ax = b is equivalent to minimizing the quadratic form 1 f (x) = xT Ax  bT x. 2 x, yA xT Ay defines an inner product on Rn . If {p1 , , pn } is an orthonormal basis for Rn equipped with , A , then the solution x can be expanded as a linear combination x=
n i=1 c k pk Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 33 / 39 Molecular Mechanics Optimization Conjugate Gradient Methods The linear CG algorithm can be implemented in the following way. Algorithm: Take p0 = r0 = b  Ax0 and then
1 2 3 4 Set pk = rk  k1 pk1 , where k1 = rk , rk /rk1 , rk1 . Set xk+1 = xk + k pk , where k = pk , rk /pk , pk A . Calculate the residual rk+1 = b  Axk+1 . Return to (1) if rk+1  > . In theory, this procedure will terminate in at most n steps, but roundoff error can cause the number of iterations to be larger. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 34 / 39 Molecular Mechanics Optimization Conjugate Gradient Methods for Optimization CG methods can also be used to minimize nonquadratic functions. Algorithm: Take p0 = g0 and then
1 2 3 4 5 6 Calculate the gradient gk = V (xk ). Calculate k (methoddependent: see below). Calculate the new search direction: pk = gk + k pk1 . Use a line search (with k ) to minimize V (xk + k pk ). Set xk+1 = xk + k pk . Test xk+1 for convergence and then return to (1). CG works best when V is approximately quadratic near the minimum. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 35 / 39 Molecular Mechanics Optimization Conjugate Gradient Methods for Optimization The k 's should be chosen so that the step vectors p1 , , pn are conjugate when V is a quadratic function. This can be done in many ways, but the following prescriptions are popular: k k k = gT gk /gT gk1 k k1 = = gT yk /gT gk1 k k1 T gk yk /pT yk k1 (FletcherReeves) (PolakRibi`re) e (HestenesStiefel) where yk = gk  gk1 . It is also common practice to set k = 0 whenever the above prescriptions give a negative value. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 36 / 39 Molecular Mechanics Optimization TruncatedNewton Method TruncatedNewton methods modify Newton's method by using step vectors that only approximately satisfy the Newton condition: rk = Hk pk + gk  k gk  () The forcing sequence k can be chosen so that the method has quadratic convergence. One choice is k = min{c/k, gk }. Step vectors pk satisfying (*) can be found iteratively, e.g., by using the linear conjugate gradient method. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 37 / 39 Molecular Mechanics Optimization Implementations AMBER offers SD, NM, and CG; CHARMM offers SD, NM, CG, and TN; GROMOS offers SD and CG. Since these methods are only guaranteed to find local minima, it is usually good practice to repeat the calculations using different initial values. The eigenvalues of the Hessian at critical points should also be checked to rule out convergence to a saddle point on the potential energy surface. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 38 / 39 Molecular Mechanics References References Field, M. J. (2007) A Practical Introduction to the Simulation of Molecular Systems. Second Edition, Cambridge University Press. Rapp, A. K. and Casewit, C. J. (1997) Molecular Mechanics across e Chemistry. University Science Books. Schlick, T. (2006) Molecular Modeling and Simulation. Springer. Jay Taylor (ASU) APM 530  Lecture 4 Fall 2010 39 / 39 ...
View
Full Document
 Fall '10
 Staff
 Energy, Optimization, Molecule, molecular mechanics, Computational chemistry, JAY TAYLOR

Click to edit the document details