Molecular Mechanics Methods

DESCRIPTION

There are three Molecular Mechanics methods available in Gaussian. They were implemented for use in ONIOM calculations, but they are also available as independent methods. No basis set keyword should be specified with these keywords.

The following force fields are available:

Amber: The Amber force field as described in [Cornell95]. The actual parameters (parm96.dat) have been updated slightly since the publication of this paper. We use this current version from the Amber web site (ambermd.org).

Dreiding: The Dreiding force field as described in [Mayo90].

UFF: The UFF force field as described in [Rappe92].

MOLECULAR MECHANICS MOLECULE SPECIFICATIONS

Molecular Mechanics calculations use atom types to determine the functions and parameters which make up the force field. For a single element, such as carbon, there can be many different MM atom types. Which one to use depends on aspects such as the hybridization, chemical environment, and so on.

Gaussian will attempt to assign atom types automatically for UFF and Dreiding calculations. However, Amber calculations require that all atom types be explicitly specified within the molecule specification section, as in these examples:

C-CT      Specifies an SP3 aliphatic carbon atom, denoted by the Amber CT keyword.
C-CT-0.32 Specifies an SP3 aliphatic carbon atom with a partial charge of +0.32.
O-O--0.5  Specifies a carbonyl group oxygen atom (type O) with a partial charge of -0.5.

The atom type keyword is specified following the element symbol and a separating hyphen. Consult the Amber paper [Cornell95] for definitions of atom types and their associated keywords.

As the second and third lines in the example illustrate, partial charges may also be specified, as the third component of the atom description, preceded again by a hyphen separator. Partial atomic (point) charges are used to compute the electrostatic interactions. Gaussian can assign these charges automatically using the QEq formalism [Rappe91]. Request this by specifying the QEq option to the MM keyword: e.g., Dreiding=QEq. Note that these charges depend on the geometry and are computed at the beginning of the calculation; however, they are not updated during the course of an optimization or other process that changes the geometry.

Atom types and charges may also be provided for UFF and Dreiding calculations if desired. Be aware that the automatic assignment of UFF and Dreiding atom types is only reliable for well-defined systems, and it is often safer to specify them explicitly.

When using MM methods, and especially when you are modifying force field terms or defining new ones (as discussed below), it is often necessary to specify the connectivity between atoms explicitly, using Geom=Connectivity. Note that GaussView includes this input section automatically when it constructs input files.

CHARGE ASSIGNMENT-RELATED OPTIONS

Unless set in the molecule specification input, no charges are assigned to atoms by default when using any molecular mechanics force field. Options are available to estimate charges at the initial point using the QEq algorithm [Rappe91], under control of the following options, for any of the mechanics keywords:

QEq
Assign charges to all atoms using the QEq method [Rappe07]. OldQEq requests the older version which was the default in Revision C [Rappe91]. QEq=Uncharged sets the charges only for atoms which have MM charge zero, constraining the other atoms to retain their specified charges.

UnTyped
Assign QEq charges only to those atoms for which the user did not specify a particular type in the input.

UnCharged
Assign QEq charges for all atoms which have charge zero (i.e., all atoms which were untyped or which were given a type but not a charge in the input).

PARAMETER PRECEDENCE OPTIONS

Terminology: Gaussian contains built-in parameter sets for the built-in force fields listed above; these are referred to as hard-wired parameters. Soft parameters are ones specified by the user in the input stream for the current job (or a previous job when reading parameters from the checkpoint file). By default, when no relevant option is given, the hard-wired parameters are the only ones used. This topic is discussed in detail in “General MM Force Field Specifications” below.

HardFirst
Read additional parameters from the input stream, with hard-wired parameters having priority over the read-in, soft ones. Hence, read-in parameters are used only if there is no corresponding hard-wired value. Note that wildcard matches within the hard-wired parameter set take precedence over soft parameters, even when the latter contains an exact match for the same item. Use SoftFirst if you want to override hard-wired parameter matches.

SoftFirst
Read additional parameters from the input stream, with soft (read-in) parameters having priority over the hard-wired values.

SoftOnly
Read parameters from the input stream and use only them, ignoring hard-wired parameters.

HANDLING MULTIPLE PARAMETER SPECIFICATION MATCHES

Since parameters can be specified using wildcards (see below), it is possible for more than one parameter specification to match a given structure. The default is to abort if there are any ambiguities in the force field. The following options specify other ways of dealing with multiple matches.

FirstEquiv
If there are equivalent matches for a required parameter, use the first one found.

LastEquiv
If there are equivalent matches for a required parameter, use the last one found.

READING MM PARAMETERS FROM THE CHECKPOINT FILE

By default, when the geometry is read-in from the checkpoint file with Geom=Check or Geom=AllCheck, any non-standard (soft) parameters present in the checkpoint file are also retrieved. These read-in parameters are used with higher priority than corresponding hard-wired parameters, unless HardFirst is also specified. The following options can be used to modify this default behavior.

ChkParameters
Read only the MM parameters from the checkpoint file. Any non-standard (soft) parameters present in the checkpoint file are used with higher priority than corresponding hard-wired parameters, unless HardFirst is also specified. Not valid with Geom=Check or Geom=AllCheck.

NewParameters
Ignore any parameters in the checkpoint file when reading the geometry.

Modify
Read additional parameter specifications and modifications after combining any soft parameters from the checkpoint file and the hard-wired parameters (as controlled by the precedence options). This option works analogously to Geom=Modify.

PRINTING MM PARAMETERS

Print
By default, only the contributions to the energy are printed (i.e., from bond stretches, bends, electrostatics, etc.) when #P is specified. The Print option will also print the energy contributions and the force field parameters at the first geometry, and the energy contributions at later ones (since the parameters don’t change). It will print two sets of parameters for ONIOM(MO:MM) since different parameters are assigned for the model and real systems (i.e., since H may replace a C link atom).

Note that the parameters are printed for all atoms in the system, but not in a form suitable for input. The internally stored parameters are already provided in the form of input files in the main Gaussian directory (amber.prm, uff.prm, etc.), along with files for some force fields which are not stored internally (amber98.prm and oplsaa.prm). Printing for all atoms in the current job is still useful because complex rules are involved in deciding what parameters are actually applied to a given bond, dihedral, etc.

CONTROL OF SWITCHING FUNCTIONS

VRange=N
Set the van der Waals cutoff to N Angstroms. This is ROff if soft cutoffs are used, which is the default. N<0 for hard cutoffs.

CRange=N
Set the Coulomb cutoff to N Angstroms. Negative for hard cutoffs. The potential is (QiQj(1-r2/R2)2)/R.

CutRad=M
Use soft cutoffs with radius M (in Angstroms), with either van der Waals or Coulomb if VRange or CRange was set. Thus with VRange=N, CutRad=M the soft cutoffs have ROn=N-M, ROff=N.

Switch=value
Switch=YK means use the York-Karplus function, which is the default. Switch==CHARMM means use CHARMM-style switching for van der Waals. Switch=CHARMMSq means to use CHARMM-style expression but with R2 instead of R for van der Waals switching.Switch=N means use switching function number N.

AVAILABILITY

Analytic energies, gradients, and frequencies.

RELATED KEYWORDS

ONIOM, Geom=Connectivity

GENERAL MOLECULAR MECHANICS FORCE FIELD SPECIFICATIONS

Gaussian offers a tremendous amount of flexibility to users that want to modify or otherwise customize MM force fields. This section discusses the specifics of doing so, and it begins with some required background information necessary for understanding the function and parameter input.

A basic MM potential function;for MM is usually of the form:

Molecular Mechanics Potential Function
Molecular Mechanics Potential Function

As noted, the terms of the expansion are the stretch terms (bonds), bend terms (bond angles), torsional terms (dihedral angles), and non-bonded interactions. The particular forms of the individual functions in this equation are from the Amber force field, which uses simple harmonic functions for stretches and bends, a cosine function for torsions, and the standard functions for electrostatic and van der Waals interactions. Other force fields use different functions, and some include additional types of terms not included by Amber, such as dipole interactions or stretch-bend couplings.

In order to evaluate the MM potential function, Gaussian needs to know which structures—stretches, bends, and torsion—are present in the system, as well as the functions and parameters to be used for them. The structures can be identified from the molecular connectivity. By default, Gaussian determines which atoms are bonded and the corresponding bond orders (single, double, and so on). It does a good job when the input geometry is reasonable and the bond orders (single, double, etc.) are all well defined. However, if the calculation is started with a more approximate geometry, or if there are bonds whose orders are not easy to identify, it is safer to specify the connectivity list explicitly in the input stream, using Geom=Connectivity keyword (see the discussion of Geom for details).

The remainder of this section discusses MM functions and their parameters in general, and the following section documents all of the available functions.

Viewing the MM functions. A useful way to begin consideration of this topic is to examine the output from a Gaussian MM job. Here is the input file we used for an Amber calculation on methane:

#P Amber Geom=Connectivity IOp(4/33=3)
 
Methane
 
0 1
C-CT--0.4    -0.85  0.42  0.00
H-HC-0.1     -0.50 -0.57  0.00
H-HC-0.1     -0.50  0.93  0.87
H-HC-0.1     -0.50  0.93 -0.87
H-HC-0.1     -1.92  0.42  0.00

Note that the atoms types are CT for the carbon atom and HC for the hydrogen atoms. We have also assigned partial charges to each atom.

Here is the portion of the output produced by IOp(4/33=3) in the job’s route section:

 Atomic parameters:                         MM functions for single atoms.
 Center  VDW
      1 1.9080  .1094000
      2 1.4870  .0157000
      3 1.4870  .0157000
      4 1.4870  .0157000
      5 1.4870  .0157000
 Molecular mechanics terms:                 MM functions of 2 of more atoms.
 NBDir 3 1      0      0                    Non-bonded interactions master function.
 HrmStr1 1-2  340.00 1.0900                 Stretch terms.
 HrmStr1 1-3  340.00 1.0900
 HrmStr1 1-4  340.00 1.0900
 HrmStr1 1-5  340.00 1.0900
 HrmBnd1 2-1-3  35.00 109.50                Bend terms.
 HrmBnd1 2-1-4  35.00 109.50
 HrmBnd1 2-1-5  35.00 109.50
 HrmBnd1 3-1-4  35.00 109.50
 HrmBnd1 3-1-5  35.00 109.50
 HrmBnd1 4-1-5  35.00 109.50 
 NBPair 2-1 3 1      0      0 -1.000 -1.000 Non-bonded interactions
 NBPair 3-1 3 1      0      0 -1.000 -1.000 of close neighbors.
 NBPair 3-2 3 1      0      0 -1.000 -1.000
 NBPair 4-1 3 1      0      0 -1.000 -1.000
 NBPair 4-2 3 1      0      0 -1.000 -1.000
 NBPair 4-3 3 1      0      0 -1.000 -1.000
 NBPair 5-1 3 1      0      0 -1.000 -1.000
 NBPair 5-2 3 1      0      0 -1.000 -1.000
 NBPair 5-3 3 1      0      0 -1.000 -1.000
 NBPair 5-4 3 1      0      0 -1.000 -1.000

Atomic parameters. Some MM functions depend only on the atom type of the atom in question. In our example, Amber includes the van der Waals interaction for each atom, and the resulting values are listed in the preceding output under Atomic parameters. The center number corresponds to the atom position within the molecule specification. The DREIDING and UFF force fields contain only terms of this type.

Molecular mechanics terms. The Molecular mechanics terms section list terms within the total potential that describe the interactions of multiple atoms. For example, there are stretch terms involving each bonded pair of atoms (1-2, 1-3 and 1-4), computed via the function called HrmStr1, using a force constant value of 340 kcal/(mol-Å2) and an equilibrium bond length of 1.09 Å. These parameter values are determined from the atom types of centers 1 and 2: CT and HC (respectively).

The particular function and parameters are chosen as follows: from the connectivity, the program knows that a bond exists between centers 1 and 2. The Amber force field calls for a stretch term between all bonded pairs of atoms, and the various functions and parameters are stored in internal tables within the program. In the simplest case, Gaussian uses the two atom types to look up functions and parameters to use. In this case, the entry

HrmStr1 CT HC 340 1.09

is used since it corresponds to the two atom types in the C–H bond in methane. The force constant and equilibrium bond length specified in the entry are then used within the computation. When we look up this function in the reference section below, we see that the actual term will be computed as ForceC*(R–Req)2 where ForceC is the force constant, Req is the equilibrium bond length and R is the bond length; in this case, this becomes 340.0*(R–1.09)2.

The list of terms includes stretch terms for all bonded pairs of atoms and bend terms for all corresponding bond angles.

Non-bonded interaction terms. The remaining terms in the list are those corresponding to non-bonded interactions: contributions to the potential arising from interactions between atoms which are not bonded. Before considering them specifically, we need to look at how they are computed.

In general, the list of non-bonded terms is formed from all possible pairs of atoms, and thus the number of terms scales quadratically with the number of atoms. However, interactions between atoms that are close together, based on the number of bonds that separate them, are usually scaled down. Typically, one- and two-bond separated interactions are scaled to zero (since these interactions are described by stretch and bend terms). Similarly, three-bond separated interactions are scaled with a factor between zero and one, depending on the force field. It is clear that for larger systems, the list of non-bonded terms can become very long. In most MM programs, cutoffs are used the keep the list manageable. In Gaussian, we use a less common approach to efficiently deal with the non-bonded interactions.

The total non-bonded energy in the system can be written as follows:

Total Non-Bonded Energy Expression
Total Non-Bonded Energy Expression

where EvdW and EQ are the van der Waals and electrostatic interactions between the two centers (respectively), and the corresponding s values are the associated scaling factors.

The first equation is expressed in the usual form for the non-bonded interactions, and the second one is an equivalent expression where all possible non-bonded interactions (regardless of distance) are evaluated without scaling (first term) followed by subtraction of the overcounted pairs (second term).

This reformulation has significant computational advantages. The first term can be evaluated in an efficient manner, using linear scaling and other techniques. In addition, most of the pairs in the second term give a zero interaction (provided the system is not too small), because the scaling factor sij is 1 for most pairs. The program simply keeps a list of those pairs that give non-zero interactions to compute the second term. The result is an efficient algorithm that does not use cutoffs, and requires only a list whose length scales linearly with the size of the system.

The functions NBDir and NBTerm are used when forming these two non-bonded interaction terms. As the previous methane job output illustrates, a single NBDir function describes the complete set of interactions between all pairs of atoms, i.e., the first term of the rewritten equation. The various NBTerm entries comprise the second term in the equation. The final two parameters to these functions are the scale factors for van der Waals and Coulomb interactions (respectively). In the methane example, both scale factors are 1 for all terms since no atoms are separated by more than two bonds in methane, and the final non-bonded interaction is accordingly 0.

Specifying missing parameters. Consider the following molecule: 2-methylpropene. We’ve set up an Amber calculation for this system, specifying the SP2 hybridized carbons as atom type CM, and the SP3 hybridized carbons as CT.

2-MethylPropene
2-MethylPropene

# Amber

2-methylpropene

0 1
 C-CM       -2.53    0.19    0.00
 H-HA       -2.00   -0.73    0.00
 H-HA       -3.60    0.19    0.00
 C-CM       -1.86    1.37    0.00
 C-CT       -2.63    2.70    0.00
 H-HC       -1.93    3.51    0.07
 H-HC       -3.24    2.76   -0.88
 H-HC       -3.24    2.76    0.85
 C-CT       -0.32    1.37    0.00
 H-HC        0.03    1.87   -0.83
 H-HC        0.03    1.87    0.81
 H-HC        0.03    0.36   -0.00

When we run this job, the program terminates with the following error message:

 Angle bend undefined between atoms    2    1    3
 Angle bend undefined between atoms    5    4    9
 MM function not complete

Atoms 2-1-3 correspond to types HA-CM-HA, and 5-4-9 correspond to CT-CM-CT. Even though the atom types and the sequence are chemically reasonable, there are no entries for these bending terms in the hard-wired tables. The reason is that the Amber force field has been developed specifically for biochemical systems, in which these particular sequences of atom types do not occur. Hence, the parameters have not been defined by Amber. Note that Gaussian only checks if functions are found for all the stretches and bends; torsional and other terms are simply left out of the total potential function if no entries exist in the tables.

Adding and replacing parameters. When hard-wired parameters are not available, either from incomplete tables or the introduction of new atom types, they need to be specified in the input stream. Determining the appropriate parameters to use can be challenging. They can often be found in the literature; otherwise they need to be derived from experimental or accurate computational data (this topic is beyond the scope of our discussion here).

To specify new parameters, use the Amber=HardFirst keyword, and then add them to the input file in a separate section. In this case, we define HrmBnd1 functions for the missing atom type sequences. The 2-methylpropene job then looks like:

# Amber=HardFirst

2-methylpropene

0 1
 C-CM       -2.53    0.19    0.00
 molecule specification continues…

HrmBnd1 CT CM CT 70.0 120.0
HrmBnd1 HA CM HA 40.0 120.0

The parameter values used in this example are a rough guess based on existing Amber parameters. The HardFirst option indicates that the usual Amber hard-wired parameters tables are checked first for functions—referred to as the hard parameters—and that functions defined in the input file are to be used only when no match is present. With the SoftFirst option, one can replace parameters from the hard-wired set. Note that the ordering specified by HardFirst and SoftFirst only matters when one or more functions is defined in both parameter lists. Finally, the SoftOnly option indicates that the complete force field definition will be provided in the input file, ignoring the hard-wired set completely.

Using external parameter files. It is possible to include MM function entries within an external file by replacing the entries within the input file with the usual Gaussian include file mechanism: e.g., @oplsaa.prm. Parameter files for various force fields are included within the main Gaussian directory (i.e., $g09root/g09). If you want to specify the entire force field via an external file, use an input file like the following:

# UFF=SoftOnly

Read-in force field example

molecule specification

@$g09root/g09/oplsaa.prm

Note that the choice of a molecular mechanics keyword is arbitrary as the entire force field depends on the read-in functions and parameters and not on any internal ones (in other words, specifying Amber would yield the same result as UFF).

Wildcards within function specifications. Wildcards can be used in function specifications given in the input file. For example, this Amber bend function specifies that all the bends with a central carbon atom of type CM will have the same parameter values:

HrmBnd1 * CM * 70.0 120.0

One could also include a generic entry using wildcards as well as other more specific entries:

HrmBnd1 *  CM *  70.0 120.0
HrmBnd1 HA CM HA 40.0 120.0

By default, the program uses the most specific entry that applies to a specific structure (here, a bond angle centered on a CM carbon). Relative specificity is determined by the number of wildcards within an entry. If there happen to be different entries that are equally specific, the program terminates with a multiple matching entries error message. In such cases, you can direct the program to use either the first or last match with the FirstEquiv or LastEquiv options, respectively. Be aware, however, that the multiple matching entries message typically indicates some inconsistency in the soft parameter set which needs to be investigated and resolved.

Note that wildcards/entry specificity do not affect whether an entry from the hard-wired or soft set is used. Thus, when SoftFirst is specified, an entry in the input file with three wildcards will be used even if there is an entry containing the exact atom type sequence in question within the hard-wired set.

Substructures: using bond orders and other structural features to select parameters. Substructures provide a mechanism for using additional structural information like bond orders and bond angle values to determine which parameters are used. Consider butadiene (illustrated below). The carbon atoms are all of type CM, but two bonds are formally double, and one is single. By default, Amber uses the same force constant and equilibrium bond length for the stretching term for each defined atom type combination. However, substructures can be used to assign different values for the different types of bonds.

Substructure numbers are appended to the function name, separated by a hyphen: e.g., HrmStr-1, HrmStr-2, and so on. In this case, the numbers refer to the bond order (see substructure definitions below). Sometimes, multiple substructure suffixes will be used (e.g., AmbTrs-1-2).

The butadiene calculations input file is shown below the molecular structure:

Butadiene Molecule
Butadiene Molecule

# Amber=SoftFirst Geom=Connectivity

Butadiene

0 1
 C-CM    -2.49  -0.07   0.00
 H-HA    -1.96  -1.00   0.00
 H-HA    -3.56  -0.07   0.00
 C-CM    -1.82   1.09   0.00
 H-HA    -2.35   2.02   0.00
 C-CM    -0.28   1.09   0.00
 H-HA     0.24   0.16   0.00
 C-CM     0.39   2.26  -0.00
 H-HA    -0.13   3.19   0.00
 H-HA     1.46   2.26  -0.00

 1 2 1.0 3 1.0 4 2.0
 2
 3
 4 5 1.0 6 1.0
 5
 6 7 1.0 8 2.0
 7
 8 9 1.0 10 1.0
 9
 10

HrmBnd1 HA CM HA 40.0 120.0   Parameters values for illustration purposes only!
HrmBnd1 CM CM CM 70.0 120.0
HrmStr1-1 CM CM 350.0 1.51    Substructures: bond order-specific stretch parameters.
HrmStr1-2 CM CM 500.0 1.37

This example specifies a larger force constant value and a smaller equilibrium bond length for double bonds for use with the HrmStr1 function.

Note that the values used for the parameters in this example are only for illustration purposes only. In fact, for a production calculation, we would also need to define torsional parameters that are specific for single or double central bonds, using the same substructures mechanism.

Besides the bond order, other examples of substructures are whether the term lies in a cyclic structure, how many hydrogen atoms are connected to it, and the like. Note that multiple substructure suffixes can be used with many functions: e.g., HrmStr1-1-0-4.

Substructure slots have relatively consistent meanings in the contexts of different functions. The first substructure suffix defines the bond order or angle range, the second specifies the number of atoms bonded to the current atom, and the third describes the ring environment (if any). A 0 value in any substructure slot acts are a wildcard/placeholder and means that the particular substructure is to be ignored for that entry. Unneeded terminal substructures can be similarly set to 0 or be simply omitted.

The following substructure applies to functions that depend only on the atom type:

Third substructure: ring context (the first and second substructures must be 0):
-0-0-1None of the following ring contexts apply
-0-0-3Atom is in a 3-membered ring
-0-0-4Atom is in a 4-membered ring
-0-0-5  Atom is in a 5-membered ring

The following substructures apply to functions related to bond stretches:

First substructure: bond order
-0Ignore this substructure (substructure “wildcard”)
-1Single bond: 0.00 ≤ bond order < 1.50
-2Double bond: 1.50 ≤ bond order < 2.50
-3Triple bond: bond order ≥ 2.50
Third substructure: ring context (if used, the second substructure must be 0):
-x-0-1None of the following ring contexts apply
-x-0-3Bond is in a 3-membered ring
-x-0-4Bond is in a 4-membered ring
-x-0-5  Bond is in a 5-membered ring

The following substructures apply to functions for bond angles (angle values are in degrees):

First substructure: angle value range
-10° ≤ θ ≤ 45°
-245° < θ ≤ 135°
-3135° < θ ≤ 180°
Second substructure: numbers of bonded atoms
-x-0 Ignore this substructure (substructure “wildcard”)
-x-n n is the number of atoms bonded to the central atom.
Third substructure: ring context
-x-y-0 Ignore this substructure (substructure “wildcard”)
-x-y-1None of the following ring contexts apply
-x-y-3Bend is in a 3-membered ring
-x-y-4Bend is in a 4-membered ring
-x-y-5  Bend is in a 5-membered ring
Fourth substructure: number of bonded hydrogens
-x-y-z-n  There are n+1 hydrogen atoms bonded to the central atom (other than ones comprising the bond angle itself).

The following substructures apply to functions involving dihedral angles.

First substructure: bond order
-0Skip this substructure (substructure “wildcard”)
-1Single central bond: 0.00 ≤ bond order < 1.50
-2Double central bond: 1.50 ≤ bond order < 2.50
-3Triple central bond: bond order ≥ 2.50
Second substructure: bond type of the central bond
-x-0 Skip this substructure (substructure “wildcard”)
-x-1 Resonance central bond (1.30 ≤ bond order ≤ 1.70)
-x-2 Amide central bond (priority over resonance)
-x-3 Other central bond type
Third substructure: ring context
-x-y-1None of the following ring contexts apply
-x-y-4Dihedral is in a 4-membered ring
-x-y-5  Dihedral is in a 5-membered ring

For example, HrmStr1-2-0-4 specifies that the HrmStr1 term applies only to a stretch that involves a double bond that lies in a four-membered ring (compare to HrmStr1-2 which applies to all double bonds). A substructure value of zero is similar to a wildcard: HrnStr1-0-0-4 specifies that the term applies to a stretch of any bond order, but it must be in a four-membered ring.

Note that substructures are not always needed. Some functions explicitly include bond orders and similar structural information within their definitions. For example, the HrmStr3 stretch term from the UFF force field includes the bond order as its third parameter. This parameter can be passed to this function during a UFF calculation for each bond by specifying –1.0 as the value, using a specification like the following in the input file (and the appropriate option to the UFF keyword):

HrmStr3 * * -1.0 0.1332

Defining new atom types. Gaussian understands the standard MM atom types from the Dreiding, UFF, and Amber force fields. However, when a system is studied that these force fields are not parametrized for, none of the available atom types may be suitable. In that case, one can define a new atom type simply by using a new atom type label within the molecule specification. Naturally, since any new atom types will not exist in the hard-wired tables, you will need to also define all required functions using them within the input file.

The non-bonded interaction master function. As we’ve seen, specifying non-bonded interaction terms can be quite lengthy, requiring many NBTerm terms. When non-bonded interactions are specified in the input (using one of the options to the MM keyword), they can be specified using the NonBon function. This function uses a single entry to define the non-bonded interaction, and it is expanded automatically by the program into the appropriate NBDir and NBTerm entries.

This function has the following general form:

NonBon V-Type C-Type V-Cutoff C-Cutoff VScale1 VScale2 VScale3 CScale1 CScale2 CScale3

V-Type and C-Type are integers indicating the kind of van der Waals and Coulomb interactions to be used (respectively). The remaining terms are cutoff values and three scale factors for each of the two interaction types; the latter are scale factors for atoms separated by 1, 2 and 3 or more bonds (respectively). When any scale factor is <0, a scale factor of 1.0/|scale-factor| is used.

As an example, here is the function used for the standard Amber force field:

NonBon 3 1 0 0 0.0 0.0 0.5 0.0 0.0 -1.2

This function specifies the Amber arithmetic van der Waals interaction, 1/R Coulomb interaction, and no cutoffs. Van der Waals interactions are scaled by 0.5 for atoms more than 2 bonds distant and are removed for nearer pairs. Coulomb interactions are scaled by 1.0/1.2 for atoms more than 2 bonds distant and are also removed for nearer pairs

Global parameters. In some cases, a parameter needed for the potential function is not related to a specific atom type at all. One example is the dielectric constant. For such cases, we use global parameters: defined values which may be used as parameters within potential functions. Any defined global parameters are listed in the MM force field definition section of the output (requested with IOp(4/33=3)). Here is an example definition of the dielectric constant:

Dielc  78.39

Global parameters also can be specified in vector or matrix format. The following extract from the MM2 force field definition illustrates the use of a matrix in functions which specify the shift of the hydrogen (i.e., interatomic distance scaling) when calculating van der Waals interactions:

VDWShf1 * 1 0.0      By default, use index value 1 with VDWShf2
VDWShf1 MM5 2 -1.0   Use index value 2 for atom type MM5
VDWShf1 MM36 3 -1.0  Use index value 3 for atom type MM36
VDWShf2 1 2 0.915    Scale factor: matrix element (1,2)
VDWShf2 1 3 0.915    Scale factor: matrix element (1,3)

We’ll consider these two functions in reverse order. VDWShf2 sets up a lower triangular matrix of scale factors for the interatomic distance. The elements (1,2) and (1,3) are both set to 0.915; the other elements—(1,1), (2,2), (2,3) and (3,3)—are left as the default value, meaning that the distance is unscaled in the van der Waals interaction computation. Note that this formulation is independent of atom type.

The VDWShf1 function defines which matrix elements to use for various atom types. The first entry specifies the default index value as 1, and the remaining two entries specify index values 2 and 3 for atom types MM5 and MM36. The third parameter to this function specifies the second center involved in the interaction. In the first entry, the 0.0 value functions as a wildcard, while in the other two entries, the value of –1 tells the program to enter the appropriate center automatically.

These functions result in only bonds involving an atom of type MM5 (H) or MM36 (D) and another of a different type being shifted. When these functions are evaluated, interactions of type MM5-X use element (1,2) and ones of type MM36-X use (1,3). Interactions involving other types use (1,1), while MM5-MM5, MM5-MM36 and MM36-MM36 use (2,2), (2,3) and (3,3), respectively—all of which use unscaled interatomic distances.

Using this mechanism, parameter specification is much clearer and compact. Without it, an entry for every pair of atoms would have been needed to achieve the same goal.

Step-down tables. Wildcards provide a mechanism for specifying generic parameters/defaults that do not explicitly depend on each of the atom types involved in a particular term. Step-down tables have the same purpose, but are more sophisticated. Consider the following example, which is taken from the UFF force field. In these entries, the UseTable and StepDown entries are not potential functions and do not contain parameters, but instead describe the step-down mechanism that is used to process entries that do refer to potential functions and contain parameters.

In the following example, the first set of UFFBnd3 lines define potential functions, while the UFFBnd2 lines and the second set of UFFBnd3 lines use the step-down mechanism:

UseTable UFFBnd2 #3                                  Specify step-down table used for specified functions
UseTable UFFBnd3 #3
StepDown UFFBnd2 0 1 0                               Matching rules for specified functions
StepDown UFFBnd2 0 2 0
StepDown UFFBnd3 0 1 0    
StepDown UFFBnd3 0 2 0    
Table #3 C_R  Trig                                   Define alternate atoms types
Table #3 N_R  Trig      
Table #3 H_   Lin       
Table #3 Li   Lin       
UFFBnd3   0 C_3   0 109.471   -1.0 -1.0 0.1332       Functions defined with
UFFBnd3   0 N_3   0 106.700   -1.0 -1.0 0.1332       explicit atom types
UFFBnd3   0 N_2   0 111.200   -1.0 -1.0 0.1332    
UFFBnd3   0 O_3   0 104.510   -1.0 -1.0 0.1332    
UFFBnd2-0-2   0 Lin   0 2     -1.0 -1.0 0.1332       Functions defined with step-
UFFBnd2-0-3   0 Trig  0 3     -1.0 -1.0 0.1332       down types & substructures
UFFBnd3-0-3   0 Lin   0 180.0 -1.0 -1.0 0.1332
UFFBnd3-0-4   0 Lin   0 180.0 -1.0 -1.0 0.1332
UFFBnd3-0-5   0 Lin   0 180.0 -1.0 -1.0 0.1332
UFFBnd3-0-6   0 Lin   0 180.0 -1.0 -1.0 0.1332
UFFBnd3-0-2   0 Trig  0 120.0 -1.0 -1.0 0.1332
UFFBnd3-0-4   0 Trig  0 120.0 -1.0 -1.0 0.1332
UFFBnd3-0-5   0 Trig  0 120.0 -1.0 -1.0 0.1332
UFFBnd3-0-6   0 Trig  0 120.0 -1.0 -1.0 0.1332

H_ and N_R are UFF atom types in the UFF force field. If there is a bond angle involving such atoms in the molecule—H_–N_R–H_—the program needs to determine the appropriate function and parameters to use. First, it determines the appropriate step-down tables for the available bend function. In this example, table #3 is used for both available bend functions: UFFBnd2 and UFFBnd3.

The Table entries define alternate types for specified atom types. Typically, they map many specific atom types to a more generic pseudo-type (e.g., Lin and Trig in the example). These alternate types allow wide-ranging default values to be specified for function parameters, effectively enabling undefined parameters to be automatically estimated.

The StepDown entries specify matching rules for selecting specific function entries. In this example, both bend functions use the same two entries: 0 1 0 and 0 2 0. They specify that first the program should look for an entry specifying the original atom type as the center atom (zero values function as wild cards). If no such entry is present, then an entry specifying the first alternate atom type as the central atom should be used.

In the example, the program would look for bend functions specifying N_R as the central atom, but none is defined. Next, it would look for functions with Trig as the type of the central atom. There are four such functions: UFFBnd2-0-3 and UFFBnd3-0 with second substructures 2, 4 and 5. Which one is used will depend on the number of atoms bonded to the central nitrogen atom.

Note that in practice, step-down tables can specify multiple alternates for the various atom types and also include many more types mapped to each listed alternate.

When a step-down table is in use, then any wildcards within other function specifications are ignored. In other words, of you choose to use a step-down table, then wildcards cannot be used anywhere else.

AVAILABLE FORCE FIELD TERMS (FUNCTIONS)

Some preliminary notes:


VDW: van der Waals parameters, used for NBDir and NBTerm (See MMFF94 below for MMFF94-type van der Waals parameters)

VDW Atom-type Radius Well-depth


VDW94: MMFF94 type van der Waals parameters (used for NBDir and NBTerm)

VDW94 Atom-type Atomic-pol NE Scale1 Scale2 DFlag

Atomic-pol Atomic polarizability (Angstrom3).
NE Slater-Kirkwood effective number of valence electrons (dimensionless).
Scale1 Scale factor (Angstrom1/4).
Scale2 Scale factor (dimensionless).
DFlag 1.0 for donor type atom, 2.0 for acceptor type, otherwise 0.0.

Buf94: MMFF94 electrostatic buffering

Buf94 Atom-type Value


NonBon: Non-bonded interaction master function. This function will be expanded into pairs and a direct function (NBDir and NBTerm) before evaluation of the MM energy.

NonBon V-Type C-Type, V-Cutoff C-Cutoff VScale1 VScale2 VScale3 CScale1 CScale2 CScale3

V-Type is the van der Waals type:

0 No van der Waals
1 Arithmetic (as for Dreiding)
2 Geometric (as for UFF)
3 Arithmetic (as for Amber)
4 MMFF94-type van der Waals
5 MM2-type van der Waals
6 OPLS-type van der Waals

C-Type is the Coulomb type:

0 No Coulomb
1 1/R
2 1/R2
3 1/R buffered (MMFF94)
10 Form dipole-dipole interactions without cutoffs

V-Cutoff and C-Cutoff are the van der Waals and Coulomb cutoffs (respectively):

0 No cutoff
>0 Hard cutoff
<0 Soft cutoff

Vscale1-3 are van der Waals scale factors for 1- to 3-bond separated pairs. CScale1-3 are Coulomb scale factors for 1- to 3-bond separated pairs. If any scale factor < 0.0, then 1.0/|scale-factor| scaling is used (e.g., a scale factor specified as -1.2 results in scaling of 1.0/1.2).


NBDir: Coulomb and van der Waals direct (evaluated for all atom pairs)

NBDir V-Type C-Type V-Cutoff C-Cutoff

V-Type, C-Type, V-Cutoff, and C-Cutoff as above.


NBTerm: Coulomb and van der Waals single term cutoffs

NBTerm Atom-type1 Atom-type2 V-Type C-Type V-Cutoff C-Cutoff V-Scale C-Scale

V-Type, C-Type, V-Cutoff, C-Cutoff, V-Scale, and C-Scale as above.


AtRad: Atomic single bond radius

AtRad Atom-type Radius


EffChg: Effective charge (UFF)

EffChg Atom-type Charge


EleNeg: GMP Electronegativity (UFF)

EleNeg Atom-type Value


Table: Step down table

Table Original-atom-type Stepping-down-type(s)


HrmStr1: Harmonic stretch I (Amber 1): ForceC*(R–Req)2

HrmStr1 Atom-type1 Atom-type2 ForceC Req

ForceC Force constant
Req Equilibrium bond length

HrmStr2: Harmonic stretch II (Dreiding 4a): ForceC*[R–(Ri+RjDelta)]2

HrmStr2 Atom-type1 Atom-type2 ForceC Delta

ForceC Force constant
Delta Delta

Ri and Rj are atomic bond radii specified with AtRad.


HrmStr3: Harmonic stretch III (UFF 1a): k*(R–Rij)2

Equilibrium bond length Rij = (1 – PropC*lnBO)*(Ri + Rj) – Ren
Force constant: k = 664.12*Zi*Zj/(Rij3)
Electronegativity correction: Ren = Ri*Rj*[SQRT(Xi) – SQRT(Xj)]2/(Xi*Ri + Xj*Rj)

HrmStr3 Atom-type1 Atom-type2 BO PropC

BO Bond order (if <0, it is determined on-the-fly)
PropC Proportionality constant

Ri and Rj are atomic bond radii defined with AtRad. Xi and Xj are GMP electronegativity values defined with EleNeg. Zi and Zj are the effective atomic charges defined with EffChg.


MrsStr1: Morse stretch I (Amber): DLim*(e-a(R–Req)–1)2 where a = SQRT(ForceC/DLim)

MrsStr1 Atom-type1 Atom-type2  ForceC  Req  DLim

ForceC Force constant
Req Equilibrium bond length
DLim Dissociation limit

MrsStr2: Morse stretch II (Dreiding 5a): DLim*(e-a(Ri+Rj–Delta)–1)2 where a = SQRT(ForceC/DLim)

MrsStr2 Atom-type1 Atom-type2 ForceC Delta DLim

ForceC Force constant
Delta Delta
DLim Dissociation limit

Ri and Rj are atomic bond radii defined with AtRad.


MrsStr3: Morse stretch III (UFF 1b): BO*DLim*(e-a(R–Rij)–1)2 where a = SQRT(k/[BO*DLim])

Equilibrium bond length Rij = (1 – PropC*lnBO)*(Ri + Rj) – Ren
Force constant: k = 664.12*Zi*Zj/(Rij3)
Electronegativity correction: Ren = Ri*Rj*[SQRT(Xi) – SQRT(Xj)]2/(Xi*Ri + Xj*Rj)

MrsStr3 Atom-type1 Atom-type2 BO PropC DLim

BO Bond order (if <0, it is determined on-the-fly)
PropC Proportionality constant
DLim Dissociation limit

Ri and Rj are atomic bond radii defined with AtRad. Xi and Xj are GMP electronegativity values defined with EleNeg. Zi and Zj are the effective atomic charges defined with EffChg.


QStr1: Quartic stretch I (MMFF94 2):

(Req/2)*(R–ForceC)2*[1+CStr*(R–ForceC+(7/12)*CStr2*(R–ForceC)2]

QStr1 Atom-type1 Atom-type2 ForceC Req CStr

ForceC Force constant (md-Angstrom-1)
Req Equilibrium bond length
CStr Cubic stretch constant (Angstrom-1)

UFFVOx: Atomic torsional barrier for the oxygen column (UFF 16)

UFFVOx Atom-type Barrier


UFFsp3: Atomic SP3 torsional barrier (UFF 16)

UFFVsp3 Atom-type Barrier


UFFVsp2: Atomic SP2 torsional barrier (UFF 17)

UFFVsp2 Atom-type Barrier


HrmBnd1: Harmonic bend (Amber 1): ForceC*(θ–θeq)2

HrmBnd1 Atom-type1 Atom-type2 Atom-type3  ForceC  θeq

ForceC Force constant (in kcal/(mol*rad2))
θeq Equilibrium angle

HrmBdn2: Harmonic Bend (Dreiding 10a): [ForceC/sin(θeq2)]*(cos(θ)–cos(θeq))2

HrmBnd2 Atom-type1 Atom-type2 Atom-type3  ForceC  θeq

ForceC Force constant
θeq Equilibrium angle

LinBnd1: Dreiding Linear Bend (Dreiding 10c): ForceC*(1+cos(θ))

LinBnd1 Atom-type1 Atom-type2 Atom-type3 ForceC

ForceC Force constant

UFFBnd3: UFF 3-term bend (UFF 11): k*(C0 + C1*cos(θ))+C2*cos(2θ) where C2=1/(4 * sin(θeq2)),
C1 = –4*C2*cos(θeq) and C0=C2*(2*cos(θeq2)+1)

Force constant: k = 664.12*Zi*Zk*(3*Rij*Rjk*(1–cos(θeq2))–cos(θeq)*Rik2)/Rik5

UFFBnd3 Atom-type1 Atom-type2 Atom-type3  θeq  BO12  BO23  PropC

θeq Equilibrium angle
BO12 Bond order for Atom-type1–Atom-type2 (when <0, it is determined on-the-fly)
BO23 Bond order for Atom-type2–Atom-type3 (when <0, it is determined on-the-fly)
PropC Proportionality constant

Ri, Rj and Rk are atomic bond radii defined with AtRad. Xi, Xj and Xk are GMP electronegativity defined with EleNeg. Zi, Zj and Zk are effective atomic charges defined with EffChg.


UFFBnd2: UFF 2-term bend (UFF 10): [k/(Per2)]*[1–cos(Per*θ)]

Force constant: k = 664.12*Zi*Zk*(3*Rij*Rjk*(1–cos(Per2))–cos(Per)*Rik2)/Rik5

UFFBnd2 Atom-type1 Atom-type2 Atom-type3  Per  BO12  BO23  PropC

Per Periodicity: 2 for linear, 3 for trigonal, 4 for square-planar.
BO12 Bond order for Atom-type1–Atom-type2 (when <0, it is determined on-the-fly)
BO23 Bond order for Atom-type2–Atom-type3 (when <0, it is determined on-the-fly)
PropC Proportionality constant

Ri, Rj and Rk are atomic bond radii defined with AtRad. Xi, Xj and Xk are GMP electronegativity defined with EleNeg. Zi, Zj and Zk are effective atomic charges defined with EffChg.


ZeroBnd: Zero bend term: used in rare cases where a bend is zero. This term is needed for the program not to protest about undefined angles.

ZeroBnd Atom-type1 Atom-type2 Atom-type3


CubBnd: Cubic bend (MMFF94 3): (ForceC/2)*(1+CBend*(θ–θeq))*(θ–θeq)2

CubBnd Atom-type1 Atom-type2 Atom-type3  ForceC  θeq  CBend

ForceC Force constant (in md*Angstrom/rad2)
θeq Equilibrium angle
CBend Cubic Bend constant (in deg-1)

LinBnd2: MMFF94 Linear Bend (MMFF94 4): ForceC*(1+cos(θ))

LinBnd2 Atom-type1 Atom-type2 Atom-type3 ForceC

ForceC Force constant (md)

AmbTrs: Amber torsion (Amber 1): Σi=1,4 (Magi*[1+cos(i*θ–I(i+4))])/NPaths

AmbTrs Atom-type1 A-type2 A-type3 A-type4 PO1  PO2  PO3  PO4  Mag1  Mag2  Mag3  Mag4  NPaths

PO1–PO4 Phase offsets
Mag1–Mag4 V/2 magnitudes
NPaths Number of paths. When zero or less, determined on-the-fly.

DreiTrs: Dreiding torsion (Dreiding 13): V*[1–cos(Period*(θ–PO))]/(2*NPaths)

DreiTrs Atom-type1 Atom-type2 Atom-type3 Atom-type4 V PO Period NPaths

V Barrier height V
PO Phase offset
Period Periodicity
NPaths Number of paths. When zero or less, determined on-the-fly.

UFFTorC: UFF torsion with constant barrier height (UFF 15): [V/2]*[1–cos(Period*PO)*cos(V*θ)]/NPaths

UFFTorC Atom-type1 Atom-type2 Atom-type3 Atom-type4 Period PO V NPaths

Period Periodicity
PO Phase offset
V Barrier height V
NPaths Number of paths. When zero or less, determined on-the-fly.

UFFTorB: UFF torsion with bond order based barrier height (UFF 17): [V/2]*[1–cos(Period*PO)* cos(Period*θ)]/NPaths where V = 5*SQRT(Uj*Uk)*[1+4.18*Log(BO12)]

UFFTorB Atom-type1 Atom-type2 Atom-type3 Atom-type4  Period  PO  BO12  NPaths

Period Periodicity
PO Phase offset
BO12 Bond order for Atom-type1–Atom-type2 (when <0, it is determined on-the-fly)
NPaths Number of paths. When zero or less, determined on-the-fly.

Uj and Uk are atomic constants defined with UFFVsp2.


UFFTor1: UFF torsion with atom type-based barrier height (UFF 16): [V/2]*[1–cos(Period*PO)* cos(Period*θ)]/NPaths where V=SQRT(Vj*Vk)

UFFTor1 Atom-type1 Atom-type2 Atom-type3 Atom-type4 Period PO NPaths

Period Periodicity
PO Phase offset
NPaths Number of paths. When zero or less, determined on-the-fly.

Vj and Vk are atomic constants defined with UFFVsp3.


UFFTor2: UFF torsion with atom type based barrier height (UFF 16) (differs from UFFTor1 in the atomic parameter that is used): [V/2]*[1–cos(Period*PO)*cos(Period*θ)]/NPAths where V=SQRT(Vj*Vk)

UFFTor2 Atom-type1 Atom-type2 Atom-type3 Atom-type4 Period PO NPaths

Period Periodicity
PO Phase offset
NPaths Number of paths. When zero or less, determined on-the-fly.

Vj and Vk are atomic constants from UFFVOx.


VDWDreiTRS: Dreiding special torsion for compatibility with Gaussian 98 code. During processing, it is replaced with DreiTRS, with the following parameters:

OldTor Atom-type1 Atom-type2 Atom-type3 Atom-type4


ImpTrs: Improper torsion (Amber 1): Mag*[1+cos(Period*(θ–PO))]

ImpTrs Atom-type1 Atom-type2 Atom-type3 Atom-type4 Mag PO Period

Mag V/2 Magnitude
PO Phase offset
Period Periodicity

Wilson: Three term Wilson angle (Dreiding 28c, UFF 19): ForceC*(C1 + C2*cos(θ) + C3*cos(2θ)) averaged over all three Wilson angles θ

Wilson Atom-type1 Atom-type2 Atom-type3 Atom-type4 ForceC C1 C2 C3

ForceC Force constant
C1, C2, C3 Coefficients

HrmWil1: Harmonic Wilson angle I (MMFF94 6): (ForceC/2)*(θ2) summed over all 3 Wilson angles θ

HrmWil1 Atom-type1 Atom-type2 Atom-type3 Atom-type4 ForceC

ForceC Force constant

MM2Wil: MM2 Wilson sixth order bend (MM2): Σi=1,2,3(ForceCi/2)*(θi2)*[1+6Bend*(θi4)]

MM2Wil Atom-type1 Atom-type2 Atom-type3 Atom-type4 ForceC1 ForceC2 ForceC3 6Bend

ForceC1–ForceC3 Force constants
6Bend Sixth order bend constant (in deg-4)

StrBnd: Stretch-bend (MMFF94 5): (ForceC1*(R12Req12)+ForceC2*(R23Req23))*(θ–θeq)

StrBnd Atom-type1 Atom-type2 Atom-type3  ForceC1  ForceC2  Req12  Req23  θeq

ForceC1, ForceC2 Force constants
Req12 , Req23 Equilibrium bond lengths (if <0, retrieved from the appropriate stretch terms)
θeq Equilibrium angle

CubStr1: Cubic Stretch 1 (MM2): (ForceC/2)*(1–CStr*(R–Req))*(R–Req)2

CubStr1 Atom-type1 Atom-type2  ForceC  Req  CStr

ForceC Force constant
Req Equilibrium bond length
CStr Cubic stretch constant (Angstrom -1)

SixBnd1: Sixth order bend (MM2): (ForceC/2)*(θ–θeq)2(1+6Bend*(θ–θeq)4)

SixBnd1 Atom-type1 Atom-type2 Atom-type3  ForceC  θeq  6Bend

ForceC Force constant
θeq Equilibrium angle
6Bend Sixth order bend constant (in deg-4)

MM2Tors: MM2 Torsional (MM2): En1/2(1+cosθ)+En2/2(1–cos2θ)+En3/2(1+cos3θ)

MM2Tors Atom-type1 Atom-type2 Atom-type3 Atom-type4 En1 En2 En3

En1–En3 Energies

MM2VDW: Parameters for the MM2 van der Waals function (MM2)

MM2VDW Par1 Par1 Par3 Par4 Par5

Par1–Par5 Parameters

VDWX: Matrix with van der Waals parameters (MM2). These values do not depend on atom types, and can be used as alternatives for the atomic parameters, as needed in the MM2 force field (see discussion above).

VDWX Index1 Index2 Radius Well-depth IsDef

IsDef Indicates whether the pair (1.0) has or (0.0) has not been defined.

VDWI: Atomic indices, used to refer to the matrix in VDWXM (MM2)

VDWI Atom-type Index


VDWShf2: Matrix with parameters for atom shifting in the van der Waals calculation (MM2)

VDWShf2 Index1 Index2 Shift


VDWShf1: Atomic indices into the matrix in VDWShf2 (MM2)

VDWShf1 Atom-type Index BondTo

BondTo Center bonded to (if –1.0, determined by program).

SixBndI: Sixth order bend (MM2): (ForceC/2)*(θ–θeq)2(1+6Bend*(θ–θeq)4)

SixBndI Atom-type1 Atom-type2 Atom-type3  ForceC  θeq  6Bend Flag

ForceC Force constant
θeq Equilibrium angle
6Bend Sixth order bend constant (in deg-4)
Flag If ≥ 0.0, suppress projection onto plane even when the central atom is trigonal.

This bend is projected onto the plane in case of a trigonal center. If the center is not trigonal, the regular bend is calculated.


Dipole: Bond dipole

Dipole Atom-type1 Atom-type2 DMom DPos

DMom Dipole moment (in Debye)
DPos Position along the Atom1-Atom2 bond of the dipole

DielC: Dielectric constant. This allows the dielectric constant to be specified via the parameter file; the default is 1.

DielC DielConst

DielConst Dielectric constant

QStr2: Quartic stretch 2 (MM2/Tinker): (ForceC>/2)*(1+CStr*(R–Req)+QStr*(R–Req)2)(R–Req)2)

QStr2 Atom-type1 Atom-type2 ForceC  Req  CStr  QStr

ForceC Force constant
Req Equilibrium bond length
CStr Cubic stretch constant (Angstrom-1)
QStr Quartic stretch constant (Angstrom-2)

CubStr2: Cubic stretch 2 (for compatibility with old MMVB): (ForceC>/2)*(1–CStr*(R–Req))*(R–Req)2)

CubStr2 Atom-type1 Atom-type2  ForceC  Req  CStr

ForceC Force constant
Req Equilibrium bond length
CStr Cubic stretch constant (Angstrom-1; set to 0 when R > 15 Å)

NBonds: Formal number of bonds, based on atom type for this center

NBonds Atom-type NumBnd

NumBnd Formal number of bonds

StrUnit: Units for the force constant in stretching functions

StrUnit Flag

Flag 0 means units of kcal/(mol*Angstrom2); 1 means units of md/Angstrom

BndUnit: Units for the force constant in bending functions

BndUnit Flag

Flag 0 means units of kcal/(mol*rad2); 1 means units of md*Angstrom/rad2

TorUnit: Units for the barrier in torsion functions

TorUnit Flag

Flag 0 means units of kcal/mol; 1 means units of md*Angstrom

OOPUnit: Units for the force constant in out-of-plane bending functions

OOPUnit Flag

Flag 0 means units of kcal/(mol*rad2); 1 means units of md*Angstrom/rad2

SBUnit: Units for the force constant in stretch-bend functions

SBUnit Flag

Flag 0 means units of kcal/(mol*Angstrom*rad); 1 means units of md/rad

StrFact: Factor for the stretching functions

StrFact Value


BndFact: Factor for the bending functions

BndFact Value


TorFact: Factor for the torsion functions

TorFact Value


OOPFact: Factor for the out-of-plane functions

OOPFact Value


SBFact: Factor for the stretch-bend functions

SBFact Value

 


Last update: 5 June 2013