Derivation of Class I1 Force Fields. I.

Methodology and Quantum Force Field

for the Alkyl Functional Group and
Alkane Molecules
J. R. MAPLE, M.-J. HWANG, T. P. STOCKFISCH, U. DINUR,t
M. WALDMAN, C. S. EWIG, and A. T. HAGLER"
Biosym Technologies, lnc., 9685 Scranton Road, San Diego, California 92121

Received 29 July 1993; accepted 19 August 1993

ABSTRACT
A new method for deriving force fields for molecular simulations has been
developed. It is based on the derivation and parameterization of analytic
representations of the ab initio potential energy surfaces. The general method
is presented here and used to derive a quantum mechanical force field
(QMFF) for alkanes. It is based on sampling the energy surfaces of 16
representative alkane species. For hydrocarbons, this force field contains 66
force constants and reference values. These were fit to 128,376 quantum
mechanical energies and energy derivatives describing the energy surface.
The detailed form of the analytic force field expression and the values of all
resulting parameters are given. A series of computations is then performed to
test the ability of this force field to reproduce the features of the ab initio
energy surface in terms of energies as well as the first and second derivatives
of the energies with respect to molecular deformations. The fit is shown to be
good, with rms energy deviations of less than 7% for all molecules. Also,
although only two atom types are employed, the force field accounts for the
properties of both highly strained species, such as cyclopropane and
methylcyclopropanes, as well as unstrained systems. The information
contained in the quantum energy surface indicates that it is significantly
anharmonic and that important intramolecular coupling interactions exist
between internals. The representation of the nature of these interactions, not
present in diagonal, quadratic force fields (Class I force fields), is shown to be
important in accounting accurately for molecular energy surfaces. The Class I1
force field derived from the quantum energy surface is characterized by
accounting for these important intramolecular forces. The importance of each
'Author to whom all correspondence should be addressed.
tl'resent address: Dept. of Chemistry, Ben Gurion Univ. of
the Negev, Beersheva 84105, Israel.
Journal of Computational Chemistry, Vol. 15, No. 2, 162-182 (1994)
0 1994 by John Wiley & Sons, Inc. CCC 0192-8651l941020162-21
 DERIVATION OF CLASS II FORCE FIELDS

of the interaction terms of the potential energy function has also been
assessed. Bond anharmonicity, angle anharmonicity, and bond/angle, bond/
torsion, and angle/angle/ torsion cross-term interactions result in the most
significant overall improvement in distorted structure energies and energy
derivatives. The implications of each energy term for the development of
advanced force fields is discussed. Finally, it is shown that the techniques
introduced here for exploring the quantum energy surface can be used to
determine the extent of transferability and range of validity of the force field.
The latter is of crucial importance in meeting the objective of deriving a force
field for use in molecular mechanics and dynamics calculations of a wide
range of molecules often containing functional groups in novel
environments. 0 1994 by John Wiley & Sons, Inc.

of this article and others that will follow in this

Introduction series6
Procedurally, although force fields have tradi-
tionally been derived almost exclusively from ex-

T he rapid pace of experimental research in

complex molecular systems, particularly in
biochemistry and biophysics, has led to greatly in-
perimental data,7-25 including thermodynamic
properties, vibrational frequencies, gas-phase mo-
lecular structures, and crystal structures, a lack of
creased need for accurate molecular simulation such experimental data makes it difficult and in
techniques. The most powerful techniques at pres- some cases impossible to parameterize and test ac-
ent are molecular mechanics and dynamics, which curate potential energy functions for use in molec-
presuppose an analytic energy function of the ular mechanics and For example,
atomic coordinates, the force field, that then gov- there is little experimental information on the gas-
erns all the relevant molecular energetic, structural, phase structures, vibrational frequencies, and en-
and dynamic properties. Although in some appli- ergetic properties of amino acids, ammonium cat-
cations only a qualitative picture derived from the ions, and carboxylate anions. Thermodynamic
simulation is sufficient, in the great majority of data, such as conformational energy differences,
cases an increasingly high level of accuracy is crit- rotational barriers, and sublimation energies, are
ical for obtaining meaningful results. Examples in- especially useful for developing force fields and are
clude refinement of X-ray crystallographic' and nu- also rare for many functional groups. In addition,
clear magnetic resonance (NMR)2structural models experimental structures and spectroscopic infor-
of biological species, simulation of enzymatic re- mation on higher-energy conformers or transition
action rates3 and binding constant^,^ and predic- states are relatively nonexistent.
tion of protein structures and interaction e n e r g i e ~ . ~ Further, because our major objective in carrying
A number of factors affect the accuracy of mo- out theoretical molecular simulations of biological
lecular simulations as compared to practical exper- and other organic systems is to study molecules
imental results, including the treatment of solvent, for which detailed structural or thermodynamic in-
counterions, and other species that may be present formation is not available from experiment, we de-
in the modeled environment, as well as other de- sire a predictive tool, and hence force-field func-
tails of the simulation techniques. However, no tional forms and parameters consequently must be
advance in these techniques can compensate for transferable. To derive a transferable force field, it
inadequacies in the underlying force field. Further, is generally necessary to fit substantial amounts of
because in practice the final results will generally data with the number of independent experimental
depend on a delicate balance of many differing observables far exceeding the number of parame-
types of terms in the model, inadequacies in a force ters, thereby enabling extrapolation or transfera-
field may be manifested by gross qualitative errors bility to molecular environments not included in
in predicted results. It is therefore critical that there the set of molecules from which the force constants
exist an accurate and systematic procedure for de- are determined. Unfortunately, because for most
riving the force fields needed to model any molec- functional groups a sufficient amount of this ex-
ular species or environment of interest. Demon- perimental data is not available, the ratio between
strating and testing such a procedure is the topic the numbers of observables and parameters can be

JOURNAL OF COMPUTATIONAL CHEMISTRY 163

MAPLE ET AL.

increased only by reducing the number of param- types of terms in the energy expression may be
eters. This essentially means using a smaller num- added, deleted, or reformulated using the fit to the
ber of force constants from a simplified and cor- theoretical energy surface and its derivatives as cri-
respondingly less accurate and less transferable teria of physical importance.26Thus, the complete
potential energy function. Thus, the problem of energy surface provides in a sense a laboratory for
limited experimental data has been one of the fun- developing increasingly more realistic and pow-
damental limitations both in deriving parameters erful types of force fields.
and in developing improved potential energy func- This article describes the derivation of a quan-
tional forms. tum mechanical force field (QMFF) that approxi-
Recently, a method for exploiting the vast mates the Hartree-Fock (HF) 6-31G" potential en-
amount of information that can be obtained from ergy surface of alkanes, including highly strained
ab initio quantum mechanical calculations to, in a small ring molecules (e.g., cyclopropane and cy-
way, amplify the experimental data was intro- clobutane). The methodology (previously used to
d u ~ e d .It~not~ ,only
~ ~ results in a considerable in- derive a quantum mechanical force field for the
crease in the ratio of observables to parameters, formate a n i ~ n ~ ~is, ~extended
') to fit the ab initio
but it also provides an objective criteria for the test- potential energy surfaces of a representative set of
ing of alternative functional forms and the trans- alkanes, thereby maximizing the transferability of
ferability of force constants. The method makes use the resulting force field to all saturated hydrocar-
of distorted structures of molecules created by de- bons, regardless of the nature of branching or cy-
formations of all internal coordinates. The energy clization or the degree of strain energy. The ac-
and first and second derivatives of the energy with curacy of the QMFF potential energy function is
respect to the atomic Cartesian coordinates are cal- assessed with regard to the fit of the quantum me-
culated for each of these distorted structures by ab chanical potential energy surface. In particular, the
initio molecular orbital methods. These derivatives importance of including anharmonicity and cross-
provide literally thousands of data that describe the terms, as reflected by the energy surfaces of organic
quantum mechanical energy surface of the mole- molecules, is reexamined for the case of alkanes.
cule. For M structures of a molecule containing N The importance of these terms has been a subject
atoms, the resulting ab initio data include M - 1 of several previous s t ~ d i e s . ~ ~ - ~ O
relative energies, M(3N - 6) independent first de- The first section of this article describes a method
rivatives, and M(0.5)(3N - 6)(3N - 5) second de- for sampling the ab initio potential energy surface
rivatives of the energy with respect to a well-de- with distorted structures generated by displacing
fined set of internal coordinates. For example, only atoms in the directions of the normal modes of
five structures of n-butane, which has 14 atoms, vibration. The next section defines the character-
yield a total of 3514 independent quantum me- istics of the functional form that emerge from the
chanical "observables." This is only a single, rel- fit to this energy surface. Next, a least-squares
atively small, molecule. Also, when parameteriz- method for parameterizing potential energy func-
ing a quantum energy surface numerous differing tions to fit ab initio energy and energy derivative
molecules are used, as shown below, thereby pro- data is described, along with the derivation of the
viding more than enough information for para- force field. Then, the quality of fit to the ab initio
meterizing complex potential energy functions. energy surface is presented. Finally, the impor-
Further, to account for the intricacies of the energy tance of anharmonic terms and cross-terms derived
surfaces of molecules, this data contains a wealth from the detailed form of the ab irzitio potential en-
of information on the energy of distorted structures ergy surface of alkanes is also examined.
(i.e., anharmonicity), transition states, and the en- Force fields derived by this method provide a
ergy changes accompanying displacements of par- good starting point for applications to large mol-
ticular internal coordinates, information that is vir- ecules when experimental data for (relatively small)
tually inaccessible to experiment. representative molecules of a given functional
It is important to note that parameterizing the group is too limited or even entirely unavailable.
energy of a molecule employing the quantum me- This method is also useful when one needs a
chanical energy surface not only provides a means prompt derivation of necessary terms in the force
for determining the numerical values of individual field and time is not available for fits of experi-
terms, but also permits development of the analytic mental data. The resulting force fields should pro-
form of the energy expression as Differing vide similar quality to that which would be ob-

164 VOL. 15, NO. 2

 DERIVATION OF CLASS II FORCE FIELDS

tained if quantum mechanical calculations could be mental data.6 Of course, this two-step methodol-
used on the complete system of interest. However, ogy should be most valuable for those functional
because systematic errors are introduced into the groups for which only a limited amount of useful
force field by the quantum mechanical data (e.g., experimental data is available (i.e., the great ma-
stretching force constants from HF are known to jority of functional groups). Nevertheless, com-
be -15% too a method for refining the parisons with alkane force fields previously de-
potential energy function to fit experimental data rived exclusively from experimental data7-23assist
is needed and is described in a following article.6 in the evaluation of the methodology and the es-
This article describes only the first part of a two- tablishment of a benchmark for the evaluation of
part method for deriving force fields from, first, alternative functional forms, as well as of force
quantum mechanical data and, second, experi- fields for other functional groups. This two-step

(A) Straight Chain Hydrocarbons

J
" 6
methane ethane propane

n-butane n-pentane

(8) Branched Hydrocarbons

isobutane isopentane neopentane

(C) Cyclic Hydrocarbons

cyclopentane cyclohexane

FIGURE 1. (A) Straight-chain, (B) branched, (C) cyclic, (D) four-membered ring, and (E) three-membered ring
molecules used to comprise the training set.

JOURNAL OF COMPUTATIONAL CHEMISTRY 165

MAPLE ET AL.

methodology for deriving force fields can be con- curvature of the ab initio potential energy surface
sidered a logical extension of Lifson's consistent are computed for the distorted structures.
force field (CFF) method"-'5 for deriving force Molecules chosen for the training set are de-
fields from as wide a variety of experimental data picted in Figure 1. Methane, ethane, propane, n-
as possible because it first takes advantage of the butane, and n-pentane were selected to represent
vast amount of data available from quantum me- straight-chain n-alkanes, while isobutane, isopen-
chanics and then uses the same wide variety of tane, and neopentane were chosen to give infor-
experimentally available data, as proposed in the mation on branched-chain hydrocarbons. Cycloal-
CFF methodology. kanes were sampled by including cyclopentane
and cyclohexane, and cyclobutane, methylcyclo-
butane, cyclopropane, and methyl cyclopropane
Potential Energy Surface Sampling derivatives were included to ensure the sampling
of highly strained interactions, as discussed below.
Several steps are involved in sampling the po- Equilibrium structures and normal modes of vi-
tential energy surface. First, a training set of mol- bration of the alkanes in the training set were gen-
ecules must be selected for the parameterization. erated by DISCOVER (Biosym Technologies), a mo-
Then, the normal modes of vibrational dispiace- lecular mechanics and dynamics program
ments are computed from an equilibrium structure employing the CVFF empirical force field,'3-'5 or
derived by energy minimization. Next, the atoms with the AM1 Hamiltonian in the semiempirical
are displaced in the directions of the normal mode programs MOPAC or AMPAC (QCPE). Gaussian
vector components. Finally, the energy, forces, and 8833(Gaussian, Inc.), an ab initio quantum mechan-

CD) Four-membered Rings

cyclobutane methylcyclobutane

(E) Three-membered Rings

)Jk-f$ H H

0
cydopropane methycyclopropane

1,2-dimethylcyclopropane 1.1-dimethylcyclopropane

FIGURE 1. (Continued)

166 VOL. 15, NO. 2

 DERIVATION OF CLASS I I FORCE FIELDS

ics program, was also used to optimize molecular 1.169 A and H-C-H angles ranging from 92.8 to
structures. For distorted structures, GRADSCF 127.8'. The first structure in this figure corresponds
(Polyatomics, Inc.) was employed to compute the to the ab initio optimized geometry. All others have
energies and energy derivatives. Because the goal
was to create random distorted structures that ex-
plore the quantum mechanical energy surface, it
does not matter in principle which of these pro-
grams is used to calculate the starting geometry
and normal modes. However, ab initio calculations
are of course computationally more expensive.
To determine the amplitude of atom displace-
1093 104.r
ment in the directions of the normal modes of vi-
bration, the energy of each individual normal mode E = 0.00 kcal/mol E = 7.52 kcal/mol
was randomly selected within a range of 0-2 kcal/
mol, generally the range of interest in molecular
mechanics and molecular dynamics. Then, the dis-
placement amplitude (AQJ in mass-weighted co-
ordinates was computed with the harmonic
appro~imation~~

E, O.5w2(AQ,)' (1)
~~

= 126.2 102.2

from the vibrational frequency ( w ) in rad/s and the E = 20.47 kcal/mol E = 8.62 kcal/mol
randomly selected energy (E,) of each normal
mode. Next, the representation of the orthogon-
alized displacement vectors was transformed from
Cartesian to internal coordinates, and each com-
ponent of these normalized displacement vectors
was multiplied by the random amplitude (AQ,) of
the mode. Each component of the resulting linearly 120.1
independent displacement vectors displaced an in-
ternal coordinate, and the corresponding compo- E = 16.33 kcal/mol E = 7.65 kcal/mol
nents for each mode were simultaneously added
to the internal coordinates of the equilibrium struc-
ture. Finally, Cartesian coordinates were deter-
mined from the internal coordinates of the result-
ing distorted structure. Usually, 5-10 distorted
structures were created per molecule. To improve
the sampling of torsion potentials, additional struc-
tures of ethane, propane, and n-butane were cre- 117.2

ated by rigid rotation about the H-C-C-H, E = 7.52 kal/mol E = 3.25 kcal/mol
H-C-4-C, and C--C-C--C torsions, respec-
tively, starting from the optimized (global mini- FIGURE 2. Examples of distorted structures used in
mum) structure. In addition, ab initio optimized the derivation of the quantum mechanical energy
structures, which were needed for the purpose of surface of hydrocarbons. Here, we show several
testing the force field, were also included in the distorted structures of methane from which we calculate
relative energies and first and second derivatives. As
potential energy surface sample.
seen from these structures, the energies range from the
Figure 2 gives a set of examples of internal co- minimum energy structure up to 20 kcalimol. The
ordinate distortions in methane as generated by H-C-H angle is sampled from 93 up to 128" and a
this procedure, as well as the resultant energies. large variety of combinations of these angles (to enable
The range of distortion is exemplified by these sampling of different couplings) are also probed. This is
structures. Thus, the sampling for methane in- representative of the types of distortions that are
cludes C-H bond lengths ranging from 1.031 to generated in the training set of alkane molecules.

JOURNAL OF COMPUTATIONAL CHEMISTRY 167

MAPLE ET AL.

C-C-C-H 1
C-H bond lengths, angstroms C-C bond lengths, angstroms

FIGURE 3. Distribution of (left) C-H and (right) G - C

bond lengths in structures used to sample the ab initio
potential energy surface. For the case of G - H bond
lengths in the 1.07-1.09 A, the histogram bar was
truncated because the 647 G - H bonds in this range
exceeded the dimensions of the ordinate. Two c-c-c.c 1
histogram bars for C-C bond lengths were similarly
truncated.

been distorted and have ab initio energies ranging

up to 20.47 kcalimol relative to the minimum.
The final distorted structures were checked to
verify that the internal coordinates were well sam- Totslon angler.

pled within a wide range to assure that an adequate

FIGURE 5. Distribution of (top) H-G-C-H, (middle)
portion of the energy surface is explored for pre-
G-C-G-H, and (bottom) C - G - C - 4 torsion angles
dictive purposes. The degree of sampling is de- in structures used to sample the ab initio potential
picted in terms of histograms in Figures 3 through energy surface.
6. From these figures, we see that indeed a large
range of distorted geometries are included. Cer-
tainly, any distortion we would need to calculate
the energetics for a molecular dynamics or me-
chanics calculation is spanned by these distorted
configurations. The C-H and C--C bond lengths
are sampled over a range of 0.96 to roughly 1.24 A
for C-H bond lengths and 1.38 to 1.66 A for C--C

50
t?c

40 300

6 250
30 C
8 200
20
8 150
"-
~

10 :'
0 0
r r - r r r r
r.-rrrr rrrrr

HC-H angles. C-C-H angles, C-C-C angles, Nonbond dislaancas, angstroms

FIGURE 4. Distribution of (left) H-C-H, (middle) FIGURE 6. Distribution of (top) H...H, (lower left)
G - G - H , and (right) C--C-C bond angles in C...H, and (lower right) C..C distances between
structures used to sample the ab inifio potential energy nonbonded atoms in structures used to sample the ab
surface. initio potential energy surface.

168 VOL. 15, NO. 2

 DERIVATION OF CLASS II FORCE FIELDS

bond lengths. Likewise, the H-C-H angle is sam- structure energy divided by the number of degrees
pled over a range of almost 50°, and the C-C-H of freedom (n,) is shown in the third column. Both
angle from 84 to 132". One can also see that these quantities depend on the extent of distortion of
are roughly normal curves and most of the inter- the structure, and the latter quantity is independent
nals are sampled at about what would be expected of molecular size. The average maximum energy
in their equilibrium nonperturbed structures. Fig- per degree of freedom is about 1.3 kcalimol for
ure 3, for example, shows that the final set of dis- these molecules. The highest value is 2.3 kcal/mol
torted structures for 16 molecules contains 647 for methane.
C-H bond lengths in the range of 1.07-1.09 A,
while approximately 270 C-H bond lengths were
in the range 1.09-1.11 A within this normal dis- QMFF Functional Form
tribution. We note that the torsion angles were
sampled in a more uniform distribution to make The QMFF functional form was derived by mod-
sure the entire 360" range was covered and that the ifying individual terms in the potential energy
C-C-C angles have peaks corresponding to the function and observing the effect on the fit of quan-
cyclopropane and cyclobutane structures. The tum mechanical data and the accuracy with which
small, strained cyclic alkanes were added to the ab initio equilibrium structures, energies, and en-
ergy derivatives were reproduced. The modifica-
sample set to enable derivation of a single force
tions in the analytic representations of the intra-
field that would predict the properties of both un-
molecular force field involved changes in the form
strained and highly strained alkanes with the same
of individual energy terms (e.g., harmonic vs. an-
set of parameters and functional form. As Figure
harmonic bond stretching) or the addition of new
5 demonstrates, torsion angles were well sampled
types of couplings or cross-terms (e.g., angleitor-
throughout the entire 360" range. Especially in the
sion coupling interactions). The effects of these
case of the C-C-C-C torsion angles, rigid ro-
functional form modifications were observed for
tations and energy minimizations at constrained
fits of quantum mechanical data, and the functional
torsion angles were used to improve the sampling
form alterations that result in the greatest im-
of torsion coordinates. Finally, the distances of
provements in the fits of the ab initio potential en-
closest approach between nonbonded atoms (at-
ergy surface were retained in the QMFF functional
oms separated by at least three bonds) were limited
form.
to 1.6, 2.0, and 2.6 for H...H, H...C, and C...C
atom pairs, respectively.
Potential energy surface sampling was com- TABLE 1.
Maximum Relative Energies (kcal/mol) of Distorted
pleted by calculating the ab initio energy and the Structures (Emax)and the Maximum Relative
first and second derivatives of the Hartree-Fock 6- Energies Divided by the Number of Degrees of
31G* energy with respect to the Cartesian coordi- Freedom (nu)Describing the Hartree-Fock Energy
nates. GRADSCF was used for the energy deriv- Surface of Hydrocarbons.
ative calculations on distorted structures of meth-
Molecule Emax 6naX/n
ane (81, ethane (421, propane (21), n-butane (36), ~ ~~

isobutane (171, n-pentane (12), isopentane (8), Methane 20.5 2.3

neopentane (6), cyclopropane (8), cyclobutane Ethane 31.4 1.7
(lo), cyclopentane (14), cyclohexane (9), methyl- Propane 55.4 2.1
cyclobutane (7), methylcyclopropane (7),1,2-di- n-Butane 67.3 1.9
methylcyclopropane (6), and 1,l-dimethylcyclo- n-Pentane 28.0 0.6
lsobutane 68.0 1.9
propane ( 3 ) . The numbers in parentheses indicate
lsopentane 20.6 0.5
the total numbers of structures, including ab initio Neopentane 21.4 0.5
equilibrium and transition state structures. The re- Cyclopentane 70.1 1.8
sulting total number of independent quantum me- Cyclohexane 35.9 0.7
chanical observables, including relative energies Cyclobutane 21.8 0.7
and first and second derivatives, is 128,376. Methylcyclobutane 61.5 1.6
The maximum energies, relative to the ab initio Cyclopropane 21.3 1 .o
global minimum, of the distorted structures of the Methylcyclopropane 41.6 1.4
molecules in the training set are listed in the second 1,2-Dimethylcyclopropane 52.7 1.1
column of Table I, and the maximum distorted- 1,l -Dimethylcyclopropane 48.2 1 .o

JOURNAL OF COMPUTATIONAL CHEMISTRY 169

MAPLE ET AL.

The individual terms of the QMFF functional As noted in the following sections, the cubic and
form are described below and compared, where quartic term in eq. (3) represent anharmonicity in
appropriate, to previous force fields often em- the angle bending. This anharmonicity was man-
ployed in molecular mechanics and dynamics, par- ifested by the energy surface of the highly strained
ticularly AMBER,24CHARMM,25and CVFF,12-15as molecules with three- and four-memberd rings.
well as MM219,20,35 and MM3.21,35 For this compar-
ison, it is useful to classify these various force fields TOHSION
in terms of their functional forms. Thus, AMBER
and CHARMM, which are simple diagonal quad- The potential energy function for torsion angles
ratic force fields, are termed ”Class I.” CVFF and 4 is given by the three-term Fourier expansion
MM2 contain some limited coupling terms and an- €6 = 1K4(l - cos 4) + 2K,(1 - cos 24)
harmonicity but are still necessarily restricted in
their functional forms. MM3 contains extensive an- + 2K6(1 - cos 24) (4)
harmonicity and coupling. MM3 and CFF936(based All three cosine terms are needed to fit the con-
on QMFF), which contain a comprehensive set of formational energies of hydrocarbon^.^^,^^ The pe-
anharmonic and coupling terms, will be referred riodic cosine function is employed by essentially
to as “Class 11.” The detailed definition of Class I1 all common types of force fields, although in most
force fields is given in Hwang et aL6 of them two of the three terms in the expansion
are generally omitted. Although the phases in the
B O N D STKtiTCHING first and third terms in eq. (4) differ from the MM3
formulation,21both energy functions are equivalent
In the resulting QMFF force field, the diagonal and differ only by a constant term that contributes
potential energy function for bonds is the fourth- to absolute, but not relative, energies of a given
degree polynomial molecule.
El, = 2&(b - b0)2 + 3Kb(b - b ~ +) 4&,(b
~ bo)4
-
(2) NONBOND
where b is the bond length and bo is the reference For nonbond interactions between atoms i and
value. By contrast, the AMBER and CHARMM j (separated by at least three bonds and by a dis-
force fields employ only the first, the simple har- tance Y),the QMFF nonbond potential energy func-
monic, term, and MM2 only the first two. CVFF tion is an inverse-power expression:
uses a Morse potential, which is roughly equivalent
but computationally less efficient. It is well known E,, = 9,9,/r + ~ [ 2 ( r * / r-) ~3(r*/r)6] (5)
that bond energies are highly a n h a r m ~ n i c ’and
~,~~ An inverse 9 power was found to be necessary
that anharmonicity is well represented by cubic and to fit both intramolecular and intermolecular prop-
quartic term^.^^,^^ The quartic term, which is also erties simultaneously including rotational barriers
employed in the MM3 force field, prevents bonds and fits to crystal structures. A similar suggestion
from breaking during molecular mechanics or dy- was hypothesized by Lifson some years ago, sug-
namics calculations, which can occur when only gesting that the 12th-power repulsion was too
the cubic term is used. ”hard.” More recently, Allinger has also concluded
that the simple 12th-power repulsion is counter-
ANGLti BtiNDING indicated by the combination of intra- and inter-
molecular properties and has used a softer expo-
A fourth-degree polynomial
nential term in MM3. Class I force fields including
E , = 2Ko(8 - 80)2 + 3Ko(8 - 80)3 CVFF, AMBER, and CHARMM, of course, all use
inverse 12th-power repulsions. We also employ a
+ 4K,(8 - 80)4 (3)
recently introduced set of combination rules for
is also used to describe the anharmonicity in the relating the individual atomic parameters to the
bond angle (8) energy function. Again, AMBER24 pairwise potential3?
and CHAIUVIM~~ rely on only a simple harmonic
approximation, as does CVFF, while MM2 has har- Y* = [(rT6 + (6)
monic and sixth-order terms. This functional form and
is still substantially less elaborate than the sixth-
degree polynomial used in the MM3 formulation. E = + r:6]
(~~&,)”~2(r:r:)~/[r:~ (7)

170 VOL. 15, NO. 2

 DERIVATION OF CLASS I1 FORCE FIELDS

Here, q,, F ; , and r: are the partial charge, well BOND/BOND COUPLING
depth, and diameter of atom i. Thus, the potential
function in eq. (5) uses Coulombs law to approx- The QMFF functional form for bondibond cross-
imate the electrostatic interaction, an inverse 9th- terms is
power repulsive interaction between atoms, and Kbi,'(b - bo)(b' - bi)
Ebb' = (10)
an attractive dispersion interaction. The new com-
bination rules in eqs. (6)and (7) have been shown and is identical in form to the cross-term used by
to reproduce experimental rare gas data to a much CVFF.14,15The inclusion of this bond/bond cou-
better extent than more commonly used combi- pling interaction, which is neglected by MM3, is
nation rules (e.g., arithmetic and geometric means) necessary for a force field to reproduce vibrational
used in previous force fields.38 frequencies with high accuracy. 15,26*27)37Bondlbond
coupling was revealed from fits of ab initio second
BOND INCKEMENTS derivatives of the energy with respect to Cartesian
coordinates, as well as from calculations of the sec-
The partial atomic charges in eq. (5) are com- ond derivatives of the energy with respect to bond
puted from the bond incremental charges or bond lengths. Although this coupling interaction is sig-
increments. The bond increment 6,, is defined as nificant in alkanes, the second derivatives in other
the partial charge contributed by atom j to atom species reveal that the coupling is particularly
i.39 In general, each of the n atoms bonded to an strong for bonds in systems with delocalized elec-
atom contributes to the atomic charge. By com- trons (e.g., the CO/CN bond coupling in amides).
puting a partial atomic charge for atom i from the Following the convention used by CVFF,'4f's
summation over the bond increments only bonds with an atom in common are coupled,
c 4,
I1
as depicted in Figure 7, which illustrates bond/
91 = (8) bond interactions (Fig. 7a), as well as the other
/=1
cross-term interactions used in the QMFF force
which depend on the relative electronegativity of field. The neglect of coupling of bonds without an
the bonded atoms, the dependence of the partial atom in common is ordinarily a good approxima-
charge on the bonding environment is automati- tion.*" During the derivation of QMFF for other
cally taken into account. Consequently, bond in- species from fits of ab inifio potential energy sur-
crements can be more transferable than partial faces, the only significant exception found for this
atomic charges.39An additional advantage of using approximation was for bonds in aromatic ring mol-
bond increments is that charge conservation is au- ecules.
tomatic. When charges are directly transferred to
another molecule, the sum of partial atomic charges
BOND/ANGLE COUPLING
in that molecule is, in general, not equal to the
charge of the molecule (i.e., zero). However, this As depicted in Figure 7b, angles are coupled to
problem does not occur when bond increments are the two bonds that comprise the bond angle by the
used to compute the atomic charges. energy function
The first atom in the bond increment subscript
is defined as the charge acceptor and the second Et,, = K d b - bo)(e - 00) (11)
as the charge donor. This convention is expressed There are two of these terms for each angle (i.e.,
by the relationship one for each bond). Bondiangle coupling is nec-
s,, = - 6I' (9) essary to reproduce f r e q u e n c i e ~ ' ~ and
, ~ ~several
,~~,~~
between all bonded atom pairs i and j . Equation structural features. For example, in tri-tert-butyl-
(9) states that the charge donated by atom j to atom methane,14 as well as in other molecules, bond/
i is opposite in sign to that donated by i to j . Con- angle coupling is needed to reproduce the bond
sequently, only a single bond increment parameter lengthening that is observed when the bond angle
is needed to describe each unique pair of bonded is compressed in strained molecules, as well as to
atoms. For alkanes, this means that the ScH bond account for the concomitant shifts in vibrational
increment is the only adjustable electrostatic pa- frequencies.
rameter. The SCc bond increment is fixed at zero The neglect of coupling of an angle to other
because two carbon atoms bonded to each other bonds (i.e., bonds other than the two that define
have the same tendency to donate or accept elec- the bond angle) is ordinarily a good approxima-
trons. tion,40and this has been verified by fits of ab inifio

JOURNAL OF COMPUTATIONAL CHEMISTRY 171

MAPLE ET AL.

imental vibrational frequencies for coupled bend-

a. bond/bond
ing modes. However, it has been suggested that
this term is relatively unimportant for bond angles
in tetracoordinate centers.40The results obtained
from the fit of the quantum energy surface, es-
b. bond/angle pecially for the relative energies and forces in
molecules with large amounts of angle strain
(i.e., cyclopropane and methyl derivatives of
cyclopropane), support the importance of the an-
t
gle/angle coupling interactions.
c. angle/angle

ANGLE/ANGLti/TOKSION COUPLING
The QMFF functional form for the angle/angle/
torsion cross-term is
d. angle/angle/torsiOn
E,,~, = K,,,~,(B - e,)(e' - e;)cos 4
(13)
As indicated in Figure 7d, this term describes the
coupling between a torsion angle A-B-C-D and the
two vicinal bond angles A-B-C and B-C-D that share
e. angle/ torsion the common bond B-C. This coupling term,
which was first introduced by Lifson and War-
~ h e l , is
~ *used by CVFF to more accurately repro-
duce vibrational frequen~ies,'~-'~ but it is not used
f. bond/torsion (type1) by MM3. It is found that this coupling interaction
makes a large contribution to the ab initio relative
energies, forces, and coupling derivatives in all al-
kanes (except methane) examined here (see Table
I). It is most notable in cyclobutanes, wherein it
accounts for a significant proportion of the ring
pucker and inversion barrier. Similar results were
noted from fits of ab initio data of all other functional
FIGURE 7. Representationof cross-terms (coupling
groups that have been examined. Thus, informa-
interactions) in QMFF.
tion implicit in the quantum energy surface sup-
ports Lifson and Warshel's conclusion as to the
potential energy surfaces of a variety of functional importance of this term.
groups. This approximation, as well as the func-
tional form in eq. (ll),were also used in CVFF'4,'5
BOND/'I'OKSION AND
as well as MM2 and MM3.21,37 None of the cross- ANGLE/IOHSION COUPLING
terms were present in the CHARMM and AMBER
force fields, which are diagonal quadratic func- The bond/ torsion and angle/torsion interactions
tional forms. Thus, consideration of these force shown in Figures 7e-7g are expressed as the dis-
fields is omitted in discussing all subsequent cou- placement of a bond length or bond angle multi-
pling interactions. plied by a Fourier expansion in the torsion coor-
dinate, as indicated in eqs. (14) and (15).
ANGLti/ANGLE COUPLING E, = (b - bO)['K, cos + ' K d , cos 24
Bond angles are coupled to other bond angles + 3 K d , cos 341 (14)
that share a common bond, as illustrated in Figure Em, = ( 0 - OO)[lK,, cos 4 + 'K,, cos 24
7c, and the QMFF functional form is given by
+ 3Kd, cos 341 (15)
E,"~ = K H H I (-e eo)(ef - e;)
(12) As shown in Figure 7e, there are two angle/torsion
This interaction has been used in the CVFF14,'5and interactions for every torsion angle 4. There are
MM3'I force fields to accurately reproduce exper- also two types of bond/torsion couplings. The first

172 VOL. 15, NO. 2

 DERIVATION OF CLASS It FORCE FIELDS

type (in Fig. 7f) involves the coupling of the central The restoring force on the bond is therefore given
bond to the torsion angle, while the second type by
(in Fig. 7g) couples the peripheral bond to the tor-
- d E l d b = - k - k cos 34 (19)
sion angle, There are two bonditorsion interactions
of the second type for every torsion angle because Thus, the first term in eq. (18) is a linear term that
there are two peripheral bonds. The CVFF force is effectively added to the polynomial in eq. (2).
field does not employ angle/torsion coupling or This linear term introduces a force [i.e., the first
bonditorsion coupling, and MM3 uses only the term in eq. 19)] on the bond length that is indepen-
first of these bond/ torsion interactions. dent of the torsion angle, thereby driving the bond
Differentiation of eq. (14) with respect to the length away from the idealized (strain free) bond
bond length gives length (i.e., bo). The net result of the use of eq. (17)
is that k and bo are correlated, and the transfera-
dEd,ldb = ' K d , cos 4 + 2 K d , cos 2 4 bility of these two parameters is reduced. This
+ 3 K d , cos 34 (16) problem is eliminated by the use of the functional
which shows that this interaction describes the de- form in eqs. (14) and (15).
pendence of the bond length restoring force on the The full functional form that describes the quan-
torsion coordinate 4. Similarly, the angle/ torsion tum energy surface is given by the sum of the in-
coupling term is used to determine the torsion an- dividual terms in eqs. (2)-(15). The result is
gle dependence of the bond angle restoring force, E = 2 [2Kr,(b - b0)2 + 3Kb(b b0)3-
that is, these interactions are important for describ- 1,
ing the structure of molecules in which different + '&(b - bo)'] + 2 [ 2 K 8 ( 8 0")'
-
H
conformers exhibit significant differences between
bond lengths or bond angles. The bonditorsion + 3Kn(8 - 80)3+ 4K,(8 - 80)4]
interaction is reflected in the bond stretching that + 2 ['K,(l - cos 4) + 2Kh(1
h
occurs with eclipsed bonds (i.e., cis torsion angles)
in structures such as eclipsed ethane. The fits of
the alkane nb initio energy surface (particularly the
energies and forces) reveals the existence of the
bond/ torsion and angle/ torsion coupling interac-
tions, both for alkanes as well as acetals, carbo-
hydrates, amides, and other functional groups.
The bond/ torsion and angle/torsion coupling
terms are found to reproduce the anomeric effect
in acetals and carbohydrates or other functional
groups in which a carbon atom separates two elec-
tronegative atoms.42It is intriguing that the same
bond i torsion and angle/ torsion coupling interac-
tions revealed in the hydrocarbon energy surface
may explain (or cause) the anomeric effect in car-
bohydrates. The implications of this are the subject
of further study.
The bonditorsion functional form used by MM3
is given by
E = k(b - bo)[l + cos 341 (17)
and differs from eq. (14). One difference is that the
MM3 form lacks the first two cosine terms in the
Fourier expansion. Another important difference is
the linear term in the MM3 form in eq. (17). This TKANSFEHASILI'l'Y OF PAHAME'I'ERS
term, which is independent of the torsion angle,
is more clearly seen by rewriting eq. (17) as The QMFF force field uses the conventional ap-
proach of assigning parameters based on atom
E = k(b - bo) + k(b - bo)COS 34 (18) types. The parameters for two internal coordinates

JOURN,AL OF COMPUTATIONAL CHEMISTRY 173

MAPLE ET AL.

or cross-terms are distinguished if the atom types ergy function and from quantum mechanics (i.e.,
that define the internal coordinate or coupling in- the “observables”), respectively. The relative ener-
teraction differ. For alkanes, all hydrogen and car- gies for each distorted structure of a molecule, b,
bon atoms have atom types of ” h ’ and “c,” re- are determined relative to the structure with the
spectively, irrespective of whether they belong to lowest energy. For the summation over energy de-
methyl, methylene, or methine groups. Parameters rivatives in eq. (21), Nh is the number of atoms in
for methyl, methylene, and methine groups are the molecule b, and xi is one of the three Nf,Cartesian
same. This is also the case for the CVFF,14,15 AM- atomic coordinates.
BER,24 and CHARMMZ5force fields. In contrast,
the MM3*I force field uses different bond angle WEIGHTING FACTOKS
reference values for these groups. Further, the
The overall weighting factor WGf,in eq. (21) al-
same atom types were used in three-, four-, and
lows all data for configuration u of molecule b to
five-membered rings, unlike the MM2 and MM3
be weighted independently. One might, for ex-
force which employ four different sets
ample, want to weight minima and barriers more
of atom types for alkanes (i.e., one each for three-, heavily than arbitrarily distorted configurations in
four-, and five-membered rings plus one more for
some cases. For the work presented here, this
all other alkanes in addition to the above-men- weighting factor was set to 1.
tioned distinction between methyl, methylene, and The weighting factor for the energies Wail,was
methine groups). QMFF has been derived with the
set equal to a value of 15,000. Larger values gave
goal of obtaining the most accurate and physically
significantly improved fits of the distorted struc-
reasonable energy function and, hence, of maxi-
ture energies. However, the predictions of the en-
mizing the force-field transferability. This was ergy of distorted structures not included in the fit
found to result in a minimal number of atom types often degraded when larger energy weighting fac-
(i.e., one set for all strained and unstrained al-
tors were employed. Thus, a value of 15,000 was
kanes). about the maximum energy weighting factor that
could be employed without overfitting the energies
Least-Squares Fit of the of the distorted structures and thereby reducing
Energy Surface the transferability of the parameters.
The weighting factor for the first derivatives,
The force constants in eq. (20) were determined Wlnfl,in eq. (21), was set to 100. Parameterization
with PROBE (Biosym Technologies), a program with values of Wlnflranging from 50 to 200 gave
that uses a least-squares method (with a Leven- comparable predictions of the ab initio equilibrium
berg-Marquardt minimizer43)to adjust the poten- structures and vibrational frequencies. Wlof,values
tial energy function parameters to fit the energy less than 50 resulted in slight improvements in the
and the first and second derivatives of the energy fit of the curvature of the quantum energy surface
(with respect to Cartesian atomic coordinates) of (d2Eldxrdx,)and vibrational frequencies, but were
one or more molecules and structures. The accompanied by a significant degradation in first
weighted sum of squared deviations (i.e., the ob- derivative and structure predictions. For values of
ject function to be minimized) was defined as the weighting factor (W,,,,,) greater than 200, minor
improvements in the fits of first derivatives were
accompanied by large errors in the fit of second
derivatives of the energy and, consequently, vi-
brational frequency predictions.
The weighting factor for the second derivatives
of the energy is W2f1ah. WZllnhwas set to 1 for all
second derivatives of the energy, with the excep-
tion of 1-4 second derivatives, for which a value
of 75 was used. This large value was used to com-
where Eat, is the potential energy calculated from pensate for the small size of 1-4 second derivatives
eq. (20), and Enhis the corresponding quantum me- relative to 1-2 and 1-3 second derivatives, which
chanical energy of the distorted configuration, a, are defined as second derivatives with respect to
of molecule b. Eib and ELb are the relative energies the Cartesian coordinates of atom pairs separated
of configurations determined from the potential en- by one or two bonds, respectively. In experiments

174 VOL. 15, NO. 2

 DERIVATION OF CLASS II FORCE FIELDS

experiments in which ab initio data for ethane, pro- have the same energy and the energy function
pane, and n-butane were fit separately, substan- must be symmetric about these 0 and 180" values
tially different results for the torsion parameters in to reflect this (see Scheme 1). As a result, there
eq. (5) were determined when WZlrnb was less than were only two independently adjustable parame-
75, but they remained constant when W2,,"h was ters (i.e., a reference value and a quadratic force
increased above a value of 75. For small-ring mol- constant) for each of the H-C-H, H-C-C, and
ecules, some 1-4 second derivatives are simulta- C-C-C diagonal energy functions. This procedure
neously 1-3 (e.g., some carbodhydrogen atom also ensures that the polynomial in eq. (3) has only
pairs in cyclopropane) or 1-2 (e.g., bonded carbon a single minimum in the 0-180" range, as required
atom pairs in cyclobutane are the 1-4 terminal at- to prevent unacceptably large or small bond angles
oms of a C-C-C-C torsion). These second deriva- during a structure optimization or molecular dy-
tives were defined by PROBE as 1-2 or 1-3 deriv- namics simulation. Quadratic force fields (e.g.,
atives rather than 1-4, and were given a weighting CVFF, CHARMM, and AMBER) do not have this
factor of 1because their values were already large. correct symmetry-required limiting behavior in the
bond angle energy function, nor is this limiting
behavior obeyed by the MM3 bending energy func-
OKIGIk OF NONBONI) PAKAMETEKS
tion.
Although charges (or bond increments) and re-
pulsive constants could in principle be determined CONSTKAINI'S ON TOKSION PAKAMtiI'tiKS
from fits of the nb initio data, a pragmatic approach
that served to enhance transferability was em- In addition, we tested whether the ' K , term for
ployed in which these parameters were obtained the torsion potential in eq. (4)could be constrained
from experimental hydrocarbon crystal structures to be the same for the H-C-C-H, H-C-C-C, and
and sublimation energies.6These parameters were C-C-C-C torsions. Constraints on the 2 K , and 3 K ,
transferred (see Table 11) and held fixed in all least- force constants were similarly tested. Thus, only
square fits of the nb initio data. The derivation of three adjustable torsion parameters were used in
these parameters will be the subject of a later study. the fit of the quantum mechanical data. The ra-
tionale for testing this constraint is to ask whether,
to a first approximation, the potential for twisting
CONS'I'RNN'IS ON BOND ANGLE
about a bond depends primarily on the identity of
ANHAKMONIC FORCE CONSTANTS
the two atoms in the bond rather than on the pe-
The cubic and quartic force constants for the ripheral 1-4 atoms. This would be desirable as it
bond angle energy function were constrained so would maximize the transferability of the torsion
that the function in eq. (3) has a slope of 0 at angles parameters so that they could be applied to the
of 0 and 180".This is, in fact, required by symmetry. large number of possible X-C-C-Y torsions, where
This can be seen by considering either a 0 or 180" X and Y are any atoms. Only a small degradation
angle. Deformations in either direction by say 25" in the quality of fit accompanied this reduction in
produce the same angle between the bonds, i.e., the number of torsion parameters from 9 to 3, and
175 or 185" are equivalent (both 175"). Thus, they thus we adopted this approximation.

SCHEME 1. Symmetry of angle function about 180".

JOURNAL OF COMPUTATIONAL CHEMISTRY 175

MAPLE ET AL.

lists the relative deviations of energes and energy

Parameterization derivatives between those calculated by the QMFF
A two-step procedure for parameterizing the ab and by ab initio. For comparison, Table I11 also gives
initio data was employed. In the first step, data deviations between ab initio values and those ob-
from molecules with a three- or four-membered tained with a harmonic diagonal (Class I) force field,
ring was excluded from the fitting procedure. In which was derived in a manner similar to that of
the second step, the data from the three- and four- QMFF (i.e., by eliminating all anharmonic and cross-
membered rings was added to the parameteriza- terms and refitting the harmonic diagonal force con-
tion data set, and the parameters in the diagonal stants to the same ab initio data). As noted in the
bond, bond angle, and torsion angle energy func- discussion of the interaction terms, the quadratic or
tions were frozen at the values determined from harmonic diagonal force field is representative of the
the first step, while all other parameters were re- standard force fields used in the simulation of large
determined. (This procedure reduced systematic organic and biomolecular systems such as AMBER,
errors in predicted vibrational frequencies.) CHARMM, etc. CVFF, as noted above, contains
Two statistical tests of the parameters were ap- some coupling terms and anharmonicity, and lies
plied. Parameter uncertainties were determined from midway between Class I and Class I1 functional
standard deviations, which were calculated from the forms in this respect. The second and third columns
diagonal elements of the variance-covariance ma- of Table 111 give the rms deviations for the relative
trix.",@ In addition, parameters were subjected to a energies of the harmonic and QMFF force fields. For
t-test to check the hypothesis that they were statis- the Class I function, the deviations range from 0.37
tically identical to zero.@If the probability that this kcal/mol for n-pentane to 6.33 kcal/mol for 1,l-di-
hypothesis was true exceeded 0.05 (at the 95% con- methylcyclopropane, while for QMFF the deviations
fidence level), then the hypothesis was accepted and range from 0.26 kcalimol for the distorted structures
the parameter was fixed at a value of zero during a of methane to 2.87 kcal/mol for methylcyclopropane.
subsequent parameter refinement. It turned out that The average rms energy deviation, which is given in
this hypothesis was rejected for all adjustable param- the last row in Table 111, is 0.87 and 2.69 kcal/mol for
eters, and so no force constants were set equal to QMFF and the harmonic diagonal force fields, re-
zero in the QMFF alkane force field. spectively. This demonstrates that on average the
To summarize the fit of the alkane quantum me- inclusion of anharmonicity and cross-term interac-
chanical data, after excluding the 4 nonbond param- tions reduces the errors in calculated energies by 1.8
eters (determined separately from fits of experimen- kcalimol. The fourth and fifth columns list the same
tal crystal structures and sublimation energies), the rms deviations as percentages of the maximum rel-
78 remaining force constants and reference values ative energy of the distorted structures (listed in Ta-
listed in Table I1 were determined. Only 66 of these ble I). The percentage error in the harmonic force
parameters were independently adjustable, as de- field energes range from 1.3% for n-pentane to
scribed above, because the values of 6 torsion pa- 17.0%for cyclobutane, while the percentage error for
rameters and 6 anharmonic bond angle parameters QMFF energies ranges from 0.8% for n-butane to
were determined from symmetry relationships or as- 6.9% for methylcyclopropane. The average percent-
sumed transferabdity. Thus, because a total of age error over all structures is 2.2 and 6.9% for QMFF
128,376 independent ab initio observables (i.e., rela- and the harmonic diagonal force field, respectively.
tive energies and the first and second derivatives of The percent rms deviations between the force-
the energies) were used an observable-to-parameter field and ab initio first and second energy derivatives
ratio of 1945:1 was achieved for the fit of the data. were calculated from the ratio between the sum of
squared deviations and the sum of squared ab initio
derivatives. The first derivative errors range from
6.2% for n-butane to 19.5% for methylcyclopropane
Quality of Fit to Ab fnitio Potential for QMFF while the corresponding range is 15.2%
Energy Surface for methane to 65.7% for cyclopropane for the har-
monic force field. Similarly, for QMFF the second
RELATIVE ENERGIES AND derivative errors range from 2.6% for ethane to 7.4%
ENERGY DERIVATIVES
for methylcyclopropane, while for the harmonic
The quality of fit of the QMFF to the ab initio po- force field these deviations range from 12.8% for
tential energy surface is indicated in Table 111, which ethane to 24.9% for methylcyclobutane. On the

VOL. 15, NO. 2

176
 DERIVATION OF CLASS II FORCE FIELDS

TABLE ti.
Force Constants and Reference Values in QMFF for Alkanes.

'Kb 3Kb 4Kb

Bond bo (4 (kcal mol-' A-*) (kcal mol-' k3) (kcal mol A 4,

H-C 1.0845 417.4 - 851.3 1040.0

c-c 1.5297 340.2 - 586.1 758.0

8. Bond Angle: E = 'K,(O - 8,)' + 3K,(8 - 8,)' + 4K,(0 - 80)4

60 'Kn 3Kn 4K,

Angle (") (kcal mot-' rad-') (kcal mot rad-3) (kcal mol rad ' 4,

H-C-H 107.7 51.5 - 8.9 - 10.8

H-C-C 110.8 52.7 - 10.9 -11.3
c-c-c 112.9 52.2 -12.1 -11.3

C. Torsion Angle: E = 'Kd'(l - cos 4) + 'K,(l - cOS 24) + 3K,(1 - COS 34)

Kd' 2K*, 3K,,,

Torsion (kcal mol-') (kcal mol - I ) (kcal mol - l )
H-C-C-H -1.152 0.012 -0.179
H-C-C-C -1.152 0.012 - 0.179
c-c-c-c -1.152 0.012 -0.179

D. van der Waals Interaction: E = ~ [ 2 ( r ' / r )-~ 3(r*/r)6] where r* = [(r,*6 + r,%)/2]1'6and E = (&,€,)I 22(r,*r;)3/(r,*6 +
r/*6)

6
Atom (kcal mot-')

H 2.995 0.020
C 4.010 0.054
E. Bond Increment: E = 332.0 qjqj/rij where qi = &6,k

k
Bond (electrons)
H-C 0.053
c-c 0.000

~ ~

Kbb'
BondIBond (kcal m o l ~A - ')

H-CIH-C 11.769
H-CIC-C 11.310
c-CIC-c 10.062
G. BondIAngle: E = Kb,(b - bo)(O - 00)
Kbn
Bond/Angle (kcal mol-' A-I rad-')
~ ~

H-C/H-C-H 22.7
H-CIH-C-C 13.9
C-CIH-C-C 34.5
c-CIC-c-c 18.3

JOURNAL OF COMPUTATIONAL CHEMISTRY 177

MAPLE ET AL.

TABLE II.
H. AngleIAngle: E = K,,,(B - 0,)(0' - 0;)

KMil
AngleIAngle Common Bond (kcal mol ~ rad ~ 2,

H-C-HIH-C-H H-C 0.93

H-C-HIH-C-C H-C 1.57
H-C-CIH-C-C H-C 3.12
H-C-CI H-C-C c-c -2.18
H-C-CIC-C-C c-c - 4.62
c-c-cIc-c-c c-c - 8.97

I. AngleIAngleITorsion: E = K,,,,(8 - 0,)(0' - 0;)cos 4

KUHr6

AngleIAngleITorsion (kcal mol-I rad-*)

H-C-CIH-C-CIH-C-C-H - 14.8
H-C-CIC-C-CIH-C-C-C - 19.2
c-c-cIc-c-cIc-c-c-c - 35.6

J. BondITorsion (type 1): E = (b - b#K& cos 4 + 2K,,b cos 24 + 3K& cos 341
K& 2K& 3Kd*
BondITorsion (kcal mol-l k') (kcal mol-I k ' ) (kcal mol- 8, ~ I)

C-CIH-C-C-H - 52.239 - 0.784 - 0.901

C-CIH-C-C-C - 49.1 11 -3.173 -0.393
c-CIC-c-c-c - 47.281 - 5.641 0.315

K. BondITorsion (type 2): E = (b' - bA)['Kdb,cos 4 + 2Kdb,cos 2 4 + 3K&, cos 341

2K&, 3K&,
Bond/Torsion (kcal mo1-I k1) (kcal mol--Ik l )

H-CIH-C-C-H 0.778 0.476 0.079

H-CIH-C-C-C 0.779 0.288 0.398
C-CIH-C-C-C 1.665 0.583 - 0.342
c-CIC-c-c-c 2.294 0.765 0.357

L. AngleITorsion: E = (0 - O,)[IK,,,, cos 4 + 2K6Hcos 2 4 + 3K6Hcos 341

KdH K" 3K6u
AngleITorsion (kcal mol-' rad-') (kcal mol-I rad-') (kcal mol ' rad-')
H-C-CIH-C-C-H - 1.622 0.739 - 0.449
H-C-CIH-C-C-C 0.260 0.702 -0.163
C-C-CIH-C-C-C - 1.886 0.101 - 0.295
c-c-cIc-c-c-c - 0.030 - 0.015 - 0.003

whole, the largest errors in the energy and in the in calculated forces and curvatures by factors of 3.3
first and second derivatives are found for methyl- and 4.0, respectively.
cyclopropane. The average first derivative deviations
for a11 molecules are 10.1 and 33.0% for QMFF and
tiFFti:CT OF SIKAIN ON QUALILY OF F l l
the harmonic force fields, while the average second
derivative deviations are 4.6 and 18.4% for-the same These results dramatically emphasize the large
two force fields. Thus, neglecting anharmonicity and degree of anharmonicity and intramolecular cou-
cross-term interactions causes an increase in the error pling interactions that exist among even the alkane

178 VOL. 15, NO. 2

 DERIVATION OF CLASS I t FORCE FIELDS

TABLE 111.
Ability of Diagonal Harmonic (Class I) Force Field and the Class II Functional Form to Account for the
Quantum Mechanical Energy Surface of Alkanes: rms Deviations from Ab Initio Energies and First and Second
Derivatives of Distorted Structures of Alkanes for a Harmonic Diagonal Force Field and for QMFF.
(kcal/mol)
AErmSa %AErmSb %AE:,,C ‘/‘A E:,,,,,d
Molecule Harmonic QMFF Harmonic QMFF Harmonic QMFF Harmonic QMFF

Methane 1.09 0.26 5.3 1.3 15.2 7.6 16.2 3.9

Ethane 1.24 0.40 3.9 1.3 26.8 7.7 12.8 2.6
Propane 2.47 1.54 4.5 2.8 26.0 13.5 14.3 3.6
n-Butane 2.55 0.56 3.8 0.8 28.6 6.2 22.0 3.2
n-Pentane 0.37 0.64 1.3 2.3 26.8 10.2 14.1 2.9
lsobutane 1.78 0.62 2.6 0.9 26.5 16.1 18.1 4.7
lsopentane 0.83 0.73 4.0 3.5 25.7 6.7 15.9 3.3
Neopentane 1.40 0.46 6.5 2.1 29.5 7.1 17.1 4.1
Cyclopentane 2.85 1.61 4.1 2.3 30.0 10.5 18.9 4.7
Cyclohexane 1.60 0.74 4.5 2.1 27.2 8.2 13.1 2.9
Cyclobutane 3.71 0.43 17.0 2.0 19.5 7.0 19.6 4.9
Methylcyclobutane 5.60 0.78 9.1 1.3 25.9 8.2 24.9 5.0
Cyclopropane 2.39 0.34 11.2 1.6 65.7 12.4 22.7 6.2
Methylcyclopropane 5.96 2.87 14.3 6.9 64.8 19.5 19.3 7.4
1,2-DimethyIcyclopropane 2.85 0.86 5.4 1.6 56.0 11.6 21.5 7.1
1 , l -Dimethylcyclopropane 6.33 1.09 13.1 2.3 43.5 9.7 24.5 7.3
Average 2.69 0.87 6.9 2.2 33.0 10.1 18.4 4.6
arms deviation of force-field relative energies from ab initio relative energies.
bPercent (rms) deviation of force-field relative energies, calculated by dividing A€,,, by the maximum relative energy of distortion
(listed in Table I).
“Percent (rms) deviation of force-field first derivatives of the energy.
‘Percent (rms) deviation of force-field second derivatives of the energy.

molecules, and in the alkane energy surface. The gree to which these terms contribute to the fit of
inclusion of terms to represent this anharmonicity the relative energies and the first and second de-
and the coupling interactions greatly increases the rivatives of the energy in the sum of squared de-
ability of the analytic form to account for the energy viations in eq. (21) (Table IV). This assessment is
surface of these alkanes and the physical and dy- based on the extent of increases in deviations from
namic properties that we wish to calculate. It be- ab initio energies and energy derivatives as selected
comes clear that this allows us to achieve a signif- terms in the energy function are omitted. The sec-
icant improvement in our ability to simulate ond column in Table IV shows the reduction in the
molecular properties of interest as they all derive number of adjustable parameters that accompanies
from the molecular energy surface. leaving out the interactions listed in the first col-
umn. The third column presents the rms deviations
in relative energies for all 16 molecules used for
Use of the Modeled Bnerm Surface the fit. Similarly, the fourth and fifth columns are
to Assess Importance of the relative rms deviations in first and second de-
Individual Terms rivatives, respectively. In all cases, the relative de-
viations were computed from the ratio of squared
It is of interest to assess the importance of deviations to squared energy or derivative values.
individual coupling interactions and of anhar- The importance of individual terms varies
monicity in determining the overall shape of the widely between differing molecules, with specific
quantum energy surface. This can be done terms such as cross-terms often being important.
quantitatively by making use of the analytic energy Thus, for example, the angle-angle-torsion cross-
expression (force field) derived from this surface. term, which is neglected in other force fields, was
The importance of the individual energy func- shown by our analysis of the quantum mechanical
tion terms in eqs. (2)-(15) are assessed by the de- energy surface to reduce the maximum error in the

JOURNAL OF COMPUTATIONAL CHEMISTRY 179

MAPLE ET AL.

TABLE IV.
Percent rms Deviationsaof DerivativesbCalculated with Various Terms Omitted from the QMFF Functional
Form.
Missing Reduction of
Energy Terms Parameters AEC %AE' Yo A €
~~ ~

None 0 0.96 10.0 4.2

Anharmonic bonds 4 1.36 22.7 18.2
Anharmonic angles 6 0.95 10.5 4.8
Bondlbond 3 0.96 10.1 4.5
Bond/angle 4 1.11 21.4 12.3
Angleiangle 6 1.17 12.0 4.3
Torsion 3 3.50 10.2 4.2
Nonbonds 2 3.20 11.5 4.2
Angleiangleitorsion 3 1.85 16.8 6.4
Bonditorsion 1 9 1.69 30.7 4.1
Bondltorsion 2 12 0.93 11.4 4.4
Angle/ torsion 12 1.17 11.7 4.4
arms deviations averaged for all 16 molecules in the fit.
bE' = first derivative and €" = second derivatives; %A is the percent deviation.
CAbsoluterms energy deviations in kcalimol.

sampled energies from 5.3 kcalimol to 1.3 kcalimol SECOND DtiKIVA'I'IVES

in isobutane.
The most important cross-terms for fitting the
In the following, we will analyze the importance curvature of the energy surface are the bondiangle
of individual terms in the force field both as re- and angleiangleltorsion terms. From the last col-
vealed by their contributions to the first and second umn in Table IV, it can be seen that neglecting these
derivatives of the energy. two cross-terms causes the second derivative de-
FIRST DEKIVATIVES viations to increase from 4.2 to 12.3 and 6.4%, re-
spectively. For fitting ab initio second derivatives,
An examination of the first derivative deviations these two cross-terms are evidently less important
in column 4 of Table IV reveals the terms that make than bond anharmonicity but much more impor-
the largest contributions to intramolecular forces in tant than bond angle anharmonicity or other cou-
alkanes. The most important term for reducing the pling interactions. This suggests that these two
deviations in forces is the bonditorsion 1 term [see cross-term types should be particularly important
eq. (14) and Fig. 6f]. An increase from 10.0 to 30.7% for accurate calculations of vibrational frequencies.
in the first derivative deviations results when this
term is neglected in the functional form. The bond/
BOND ANHAKMONlCllY
torsion coupling interaction is used to model the
torsion angle dependence of forces on the bond The importance of anharmonicity in the bond
length, as discussed above [see eq. (16)]. The im- energy function is indicated by the results in the
portance of this coupling term for calculations of second row of Table IV. These results were ob-
bond length forces is consistent with the large in- tained by omitting the cubic and quartic bond en-
creases (0.014 and 0.019 A) in bond lengths that ergy terms from eq. (2) and calculating the energy
are observed in the eclipsed ethane and neopen- and derivatives, along with the rms deviations
tane structures. Table IV shows that the anhar- from the ab initio results, with the remaining QMFF
monic bond, bond/angle, and angle/angle/torsion energy function terms and parameters in Table 11.
terms are the three other most important terms for The results in Table I show that elimination of bond
reproducing intramolecular forces in alkanes. anharmonicity costs 0.4 kcalimol in the fit of ab
When these terms are neglected, the first derivative initio energies, that is, the rms energy deviation
deviations increase from 10.0 to 22.7, 21.4, and increased from 0.96 to 1.36 kcal/mol. Errors in the
16.8%, respectively. This indicates that these three first derivatives more than double (from 10.0 to
terms are also important for calculations of struc- 22.7%) when bond anharmonicity is neglected,
tural changes in strained alkanes. while errors in the second derivatives increase from

180 VOL. 15, NO. 2

 DERIVATION OF CLASS II FORCE FIELDS

4.2 to 18.2%. This fourfold increase in the second tocol described here can be used to derive a force
derivative error dramatically demonstrates the im- field for any functional group for which energies
portance of bond anharmonicity in the ab initio po- and first and second derivatives can be obtained
tential energy surface. from quantum mechanics. The result will be an
analytic representation that, if used in a simulation,
will give essentially comparable results to that
Conclusions which would have been obtained had one used the
quantum energy surface itself to within the range
A general method for determining force fields of deviations reported here.
for potential energy functions from quantum me-
chanical energy surfaces has been presented. The
functional form of the resulting force field, a Class Acknowledgment
I1 force field, was given and the method was used
to derive an alkane quantum mechanical force field, The authors gratefully acknowledge the support
that is, a force field that fits the quantum mechan- of the Consortium for Research and Development
ical surface. All force constants and reference val- of Potential Energy Functions. Former and present
ues were determined. The force field was tested in members of the Consortium include Abbott Lab-
terms of the fit to the ab initio energy surface by oratories; Battelle Molecular Science Research Cen-
comparing energies and first and second deriva- ter; Bayer AG; Bristol-Myers Squibb Pharmaceut-
tives. The results showed that the quality of the fit ical Group; Convex Computer Corp.; Cray
is generally good. The quantum mechanical energy Research, Inc.; Deutsches Kunstoff Institute; Dow
surface was found to contain significant anhar- Chemical Co.; DuPont Merck Pharmaceutical Co.;
monicity and intramolecular coupling interactions. E. I. DuPont de Nemours & Co.; Eastman Kodak
This was demonstrated by comparing the Class I1 Co.; Esteve SA, Laboratorios de Dr; ETA Systems;
force field derived here including analytic repre- Farmitalia Carlo Erba; Fujitsu Systems Business of
sentations of these effects with a diagonal quadratic America; Gessellschaft fur Biotechnologische For-
Class I force field, that is, the functional form used schung mbH; Glaxo Group Research Ltd.; Hitachi,
commonly in potential functions such as AMBER Ltd.; Hybritech, Inc.; IBM Almaden Research Cen-
and CHARMM for the simulation of organic and ter; ICI Americas, Inc.; Immunex Research & De-
biomolecular systems. (It was pointed out that velopment Corp.; Institut Francais du Petrole; Jans-
CVFF, which contains some coupling terms, is sen Research Foundation; Ludwig Institute;
somewhere between a Class I and Class I1 force Mitsubishi Petrochemical; Marion Merrell Dow Re-
field.) search Institute; Merck, Sharp & Dohme Research;
Finally, the importance of the various energy Monsanto Co.; National Cancer Institute; North
terms was determined quantitatively by assessing Carolina Supercomputing Center; Protein Engi-
the effect on the fit to the energy, first, and second neering Research Institute; Rhbne-Poulenc; Rohm
derivatives from omitting the individual terms. In & Haas Co.; Sandia National Laboratories; Sandoz,
a subsequent article, the quantum force field de- Ltd.; Searle Research & Development; Shell Re-
rived here will be scaled to a fit of experimental search, Ltd.; SINTEF; SRI Corp.; Synergie Tech-
data and an experimental force field will be de- nologies, Ltd.; Takeda Chemical Industries, Ltd.;
rived. The quantum force field presents enough Universite De Lille I Cerim; Universite De Nancy
information on the shape and relative values of the I; U.S. Army, Ballistic Research Laboratory and Vi-
force constants such that only a limited number of rology Division; and Upjohn Laboratories.
experimental data need to be used and only a few
scaling factors and reference values need to be de-
termined, thus greatly reducing the amount of in- References
formation required from experimental data.
1. A. T. Brunger, J. Kuriyan, and M. Karplus, Scierzcc, 235,458
Finally, it was pointed out that this methodology (1987).
should provide a useful starting point for use in 2. G. M. Clore, M. Nilges, D. K. Sukumaran, A. T. Brunger,
simulations for functional groups where insuffi- M. Karplus, and A. M. Gronenborn, EMBO I . , 5, 2729
cient amounts of experimental data are available or (1986).
where, for pragmatic reasons, there is a need to 3. P. A. Bash, U. C. Singh, F. K. Brown, R. Langridge, and
proceed without fitting experimental data. The pro- P. A. Kollman, Science, 235, 574 (1987).

JOURNAL OF COMPUTATIONAL CHEMISTRY 181

MAPLE ET AL.

4. A. Warshel and F. Sussman, Proc. Natl. Acad. Sci. U S A , 83, and U. Dinur, In Computer Simulation of Bioriiolecular Sys-
3806 (1986). tems-Theoretical and Experimental Applications, van Gun-
5. F. Avbelj, J. Moult, D. H. Kitson, M. N. G. James, and steren and Weiner, Eds., ESCOM, Leiden, 1989, pp. 149-
A. T. Hagler, Biocheinistry, 29, 8658 (1990). 167.
6. M.-J. Hwang, T. P. Stockfisch, and A. T. Hagler, submitted. 28. U. Dinur and A. T. Hagler, In Revierus in Computational Chem-
7. L. S. Bartell, 1. Cheni. Phys., 32, 827 (1960); L. S. Bartell, istry, Vol. 2, K. 8. Lipkowitz and D. B. Boyd, Eds., VCH
Tetrahedron, 17, 177 (1962). Publishers, New York, 1991, chap. 4, pp. 99-164.
8. S. Fitzwater and L. S. Bartell, 1. A m . Chem. Soc., 98, 5107 29. A. T. Hagler, In The Peptides, Vol. 7, V. J. Hruby and
(1976). J. Meienhofer, Eds., Academic, New York, 1985, pp. 213-
299.
9. J. B. Hendrickson, 1. A m . Cheni. Soc., 83, 4537 (1961); J. B.
Hendrickson, 1. A m . Cheni. SOC.,86, 4854 (1964). 30. S. Dasgupta and W. A. Goddard, 1. Chem. Pkys., 90, 7207
10. T. Simanouti, I. Ckerii. Phys., 17, 734 (1949); T. Simanouti (1989).
and S. Mizushima, 1. Chern. Phys., 17, 1102 (1949). 31. W. J. Hehre, L. Radom, P. v. R. Schleyer, and J. A. Pople,
11. J. H. Schachtschneider and R. G. Snyder, Spectrochim. Acta, Ab Initio Molecular Orbital Theory, John Wiley & Sons, New
19, 117 (1963); R. G. Snyder and J. H. Schachtschneider, York, 1986.
Spectrocllirn. Acta, 21, 169 (1965). 32. C. E. Blom and C. Altona, Mol. Phys., 31, 1377 (1976).
12. S. Lifson and A. Warshel, 1. Chew. Phys., 49, 5116 (1968). 33. M. J. Frisch, M. Head-Gordon, H. 8. Schlegel, K. Raghav-
13. A. Warshel and S. Lifson, 1. Chem. Phys., 53, 582 (1970). achari, J. S. Binkley, C. Gonzalez, D. J. Defrees, D. J. Fox,
14. A. T. Hagler, P. S. Stern, S. Lifson, and S. Ariel, 1. A m . R. A. Whiteside, R. Seeger, C. F. Melius, J. Baker, L. R.
Kahn, J. J. P. Stewart, E. M. Fluder, S. Topiol, and J. A.
Cllcvi. Soc., 101, 813 (1979).
Pople, Gaussian 88, Gaussian, Inc., Pittsburgh, PA, 1988.
15. S. Lifson and P. S. Stern, 1. Cherii. Phys., 77, 4542 (1982).
34. E. B. Wilson, Jr., J. C. Decius, and P. C. Cross, Molecirlar
16. R. H. Boyd, 1. Chern. Pliys., 49, 2574 (1968). Vibrations: The Theory of Infrared and Rariiari Vibratiorial Spec-
17. C. Shieh, D. McNally, and R. H. Boyd, Tetrahedron, 25,3653 tra, Dover Publications, New York, 1955.
(1969).
35. J. P. Bowen and N. L. Allinger, In Rez1iet.c~in Coiiipirtatiorial
18. S. Chang, D. McNally, S. Shary-Tehrany, S. M. J. Hickey, Cheniistry, Vol. 2, K. B. Lipkowitz and D. B. Boyd, Eds.,
and R. H. Boyd, 1. A ~ I IChevi.
. Soc., 92, 3109 (1970). VCH Publishers, New York, 1991, chap. 3, pp. 81-97.
19. N. L. Allinger, M. A. Miller, F. A. Van Catledge, and J. A. 36. P. M. Morse, Phys. Rev., 34, 57 (1929).
Hirsch, 1. A m . Clleirr. Soc., 89, 4345 (1967); N. L. Allinger,
J. A. Hirsch, M. A. Miller, I. J. Tyminski, and F. A. Van 37. U. Burkert and N. L. Allinger, Molecidnr Meclinnics, Amer-
Catledge, 1. A m . Cheni. Soc., 90, 1199 (1968). ican Chemical Society, Washington, DC, 1982.
20. N. L. Allinger, 1. Atti. Cheiii. Soc., 99, 8127 (1977). 38. M. Waldman and A. T. Hagler, 1. C o r y . C h i . , 14, 1077
(1993).
21. N . L. Allinger, Y. H. Yuh, and J.-H. Lii, 1. A m . Chein. Soc.,
111, 8551 (1989). 39. T. Oie, G. M. Maggiora, R. E. Christoffersen, and D. J.
22. J.-H. Lii and N. L. Allinger, 1. A m . Chenl. Soc., 111, 8566 Duchamp, Int. 1. Quantuiii Chem. Qiiantiirfi B i d . S y i i p . , 8, 1
(1989). (1981).
23. J.-H. Lii and N. L. Allinger, 1. A m . Chern. Soc., 111, 8576 40. T. A. Halgren, 1. A m . Chefif.Soc., 112, 4710 (1990).
(1989). 41. A. Warshel and S. Lifson, Chem. Phys. Lett., 4, 255 (1969).
24. S. J . Weiner, P. A. Kollman, D. T. Nguyen, and D. A. Case, 42. L. Norskov-Lauritsen and N. L. Allinger, 1. Cotrip. Cherii.,
1. Corrip. Cherri., 7, 230 (1986). 5, 326 (1984).
25. L. Nilsson and M. Karplus, 1. Coaip. Cheni., 7, 591 (1986). 43. R. Fletcher, Practical Methods of Oytirnization, Vol. 1, Wiley,
26. J . R. Maple, U. Dinur, and A. T. Hagler, Proc. Natl. Acad. New York, 1980.
Sci. U S A , 69, 5350 (1988). 44. N. Draper and H. Smith, Applied Regression Aiinlysis, 2nd
27. A. T. Hagler, J. R. Maple, T. S. Thacher, G. B. Fitzgerald, ed., Wiley, New York, 1981.

182 VOL. 15, NO. 2