Colby College Molecular Mechanics Tutorial
Thomas W. Shattuck
Department of Chemistry
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 e-mail at firstname.lastname@example.org
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 Mechanics
2: Conformational Preference of Methylcyclohexane
3: Geometry (or How Does Molecular Mechanics Measure Up?)
4: Building More Complex Structures: 1-Methyl-trans-Decalin
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 Gramicidin-S
14. Solvation and -Cyclodextrin
15. Rulers in QUANTA and -Cyclodextrin
16. Docking: -Cyclodextrin and -Naphthol
Introduction to Molecular Mechanics
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.
A molecule can possess different kinds of energy such as bond and thermal energy. Molecular mechanics calculates the steric energy of a molecule--the 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:
Esteric energy = Estr + Ebend + Estr-bend + EVdW + Etor + Eqq (1)
The bond stretching, bending, and stretch-bend 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 non-bonded atoms.
Estr 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, ro, and the energy required to stretch or compress it can be approximated by the Hookian potential for an ideal spring:
Estr = 1/2 ks,ij ( rij - ro )2 (2)
where ks,ij is the stretching force constant for the bond and rij is the distance between the two atoms, Figure 1.
Ebend 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:
Ebend = 1/2 kb,ij ( ij - )2 (3)
where kb,ij is the bending force constant and ij is the instantaneous bond angle (Figure 2).
Figure 2. Bond Bending
Estr-bend is the stretch-bend 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:
Estr-bend = 1/2 ksb,ijk ( rij - ro ) (ik - o ) (4)
where ksb,ijk is the stretch-bend force constant for the bond between atoms i and jwith the bend between atoms i, j, and k.
Figure 3. Stretch-Bend 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.
Van der Waals interactions, which are responsible for the liquefaction of non-polar gases like O2 and N2, also govern the energy of interaction of non-bonded 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 three-dimensional 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:
EVdW,ij = - + (5)
where A and B are constants dependent upon the identities of the two atoms involved and rij is the distance, in Angstroms, separating the two nuclei. This equation is also called the Lennard-Jones potential. Since, by definition, lower energy is more favorable, the - A/r6 part is the attractive part and the + B/r12 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 H2O2 or CH3-CH3
Torsional Interactions: Etor 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:
Etor = 1/2 ktor,1 (1 - cos ) +1/2 ktor,2 (1 - cos 2 ) + 1/2 ktor,3 ( 1 - cos 3 ) (6)
The angle is the dihedral angle about the bond. The term in 3is important for sp3 hybridized systems ( Figure 5a and b ) and the term in just is useful for the central bond in molecules such as butane that have C-C-C-C frameworks(Figure 5c). The constants ktor,1, ktor,2 and ktor,3 are the torsional constants for one-fold, two-fold and three-fold rotational barriers, respectively.
a. b. c.
Figure 5. Torsional Interactions, (a) dihedral angle in sp3systems. (b) three-fold, 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., C-C-C-C 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:
Eqq,ij = (7)
The Qi and Qj are the partial atomic charges for atoms i and j separated by a distance rij. 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 H2O2 is shown in Figure 6b.
igure 6. (a) Coulomb attraction of a positive and a negative charge. (b) Coulomb repulsion of the two hydrogens in H2O2, with the charge on each hydrogen as Q1 = Q2 = 0.210.
The bond stretching, bond bending, stretch-bend, 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 non-bonded 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 oxygen-containing compounds and even less reliable for nitrogen and sulfur species. This is because there aren't as many hetero-atom containing compounds in the data base for MM2 and hydrocarbons are a more homogeneous class of compounds than substances with hetero-atoms. 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:
Estr = 1/2 ks,ij ( rij - ro )2 + 1/2 ks,ij ij ( rij - ro )3 (8)
where ij is the anharmonicity constant. For example, for a C(sp3)-C(sp3) bond the anharmonicity is -2, see Figure 7b:
Estr = [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 C-C bond with only the (r-ro)2 harmonic term., Eq. 2
(b), Comparison of the harmonic term with Eq. 8, which includes the (r-ro)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 non-linear 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 non-methyl 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, CH3-CH2-CH2-CH3 , 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,
fH° = 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.
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 rH for a reaction is given from H°(bonds broken)- H°(bonds formed).
Table I. Bond Enthalpies, H°(A-B) (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 H2 (g) + 1/2 O2 (g) -> CH3-CH=O (g)
# Bonds Broken - # Bonds Formed
2 C (graph) 2 (716.7 kJ/mol) 1 C=O 743 kJ/mol
2 H-H 2 (436 kJ/mol) 4 C-H 4 (412 kJ/mol)
1/2 O=O 1/2 (497 kJ/mol) 1 C-C 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 C-H ALDEHYDE -2.500 -2.500 26.800
1 C-C SP3-SP2 C=O -3.000 -3.000 -0.600
1 ME-CARBONYL -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 fH° = -169.33 kJ/mol, which is a significant improvement over the bond energy calculation from Table I of -185.1 kJ.