Colby College Molecular Mechanics Tutorial September 2001 Thomas W. Shattuck Department of Chemistry Colby College Waterville, Maine 04901 Please, feel free to use this tutorial in any way you wish , provided that you acknowledge the source and you notify us of your usage. Please notify us by email at twshattu@colby.edu or at the above address. This material is supplied as is, with no guarantee of correctness. If you find any errors, please send us a note. Table of Contents Introduction to Molecular Mechanics1: QUANTA/CHARMm 2: Conformational Preference of Methylcyclohexane 3: Geometry (or How Does Molecular Mechanics Measure Up?) 4: Building More Complex Structures: 1MethyltransDecalin 5: Conformational Preference for Butane 6: Working With MM2 from QUANTA 7: Comparing Structures 8: Plotting Structures 9: Conformational Preference of Small Peptides 10: Molecular Dynamics 11: Dynamics in Small Peptides 12: Free Energy Perturbation Theory and Henry's Law Constants and Gibb's Free Energy of Solvation 13. Protein Structure and GramicidinS 14. Solvation and Cyclodextrin 15. Rulers in QUANTA and Cyclodextrin 16. Docking: Cyclodextrin and Naphthol Introduction to Molecular Mechanics Summary The goal of molecular mechanics is to predict the detailed structure and physical properties of molecules. Examples of physical properties that can be calculated include enthalpies of formation, entropies, dipole moments, and strain energies. Molecular mechanics calculates the energy of a molecule and then adjusts the energy through changes in bond lengths and angles to obtain the minimum energy structure. Steric Energy A molecule can possess different kinds of energy such as bond and thermal energy. Molecular mechanics calculates the steric energy of a moleculethe energy due to the geometry or conformation of a molecule. Energy is minimized in nature, and the conformation of a molecule that is favored is the lowest energy conformation. Knowledge of the conformation of a molecule is important because the structure of a molecule often has a great effect on its reactivity. The effect of structure on reactivity is important for large molecules like proteins. Studies of the conformation of proteins are difficult and therefore interesting, because their size makes many different conformations possible. Molecular mechanics assumes the steric energy of a molecule to arise from a few, specific interactions within a molecule. These interactions include the stretching or compressing of bonds beyond their equilibrium lengths and angles, the Van der Waals attractions or repulsions of atoms that come close together, the electrostatic interactions between partial charges in a molecule due to polar bonds, and torsional effects of twisting about single bonds. To quantify the contribution of each, these interactions can be modeled by a potential function that gives the energy of the interaction as a function of distance, angle, or charge. The total steric energy of a molecule can be written as a sum of the energies of the interactions: E_{steric energy} = E_{str} + E_{bend} + E_{strbend} + E_{VdW} + E_{tor} + E_{qq} (1) The bond stretching, bending, and stretchbend interactions are called bonded interactions because the atoms involved must be directly bonded or bonded to a common atom. The Van der Waals, torsional, and electrostatic (qq) interactions are between nonbonded atoms. Bonded Interactions E_{str} represents the energy required to stretch or compress a bond between two atoms, Figure 1. Figure 1. Bond Stretching A bond can be thought of as a spring having its own equilibrium length, r_{o}, and the energy required to stretch or compress it can be approximated by the Hookian potential for an ideal spring: E_{str} = 1/2 k_{s,ij} ( r_{ij}  r_{o} )^{2} (2) where k_{s,ij} is the stretching force constant for the bond and r_{ij }is the distance between the two atoms, Figure 1. E_{bend} is the energy required to bend a bond from its equilibrium angle, _{o}. Again this system can be modeled by a spring, and the energy is given by the Hookian potential with respect to angle: E_{bend} = 1/2 k_{b,ij} ( _{ij}  _{} )^{2} (3) where k_{b,ij} is the bending force constant and _{ij} is the instantaneous bond angle (Figure 2). Figure 2. Bond Bending E_{strbend}_{ }is the stretchbend interaction energy that takes into account the observation that when a bond is bent, the two associated bond lengths increase (Figure 3). The potential function that can model this interaction is: E_{strbend} = 1/2 k_{sb,ijk} ( r_{ij}  r_{o} ) (_{ik}  _{o} ) (4) where k_{sb,ijk} is the stretchbend force constant for the bond between atoms i and jwith the bend between atoms i, j, and k. Figure 3. StretchBend Interaction Therefore, when intramolecular interactions stretch, compress, or bend a bond from its equilibrium length and angle, it resists these changes with an energy given by the above equations. When the bonds cannot relax back to their equilibrium positions, this energy raises the steric energy of the entire molecule. Nonbonded Interactions Van der Waals interactions, which are responsible for the liquefaction of nonpolar gases like O_{2} and N_{2}, also govern the energy of interaction of nonbonded atoms within a molecule. These interactions contribute to the steric interactions in molecules and are often the most important factors in determining the overall molecular conformation (shape). Such interactions are extremely important in determining the threedimensional structure of many biomolecules, especially proteins. A plot of the Van der Waals energy as a function of distance between two hydrogen atoms is shown in Figure 4. When two atoms are far apart, an attraction is felt. When two atoms are very close together, a strong repulsion is present. Although both attractive and repulsive forces exist, the repulsions are often the most important for determining the shapes of molecules. A measure of the size of an atom is its Van der Waals radius. The distance that gives the lowest, most favorable energy of interaction between two atoms is the sum of their Van der Waals radii. The lowest point on the curve in Figure 4 is this point. Interactions of two nuclei separated by more than the minimum energy distance are governed by the attractive forces between the atoms. At distances smaller than the minimum energy distance, repulsions dominate the interaction. The formula for the Van der Waals energy is: E_{VdW,ij} =  + (5) where A and B are constants dependent upon the identities of the two atoms involved and r_{ij} is the distance, in Angstroms, separating the two nuclei. This equation is also called the LennardJones potential. Since, by definition, lower energy is more favorable, the  A/r^{6} part is the attractive part and the + B/r^{12} part is the repulsive part of the interaction. For two hydrogen atoms in a molecule: A = 70.38 kcal Å^{6} B = 6286. kcal Å^{12} Figure 4: Van der Waals interactions between two hydrogen atoms in a molecule, such as H_{2}O_{2} or CH_{3}CH_{3} Torsional Interactions: E_{tor} is the energy of torsion needed to rotate about bonds. Torsional energies are usually important only for single bonds because double and triple bonds are too rigid to permit rotation. Torsional interactions are modeled by the potential: E_{tor} = 1/2 k_{tor,1} (1  cos ) +1/2 k_{tor,2} (1  cos 2 ) + 1/2 k_{tor,3} ( 1  cos 3 ) (6) The angle is the dihedral angle about the bond. The term in 3is important for sp^{3} hybridized systems ( Figure 5a and b ) and the term in just is useful for the central bond in molecules such as butane that have CCCC frameworks(Figure 5c). The constants k_{tor,1, }k_{tor,2} and k_{tor,3} are the torsional constants for onefold, twofold and threefold rotational barriers, respectively. a. b. c. Figure 5. Torsional Interactions, (a) dihedral angle in sp^{3}systems. (b) threefold, 3rotational energy barrier in ethane. (c) butane, which also has a contribution of a one fold, barrier. The origin of the torsional interaction is not well understood. Torsional energies are rationalized by some authors as a repulsion between the bonds of groups attached to a central, rotating bond ( i.e., CCCC frameworks). Torsional terms were originally used as a fudge factor to correct for the other energy terms when they did not accurately predict steric energies for bond twisting. For example, the interactions of the methyl groups and hydrogens on the "front" and "back" carbons in butane were thought to be Van der Waals in nature (Figure 5). However, the Van der Waals function alone gives an inaccurate value for the steric energy. ^ If bonds in the molecule are polar, partial electrostatic charges will reside on the atoms. The electrostatic interactions are represented with a Coulombic potential function: E_{qq,ij } = (7) The Q_{i} and Q_{j} are the partial atomic charges for atoms i and j separated by a distance r_{ij}. _{}is the molecular dielectric constant (normally 1.0), which approximates the dielectric effect of intervening solute or solvent atoms. k is a units conversion constant; for kcal/mol, k=2086.4. Like charges will raise the steric energy, while opposite charges will lower the energy. The Del Re method is often used for estimating partial charges. The Coulomb potential for a unit positive and negative charge is shown in Figure 6a and the Coulomb potential for the hydrogens in H_{2}O_{2} is shown in Figure 6b. F igure 6. (a) Coulomb attraction of a positive and a negative charge. (b) Coulomb repulsion of the two hydrogens in H_{2}O_{2}, with the charge on each hydrogen as Q_{1 }= Q_{2 }= 0.210. The bond stretching, bond bending, stretchbend, Van der Waals, torsion, and electrostatic interactions are said to make up a force field. Each interaction causes a steric force that the molecule must adjust to. The above potential functions represent the various nonbonded interactions that can occur between two atoms i and j. A full force field determines the steric energy by summing these potentials over all pairs of atoms in the molecule. Empirical Force Fields All the potential functions above involve some force constant or interaction constant. Theoretically, these constants should be available from quantum mechanical calculations. In practice, however, it is necessary to derive them empirically. That is, the constants are adjusted so that the detailed geometry is properly predicted for a number of well known compounds. These constants are then used to calculate the structures of new compounds. The accuracy of these constants is critical to molecular mechanics calculations. Unfortunately, no single best set of force constants is available because of the diversity of types of compounds. For example, the MM2 force field works best on hydrocarbons because most of the known compounds used in deriving the force field were hydrocarbons. MM2 is less accurate for oxygencontaining compounds and even less reliable for nitrogen and sulfur species. This is because there aren't as many heteroatom containing compounds in the data base for MM2 and hydrocarbons are a more homogeneous class of compounds than substances with heteroatoms. However, the MM2 force field is one of the best available and the most widely accepted force field for use with organic compounds. MM2 was specifically parameterized to reproduce experimental enthalpies of formation. It is important to realize that the force field is not absolute, in that not all the interactions listed in Equation 1 may be necessary to accurately predict the steric energy of a molecule. On the other hand, many force fields use additional terms. For example, MM2 adds terms to the bonded interactions to better approximate the real potential function of a chemical bond. These additional terms take into account anharmonicity, which is a result of the fact that given enough vibrational energy, bonds will break. Purely quadratic potentials have steep "walls" that prevent bond dissociation (Figure 7a). Cubic terms are added to Equation 2 to adjust for this: E_{str} = 1/2 k_{s,ij} ( r_{ij}  r_{o} )^{2} + 1/2 k_{s,ij} _{ij} ( r_{ij}  r_{o} )^{3} (8) where _{ij} is the anharmonicity constant. For example, for a C(sp^{3})C(sp^{3}) bond the anharmonicity is 2, see Figure 7b: E_{str} = [4.40 mdynes/Å] ( r  1.532Å )^{2} + [4.40 mdynes/Å] [2.00] ( r  1.532Å )^{3} (9) The addition of the cubic term makes the small r portion steeper or more repulsive. This is realistic for real bonds. At larger r the curve is less steep, as desired. For r very large (r> 3Å) the energy decreases, which is unphysical; the curve should approach a constant value. Even though the large r behavior is incorrect, the bond length in compounds remains less than this value, so this region is unimportant under normal conditions.
Figure 7. (a). Energy for the stretching of a CC bond with only the (rr_{o})^{2} harmonic term., Eq. 2 (b), Comparison of the harmonic term with Eq. 8, which includes the (rr_{o})^{3} term for anharmonicity. Enthalpy of Formation The steric energy of a molecule can be used to calculate the enthalpy of formation. First, a bond energy calculation is done using standard tabular values. It is customary to use bond increments rather than the bond energy calculations that you did in General Chemistry. However, the principle is the same. Thermal energy terms must then be added to account for the energy of translation and rotation of the molecule. The energy of translation (x, y, z motion of the center of mass of the molecule) is 3/2RT. The rotational energy of a nonlinear molecule is also 3/2RT ( 1/2RT for each rotational axis). The steric energy calculation in molecular mechanics corresponds to an internal energy calculation. Since H=U+(PV), and PV=nRT for an ideal gas, we must also add RT to convert from internal energy to enthalpy. We have not yet considered molecular vibrations, especially internal rotations. In principle, every vibration, including internal rotations, contributes to the enthalpy. However, the contribution of vibrations is difficult to calculate. In practice the contributions are often small so they can be ignored. However, the internal rotation of the methyl group is always included; in fact the effect is included in the bond increment calculation. Extra terms must also be added for nonmethyl free internal rotations. This contribution, which is called the torsional increment, is estimated as 0.36 kcal/mol or 1.51 kJ mol^{1 }for each internal rotation. For example, butane, CH_{3}CH_{2}CH_{2}CH_{3} , has one additional internal rotation, other than the methyl group rotations; so the torsional increment for butane would be 0.36 kcal/mol. In summary the enthalpy of formation is then, _{f}H° = 3/2RT + 3/2RT + RT + bond energy + steric energy + torsional increments (10) This formula also assumes that there is only one low energy conformation of the molecule. If there are several low energy conformations, each must be accounted for in Equation 10. Bond Energy You are familiar with bond energy calculations from General Chemistry. The energy of a molecule is assumed to be an additive function of the energy of individual bonds (Table I). The _{r}H for a reaction is given from H°(bonds broken) H°(bonds formed). Table I. Bond Enthalpies, H°(AB) (kJ/mol)
C (graph) > C (g) H°= 716.7 kJ/mol For example, the enthalpy of formation of acetaldehyde is calculated as: 2 C(graph) + 2 H_{2} (g) + 1/2 O_{2} (g) > CH_{3}CH=O (g) # Bonds Broken  # Bonds Formed 2 C (graph) 2 (716.7 kJ/mol) 1 C=O 743 kJ/mol 2 HH 2 (436 kJ/mol) 4 CH 4 (412 kJ/mol) 1/2 O=O 1/2 (497 kJ/mol) 1 CC 348 kJ/mol total 2553.9 kJ/mol  total 2739 kJ/mol = 185.1 kJ The experimental value is 166.19 kJ, so the value derived from Table I is not very accurate. The bond energy calculations in molecular mechanics are done slightly differently, using bond increments. Again the bond energies are assumed to be additive. The contributions are taken not only from each bond, but increments are added for certain structures, such as tertiary carbon linkages. The bond energy calculation for acetaldehyde from the MMX program is given below, with energies in kcal. MMX also calculates entropies, which are also listed for your interest. # Bond or Structure Each Total Tot S contrib. ^ 1 C=O 25.00 25.00 2.300 1 CH ALDEHYDE 2.500 2.500 26.800 1 CC SP3SP2 C=O 3.000 3.000 0.600 1 MECARBONYL 2.000 2.000 ______ bond energy = 42.115 kcal S° = 62.600 cal/K The bond energy is 42.115 kcal or 176.2 kJ. However, caution should be used since these calculations are designed to be used in conjunction with steric energies in a molecular mechanics calculation and not as general bond energy values. Using Equation 10, with the steric energy calculated by molecular mechanics gives the final _{f}H° = 169.33 kJ/mol, which is a significant improvement over the bond energy calculation from Table I of 185.1 kJ.
