Computational Chemistry/Molecular mechanics

 Previous chapter -  ../

Introduction
A good introduction is molecular mechanics.

In molecular mechanics we treat a group of molecules as a  classical collection of balls and springs rather than a quantum collection of electrons and nuclei. This means we can readily make physical models and have these physical models turned into computer programs. There is a hierarchy of models, the minimal being  atoms as hard spheres of radius  equal to the covalent radius and using VSEPR (Valence Shell Electron Repulsion) for the lonepairs. Angles are approximately determined by best mutual avoidance in the hierarchy lone pairs > bond pairs. The electronegativities of atoms $$\chi$$ can be used to increase the flexibility of the model. Substituents with greater electronegativity decrease the size and repulsiveness of bond pairs. The model is incomplete because of neglect of bonding overlap and there is no description  of any $$\pi$$ bonds the molecule may have. Also we are making invalid additivity assumptions illustrated as follows.

The longest bondlengths are adjacent to electron pairs as in the classic example of F-F. ( HF is 91.7 pico-metres, FF 141.7 (covalent radii H 37.0 pm, F 71.0pm).) HF is quite clearly shorter than average and FF longer. The Fluorine covalent radius F value determined as 1/2 F-F, will be  longer than the value determined from fluoro-alkanes. Once you move to using a computer it is usual to go straight to a full molecular mechanics model with a force field such as MMx, Merck or AMBER.

The data which describes how the balls and springs behave is known as the force field. The earliest force fields were only for the atoms C, H, O and N and gradually added extra atoms, atom types and functional groups. Rappe and coworkers at CalTec, Los Angeles have produced a deductive paradigm for producing a whole periodic table force field. Everything necessary to code it is in the literature. Recent versions of CHARMM also have a whole periodic table force field. The Merck Forcefield which is currently thought to be the most accurate, has  more metal elements than the older forcefields.

Force fields
A formal description of a typical force field follows.

harmonic force constant.

$$ V_{s} = \sum_{l} \frac {1} {2} k_{l} {( l - l_{0} )}^2 $$

optional anharmonic force constant.

$$ - \sum_{l} \frac {1} {3!} k_{l}^{anh} {( l - l_{0} )}^3 $$

angle bending constants.

$$ + \sum_{\theta} \frac {1} {2} k_{\theta} \left\{ {( \theta - \theta_{0} )}^2 - k_{\theta}^{'} {( \theta - \theta_{0} )}^3 \right\} $$

Fourier representation of torsional potential.

$$ + \sum_{n} \frac {1} {2^{n}} k_{\omega} ( 1 + s \cos (n\omega) ) + \sum_{\omega} \frac {1} {2} k_{\omega} ( 1 + s \cos (\omega) ) $$

Out of plane bending for trigonal groups.

$$ + \sum_{\Gamma} \frac {1} {2} k_{\Gamma} {( \Gamma - \pi )}^2 $$

Van der Waals terms.

$$ \sum_{i<j} \left\{ \frac {A_{ij}} {r_{ij}^{12}} - \frac {B_{ij}} {r_{ij}^{6}} \right\} $$

electrostatic interactions.

$$ \sum_{i<j} \frac {q_{i} q_{j}} {r_{ij}^{2}} $$

What do the force field components mean?
Let us first consider bond lengths. The most important number is $$r_0$$ the average bondlength. Then we have the harmonic force constant, which is the rigidity of the bond, obtainable from either IR spectroscopy or an energy derivative electronic structure calculation. The anharmonic constants are always negative corresponding to a steeper slope on the repulsive, compressive side and a shallowing out on moving towards dissociation. Bending constants are complicated and not so subject to simple rules of thumb. Four body dihedral functions are of course periodic and can be almost any shape with a number of minima related to valency and lonepair content of the four entities.

Pair potentials
The pair potential term is always clearly divided into a repulsive and an attractive part. In van der Waals complexes it is this repulsive function which ultimately prevents the molecules locking together usually giving the peculiarity of a positive van der Waals energy which seems counterintuitive. (Perhaps it should be relabelled, Lennard-Jones energy, or pair energy.)

For some purposes the $$\frac 1 {r^{12}} $$ term in the pair energy is too steep and a softer repulsive potential the Buckingham potential is used:

The Buckingham potential
$$ E_{vdw} = \sum { A e ^{- B r} - C r ^{-6} } $$

Another modification to the pair potential is the addition of an inverse $$r^{10}$$ term, used in some programs to simulate the extra attractive energy of the hydrogen bond.

The Lennard-Jones potential
$$ E(r) =  4 \epsilon \left( { { \left( \frac \sigma r \right) } ^ {12} - { \left( \frac \sigma r \right) } ^ 6 } \right) $$

$$\epsilon$$ is the energy at the pair minimum. The location of the pair minimum is at $$r = 2 ^{1/6} \sigma$$ i.e. $$1.122462048 \sigma$$.

MOLWEB, the software written by MG to do interfacing and MM style property calculations, internally uses either the MM2 atom types or an atom type formed from the atomic number multiplied by 100, giving the possibility of 100 different atom types. This is not a ludicrous requirement if one thinks of transitional metal complexes where differentiation between half a dozen oxidation states, several coordination numbers and high and low spin environments might be required.

The $$\frac {A} {r^{12}} - \frac {B} {r^{6}}$$ is often expressed using the 2 parameters of an energy $$\epsilon$$ and a distance $$\sigma$$.

Table 1 - The MM2 atom types
U. Burkert  and  N.  L.  Allinger,   Molecular Mechanics,  ACS  Monograph 177, (ACS,Washington,1982). A. K. Rappé and C. J. Casewit, Molecular Mechanics across Chemistry, (University Science Books, California, 1997).

Notice the absence of metals from this table. Rappé and coworkers at CalTec, Los Angeles have produced a deductive paradigm for producing a whole periodic table force field (Rappé and Casewit, 1997).

If you do not have access to a code with this in it everything necessary to code it is in the literature. It is believed recent versions of CHARMM also have a whole periodic table force field. The new Merck Forcefield (Halgren,1996)), which is currently thought by some to be the most accurate, has more metal elements than the older forcefields.

Practical Problems
i) Conformations of cyclohexanes

Check that the available force fields predict the correct energy differences between chair and skew-boat forms of cyclohexane. Confirm that the boat is impossible to make as it is not a minimum on the potential surface. Trial boat forms of cyclohexane are rather difficult to make. One possible route is to start with benzene, put the carbons into the right place after removing the double bonds. After that it is necessary to add hydrogens to resaturate valency with H-ADD. Another way is to pull the chair form down from the menu and delete the 2 hydrogens at the foot of the chair. Move the carbon with 'MXY', and then replace the hydrogens.

Establish whether 1-chlorocyclohexane prefers the axial or equatorial positions. (Make the first conformer and optimise for energy. Then invert at the 1 centre to get the other conformer and check the energy). Try this for both MM2 and AMBER force fields.

2) Geometry of cis-Chlordane

Use your bibliographic skills to get a model of cis-Chlordane, (whatever that is). Minimize it with MM2. Calculate its energy when solvated in water and chloroform, using continuum models, (if available).