By far, the most popular application of the empirical potential energy function is to find the geometry of a molecule (or an assemblage of molecules) which corresponds to a minimum potential energy. The process of finding a minimum of an empirical potential energy function is called frequently Molecular Mechanics (MM). This process produces a motionless molecule of idealized geometry. In reality, the molecules undergo thermal motions and constantly change their geometry. Molecular Dynamics, (MD), and Monte Carlo, (MC) are based on principles of statistical mechanics and allow calculation of thermodynamic properties of molecular systems by simulating these thermal fluctuations of geometry and corresponding energies. These methods will be described briefly later in this chapter. For a review of molecular mechanics consult Burkert & Allinger, (1982).
The potential energy function refers to a ground electronic state of the molecule and does not explicitly incorporate electronic structure. Hence, it cannot be used to study electronic spectra or covalent bond formation/breaking. These phenomena are studied with quantum methods. There are methods which combine a quantum description for the small portion of the molecule (e.g., groups of substrate and enzyme interacting in the active center) with a molecular mechanics description for the rest of the molecule (see e.g., Field et al., 1990).
The minimization of the potential energy function (i.e., geometry optimization) involves a search for the minimum of a function, and, to be efficient, requires calculations of derivatives of a function (in this case, the potential energy) versus independent variables (in our case, coordinates). Most programs use cartesian coordinates as independent variables, however, in some cases, internal coordinates may be used. The derivatives of potential energy are denoted as:
where 
is the gradient (i.e., first derivative) of the potential energy
with respect to a cartesian coordinate 
of an atom
,
and 
is the second derivative of the energy with respect to
the cartesian coordinates. In most modern programs these derivatives are
calculated analytically, i.e., the appropriate mathematical formulae
for corresponding terms are incorporated into the program. Some older codes
compute derivatives numerically by approximating the slope
of an energy function (or its gradient in the case of second derivatives) from
finite differences.
The table of all possible second derivatives versus cartesian coordinates
of atoms has 3N rows and 3N columns
and is called a Hessian matrix. The derivatives are used not
only in function minimization but also yield forces acting on
atoms (from energy gradients) and normal modes of vibration (from the Hessian
matrix). There are three
major approaches to finding a minimum of a function of many variables:
(i.e., prohibitive for large macromolecules). The
BFGS algorithm
is considered the most refined one.
In general, the minimization methods are iterative. They require on input some initial estimate for the position of the minimum, and provide a better estimate for the minimum as a result. This corrected estimate is used as an input into the next cycle (i.e., iteration) and the process is continued until there is no significant improvement in the position of the minimum.
Figure 6.23: Local and global minima for a function of one
variable and
an example of a sequence of solutions: 
for a descent series minimization algorithm.
Most search methods and minimization methods using derivatives are the descent series methods, i.e., each iteration results in a solution which corresponds to a lower (or equal) value for the energy function:
As a consequence, these methods can only find the minimum closest
to the starting estimate and will never cross to a minimum (however deep)
if it is separated from the starting estimate by a maximum (however small).
This situation is schematically illustrated in Fig. 6.23 where arrows
indicate the direction of geometry optimization depending upon the
starting point. Fig. 6.23 is only a cartoon to illustrate the
behavior of descent series minimization methods. Geometry optimization of
real molecular systems (e.g. a protein with a drug molecule bound) involves
simultaneous optimization of 3N cartesian coordinates, i.e, sometimes
many thousands of variables. There is no
general way of finding a global minimum (i.e., the minimum corresponding to
the lowest possible value of the function). A different
initial geometry will usually lead to a different final minimum. Only on very
simple molecules will the single geometry optimization yield the global
minimum on the first trial. It cannot be overemphasized
that the result of a single minimization is usually a local
minimum, not a global one. To find a global minimum (or at least, to be
more confident about it) you have to perform many minimizations and use
different initial coordinates for each run.
There is another practical difficulty with contemporary molecular
mechanics programs. Most of them optimize the geometry in atomic cartesian
coordinates. To illustrate the problem let us join Alice on her journey
to the other side of the mirror. The appearance of the
potential energy function in cartesian space
is dramatically different from the one given by
internal coordinates. If we could
see a cartesian representation of a potential energy surface in, let us say,
fifty dimensions, the picture would look like the cortex of
a human brain --
lots of narrow tortuous valleys of similar depth. This is because low energy
paths for individual atoms are very narrow due to the presence of
hard bond stretching and angle bending terms. The low energy paths
correspond only to the rotation of groups or large portions of the molecule
as a result of varying torsional angles. In other words, the low energy paths
for atoms
are strongly correlated in cartesian space and atoms can move easily only
along narrow grooves. On the other hand, the
outlook of the potential energy surface is very different in the space of
internal coordinates. The surface looks like a valley surrounded by high
mountains. The high peaks correspond to stretching and bending terms,
and close van der Waals contacts, while the bottom of the valley represents
the torsional degrees of freedom. If you happen to start at the mountain tops
in the internal coordinate space, the
minimizer sees the bottom of the valley clearly from above.
If you are in the valley, you also see where the mountains are.
In cartesian space, the minimizer walks along the bottom of a narrow winding
channel which is frequently a dead-end.
In more scientific terms, all contemporary function minimization procedures
operate efficienty for functions close to quadratic (i.e., for two dimensions
a parabola, for three dimensions a boat like shape, for n dimensions:
)
and are inefficient for functions of other shapes. Assuming that bond
lengths and valence angles are already optimized and the goal of minimization
is to find the minimum corresponding to optimal values of torsional angles,
cartesian space is probably the worst possible. Look at equations
6.16, 6.27 and 6.9, which relate the torsional energy
to cartesian coordinates.
The dependance of the torsional energy on cartesian coordinates
is as far from a quadratic shape as it can be! Moreover, it has several
minima. There are also additional advantages to using minimization
in internal coordinate space. In this space there is a clear separation of
variables into the hard ones (i.e., those whose small changes produce
large changes in the function value) and soft ones (i.e., those whose
changes do not affect the function value substantially). During function
optimization in internal coordinates, the minimizer first optimizes
the hard variables and in the subsequent iterations ``cleans up'' the details
by optimizing the soft variables. In cartesian coordinate space all
variables are of the same type, i.e., hard or soft depending on
preference. Changing a cartesian coordinate of an atom results
in simultaneous change of bond lengths, valence angles and torsonal angles
associated with this atom, and distances to non-bonded atoms.
For this reason, the minimizer is only involved with details,
and therefore in the
last stages of the optimization all variables are as important as at the
beginning. The reason why most programs optimize geometry in
cartesian space is that it is not easy to derive and manage expressions for
derivatives of potential energy function in internal coordinates.
Note that the derivative of a van der Waals energy term between two atoms
versus the cartesian
coordinates of some third atom is zero. This derivative depends only
on the
cartesian coordinates of the two atoms involved. The situation is dramatically
different in internal coordinate space. Since the distance between
two atoms depends on all bond lengths, bond angles and torsional angles
on the path between two atoms, the derivative of the van der Waals term
for the two atoms will depend on all these variables. There are
matrix methods of calculating derivatives of distances and angles as
contributions from individual internal coordinates, and the final
derivative is calculated from the chain rule using these parameters. But as
you can see, the problem is far from trivial.
To avoid evaluation of derivatives, some programs use search techniques
like SIMPLEX when optimization
is requested in internal coordinate space. Since SIMPLEX is inefficient
for large numbers of variables, usually only torsional angles are
optimized by this technique. If the number of variables is too large,
then a block technique is used where only a portion of all independent
variables is optimized at a given stage. Another reason why
internal coordinates are not used routinely for geometry optimization
is that modern molecular mechanics programs are usually a portion of
a molecular mechanics/dynamics package. As will be shown later in this
chapter, molecular dynamics in coordinates
other than cartesian is a difficult problem, since the terms containing
generalized momenta and generalized coordinates in the Hamiltonian
cannot be separated. Molecular mechanics and dynamics
share large portions of the actual code, and hence there is a preference by
software authors to use the space which is common to both methods.
The potential energy of the molecule calculated from a well balanced
force field represents a strain in the molecule. Augmented
with bond/group equivalents and statistical mechanical corrections,
it can be used to estimate the heat of
formation of a compound (i.e., its molar enthalpy of formation,
). The quantity can be used
to compare the relative stability of
different compounds. In most cases however, you should assume that the
calculated potential energy incorporates some arbitrary component which
depends upon the types of atoms and covalent bonds in the molecule,
i.e., comparison of the energies calculated for
molecules with different numbers and/or types of atoms and bonds cannot be
rigorous. For this reason, potential energy will in most cases reliably
assess the difference in energy between conformers,
but will fail if you attempt to calculate the energetics of incorporating
a new functional group into the molecule. Molecular mechanics can also provide
the interaction energy, 
, of two molecules A and B as:
where 
, 
, and 
are
potential energies of the optimized complex, the optimized molecule A, and
the optimized molecule B;
respectively. Note that the type and number of atoms and covalent bonds
in the complex AB is equal to their sum in isolated molecules A and B,
and the arbitrary ``energy zero'' should cancel out in this case.
For this reason, the difference between interaction energies calculated for
different complexes, (i.e., 
) is
the preferred method over direct comparison of the energies of different
complexes (i.e., 
).
Potential energy functions can also be used to estimate contributions
from intramolecular vibrations to the vibrational free
energy and vibrational entropy. These quantities, and
contributions from translation and rotation of the molecule as a whole,
vary with temperature
and are main contributors to the thermodynamic functions such as enthalpy,
free energy, specific heat, etc., (see, e.g., Hill, 1960).
A well known approach is to use the frequencies, 
,
corresponding to normal modes within harmonic approximation (i.e, to
derive them from a mass scaled Hessian matrix at energy minimum).
The expressions for relating classical vibrational contributions to Helmholtz
free energy, 
, internal
energy, 
, heat capacity at constant volume
,
and entropy 
of the nonlinear molecule, are derived
in many standard
textbooks for statistical mechanics (see e.g., McQuarrie, 1976):
where R, T, and h are the gas constant, the
absolute temperature, and Planck's
constant respectively; and N denotes the number of atoms in the molecule.
Frequently, these values are augmented with a
correction to account for vibrations at 0 
K, which is of quantum origin,
by adding a zero-point energy (ZPE) to the
free energy and the internal energy :
The harmonic approximation is quite accurate for isolated molecules. However, for complexes of two or more molecules or systems containing the solvent, the harmonic approximation breaks down, because the motion of the individual molecules with respect to each other is diffusional, i.e., not constrained by a skeleton of rigid covalent bonds. The potential of this motion is different than the simple parabola of a harmonic oscillator and contains multiple minima and maxima. In this case, molecular dynamics or Monte Carlo approaches are more reliable for estimation of thermodynamic functions.