| http://server.ccl.net/cca/documents/molecular-modeling/node7.shtml |
|
CCL node7 |
   Next: Molecular Mechanics Up: Molecular Modeling Previous: Types of molecular images
Empirical Potential Energy FunctionMolecular modeling raison d'être is to perform computer simulations, or to do chemical experiments with a computer rather than a laboratory bench. Though the results of such simulations have to be validated by experiment, it is known that in many cases we may rely on them as much as on the true experimental values. On many occasions calculations showed flaws in original experimental data and led to corrected experiments which agreed with the theory. Simulation requires a model, a mathematical description of a system to be simulated. It is important to accept that the particular model need not be ``real'', but good models usually have a sound scientific justification to yield useful results. This section briefly outlines potential energy functions used in molecular mechanics and dynamics methods. Intelligent use of these methods requires more information than can be given here. An excellent monograph on this topic by Burkert & Allinger (1982) is available. Also the actual software instructions accompanying molecular modeling systems describe the terms in great detail. There are in general two approaches to molecular simulation: quantum methods, and empirical methods. Quantum methods are based on principles of quantum mechanics. The most rigorous implementations of quantum formalism, called ab initio methods, do not require and experimental data, beside a handful of basic physical constants like: speed of light, electron mass, elementary charge, etc. The empirical methods reviewed in this section derive their origin from classical chemical concepts and are based for the most part on experimental observations. This does not mean however that they do not use quantum results and procedures. In fact there are numerous examples of combination, quantum/empirical approaches, where some parts of a system are treated quantum mechanically and others are described by classical mechanics. The boundary between these methods becomes more and more blurred with time.
It is well known from structural chemistry and from quantum calculations
that bond lengths and
valence angles in typical units and groups are very similar even if they
appear in different molecules. The single bond between two sp
where
Merging parameters from
different force fields is therefore discouraged.
While the internal coordinates can easily be
calculated from the cartesian coordinates of the atoms
of the actual molecule, the parameters entering into the
energy function must be known in advance for all types of atoms comprising
the molecular system. As will become clear from the following discussion, many
parameters are needed to describe even simple biological molecules.
For some more exotic atom types, such parameter are still not available.
You can always put some number
in place of the missing force field parameter, however, your results will be
as reliable as these parameter values. There is a growing tendency to
estimate needed force field parameters from ab initio quantum
calculations, by performing such calculations on small, model molecules
assuming that the results are transferable to larger molecules. This approach
has the added advantage that calculations can be performed for all, even
non-existing species, while reliable experimental data are only available
for the most popular molecules.
Individual terms
The bond stretching term,
where
The form of the angle bending term is very similar. In most cases the
contribution to potential energy,
The out-of-plane term is used to account for the energy contribution from distorting an aromatic or conjugated system of bonds from planarity and is most frequently a harmonic term:
where in the case of sp The functional form of the torsional term is unique. In most cases, it depends on the cosine of the product of a torsional angle and the periodicity (see Fig. 6.21):
where
Non-bonded terms describe contributions brought by the interaction of atoms which are not covalently bonded. Moreover, atoms which are two bonds apart (1-3 interactions), are usually not included in the non-bonded interaction list, since it is assumed that their interaction is satisfactorily accounted for by the angle bending term. However, non-bonded interactions of atoms separated by 3 covalent bonds in the molecular graph (1-4 interactions) are routinely included in spite of the fact that they appear in the torsional term, though their magnitude is frequently scaled by 50%. Generally, the non-bonded interactions represent contributions to energy from van der Waals interactions, electrostatic interactions and they frequently account explicitly for hydrogen bonds. There are many expressions for these terms and only the simplest ones will be described here to illustrate the underlying ideas.
The non-bonded terms are usually represented mathematicaly
as two-body interactions. They depend on the coordinates of only two
interacting atoms, i.e., are represented as pairwise potentials.
This is an approximation,
and it is well known that the interaction between two atoms depends
significantly upon the positions of other atoms, especially those at close
proximity. However, even for pairwise potentials, calculation of
non-bonded terms is at present
a significant computational effort proportional to the square in
number of atoms. For three-body potentials (i.e., functions which
depend on positions of three non-bonded atoms), this effort is
proportional to the cube in number of atoms and unmanageable for larger
molecules. For example, for a small sizes protein containing 1000 atoms,
there are roughly
The van der Waals interaction energy,
The repulsive term,
where
The electrostatic interaction energy,
where
The electrostatics
of molecules (especially macromolecules) is an important factor in their
biological function since it is responsible
for long-range inter- and intramolecular forces. These are large forces
compared to other interactions, e.g., the interaction energy of
two point charges -0.3 a.u. and +0.3 a.u., respectively,
separated by a distance of 10 Å in vacuum is
The major ambiguity, however, is with the concept of atomic charge itself.
The atom is made up of
a nucleus and electrons orbiting around it, and its charge is zero
by definition. However, in the molecule, averaged electron trajectories are
no longer centrosymmetric around nuclei, due to the formation of chemical bonds
between atoms, and the resulting average electron density can be a very
convoluted function in three dimensions.
Depending on atom electronegativity There are several ways of deriving atomic charges depending on which molecular property is chosen to be represented best. In many popular force fields they are chosen to reproduce known experimental data, e.g.: dipole moments, geometries, vibrational spectra, etc. They are frequently derived theoretically from a variety of electronegativity equalization schemes and quantum calculations at different levels of sophistication. In the case of quantum chemical calculations, the charges may be derived either from population analysis or fitted by least squares to the electrostatic potential. The charges derived from ab initio electrostatic potential are considered superior but this method is also the most computationally demanding.
The effect of hydrogen bonding is incorporated into the
potential energy surface
of a molecular system in a variety of ways.
In general, hydrogen bonds are formed
between hydrogen donor and hydrogen acceptor groups:
--D--H...A--AA--. The simplest approach, which at the
same time performs very well, is to assume that the bond between
hydrogen donor atom D and hydrogen atom H is polarized to the
extent that the electron density in the vicinity of the proton is negligible,
i.e., that the proton does not contribute to the
van der Waals energy. In this model (introduced by Hagler and coworkers),
the proton and the acceptor atom attract each another since they possess
opposing net charges. This results in the stretching of the
D--H bond and driving atom
A towards the hydrogen. Balance is achieved by contributions from
the D--H
stretching term, and electrostatic and van der Waals repulsions of atoms
D and A. This approach does not impose any directionality
on the hydrogen bonds but the balance of all of these
forces usually yiels a geometry close to linear (i.e., when atoms
D, H, A and AA lie on the straight line).
Some more elaborate models explicitly use powers of
The last category of terms which appears in the potential energy function is constraints. Their utility depends on the intended use of the potential energy function and they will be described later in this chapter. The individual energy terms differ in the magnitude of their contribution to the total potential energy of the molecule. The bond stretching and angle bending terms are considered ``hard'' terms, since even a small distortion in the ``optimal'' value of bond length or valence angle produces a very large positive contribution to the potential energy. For this reason, the variation in bond lengths and valence angles is usually negligible in molecules without substantial strain. Larger variations are only evident in strained ring systems or crowded molecules. In fact, the overall shape of the molecule results from a balance between the softer contributions, i.e., torsional terms and non-bonded terms. It is therefore a common practice to freeze bond lengths and valence angles at their optimal values and consider only torsional and non-bonded interactions as the first approximation.
Next: Molecular Mechanics Up: Molecular Modeling Previous: Types of molecular images Computational Chemistry Wed Dec 4 17:47:07 EST 1996 |
[ CCL Home Page ]
[ About CCL ]
[ Resources ]
[ Search CCL ]
[ Announcements ]
[ Links ]
[ E-mail us ]
[ Raw Version of this page ]
| Modified: Sat May 23 16:00:00 1998 GMT |
| Page accessed 11042 times since Mon Apr 26 10:40:55 1999 GMT |
hybridized
carbon atoms is around 1.5 Å long and the valence angle for an
sp3 carbon is
usually close to 109 
. On the other hand, these values are frequently
distorted in some strained ring systems (like cyclopropane) or
crowded molecules (like tetramethylmethane). The concept of molecular
strain can be imagined as if bonds were elastic springs whose distortion
is reflected by a positive increment (i.e., thermodynamically
unfavorable) to the potential energy of the molecule.
These approaches were originally called empirical potential energy
functions but the term molecular mechanics is now in common use.
The mathematical form of this
potential energy function (also called potential energy surface) is:
represents the potential
energy of the molecular system and depends upon the cartesian coordinates
of all atoms denoted here as 
. 
's are individual energy terms
representing contributions from individual interactions which
depend upon the cartesian coordinates of partaking atoms.
Having 
, for a covalent bond between atoms
i and j, represents a contribution to the potential energy
resulting from deformation of the optimal bond length. Most frequently a simple
symmetric parabola (i.e., harmonic approximation) is used:
is a bond stretching parameter (frequently called
a bond stretching force constant) and its value represents the stiffness
of the bond (i.e., the ease with which the bond can be distorted). Large
values of 
corresponds
to the optimal (i.e., unstrained) bond length. Contribution
from this term is never negative and is equal to zero only when the actual bond
length, 
, is exactly equal to the optimal bond length 
, is represented by a
harmonic expression depending upon the bending force constant,
, the actual value of the valence angle 
and its
``natural'', optimal value, 
:
hybridized carbon, 
is an angle between
one of the bonds originating at the central atom and a plane passing through
the other two bonds.
is the height of the torsional barrier; s is -1 if a minimum
occurs for the eclipsed conformation and +1 if a minimum corresponds to
the staggered conformation; n is a periodicity, i.e., the number of maxima
per full revolution; and 
is a torsional angle.
For the torsional angle involving C--C bond in ethane, s=1 and n=3, while
for the double bond in ethylene s=-1 and n=2. More elaborate
force fields use as a torsional term the sum of several cosine functions
with different periodicities to account for smaller ``humps'' appearing
on the torsional energy dependence. On occasion, the harmonic potential
is used for torsional angles around double bonds due to their stiffness.
.
pairwise
interactions,
or 
terms. The 
factor
comes from the fact that each pair is included only once, i.e., the pair
i--j is equivalent to j--i. For the three-body potentials
we would have to calculate interaction for each possible combination of
three atoms (excluding valence angles), i.e., close to
or
terms! For that reason, the three-body potentials
are not routinely calculated. However, even two-body potentials are
computationally demanding and approximations are used.
The large number of terms is frequently
reduced by omitting the
interactions of distant atoms (i.e., skipping calculation
of the non-bonded term for atoms which are further apart than some
predetermined
cut off distance). However, checking if atoms are outside the cut off distance
is also a substantial computational effort for larger molecules.
Å, 
kcal/mol,
Å
, is composed
of two terms opposing each another, a repulsive term 
, and
an attractive term, 
:
:
is the depth of a well (i.e., the maximal attraction energy),
n and m are exponents (typically n=12 or 9, and m=6),
and 
is the distance between atoms which corresponds to a minimum.
As you can see from Fig. 
,
between two non-bonded atoms i and j is usually
represented by Coulomb's law:
, 
are the net atomic charges of atoms
i and j, respectively,
kcal/mol, i.e.,
very substantial for such a large separation between these modestly charged
atoms. For this reason,
the research in the area of adequate representation of molecular
electrostatics, yet computationally feasible, is very active.
There are methods which allow one
to incorporate polarizability, changes in dielectric constant (i.e., on the
protein/solution interface), etc., and these have proved to be very successful
in describing long range electrostatic interactions. However, these
methods are computationally expensive and for routine calculations
simplistic approximations are used.
The dielectric constant, k, is a well defined quantity only for
macroscopic systems. At the molecular level, a common approximation
is to use its value for a vacuum (i.e., k=1). To take into account
the presence of other atoms in the molecule or to incorporate the
influence of the solvent, different values are used in the range
between k=1 (vacuum) to k=80 (water). Another approach
is to assume that the
dielectric constant is proportional to the distance between
charges, usually in a linear fashion, 
. This reflects the
fact that for distant charges there is a higher probability for other
atoms or solvent molecules (i.e., polarizable entities) to enter the
space between the two charges.
and
to enforce the linearity of the observed hydrogen
bond geometry and incorporate the
12-10 Lennard-Jones potential for the distance
between H and A.