Chapter 3 Calculation Methods HyperChem uses two types of methods in calculations: molecular mechanics and quantum mechanics. The quantum mechanics methods implemented in HyperChem include semi-empirical, ab initio, and density functional quantum mechanics methods. The molecular mechanics and semi-empirical quantum mechanics methods have several advantages over ab initio and density functional methods. Most importantly, these methods are fast. While this may not be important for small molecules, it is certainly important for biomolecules. Another advantage is that for specific and well-parameterized molecular systems, these methods can calculate values that are closer to experiment than lower level ab initio and density functional techniques. The accuracy of a molecular mechanics or semi-empirical quantum mechanics method depends on the database used to parameterize the method. This is true for the type of molecules and the physical and chemical data in the database. Frequently, these methods give the best results for a limited class of molecules or phenomena. A disadvantage of these methods is that you must have parameters available before running a calculation. Developing parameters is time-consuming. The ab initio or density functional methods may overcome this problem. However they are slower than any molecular mechanics and semi-empirical methods.

Molecular Mechanics Molecular mechanical force fields use the equations of classical mechanics to describe the potential energy surfaces and physical properties of molecules. A molecule is described as a collection of atoms that interact with each other by simple analytical functions. This description is called a force field. One component of a force field is the energy arising from compression and stretching a bond. 21

Molecular Mechanics

This component is often approximated as a harmonic oscillator and can be calculated using Hooke’s law. 2 1 V spring = --- K r ( r – r 0 ) 2

(7)

The bonding between two atoms is analogous to a spring connecting two masses. Using this analogy, equation 7 gives the potential energy of the system of masses, Vspring, and the force constant of the spring, Kr. The equilibrium and displaced distances of the atoms in a bond are r0 and r. Both Kr and r0 are constants for a specific pair of atoms connected by a certain spring. Kr and r0 are force field parameters. The potential energy of a molecular system in a force field is the sum of individual components of the potential, such as bond, angle, and van der Waals potentials (equation 8). The energies of the individual bonding components (bonds, angles, and dihedrals) are functions of the deviation of a molecule from a hypothetical compound that has bonded interactions at minimum values. E Total = term 1 + term 2 + … + term n

(8)

The absolute energy of a molecule in molecular mechanics has no intrinsic physical meaning; ETotal values are useful only for comparisons between molecules. Energies from single point calculations are related to the enthalpies of the molecules. However, they are not enthalpies because thermal motion and temperaturedependent contributions are absent from the energy terms (equation 8). Unlike quantum mechanics, molecular mechanics does not treat electrons explicitly. Molecular mechanics calculations cannot describe bond formation, bond breaking, or systems in which electronic delocalization or molecular orbital interactions play a major role in determining geometry or properties. This discussion focuses on the individual components of a typical molecular mechanics force field. It illustrates the mathematical functions used, why those functions are chosen, and the circumstances under which the functions become poor approximations. Part 2 of this book, Theory and Methods, includes details on the implementation of the MM+, AMBER, BIO+, and OPLS force fields in HyperChem. 22

Chapter 3

Molecular Mechanics

Bonds and Angles HyperChem uses harmonic functions to calculate potentials for bonds and bond angles (equation 9). V stretch =

∑ K (r – r ) r

0

2

V bend =

bond

∑ K (θ – θ ) θ

0

2

(9)

angle

Example: For the AMBER force field, a carbonyl C–O bond has an equilibrium bond length of 1.229 Å and a force constant of 570 2

kcal/mol Å . The potential for an aliphatic C–C bond has a minimum at 1.526 Å. The slope of the latter potential is less steep; a C– 2

C bond has a force constant of 310 kcal/mol Å . Kr = 570 kcal/mol Å2; r0 = 1.229Å

energy (kcal/mol)

Kr = 310 kcal/mol Å2; r0 = 1.526Å

bond length (Å)

Calculation Methods

23

Molecular Mechanics

A Morse function best approximates a bond potential. One of the obvious differences between a Morse and harmonic potential is that only the Morse potential can describe a dissociating bond.

Morse

energy (kcal/mol)

harmonic

bond length (Å)

The Morse function rises more steeply than the harmonic potential at short bonding distances. This difference can be important especially during molecular dynamics simulations, where thermal energy takes a molecule away from a potential minimum. In light of the differences between a Morse and a harmonic potential, why do force fields use the harmonic potential? First, the harmonic potential is faster to compute and easier to parameterize than the Morse function. The two functions are similar at the potential minimum, so they provide similar values for equilibrium structures. As computer resources expand and as simulations of thermal motion (See “Molecular Dynamics”, page 71) become more popular, the Morse function may be used more often.

24

Chapter 3

Molecular Mechanics

Torsions In molecular mechanics, the dihedral potential function is often implemented as a truncated Fourier series. This periodic function (equation 10) is appropriate for the torsional potential. V dihedrals =

∑

dihedrals

Vn ------ 1 + cos ( n φ – φ 0 ) 2

(10)

In this representative dihedral potential, Vn is the dihedral force constant, n is the periodicity of the Fourier term, φ0 is the phase angle, and φ is the dihedral angle. Example: This example of an HN–C(O) amide torsion uses the AMBER force field. The Fourier component with a periodicity of one (n = 1) also has a phase shift of 0 degrees. This component shows a maximum at a dihedral angle of 0 degrees and minima at both –180 and 180 degrees. The potential uses another Fourier component with a periodicity of two (n = 2).

energy (kcal/mol)

sum

n=1

n=2

dihedral angle (degrees)

The relative sizes of the potential barriers indicate that the V2 force constant is larger than the V1 constant. The phase shift is 180 degrees for the Fourier component with a two-fold barrier. Minima occur at –180, 0, and 180 degrees and maxima at –90 and 90

Calculation Methods

25

Molecular Mechanics

degrees. Adding the two Fourier terms results in potential with minima at –180, 0, and 180 degrees and maxima at –90 and 90 degrees. (The “sum” potential is shifted by 2 kcal/mol to make this illustration legible.) Note that the addition of V to V shows that 1

2

the cis conformation (dihedralhnco = 0 degrees) is destabilized relative to the trans conformation (dihedralhnco = 180 degrees).

van der Waals Interactions and Hydrogen Bonding A 6–12 function (also known as a Lennard-Jones function) frequently simulates van der Waals interactions in force fields (equation 11). V VDR =

∑ i