Department of Chemistry and Chemical Biology

Atom types.

The first step in parameter development is to make a two-dimensional sketch of the fragment for which parameters are needed, and then to assign atom types to the atoms. The comments in the first section of the parm94.dat file describe the hybridization and other attributes of the atom types for the 1994 force field; [1] for the older force field, one may consult the Weiner et al. force field papers. [2] [3] This approach may be augmented by looking at the atom types in the existing residues in the files all_*94.in. For example, in pyridine the nitrogen would be assigned the same type (NC) as N1 and N3 in adenine. Note that an atom type is intrinsic to an array of distance, angle, and dihedral parameters involving the types of the neighboring atoms, as well as having its own van der Waals (VDW) parameters, and, of course, atomic mass (charge is not fixed per atom type). Therefore, if a new atom type is required, the first step is to attempt to reason by analogy and clone as many of the pre-existing parameters as possible to account for the environment of the new atom. Here it is instructive to consider the variability of the existing parameters, which tend to be duplicated over various combinations of atoms. This step may also be required if old atom types are used in a new topological relation.

For example, consider the oxygen of a sulfoxide or a sulfone:

                R                            R
                |                            |
               :S=O                        O=S=O
                |                            |
                R'                           R'

We would expect the van der Waals parameters of this oxygen to be similar to those of a carbonyl oxygen of the force field (type `O'):

                            R
                             \
                              C=O
                             /
                            R'

or to those of carboxyl or phosphate oxygens (type `O2' in both examples)

                                             R
                                             |
                    O2                       OS
                  /                          |
             R - C                      O2 - P - O2
                  \                          |
                    O2                       OS
                                             |
                                             R'

because VDW radii are dominated by the number of electrons in an atom and are not very sensitive to chemical environment. In fact, the VDW parameters for types `O' and `O2' are identical. Can one of these types be used for the oxygen in sulfoxide/sulfone? The environment must now be considered. If no suitable analogy can be found, a new atom type must be created and a complete set of bond, angle, and dihedral parameters for its neighbors added to the force field. In this case, both sulfoxide and sulfone oxygens are substituents of tetrahedral sulfurs rather than trigonal planar carbons, so the phosphate oxygen case makes an appealing analogy. We then check whether any existing bond, angle or dihedral parameters involve a S=`O2' bond, and if they do, we check that those parameters are appropriate for this case of an S=O bond. But there are no such parameters - it is therefore reasonable to extend the use of type `O2' for this case. (Had there been inappropriate S=`O2' parameters, we might want to either make a new sulfur type, or make a new oxygen type starting with the VDW parameters of `O' and `O2'.) We continue discussion of bond and angle parameters for these example fragments below.

van der Waals parameters.

What if no atom type lends itself to adaptation? When creating a new type, the first thing one must consider is VDW parameters. As illustrated above, for organic compounds these parameters may be straightforward to find by analogy based on element and bond order alone. Monoatomic ions, however, do not present such analogies in the AMBER force field and are discussed in more detail as an example.

The shape of the VDW potential for a given atom type is specified in terms of the distance between two atoms of the same type at the minimum energy point. Half the interatomic distance at that point is treated as the basic radius, or R*, parameter for that type. The form for the radial potential for two atoms is the sum of the R* values of their types. The potential well depth (`e') of the minimum energy point between two atoms of the same type is combined with the potential of another atom type by taking the root of the product. (Other parametric forms can be used which tend to have different type-type `combining rules'.)

The simplest approach to deriving VDW parameters is to match a relevant experimental determination of the size of the atom in question. One source of such measurements is diffraction data. The sum of metal and oxygen Pauling radii [4] tends to be 3% smaller than indicated by water-ion neutron and X-ray diffraction data for Li+ and Na+ ions, 2% smaller than for K+, and 1% smaller than for Rb+ and Cs+. [5] Another source of ion `size' information is crystallographic studies of ion complexes. [6] Since van der Waals parameters consist of two terms, the parameters that e.g. yield a given first peak of the radial distribution of the distance between two types of atoms are not unique. Another variety of experimental data that can contribute to parameterization is the free energy of solvation in water or another relevant solvent. However, it is still not clear whether the combination of experimental size and solvation free energy is sufficient to determine unique R* and `e' parameters for an atom in relation to an existing type. A further complication arises because an atom type may come into contact with more than one other type, and nothing in principle guarantees that VDW parameters for a group of types can be fitted to yield uniformly correct pairwise potentials. Therefore it is important to choose parameters consistent with the most significant atom types that the new type will come in contact with. To a first approximation, atom types that tend to be oppositely charged, if present, are of most interest. In a particularly important case for ions, the TIP [7] water models (as well as some other waters) have a spherical van der Waals potential centered on the water oxygen (type `OW'), which is somewhat inflated to enclose the hydrogen atoms in the molecule. Thus a cation that has been parameterized to give a correct ion-`OW' radial distance distribution function in such a water model will be ``smaller'' and come in closer contact with neighboring atoms if it is bound in a molecule consisting of AMBER atoms.

Moreover, remembering that different pairwise combining rules are in use in the modeling community, parameters from one convention must be adapted to yield the same results for a given pair of types in another scheme. Thus it was necessary to adapt the monovalent cation parameters of Aqvist [8] (found in parm94.dat and parm91.dat) for AMBER so that the ion-water combined potential gave the same optimal distance as with the combining rules used by Aqvist. Matching the ion size in the environment seems to be sufficient in the case with small monoatomic monocations; the default has traditionally [9] been to use a somewhat arbitrary well depth (epsilon) of 0.1 kcal/mol, characteristic of a rather nonpolarizable atom, and fit an R* parameter (see the Ross and Hardin reference). For the multivalent ions, different further approaches may be considered to capture the quasi-bonding electron mobility, including the use of explicit bonds or hydrogen bonding terms (see the Vedani and Huhta reference).

We have also discussed the derivation of van der Waals parameters for hydrogen in different bonding environments. [10] Using ab-initio calculations to study the interaction between water and various hydrogens, we found that a reduction in R* was required for hydrogens attached to carbons with adjacent electronegative atoms. This trend is nicely paralleled in the progessively smaller R* for hydrogens attached to carbon (HC), nitrogen (H), and water oxygen (HW).

Bond and angle parameters.

Having chosen or created one or more atom types and sets of van der Waals parameters, the bond, angle and dihedral parameters must be created. Equilibrium bond lengths and angles may be obtained from tabulations of experimental data in the literature. [11] [12] Initial bond and angle force constants may be chosen based upon analogy to similar parameters in the force field or using the method of Hopfinger and Pearlstein. [13] See [14] for an example of extending the Weiner et al. force field to guanosine triphosphate and analogs. The AMBER-1994 force field contains a limited number of unique bond and angle force constants and therefore selection by analogy is a feasible starting point. Returning to our sulfoxide/sulfone example, we find that the only existing bonds involving O2 are:

             Kbond    Rbond
      C -O2  656.0    1.250       JCC,7,(1986),230; GLU,ASP
      O2-P   525.0    1.480       JCC,7,(1986),230; NA PHOSPHATES

and it would be reasonable to use the O2-P force constant with a bond distance from the literature. Similarly, it would be reasonable to use the same angle bending parameters as for phosphates:

        Ktheta(O2-P-O2) = Ktheta(O-S-O)
        Ktheta(O2-P-OS) = Ktheta(R-S=O)
        Ktheta(OS-P-OS) = Ktheta(R-S-R)

Unless only crude parameters are desired, one should check the force constants by means of normal mode calculations if spectroscopic measurements are available for comparison; such calculations on N-methylacetamide are described in the Weiner et al. JACS 1984 paper. The bond and angle parameters are the primary determinants of the high and middle frequency vibrational modes of a molecule. For applications which require the highly accurate reproduction of vibrational frequencies, it is necessary to use a force field which includes higher order terms (anharmonicity) and cross-terms. For modeling the structures and interactions of molecules which are not highly strained, however, the simple harmonic approximation used in AMBER appears to be quite adequate.

Dihedral parameters.

The dihedral parameters, in conjunction with the atomic charges and van der Waals parameters, are the primary determinants of the relative conformational energies of a molecule. The AMBER parameters IDIVF, PK, PN, and PHASE are used to define the torsional potential energy function. Each bonded series of atoms I-J-K-L must have at least one set of these dihedral parameters in the force field (just as every bonded pair I-J or triplet I-J-K must have bond or angle parameters, except that for dihedrals multiple terms may be used). The torsional energy function formula is:

      Etors = ( PK / IDIVF ) * ( 1 + cos( PN * phi - PHASE) )

Let us look at a few examples in order to illustrate the nature of the dihedral parameters.

Figure 1

For our first example (Figure 1), if atoms J and K are sp3 carbons (type CT) as in the molecule ethane (H3C-CH3), then the intrinsic barrier to rotation about the J-K bond is on the order of 3 kcal/mol. PK is equal to one-half of the barrier magnitude [15] and would therefore be equal to 1.5 kcal/mol. The topology about the dihedral of interest has a three-fold periodicity (PN), that is, there are three potential barriers as the C-C bond is rotated -180 to 180 degrees. These barriers occur when the methyl hydrogens eclipse each other: at 0, -120, and 120 degrees. Since the dihedral formula is a Fourier series truncated to a single cosine term, no phase shift would be needed to reproduce the potential energy barriers and PHASE = 0 degrees. (PHASE = 0 degrees if an energy maximum is at 0 degrees; PHASE = 180 degrees if an energy minimum is at 0 degrees.) The final dihedral parameter that must be specified is the number of torsions associated with the central bond (IDIVF), which is the product of the number of substituents on the two central atoms. For ethane, this is 3x3 = 9. So we have:

             PK    = 3.0 kcal/mol / 2.0 = 1.5 kcal/mol
             PN    = 3
             PHASE = 0.0 degrees
             IDIVF = 9

These same torsional parameters can be used for n-butane, and the results are in good agreement with experiment and higher-level calculations for the relative energy of trans and gauche minima and cis and skew energy barriers.

Consider now the molecule ethylene, H2C=CH2, whose dihedral potential energy is shown in Figure 2.

Figure 2

The lowest-energy conformation of this molecule is planar with a two-fold (PN = 2), 60 kcal/mol (PK = 30.0 kcal/mol) barrier to rotation about the C=C bond. The barriers are found at dihedral angles of -90 and 90 degrees (energy minimum at 0 degrees), and can be reproduced by the truncated Fourier series only if a phase shift of 180 degrees (PHASE = 180.0 degrees) is used.

             PK    = 60.0 kcal/mol / 2.0 = 30.0 kcal/mol
             IDIVF = 4
             PN    = 2
             PHASE = 180.0

Finally, we examine a hypothetical molecule ZH2C-CH2Z, where Z represents an electronegative functional group. Let us imagine that we either have experimental data on the relative conformational energies or we have simulated the rotational potential of this molecule with a series of quantum mechanical calculations. In practice, this is only done for minimum and maximum energy conformations - trans, gauche+, gauche-, eclipsed, skew, etc. In our example, the energy profile shows that the trans conformation (Z-C-C-Z = 180 degrees) is about 0.5 kcal/mol less stable than the gauche.

Before fitting the torsional parameters, we must generate the energy profile for the molecular mechanical nonbonded potential as was done for the quantum potential, subtract this curve from the quantum curve, and fit the torsional potential to the difference potential.

Before these calculations can be done, atomic charges need to be calculated, also by fitting to quantum mechanical results. The difference potential is then deconvoluted into Fourier series terms (Figure 3)

Figure 3

These Fourier series terms give the force field parameters:

                        IDIVF     PK       PHASE    PN
            Z-CT-CT-Z     1      0.260       0      -3
            Z-CT-CT-Z     1      0.384       0      -2
            Z-CT-CT-Z     1      0.241       0       1

which result in the total torsional potential shown in Figure 4. (In AMBER, PN is set to less than zero when additional terms remain to be read.)

Figure 4

Care must be taken when deconvoluting the torsional potential not to introduce spurious minima or maxima into the rotational energy profile. The combined potential of the deconvoluted parameters can be plotted directly by a graphing program, or the torsional energy profile can be `empirically' generated at 20-30 degree intervals in AMBER.

In practice, such an elaborate Fourier series treatment may not be appropriate because (a) the quantum mechanical treatment may not be accurate enough to warrant it, (b) one would rather have a simpler torsional potential that is more consistent with the existing force field and (c) the electrostatic potential fitting procedure may capture the torsional energy profile well enough so that many terms are not needed. For example, in the case of 1,2-difluoroethane, the known gauche tendency of the fluorines can be simulated by adding a twofold torsion

                           IDIVF     PK       PHASE    PN
              F-CT-CT-F      1       X        0        2

with `X' adjusted to make the total molecular mechanical energy of the gauche conformation 1 kcal/mol lower than the trans conformation.

Improper torsions.

Improper torsions are so named because the atoms involved are not serially bonded; rather they are branched:

                             J
                             |
                             K
                            / \
                           I   L

                      Improper I-J-K-L

The convention is that the central atom is listed in the third position of the dihedral (`K' in the figure), and the order of the other three is determined alphabetically by atom type - and when types are the same, by atom number (order in the molecule). This convention guarantees only that programs will give the same results when given the same residue definitions; i.e. there is no physical meaning involved.

Improper dihedral potentials are sometimes necessary to reproduce out-of-plane bending frequencies, i.e. they keep four atoms properly trigonal planar for a two-fold torsional potential (PN=2). They are additionally used in the united-atom force field model when a carbon with an implicit hydrogen is a chiral center; in effect they keep the position from inverting (PN=3).

The PHASE for improper torsions is always 180 degrees. Improper torsional parameters listed in the force field file can use ``wild-card'' specifications (`X') for the non-central atoms (note that wild-card impropers must follow the explicit ones in the parm.dat force field file). When using the Prep-Link-Edit-Parm (`PLEP') programs to create coordinate and parameter input files, the improper torsions must also be specified by atom name in the Prep residue input file. In LEaP, every atom with three substituents is matched against the impropers in the force field file, and all matches are applied (discarding any wild-card terms if an explicit match is found). Thus care must be taken that a new improper does not inadvertently match other cases. In both PLEP and LEaP, an improper with no wild cards causes all wild-card-containing impropers to be ignored. Except for not mixing wild-card with explicit cases, all improper terms that match a given central atom are applied. In LEaP, if no match is found, no improper term is applied (unlike bonds, angles and `proper' torsions, for which parameters must exist). In PLEP, a parameter must be in parm.dat if there is an improper in the Prep input.

Hydrogen bonding parameters.

Unlike the previous AMBER force fields, the 1994 force field does not include a 10-12 hydrogen bonding function. This function, however, is still supported by the software. When using the hydrogen bonding function, all relevant pairs of atom types need to have parameters. Note that if a pair of atom types has H-bond (10-12) parameters, these will override any van der Waals (6-12) parameters for that pair.

Atomic charges.

To generate atomic (`point' or `partial') charges for a new fragment, the electrostatic potential around the fragment is first calculated at a grid of points using quantum or semiempirical methods, then charges at the atom centers are fit to reproduce the potential using the RESP (Restrained ElectroStatic Potential) program. The amount of effort applied can vary considerably; multiple conformations and even families of analogs can be considered, as has been done with the proteins and nucleic acids in the AMBER-1994 [16] force field. The charges for the AMBER-1994 force field have been calculated at the 6-31G* basis set level of ab initio theory and to maintain consistency, this basis set should be used when calculating charges for a new molecule. (See Appendix D regarding the charge-fitting philosophy.)

It is important to energy minimize the fragment at a suitable level of theory before calculating the electrostatic potential around it. Quantum mechanical geometry optimizations were considered to be too expensive for the amino acid charges in the AMBER-1994 force field. Rather, the RESP charges were calculated after first minimizing the geometries using the Weiner et al. force field. MM2/MM3-minimized geometries would also be reasonable for examples when a quantum mechanical optimization is deemed too expensive. In general, the charges are likely to vary more as a function of conformation than degree of optimization. This is illustrated by the results obtained for the alanyl and glycyl dipeptides (Cornell et al. paper).

  [1]   Cornell, W.D., Cieplak, P., Bayly,  C.I.,  Gould,  I.R.,
      Merz,  Jr.  K.M.,  Ferguson,  D.M.,  Spellmeyer, D.C., Fox, T.,
      Caldwell, J.W., Kollman, P.A.  J.  Am.  Chem.  Soc.  (1995,  in
      press).
  [2]   Weiner, S.J., Kollman, P.A., Case,  D.A.,  Singh,  U.C.,
      Ghio,  C.,  Alagona,  G.,  Profeta, S., Jr., Weiner, P.  J. Am.
      Chem. Soc.  106, 765-784 (1984).
  [3]   Weiner,  S.J.,  Kollman, P.A., Nguyen, D.T., Case, D.A.
      J. Comput. Chem.  7, 230-252 (1986).
  [4]   Pauling,  L.  The  Nature  of the Chemical Bond and the
      Structure of Molecules and Crystals; Cornell University  Press:
      Ithaca, New York, 1960.
  [5]   Ross, W.S. and Hardin, C.C., J.  Am.  Chem.  Soc.   116,
      6070-6080  (1994).  (This discussion of VDW parameters is based
      on that work.)
  [6]   Vedani, A., Huhta,  D.  W.,  J.  Am.  Chem.  Soc.   112,
      4759-4767 (1990) and references therein.
  [7]   Jorgensen, W.L., Chandreskhar, J., Madura, J.D.,  Impey,
      R.W., and Klein, M.L.  J. Chem. Phys.  79, 926-935 (1982).
  [8]   Aqvist, J., J. Phys. Chem. 94, 8021-8024 (1990).
  [9]   Wipff, S. J., Weiner, P., Kollman, P. A.  J.  Am.  Chem.
      Soc.  104, 3249 (1982).
  [10]  Veenstra, D. L., Ferguson, D. M., and  Kollman,  P.  A.
      J. Comput. Chem.  13, 971-978 (1992).
  [11] Allen, F. H., Kennard, O., Watson, D. G., Brammer,  L.,
      Orpen,  A.  G., and Taylor, R.  J. Chem. Soc. Perkin Trans. II,
      S1-S19 (1987).
  [12]  Harmony, M. D.,  Laurie,  R.  W.,  Kuczkowski,  R.  L.,
      Schwendemann,  R. H., Ramsay, D. A., Lovas, F. J., Lafferty, W.
      J., and Maki, A. G.  J. Phys. Chem. Ref. Data, 8, 619 (1979).
  [13]  A.J. Hopfinger and R.A. Pearlstein, J. Comp. Chem., 5,
      486 (1985).
  [14]  J.F. Cannon, J. Comp. Chem., 14, 995-1005 (1993).
  [15]  Because the ( 1 + cos() ) term ranges from 0..2, it is
      scaled to 0..1 by including the divisor  of  2  in  the  energy
      term, PK.
  [16]  Cornell,  W.D., Cieplak, P., Bayly, C.I., Gould, I.R.,
      Merz, Jr. K.M., Ferguson,  D.M.,  Spellmeyer,  D.C.,  Fox,  T.,
      Caldwell,  J.W.,  Kollman,  P.A.   J.  Am. Chem. Soc. 1995,
      117:19, 5179-5197. 
      (parameter files)

Department of Chemistry and Chemical Biology, 102 Hurtig Hall, 360 Huntington Ave., Boston, MA 02115
Tel: 617-373-2822 Fax: 617-373-8795

Send mail to e.chudin@neu.edu with questions or comments about this web site.
Last modified: Tuesday, April 01, 2003