Techniques of Molecular Modelling: Molecular Mechanics

As noted before, the origins of the molecular mechanics method can be traced to work done by eg Barton (see reprints on van der Waals, electrostatic or the joined up analysis), Bright-Wilson (Molecular Vibrations, 1938. DOI: 10.1063/1.1750232) and others in the late 1940s and early 1950s and e.g. Allinger (starting in 1959, DOI: 10.1021/ja01530a049) and the MM2/3/4 force fields and Kollman/Amber in the 1970s-80s.

1. The Mechanics Force Field

2. Using MM to optimise a molecular geometry

3. Summary of Pros and Cons of Molecular Mechanics


  1. Simple energy functions are quick to solve,
  2. can cope with large molecules,
  3. accurate for systems which are close to the models used to reproduce the force field (interpolative),
  4. The minima found are restricted to the "atom-connectivity" specified in the input file; no bonds can make/break during the optimisation process (and hence no matter how bad the initial geometry guess is, that connectivity will be preserved). Hence very useful for instantly converting sketches of molecules into tidied coordinates, even in an on-the-fly manner (elastic band optimisers)
  5. deals reasonably accurately with van der Waals interactions at non-bonded distances for larger molecules (which increase vastly in number as the molecule gets bigger and which are notoriously difficult to model using purely quantum mechanics)
  6. Mechanics, when used with high quality parameter sets, can be very reliable for calculating vibrational properties (IR frequencies, but not intensities which are an electronic property). See also the cons
  7. Mechanics is good for conformational exploration (Monte-Carlo, Genetic algorithms, etc) where the number of possible conformations can number millions. It is used in this context for NMR predictions of molecules where many possible conformations may exist.
  8. Mechanics is highly applicable to molecular dynamics, and hence derived statistical thermodynamic quantities such as free energies etc. In this context it is frequently used for analysing proteins, and other large biomolecules.


  1. No unique functional form for the energy terms (= proprietary?)
  2. Number of parameters rises rapidly for non-CHNO elements,
  3. Need to know atom/bond type for all the molecule (non-classical bonding)
  4. Does not include majority of cross terms in force field
  5. Most mechanics implementations do not compute 2nd derivatives (entropy)
  6. Cannot cope with unusual (= unpredicted) bonding or interactions
  7. Mechanics is prone to "false minima", allowed by the equations, but physically unrealistic (e.g. hemispherical carbon)
  8. Cannot be used for bond-breaking transition states and hence reactions.

Case Study 5: The Relative Stability of Carbonium Ions

Organic chemists are taught that carbonium ion stability is predominantly due to electronic factors such as whether they are primary, secondary, tertiary, allylic, benzylic etc. What is mentioned much less is that their relative stability is also very sensitive to their geometries, and in particular the angles of the three substituents at the carbon. Geometries of carbonium ions are particularly inaccessible; it is not easy at all to get crystal structures! So how can Molecular Mechanics model the following carbonium ions? They are arranged in the expected order of increasing stability. Notice that the MM2 method gradually decreases the predicted angle energies. Notice also however how it predicts that the tertiary carbonium ion is actually LESS stable than the preceeding primary ion, which is surprising to say the least. Clearly, electronic factors must be playing an important role as well. Suffice to say that even without ANY consideration of electronic factors, the MM2 method does not do too badly in its prediction of relative stabilities. This re-enforces the conclusion that bond angles in carbonium ions are just as important as substitution! See also DOI: 10.1021/ed800058c for more detail.

stretch       =   0.868  angle         = 19.173
stretch bend  =  -0.672  dihedral      =  6.497
improp torsion=   0.155  van der Waals =  2.663
electrostatics=  -1.341  hydrogen bond =  0.000
Energy of the final structure is 27.342 kcal/mol.
(QM B3LYP/cc-pVTZ → rearranges to final
product directly by breaking/forming bonds)
stretch       =   0.447 angle         =  18.758
stretch bend  =  -0.078  dihedral      =   6.916
improp torsion=   1.124  van der Waals =   1.622
electrostatics=   0.000  hydrogen bond =   0.000
Energy of the final structure is 28.788 kcal/mol.
(Relative QM B3LYP/cc-pVTZ free energy
= 25.7 Kcal/mol)
stretch       =   0.571  angle         =  5.694
stretch bend  =  -0.027  dihedral      =  7.179
improp torsion=   0.009  van der Waals =  5.678
electrostatics=  -0.545  hydrogen bond =  0.000
Energy of the final structure is 18.559 kcal/mol.
(Relative QM B3LYP/cc-pVTZ free energy
= 23.3 Kcal/mol)
stretch       =   0.512 angle         =  3.439
stretch bend  =  -0.153 dihedral      =  6.662
improp torsion=   0.002 van der Waals =  3.872
electrostatics=   0.683 hydrogen bond =  0.000
Energy of the final structure is 15.017 kcal/mol.
(Relative QM B3LYP/cc-pVTZ free energy
= 0 Kcal/mol)

Case Study 3 revisited: π-facial Hydrogen Bonding

The H-bond energy functions in the Mechanics method are normally specifically parameterised against standard models, which involve three atoms, ie X...H...Y. The model of hydrogen bonding in this molecule is unusual to say the least, since it involves a whole ring and not a single atom as the H acceptor, and hence the Mechanics method (in theory) has no functional form which can reproduce it. In particular, aromatic rings are attractive to electrophiles by virtue of their quadrupole moment, which is a higher order property of the charge distribution (thus the so-called multipole expansion goes as monopole charge, dipole, quadrupole, octapole and hexadecapole. The equations above in effect only include the first and possibly second term of this expansion).

The optimised geometry does predict what looks like H-bonding, but its directionality is not quite correct. The steric energy is -49 compared with -18 if the two components are infinitely separated. This makes the binding energy = -31 kcal/mol. In this case, it is predicted to arise almost entirely from the non-bonded van der Waals term rather than electrostatic or hydrogen bonding (dipole-dipole) terms.

One monomer Dimer
Stretch: 1.3257
Bend: 2.8090
Stretch-Bend: 0.2180
Torsion: -24.8762
Non-1,4 VDW: -3.1495
1,4 VDW: 10.7361
Dipole/Dipole: 4.1782
Total Energy: -8.7587 kcal/mol
Stretch:          2.5669
Bend:             6.6539
Stretch-Bend:     0.3319
Torsion:        -54.0415
Non-1,4 VDW:    -29.6591
1,4 VDW:         15.6989
Dipole/Dipole:    9.4549
Total Energy:   -48.9942 kcal/mol

There is one final aspect which must be considered, namely the entropy of the interaction. The equation

ΔG = ΔH - T.ΔS

implies that if entropy decreases (i.e. ΔS is -ve) as a result of a reaction (or in this case a binding, in which two free molecules are reduced to one bound complex), then ΔG must increase (two -ve signs = a positive result) as a result by the amount T.ΔS. The Mechanics method could in principle provide an estimate of this term via calculation of the normal modes of vibration (and solution of the appropriate partition function equations), but few Mechanics programs actually perform this calculation (which requires the second derivatives to be evaluated, something only available using MM4, accurate to between 25-50 cm-1). An approximate estimate of T.ΔS is about -(12-15) kcal/mol at 298K (i.e. approx 2 kcal/mol per degree of lost freedom), which has to be added (-/- = +) to the binding energy obtained above, to obtain a final value of -(15) kcal/mol. This value is certainly in the correct ballpark.

Another property of this molecule can also be studied using molecular mechanics, namely the rotation and preferred orientation of the aryl-C bond, for which mechanics is very well suited. The following "dihedral driving" calculation shows the MM2 energy as a function of the dihedral angle between the trifluoromethyl group and the anthracene ring. The barrier (about 10-11 kcal/mol) compares reasonably well with the measured value of ΔG (from NMR lineshape analysis) of about 14 kcal/mol. The discrepancy may well be due to ignoring solvation of the OH group, which increases its size, and hence the barrier.

Reference: M. L. Webb and H. S. Rzepa, Chirality, 1994, 6, 245-250. DOI: 10.1002/chir.530060406

Fast forward to next installment of this story

Case Study 6: Molecular Anvils and Tennis Balls

Now that we have a quantitative method for estimating energies of molecules, we can tackle a problem such as shown below. An application of "shape design" is the recent area of designing catalysts to promote reactions by ensnaring them in a cavity, or binding pocket. Whilst designing enzymes to do this is somewhat beyond the scope of the current course, a simple example of small molecule design can be shown here. A resorcinol-anthracene based aromatic molecule was designed in the form of spacers (the aromatic rings) and links (in the form of hydrogen-bonds) to create in the solid state a system with well defined cavities, which are large enough to accommodate small organic molecules, such as ethyl methacrylate and cyclohexadiene. The former is held within the cage by further hydrogen bonds to the cavity lattice. In particular, various changes to the structure can be made, to increase or decrease steric bulk, etc and the effects of these changes can be modelled, ie predicted.

The previous systems were modelled by drawing/creating the model using the built-in editor (ChemBio3D, etc). Here we will be importing the X-ray derived coordinates. There are some important considerations to doing so.

When a standard MM2 mechanics force field is applied to the resorcinol-anthracene anvil, the following results are obtained:

Component MM2, kcal/mol
Supermolecule -157.2
Supermolecule without substrates (Anvil) -140.1
Anvil computed as four components -92.8
Three Substrates together, without anvil +14.2
Three separated Substrates +23.4
Anvil + Substrates separately -140.1+14.2=-125.9
Binding Energy of substrates to anvil -157.2-(-125.9)=-31.3

These energies are NOT corrected for entropy (i.e. they are closer (but not quite) enthalpic, and not free energies) a significant limitation. Very likely this sort of catalyst functions by converting enthalpic energies (the binding of the catalyst components) into entropic energies (i.e. pre-organising the reactants in the cavity, which of course consumes entropy. The free energy cost of assembling the three substrates is likely to be around +24-30 kcal/mol, but this is still likely to result in a win!).

What IS useful with the mechanics approach is the ability now to make minor variations to the environment, such as changing the nature of the substrates included in the centre of the anvil to see how the binding energy responds. Its a tuning procedure!

Literature Citation. Catalysis by Organic Solids. Stereoselective Diels-Alder Reactions Promoted by Microporous Molecular Crystals Having an Extensive Hydrogen-Bonded Network, K. Endo, T. Koike, T. Sawaki, O. Hayashida, H. Masuda, and Y. Aoyama, J. Am. Chem. Soc., 1997, 4117 - 4122; DOI: 10.1021/ja964198s

An interesting alternative to the above anvil is the molecular "tennis ball" (e.g. DOI) elaborated by Julius Rebek. Quite large molecules can be enclosed in these balls, and once in, they organize themselves very specifically, a phenomenon Rebek has termed molecular "social isomerism" (DOI: 10.1002/anie.200462839)

Just as with the anvil, the tennis ball is held together with hydrogen bonds, and this can be modelled using the Mechanics method, as can the organisation of the ligands inside the cavity.

One monomer Dimer
  Stretch:          5.1147
  Bend:            49.5243
  Stretch-Bend:     0.4564
  Torsion:         -3.2591
  Non-1,4 VDW:    -20.0160
  1,4 VDW:         21.6595
  Dipole/Dipole:   25.2888
  Total Energy:    78.7685 kcal/mol
  Stretch:          9.8872
  Bend:           103.0995
  Stretch-Bend:     1.2134
  Torsion:         -3.8372
  Non-1,4 VDW:    -87.3993
  1,4 VDW:         43.0785
  Dipole/Dipole:   29.6762
  Total Energy:    95.7182 kcal/mol

Here, the energy liberated upon dimerisation is 62 kcal/mol, of which 41 is due to the non-bonded dispersion terms and 21 to hydrogen bonds. This is more than enough to overcome the entropic terms which acrue when substrates are trapped inside.

Back to scales|Back to Visualisation| Forward to MO Reactant|Forward to MO TS| Forward to MO Advanced|
(c) H. S. Rzepa 1998-2013. No reproduction rights granted to this material without permission.