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:
: The Amber force field as described in []. 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 ().
: The Dreiding force field as described in [].
: The UFF force field as described in [].
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 [] 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 []. Request this by specifying the QEq option to the MM keyword: e.g., . Note that these charges depend on the geometry and are computed at the beginnin 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 . 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 [], under control of the following options, for any of the mechanics keywords:
Assign charges to all atoms using the QEq method []. OldQEq requests the older version which was the default in Revision C []. QEq=Uncharged
sets the charges only for atoms which have MM charge zero, constraining the other atoms to retain their specified charges.
Assign QEq charges only to those atoms for which the user did not specify a particular type in the input.
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 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.
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.
Read additional parameters from the input stream, with soft (read-in) parameters having priority over the hard-wired values.
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.
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
or , 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
NewParameters
Ignore any parameters in the checkpoint file when reading the geometry.
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 .
PRINTING MM PARAMETERS
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
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.
Set the Coulomb cutoff to N Angstroms.
Negative for hard cutoffs.
The potential is (QiQj(1-r2/R2)2)/R.
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
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 MMfor MM is usually of the form:
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
keyword (see the discussion of
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)
-0.50 -0.57
0.93 -0.87
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.
Molecular mechanics terms: &&&&&&&&&&&&&&&&MM functions of 2 of more atoms.
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 -1.000 -1.000 Non-bonded interactions
NBPair 3-1 3 1
0 -1.000 -1.000 of close neighbors.
NBPair 3-2 3 1
0 -1.000 -1.000
NBPair 4-1 3 1
0 -1.000 -1.000
NBPair 4-2 3 1
0 -1.000 -1.000
NBPair 4-3 3 1
0 -1.000 -1.000
NBPair 5-1 3 1
0 -1.000 -1.000
NBPair 5-2 3 1
0 -1.000 -1.000
NBPair 5-3 3 1
0 -1.000 -1.000
NBPair 5-4 3 1
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-&A2) and an equilibrium bond length of 1.09 &A. 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 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
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
When we run this job, the program terminates with the following error message:
Angle bend undefined between atoms
Angle bend undefined between atoms
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 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 fou 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
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
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
would yield the same result as ).
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:
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
# Amber=SoftFirst Geom=Connectivity
1 2 1.0 3 1.0 4 2.0
4 5 1.0 6 1.0
6 7 1.0 8 2.0
8 9 1.0 10 1.0
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
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 t 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:
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
Table #3 H_
Table #3 Li
-1.0 -1.0 0.1332 &&&&&&Functions defined with
-1.0 -1.0 0.1332 &&&&&&explicit atom types
-1.0 -1.0 0.1332
-1.0 -1.0 0.1332
UFFBnd2-0-2
-1.0 -1.0 0.1332 &&&&&&Functions defined with step-
UFFBnd2-0-3
-1.0 -1.0 0.1332 &&&&&&down types & substructures
UFFBnd3-0-3
0 180.0 -1.0 -1.0 0.1332
UFFBnd3-0-4
0 180.0 -1.0 -1.0 0.1332
UFFBnd3-0-5
0 180.0 -1.0 -1.0 0.1332
UFFBnd3-0-6
0 180.0 -1.0 -1.0 0.1332
UFFBnd3-0-2
0 120.0 -1.0 -1.0 0.1332
UFFBnd3-0-4
0 120.0 -1.0 -1.0 0.1332
UFFBnd3-0-5
0 120.0 -1.0 -1.0 0.1332
UFFBnd3-0-6
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:
Units: Unless otherwise indicated, distances are in Angstroms, angles are in degrees, energies are in kcal/mol and charges are in atomic units.
In equations, R refers to distances and & refers to angles. Numerical and variable subscripts refer to centers, and &eq& refers to equilibrium: e.g., R12 and Rij are the bond distances between atoms 1 and 2 and atoms i and j (respectively), Ri is the bond radii for atom i, Req and &eq are the equilibrium bond length and bond angle, and Req12 is the equilibrium bond length between atoms 1 and 2 (typically used in an angle term).
Wildcards may be used in any function definition. They are indicated by a 0 or an asterisk.
Function equivalencies to those found in standard force fields are indicated in parentheses following the description.
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
Atomic-pol
Atomic polarizability (Angstrom3).
Slater-Kirkwood effective number of valence electrons (dimensionless).
Scale factor (Angstrom1/4).
Scale factor (dimensionless).
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:
No van der Waals
Arithmetic (as for Dreiding)
Geometric (as for UFF)
Arithmetic (as for Amber)
MMFF94-type van der Waals
MM2-type van der Waals
OPLS-type van der Waals
C-Type is the Coulomb type:
No Coulomb
1/R buffered (MMFF94)
Form dipole-dipole interactions without cutoffs
V-Cutoff and C-Cutoff are the van der Waals and Coulomb cutoffs (respectively):
Hard cutoff
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
Force constant
Equilibrium bond length
HrmStr2: Harmonic stretch II (Dreiding 4a): ForceC*[R&(Ri+Rj&Delta)]2
HrmStr2 Atom-type1 Atom-type2
ForceC Delta
Force constant
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
Bond order (if &0, it is determined on-the-fly)
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&
Force constant
Equilibrium bond length
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
Force constant
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
PropC DLim
Bond order (if &0, it is determined on-the-fly)
Proportionality constant
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
Force constant (md-Angstrom-1)
Equilibrium bond length
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&
Force constant (in kcal/(mol*rad2))
Equilibrium angle
HrmBdn2: Harmonic Bend (Dreiding 10a): [ForceC/sin(&eq2)]*(cos(&)&cos(&eq))2
HrmBnd2 Atom-type1 Atom-type2 Atom-type3
&ForceC &&eq
Force constant
Equilibrium angle
LinBnd1: Dreiding Linear Bend (Dreiding 10c): ForceC*(1+cos(&))
LinBnd1 Atom-type1 Atom-type2 Atom-type3 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
Equilibrium angle
Bond order for Atom-type1&Atom-type2 (when &0, it is determined on-the-fly)
Bond order for Atom-type2&Atom-type3 (when &0, it is determined on-the-fly)
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&
Periodicity: 2 for linear, 3 for trigonal, 4 for square-planar.
Bond order for Atom-type1&Atom-type2 (when &0, it is determined on-the-fly)
Bond order for Atom-type2&Atom-type3 (when &0, it is determined on-the-fly)
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
Force constant (in md*Angstrom/rad2)
Equilibrium angle
Cubic Bend constant (in deg-1)
LinBnd2: MMFF94 Linear Bend (MMFF94 4): ForceC*(1+cos(&))
LinBnd2 Atom-type1 Atom-type2 Atom-type3 ForceC
Force constant (md)
AmbTrs: Amber torsion (Amber 1): &Si=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
Phase offsets
V/2 magnitudes
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
Barrier height V
Phase offset
Periodicity
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
Periodicity
Phase offset
Barrier height V
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
Periodicity
Phase offset
Bond order for Atom-type1&Atom-type2 (when &0, it is determined on-the-fly)
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
Periodicity
Phase offset
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
Periodicity
Phase offset
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:
If there are three atoms bonded to the third center and the fourth center is H, it is removed.
If there are three atoms bonded to the third center, and at least one of them is H, but the fourth center is not H, then these values are used: V=4.0, PO=0.0, Period=3.0, and NPaths=&1.0.
Otherwise, these values are used: V=1.0, PO=0.0, Period=6.0, and NPaths=&1.0.
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
V/2 Magnitude
Phase offset
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
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
Force constant
MM2Wil: MM2 Wilson sixth order bend (MM2): &Si=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
Sixth order bend constant (in deg-4)
StrBnd: Stretch-bend (MMFF94 5): (ForceC1*(R12&Req12)+ForceC2*(R23&Req23))*(&&&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)
Equilibrium angle
CubStr1: Cubic Stretch 1 (MM2): (ForceC/2)*(1&CStr*(R&Req))*(R&Req)2
CubStr1 Atom-type1 Atom-type2 &ForceC& Req& CStr
Force constant
Equilibrium bond length
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
Force constant
Equilibrium angle
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
MM2VDW: Parameters for the MM2 van der Waals function (MM2)
MM2VDW Par1 Par1 Par3 Par4 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
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)
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
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
Force constant
Equilibrium angle
Sixth order bend constant (in deg-4)
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
Dipole moment (in Debye)
Position along the Atom1-Atom2 bond of the dipole
DielC: Dielectric constant. This allows the dielectric constant to be specified vi the default is 1.
DielC 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
Force constant
Equilibrium bond length
Cubic stretch constant (Angstrom-1)
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
Force constant
Equilibrium bond length
Cubic stretch constant (Angstrom-1; set to 0 when R & 15 &A)
NBonds: Formal number of bonds, based on atom type for this center
NBonds Atom-type NumBnd
Formal number of bonds
StrUnit: Units for the force constant in stretching functions
StrUnit 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
0 means units of kcal/(mol*rad2); 1 means units of md*Angstrom/rad2
TorUnit: Units for the barrier in torsion functions
TorUnit Flag
0 means units of kcal/ 1 means units of md*Angstrom
OOPUnit: Units for the force constant in out-of-plane bending functions
OOPUnit 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
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