Consistent Force Field (CFF)
* Menu:
* Usage:: How to use CFF with CHARMM standalone
* Status:: Current status of CFF implementation in CHARMM
* Theory:: Basis for, parameterization and performance of CFF
* Funcform:: Functional form of the CFF energy expression
* Refs:: References to papers describing CFF
In order to use CFF in CHARMM, the user has to issue the following
commands:
1. use cff
2. read cff parameter file
3. (a) read rtf name <CFF-capable rtf file>, or
(b) read psf name <file_name>
4. read sequence ! if input is via the rtf route (step 3 (a))
5. generate
6. read coord, or ic build ! if input is via the read rtf/sequence route.
When using CFF95 or later Step 3a requires a CFF-capable rtf file. This means
a file in which BOND records have been replaced by analogous DOUBLE records for
cases in which the chemical structure has a double bond. Note that CFF-capable
rtf files are *back compatible*. That is, such rtf files can equally well be
used for calculations that utilize the CHARMM force field. Thus, it is *not*
necessary to maintain two versions of the rtf files.
NOTE: 1. no binary parameter files are supported for CFF.
2. CFF is an all hydrogen force field -- i.e., extended atoms
are not supported
Examples of CFF usage in CHARMM are given in the ccfftest directory.
Status of CFF implementation into CHARMM (May 2001)
=============================================================
This implementation of CFF in CHARMM is principally due to Rick
Lapp (MSI) and William Young (MSI).
Features currently supported in CHARMM/CFF
(1) energy and first derivatives
(2) minimization
(3) dynamics
(4) most ATOM based cutoff options
Major features NOT currently implemented in CHARMM/CFF:
(1) bonds between primary atoms and image atoms.
(2) Cutoff options currently not supported are group-based cutoffs,
distance shifting and force-based switching.
(3) Fast multipoles.
Other known limitations:
(1) correlation analysis tools have not been implemented for CFF specific
energy terms -- e.g. it is not possible to calculate the correlation
function for an out-of-plane bending angle, etc ...
(2) only all-atom models (no extended atoms)
There are probably other problems/limitations/bugs. Your comments about
limitations of the current CFF implementation in CHARMM (and bugs) will be
very valuable.
Please direct comments to:
William Young, MSI
e-mail: wyoung@msi.com
phone: (619)799-5348
KNOWN BUGS:
The aim of the CFF development is a force field that is:
* broad, covering a relatively large number of differing
functional groups,
* accurate, achieved via accurate reproduction of the
quantum mechanical energy surfaces,
* consistent between differing phases and molecular environments,
* applicable to a wide range of molecular properties,
* consistent between differing types of molecules,
such as interaction of protein active sites with ligands, or
assemblies of proteins with nucleic acids or with solvent.
Quantum mechanical forcefields
The intramolecular parameters constituting the current generation of
forcefields are based on the energies and energy derivatives computed by
ab initio quantum mechanical procedures for a series of model compounds.
CFF uses quantum computations in the Hartree-Fock approximation with the
6-31G* basis set to expand the wavefunctions [1][2]. The quantum
mechanical energies and the energy first derivatives (gradients) and
second derivatives (Hessians) were computed for the equilibrium molecular
structures, at conformational energy barriers, and for a set of distorted
structures. The distorted molecular structures were generated by randomly
deforming all the internal coordinates, as well as by systematically
rotating about individual bonds. These quantum observables were fit to
the energy expression to obtain the Class II parameters [3][4]. Many of
the atomic partial charges were also determined quantum mechanically.
The intermolecular parameters of the forcefield may also, in principle,
be computed quantum mechanically [5]. The remaining CFF forcefield
intermolecular or nonbond parameters were computed by fitting to
experimental crystal lattice constants and sublimation energies of
crystals [6][7][8].
Internal energy terms
The energy of the molecule or assembly is expressed in terms of internal
coordinates such as bond lengths, bond angles, and dihedral angles. For
Class II forcefields this set of descriptors is greatly expanded by
including cross terms, that is, the interactions between bond
lengths and angles, between pairs of angles, etc. CFF contains, in all,
twelve types of energy terms: bond stretching, valence angle bending,
valence dihedral angles, out-of-plane deformation, and eight cross terms.
The cross terms extend the accuracy and range of application of the
forcefield by including the effect of neighboring atomic positions on
each of the bond lengths, valence angles, and dihedral angles.
Energy functional forms
The energy expression may be decomposed into diagonal terms that depend
on a single molecular internal coordinate such as a bond length, coupling
terms between internal coordinates, and nonbond internuclear distances.
This energy is fit to the quantum mechanical energy.
1. Bond stretching.
Ebond = K2 * (b - b0)^2 + K3 * (b - b0)^3 + K4 * (b - b0)^4 (1)
where K2, K3 and K4 are the quadratic, cubic and quartic forcefield
parameters or force constants, b is the bond length, and b0 is the
reference value of the bond length.
2. Angle bending.
Eangle = K2 * Delta^2 + K3 * Delta^3 + K4 * Delta^4 (2)
where Delta = Theta - Theta0 is the difference between the actual and
reference bond angles.
3. Out-of-plane bending.
Eoop = K * (Chi - Chi0)^2 (3)
where chi is an out-of-plane coordinate as defined by Wilson et al.[9]
4. Torsion energy, in order to reflect differing hybridizations about
the bonded atoms, must contain one-, two-, and threefold periodic terms:
Etorsion = SUM(n=1,3) { V(n) * [ 1 - cos(n*Phi - Phi0(n)) ] } (4)
where phi is a dihedral angle.
5. Stretch-Stretch interaction between two bonds in a valence angle.
Ebond-bond = K(b,b') * (b - b0) * (b' - b0') (5)
6. Stretch-Bend interaction between an angle and its bonds.
Ebond-angle = K * (b - b0) * (Theta - Theta0) (6)
7. Bend-Bend-Twist interaction between a dihedral angle and its
two valence angles.
Eangle-angle-torsion = K * (Theta - Theta0) * (Theta' - Theta0') *
(Phi - Phi1(0)) (7)
8. Stretch-Twist interaction between a dihedral angle and its end bonds.
Eend_bond-torsion = (b - b0) * SUM { V(n) * cos[n*phi] } (8)
9. Stretch-Twist interaction between a dihedral angle and its middle bond.
Emiddle_bond-torsion = (b - b0) *
{ F(1) * cos(phi) + F(2) * cos(2 * phi) + F(3) * cos(3 * phi) } (9)
10. Bend-Twist interaction between a dihedral angle and its valence angles.
Eangle-torsion = (Theta - Theta0) *
{ F(1) * cos(phi) + F(2) * cos(2 * phi) + F(3) * cos(3 * phi) } (10)
11. Bend-Bend interaction between two valence angles with a common
vertex atom.
Eangle-angle = K * (Theta - Theta0) * (Theta' - Theta0') (11)
12. Stretch-Stretch interaction between the two end bonds in a dihedral
angle.
Ebond-bond_1_3 = K(b,b') * (b - b0) * (b' - b0') (12)
Finally, the nonbond energy between atoms in different molecules or
between atoms separated by three or more bonded atoms is given by the
sum of the Coulombic electrostatic interaction and a van der Waals
energy of the 9-6 form:
13. Coulombic electrostatic interaction.
Ecoul = 332.0716*qi*qj/(D*Rij) (13)
where qi and qj are the atomic partial charges on atoms i and j,
Rij is the distance between them and D is the dielectric constant.
14. Van der Waals interaction.
Evdw = eps(ij) [2*r*(ij)/r(ij)**9 - 3*r*(ij)/r(ij)**6] (14)
where r*(ij) = [(r(i)**6 + r(j)**6))/2]**(1/6) (15)
eps(ij) = 2 sqrt(eps(i) * eps(j)) *
r(i)^3 * r(j)^3/[r(i)^6 + r(j)^6] (16)
where eps(ij) and r*(ij) are the negative of the minimum van der Waals
energy and that distance between atoms i and j where the minimum occurs,
respectively. Eps(ij) and r*ij are computed from the individual atomic
parameters eps(i), eps(j), r*i, and r*j by the Waldman-Hagler combination
rules [10].
The Hartree-Fock method, and to a lesser extent other quantum mechanical
methods, results in systematic deviations from experiment. For example,
bond lengths tend to be too short and bond-stretching vibrational
frequencies too high [11]. However, by comparison with experimental
gas-phase molecular structures and vibrational frequencies, these
deviations may be compensated for. In general, the energy expression
may be scaled using five constant factors, one for each of the classes
of energy terms: bonds, angles, torsion angles, out-of-planes and all
coupling terms [12]. The scaled energy is then:
Ediagonal = Sb * SUM{Ebond} + Stheta * SUM{Eangle} +
Sphi * SUM{Etorsion} + Schi * SUM{Eoop} (17)
Ecross = Sc * SUM{eight cross terms} (18)
The reference values b0 and q0 are also adjusted to fit experimental data.
All these values may differ among different types of bonds, bond angles,
and torsion angles. For the special case of hydrocarbons, the corrections
are especially well determined by gas-phase measurements. For hydrocarbons,
the best values of the scale factors are:
Sb(C-C) 0.88
Sb(C-H) 0.83
Stheta (all angles) 0.81
Sphi (all torsions) 0.84
Schi (all out-of-planes) 1.00
Sc (all cross terms) 0.87
The reference bond lengths for hydrocarbons were also adjusted.
Although the use of the quantum calculation greatly amplifies the
available data so that only a few such corrections are necessary for
the complete Class II forcefield, for the majority of functional groups
(molecular types) no accurate gas-phase data are available. However,
the Sb, Stheta, Sphi, and Sc constants are transferrable among different
types of bonds, bond angles, and torsion angles. Therefore, the same
scale factors are used in Eq. 17 and Eq. 18 in the final empirically
scaled forcefield. In general, the reference values b0 and theta0 are
determined from high-level quantum mechanical calculations on the model
compounds.
Validation of the CFF forcefield
Table 1 shows the accuracy of the CFF forcefield for several common
classes of molecules, compared with experimental gas-phase results.
Table 1. Summary of rms deviations between experimental and
CFF-calculated structural parameters, vibrational frequencies,
and energy differences.
bond valence torsion freq. energy
length angle angle diff.
(Ang) (deg) (deg) (cm-1) (kcal mol-1)
hydrocarbons 0.02 0.9 1.2 40 0.93
alcohols 0.02 1.7 1.7 37 0.71
aldehydes & ketones 0.01 1.1 2.3 32 0.62
amines 0.00 0.9 -- 18 0.62
carboxylic acids 0.02 1.6 1.0 34 0.78
esters 0.02 1.7 0.5 42 1.88
ethers 0.01 0.9 1.1 41 0.40
heterocycles 0.01 1.0 0.0 35 ---
sulfides 0.01 1.4 2.5 45 ---
disulfides 0.01 0.9 2.0 43 ---
thiols 0.01 1.6 1.0 -- 0.21
average 0.01 1.2 1.3 37 0.77
The frequencies are harmonic vibrational frequencies and the energy
differences include conformational energy differences and energy
barriers to internal rotation between stable conformers.
References [1] Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 28, 213-222 (1973). [2] Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. J. Chem. Phys. 77, 3654-3665 (1982). [3] Dinur, U.; Hagler, A. T. In Reviews in Computational Chemistry, Vol. 2, K. B. Lipkowitz; D. B. Boyd, Eds., VCH Publishers: New York, 99-164 (1991). [4] Maple, J. R.; Hwang, M.-J.; Stockfisch, T. P.; Dinur, U.; Waldman, M.; Ewig, C. S.; Hagler, A. T. J. Comp. Chem. 15, 162-182 (1994). [5] Dinur, U.; Hagler, A. T. J. Amer. Chem. Soc. 111, 5149-5151 (1989). [6] Hagler, A. T.; Huler, E.; Lifson, S. J. Amer. Chem. Soc. 96, 5319-5327 (1974). [7] Hagler, A. T.; Lifson, S.; Dauber, P. J. Amer. Chem. Soc. 101, 5122-5130 (1979a). [8] Hagler, A. T.; Dauber, P.; Lifson, S. J. Amer. Chem. Soc. 101, 5131-5141 (1979b). [9] Wilson, E. B., Jr; Decius, J. C.; Cross, P. C., Molecular Vibrations; Dover: New York, 1955, Chapter 4. [10] Waldman, M.; Hagler, A. T. J. Comp. Chem. 14, 1077-1084 (1993). [11] Michalska, D.; Schaad, L. J.; Carsky, P.; Hess, Jr., B. A.; Ewig, C. S. J. Comp. Chem., 9, 495 (1988). [12] Hwang, M.-J.; Stockfisch, T. P.; Hagler, A. T. J. Amer. Chem. Soc. 116, 2515-2525 (1994).
CHARMM Documentation / Rick_Venable@nih.gov