SourceForge.net Logo
MCCCS Towhee: towhee_ff

Overview
    This section describes the format that is used to specify Towhee force field files (towhee_ff) for the current version of those file (currently force field version 15, starting with code version 7.1.0). This documentation was last updated for version 8.2.3. Files for many commonly used force fields are provided with the download package and if you are using one of those you do not need to know the details of this section. However, this feature also allows the user to utilize their own force field parameters without having to modify the code themselves. The formatting for this file is extremely sensitive and the user is encouraged to look at the examples in the ForceFields directory. The names of the variables must be listed precisely as shown in this file in order to facilitate error checking of the force field file.
Differences from previous versions
  • Version 14: had a vibcoeff(0) value listed for the Square Well Bond Type.
  • Version 13: did not have the Vibration Order variable for each vibration type, the Angle Order variable for each angle type, the Torsion Order variable for each torsion type, or the Bond Increment Order for each bond increment type. It also included the deprecated Bond-Angle Strings variables for certain class 2 angle potentials.
  • Version 12: did not have the Polarizability variable for each nonbonded atom type and had a typo in the Number of One-Five types label.
  • Version 11: included a coefficient for the 'Nonbonded Potential only' Torsion style that was not used for anything.
  • Version 10: included a third term for the 'Simple Harmonic Potential' Torsion style that was not used.
  • Version 9: had a slightly different potential for the #9 bending angle style.
  • Version 8: did not allow the embedding density to be dependent upon both atoms in the pair.
  • Version 7: did not have any One-Five interaction information.
  • Version 6: had a the Cross Term Logical instead of specifying the acceptable classical_mixrule options. Also did not contain the bond increment information.
  • Version 5: used the old format of specifying the potential form as an integer instead of as a text string
  • Version 4: had a completely separate file format for the Embedded Atom Method potential. This has now been incorporated into a single version. The nbcoeffs for the implicit solvation potential used a different format in this version and is no longer compatable with current versions.
  • Version 3: allowed the specification of multiple atoms with the same nonbond parameters and did not utilize the vibration, angle, and torsion name definations. Version 3 and earlier are no longer valid inputs into Towhee.
  • Version 2: used to allow the specification of multiple force fields for a single set of parameters. This feature was made obsolete by other changes in the code and was removed in version 3.
  • Version 1: did not have the Cross Term Logical variable. Instead it was assumed from the classical_potential.
Variable explanations for towhee_ff The variable names must all be written exactly as follows. All entries provide the variable name of the descriptive character string, the FORTRAN formatting for the descriptive character string, and the FORTRAN formatting for the actual variable entries.
    'towhee_ff Version' (a17) [integer]
    • The Version number of the force field file. This documentation is accurate for version 3.
    'Number of Nonbonded Types' (a25) [integer]
    • The number of nonbonded atom types in the force field file.
    'Potential Type' (a14) [character string a30]
    • The nonbonded potential type. See the classical_potential documentation for information about how to set this variable. Only the classical_potential is needed here and not all of the variables associated with the classical_potential.
    'Classical Mixrule' (a17) [character string a30]
    • The classical mixing rule that works in combination with this forcefield. Generally this must match exactly with the classical_mixrule variable (set in towhee_input). However, certain special combination cases are also allowed. Currently the only combination case is 'LB or Geometric' which matches with a classical_mixrule of 'Lorentz-Berthelot' or 'Geometric'.
    This section is repeated once for each of the Number of Nonbonded Types.
      'Atom Type Number' (a16) [integer]
      • The number coresponding to the atom type. These numbers must run from 1 to the Number of Nonbonded Types.
      The following set of variables are repeated once for each atom type pair represented here. If the 'Classical Mixrule' is 'Explicit' then there are multiple listings (interaction of this Atom Type Number with itself and with all other atoms that have a larger Atom Type Number). For any other value of the 'Classical Mixrule' there is only a single listing (interaction of this Atom Type Number with itself).
        'eam_pair_style' (a14) [character string a20]
          This variable is only listed if Potential Type is set to 6 (Embedded Atom Method). This is the style of the EAM pair interaction for this combination of atom types. Currently supported options are listed here.
        • 'Ackland 3-part': The pair potential is a complicated 3 part function as described in Mendelev et al. 2003. If r < nbcoeff(1)
            x = r/nbcoeff(4)
            upair = ( 0.1818 * exp[-3.2 * x] + 0.5099 * exp[-0.9423 * x] + 0.2802 * exp[-0.4029 * x] + 0.02817 * exp[-0.2016 * x] ) * nbcoeff(3) / r
          elseif r < nbcoeff(2) then
            x = nbcoeff(5) + nbcoeff(6) * r + nbcoeff(7) * r2 + nbcoeff(8) * r3
            upair = nbcoeff(9) * exp[x]
          else
            upair = ∑i = 1,table_npair table_pair(2,i) * [table_pair(1,i) - r]3 * Heaviside[table_pair(1,i)-r]
        • 'Ackland Power': The pair potential is a complicated 3 part function as described in Ackland et al. 2004.
          if r < nbcoeff(1) then
            upair = ( nbcoeff(4) * exp[nbcoeff(5) * r] + nbcoeff(6) * exp[nbcoeff(7) * r] + nbcoeff(8) * exp[nbcoeff(9) * r] + nbcoeff(10) * exp[nbcoeff(11) * r] ) * nbcoeff(3) / r
          elseif r < nbcoeff(2) then
            x = nbcoeff(12) + nbcoeff(13) * r + nbcoeff(14) * r2 + nbcoeff(15) * r3
            upair = nbcoeff(16) * exp[x]
          else
            upair = ∑i = 1,table_npair table_pair(2,i) * [table_pair(1,i) - r]i+3 * Heaviside[table_pair(1,i)-r]
        • 'Belashchenko 5-8': Two part pair potential as described in Belashchenko 2013.
          if r < r1 (where r1 is nbcoeff(1)) then
            upair = nbcoeff(7) + nbcoeff(8) * (r1 - r) + nbcoeff(9)*(exp[nbcoeff(10)*(r1-r)]-1)
          else
            upair = ∑i = 1,5 ∑m=0,8 ai,m (r - ri+1)m Heaviside[ri,ri+1]
          where the ri are stored in nbcoeff(i) for i=1,6
          and ai,m are stored in nbcoeff(2 + 9*i + m) for i = 1,5 and m = 0,8
        • 'exponential': The pair potential is described by an exponential function.
          upair = nbcoeff(1) e[nbcoeff(2) rij].
        • 'morse': The pair potential is described by the Morse function.
          upair = nbcoeff(1) * { exp[-2 nbcoeff(2) (r - nbcoeff(3))] - 2 exp[-nbcoeff(2) (r-nbcoeff(3))]}
        • 'table': This is the traditional way of describing EAM potential information and presents a table of data that is interpolated to determine the energy for any distance.
        'table_pair' (a11) [integer, integer, integer]
          This variable is only listed for tabular potentials. This occurs in the following cases
          • classical_potential 'Embedded Atom Method' and eam_pair_style 'Ackland 3-part'
          • classical_potential 'Embedded Atom Method' and eam_pair_style 'Ackland Power'
          • classical_potential 'Embedded Atom Method' and eam_pair_style 'table'
          • classical_potential 'Multiwell'.
          • classical_potential 'Repulsive Multiwell'.
          • classical_potential 'Tabulated Pair'
        • Three values on a single line. First atom type, second atom type, and the number of tabular data lines of table_pair_data
        'table_pair_data' (a15) [double precision, double precision]
        • The tabular values used to determine the pair potential. This is only listed for tabular potentials. There are two entries per line (table_pair(1,n) and table_pair(2,n)) and the number of lines is equal to the value set in table_pair
        'Nonbond Coefficients' (a20) [double precision]
        • These values are listed for any non-tabular potential and they contain all of the nonbonded coefficients. These are listed with one coefficient per line. The number of nonbonded coefficients depends upon the classical_potential. The native units of the code are used (Kelvin for energy, Å for distances). If the 'Classical Mixrule' is 'Explicit' then the cross terms must be explicitly declared for all interactions of this atom type with any atom type that has a larger Atom Type Number. Otherwise the nonbonded coefficients are listed only for self-interaction and the mixing rule will set the cross terms.
      End of the section that is repeated based on the Number of Atom Type Pairs.
      The following section is only included for classical_potential 'Embedded Atom Method'. The density contributed to all of the other atom types is listed here for the current atom type.
      'eam_dens' (a8) [integer, integer, integer]
      • Three values on a single line. Atom type the density is contributed to, atom type the density is contributed from, and number of lines of density coefficients
      'eam_dens_style' (a14) [character string a20]
        The style of the EAM density interaction for this atom type. There are several options currently supported by Towhee.
      • 'Ackland cubic sum'
        ρ(rij) = ∑i = 1,eam_ndens eam_dens(2,i) * [eam_dens(1,i) - r]3 * Heaviside[eam_dens(1,i)-r]
      • 'exponential': The density is described by an exponential function.
        ρ(rij) = eam_dens(1,2) eeam_dens(2,2) rij for rij < eam_dens(2,1), and 0 otherwise.
      • 'table': This is the traditional way of describing EAM density information and presents a table of data that is interpolated to determine the density for any distance.
      'eam_dens_data' (a14) [double precision, double precision]
        The values used to determine the embedding density. The meanings are different depending on the setting of eam_dens_style above. In all cases there are a number of lines equal to the value set in eam_dens
      • 'table': Each line contains two values. The first is the distance (Å) and the second is the corresponding embedding density (arbitrary units).
      • 'exponential': The embedding density is described by an exponential function. The format is a bit clumsy as the first line contains a dummy value and then the A prefactor (arbitrary units). The second line contains the cutoff distance for the density potential (Å) and then the B exponential factor (inverse Å).
      • 'Ackland cubic sum': Each line contains two values. The first is a distance (Å) and the second is a density prefactor (arbitrary units).
      'eam_embed' (a9) [integer, integer]
      • Only listed for classical_potential 'Embedded Atom Method'. Two values on a single line. Atom type and number of lines of embedding function coefficients.
      'eam_embed_style' (a15) [character string a20]
        Only listed for classical_potential 'Embedded Atom Method'. The style of the EAM embedding function for this atom type. Currently supported options are listed below.
      • 'Belashchenko': uembed(ρ(rij))
          = eam_embed_data(4,1) * √  ρ (rij) + eam_embed_data(4,2) * ρ(rij)
            for 0 ≤ ρ(rij) < eam_embed_data(5,1)
          = eam_embed_data(1,2) + eam_embed_data(2,2) * [ρ(rij) - eam_embed_data(5,2)]2 + eam_embed_data(3,2)*[ρ(rij) - eam_embed_data(5,2)]3
            for eam_embed_data(5,1) ≤ ρ(rij)
      • 'Belashchenko 10 rho':
          belrho(0) = eam_embed_data(11,1)
          belrho(1) = eam_embed_data(11,2)
          belrho(2) = eam_embed_data(11,3)
          belrho(3) = eam_embed_data(12,1)
          belrho(4) = eam_embed_data(12,2)
          belrho(5) = eam_embed_data(12,3)
          belrho(6) = eam_embed_data(13,1)
          belrho(7) = eam_embed_data(13,2)
          belrho(8) = eam_embed_data(13,3)
          belrho(9) = eam_embed_data(14,1)
        uembed(ρ(rij))
          = [eam_embed_data(8,1) + eam_embed_data(8,2) * (ρ(rij) - belrho(7)) + eam_embed_data(8,3) * (ρ(rij) - belrho(7))2] * [2*ρ(rij)/belrho(7) - {ρ(rij)/belrho(7)}2]
            for 0 < ρ(rij) < belrho(7)
          = eam_embed_data(i,1) + (ρ(rij) - belrho(i-1)) * [eam_embed_data(i,2) + eam_embed_data(i,3)*(ρ(rij)-belrho(i-1))]
            for belrho(i) ≤ ρ(rij) < belrho(i-1) where 2 ≤ i ≤ 7
          = eam_embed_data(1,1) + eam_embed_data(3,1) * (ρ(rij) - belrho(0))2
            for belrho(1) ≤ ρ(rij) < belrho(8).
          = eam_embed_data(9,1) + eam_embed_data(9,2) * [ρ(rij) - belrho(8)] + eam_embed_data(9,3)*[ρ(rij) - belrho(8)]eam_embed_data(14,2)
            for belrho(8) ≤ ρ(rij) < belrho(9)
          = eam_embed_data(10,1) + eam_embed_data(10,2) * [ρ(rij) - belrho(9)] + eam_embed_data(10,3)*[ρ(rij) - belrho(9)]eam_embed_data(14,3)
            for belrho(9) ≤ ρ(rij)
      • 'logarithmic': uembed(ρ(rij)) = ρ(rij) * (eam_embed_data(1,2) * ln[ρ(rij)] + eam_embed_data(2,2))
      • 'power 0.5 and 2': uembed(ρ(rij)) = eam_embed_data(1,2) √  ρ (rij) + eam_embed_data(2,2) *[ρ(rij)]2
      • 'power 0.5, 2, and 4': uembed(ρ(rij)) = eam_embed_data(1,2) √  ρ (rij) + eam_embed_data(2,2) *[ρ(rij)]2 + eam_embed_data(3,2) *[ρ(rij)]4
      • 'squareroot': uebmed(ρ(rij)) = eam_embed_data(2,2) * ρ(rij) - eam_embed_data(1,2) √  ρ (rij)
      • 'table': This is the traditional way of describing EAM embedding function information and presents a table of data that is interpolated to determine the embedding function energy for any density.
      'eam_embed_data' (a15) [double precision, double precision]
        Only listed for classical_potential 'Embedded Atom Method'. The values used to determine the embedding function energy. The meanings are different depending on the setting of eam_embed_style above. In all cases there are a number of lines equal to the value set in eam_embed
      • 'table': Each line contains two values. The first is the density (arbitrary units) and the second is the corresponding embedding function energy (Kelvin/kb).
      • 'Ackland 2-term': A single line that contains a prefactor (eam_embed(1,1)) and a square factor (eam_embed(2,1) as described above.
      • 'Belashchenko': 5 lines where the first 4 lines contain the energy terms (K) and the fifth line contains the density terms (arbitrary units) as described above.
      End of the special 'Embedded Atom Method' section.
      'Mass' (a4) [double precision]
      • The atomic mass of this element. Units are grams/mol.
      'Element' (a7) [character string a2]
      • The elemental name of this Atom Type.
      'Bond Pattern' (a12) [character string a5]
      • The bond pattern for this Atom Type. The bond pattern is used with configurational-bias moves that utilize a non-uniform generation of the bond angles and dihedrals. See the configurational-bias manual for more information. Set to the word 'null' if you do not want to use this feature.
      'Base Charge' (a12) [double precision]
      • The base charge on this atom type for use with the 'bond increment' method of automatically assigning charges. This value is zero for most atom types, but can take on different values for ionic systems.
      'Polarizability' (a14) [double precision]
      • The polarizability of each atom type for use with some of the classical mixing rules. Expressed in units of cubic Å.
      'Force Field Name' (a16) [character string a10]
      • The name of the Force Field
      'Atom Names' (a10) [character string a10 on four lines]
      • Each Atom Type has a nonbonded name, a bonded name, an angle name, and a torsion name. These are listed in that order, one per line. These names are used in combination with the assemble features in Towhee to make it easier to build molecule structures in the towhee_input file.
    End of the section that is repeated based on the Number of Nonbonded Types.
    'Number of Bonded Terms' (a22) [integer]
    • The number of different Bonded terms in the force field file.
    This section is repeated once for each of the Number of Bonded Terms
      'Bond Type Number' (a16) [integer]
      • The number coresponding to the bond type. These numbers must run from 1 to the Number of Bonded Terms.
      'Bond Style' (a10) [integer]
      • The number coresponding to the bond style. Below is a list of Bond Styles, their potential form, and the required vibcoeffs. Previously, the configurational-bias algorithm required the equilibrium bond length term for any bond potential to be stored in vibcoeff(0), but this is no longer mandatory starting with version 7.1.0.
        • 1: Fixed Bond Length
            Ubond = Infinity, if ( length - vibcoeff(0) ) / vibcoeff(0) > 0.01
            Ubond = 0 otherwise.
        • 2: Standard Harmonic
            Ubond = vibcoeff(1)*[ length - vibcoeff(0) ]2
        • 3: Gromos Quartic
            Ubond = vibcoeff(1)*[ length2 - vibcoeff(0)2 ]2
        • 4: Nonlinear
            distance = [ length - vibcoeff(0) ]2
            Ubond = [ vibcoeff(1)*distance ] / [ vibcoeff(2) - distance ]
        • 5: MM2 Quadradic and Triplet
            distance = [ length - vibcoeff(0) ]
            Ubond = vibcoeff(1)*distance2 * [ 1.0 - 2.0*distance ], if distance ≤ 1/3
            Ubond = Infinity, if distance > 1/3
        • 6: Compass Quartic
            distance = [ length - vibcoeff(0) ]
            Ubond = vibcoeff(1)*distance2 + vibcoeff(2)*distance3 + vibcoeff(3)*distance4
        • 7: Nonbonded interactions used as the Bonding potential
            This has no vibcoeffs as it uses the nonbond van der Waals and coulombic potentials instead.
        • 8: No interaction
            Zero energy regardless of interatomic separation.
        • 9: Morse
            distance = [ length - vibcoeff(0) ]
            Ubond = vibcoeff(1)*[ evibcoeff(2)*distance - 1 ]2
        • 10: Infinite Square Well
            vibcoeff(0) is no longer listed for this potential starting with version 7.1.0
            Ubond = infinity if length < vibcoeff(1)
            Ubond = vibcoeff(3) if vibcoeff(1) ≤ length ≤ vibcoeff(2)
            Ubond = infinity if vibcoeff(2) < length
        • 11: Standard Harmonic plus nonbonded interactions
            Ubond = vibcoeff(1)*[ length - vibcoeff(0) ]2 + Unonbond(rij) + Uelectrostatic(rij)
        • 12: FENE plus nonbonded interactions (Reference: Kremer and Grest, 1990
            Ubond = - 0.5 * vibcoeff(2) * vibcoeff(1)2 * ln[ 1.0 - (length/(vibcoeff(1)))2 ] + Unonbond(rij) + Uelectrostatic(rij) if length < vibcoeff(1).
            Ubond = Infinity otherwise.
      'Bond Coefficients' (a17) [double precision]
      • The values of the vibcoeffs listed in order (starting with 0) with a single coefficient per line.
      'Vibration Order' (a15) [character string a10]
      • This is used to distinguish between bond types for atoms that can make more than one kind of bond (for example, sp2 carbons bond to each other with either single or double bonds). Commonly used options are listed in the advanced connectivity map manual.
      'Force Field Name' (a16) [character string a10]
      • The name of the Force Field
      'Number of Atoms with Same Parameters' (a38) [integer]
      • The number of pairs of atoms in this force field which use the same bond parameters.
      'Atom Names' (a10) [pairs of character strings a10,1x,a10]
      • All of the pairs of Atom Names for this Bond Type. These are listed with one pair of atom names per line.
    End of the section repeated based on the Number of Bonded Terms
    'Number of Angle Terms' (a21) [integer]
    • The number of different Angle terms in the force field file.
    This section is repeated based on the Number of Angle Terms
      'Angle Type Number' (a17) [integer]
      • The number coresponding to the angle type. These numbers must run from 1 to the Number of Angle Terms.
      'Angle Style' (a11) [integer]
      • The number coresponding to the angle style. Below is a list of Angle Styles,their potential form, and the bencoeffs. All of the angles are listed using degrees in the force field file, and are converted into radians for use in the code. Thus, all of the force constants need to be appropriate for use with radians, while the angles are listed in degrees. The equilibrium angle is always listed in bencoeff(0) and only this value is converted from degrees to radians.
        • 0: Rigid Angle
            Uangle(θ) = Infinity, if |(θ - bencoeff(0))| > bencoeff(1)
            Uangle(θ) = 0 otherwise.
        • 1: Standard Harmonic
            Uangle(θ) = bencoeff(1)*[ θ - bencoeff(0) ]2
        • 2: DREIDING 1 + Cos(θ)
            Uangle(θ) = bencoeff(1)*[ 1.0 + Cos(θ) ]
        • 3: Harmonic Cosine
            Uangle(θ) = bencoeff(1)*[ Cos(θ) - Cos(bencoeff(0)) ]2
        • 4: Compass Quartic Angle with autodetection
            The Bond-Angle and Bond-Bond (1-2) Class 2 cross terms are incorporated into this angle and that makes the potential assymetric with respect to the atom order. This uses the Compass atom type definitions to decide the order of the terms. Then, it finds the appropriate bond terms for use in the bond-angle and bond-bond cross terms. Two logicals control whether there are any Bond-Angle or Bond-Bond cross terms. Definitions: rij is the distance between the first and second atoms in the angle, rjk is the distance between the second and third atoms of the angle, and the appropriate vibcoeff(0) values are determined automatically.
            differ = θ - bencoeff(0)
            Uangle(θ) = bencoeff(1)*differ2 + bencoeff(2)*differ3 + bencoeff(3)*differ4 + Ubond-angle(rij,rjk,θ) + Ubond-bond(rij,rjk)
            if Bond-Angle Logical is true then
              Ubond-angle(θ,rij,rjk) = differ*[ bencoeff(4)*(rij - vibcoeff(0)ij) +bencoeff(5)*(rjk - vibcoeff(1)jk) ]
            otherwise
              Ubond-angle(rij,rjk,θ) = 0
            if Bond-Bond Logical is true then
              Ubond-bond(rij,rjk) = bencoeff(6)*[ (rij - vibcoeff(0)ij) *(rkj - vibcoeff(0)jk) ]
            otherwise
              Ubond-bond(rij,rjk) = 0
        • 5: Charmm Harmonic with Urey-Bradley
            Uangle(θ,dik) = bencoeff(1)*[ θ - bencoeff(0) ]2 + bencoeff(3)*[dik - bencoeff(2) ]2
            where dik is the distance between the first and third atoms of the angle.
        • 6: Nonbonded Terms only between the 1-3 atoms of the angle
            The angle energies are actually the nonbonded van der Waals plus the nonbonded coulomb terms. Towhee requires an angle term between all atoms connected by two bonds, so this term is used if you actually want just a nonbonded interaction (for example, with a freely jointed chain). There are no bencoeffs for this type.
        • 7: Nonbonded Terms between the 1-3 atoms plus a harmonic angle
            The angle energies are the nonbonded van der Waals, the nonbonded coulomb terms, and the following.
              Uangle(θ) = bencoeff(1)*[ θ - bencoeff(0)]2
        • 8: Compass Quartic Angle with explicit ordering of terms
            The Bond-Angle and Bond-Bond (1-2) Class 2 cross terms are incorporated into this angle and that makes the potential assymetric with respect to the atom order. This version uses a positive or negative convenction on the type numbers (assigned in towhee_input) to tell the code whether to follow the order as presented here (positive) or to use the inverse order (negative). Two logicals control whether there are any Bond-Angle or Bond-Bond cross terms.
            differ = angle-bencoeff(0)
            Uangle(θ) = bencoeff(1)*differ2 + bencoeff(2)*differ3 + bencoeff(3)*differ4
            if there are Bond-Angle cross terms then
            Ubond-angle(θ,dij,djk) = [ θ - bencoeff(0) ] * {bencoeff(4)*[ dij - bencoeff(5) ] + bencoeff(6)*[ djk - bencoeff(7) ]}
            Ubond-bond(dij,djk) = bencoeff(8) * [ dij - bencoeff(9) ] * [ djk - bencoeff(10) ]
        • 9: Fourier Expansion with constant plus single term
            This was originally implemented for use with the UFF forcefield. The prefactor is computed in a much more elaborate manner than is usual for angle force constants and only works in combination with the UFF 12-6 classical_potential. This prefactor is computed for each angle in the simulation during the structure check at the start of a simulation. Note: this potential is not well designed as it also has a minimum for a nonphysical θ value of zero.
            θ = bencoeff(0)
            rik = [vibcoeff(ij,0)2*vibcoeff(jk,0)2 - 2*vibcoeff(ij,0)*vibcoeff(jk,0)*Cos(θ) ]1/2
            benprefact(ijk) = kcaltok*664.12*nbcoeff(10,i)*nbcoeff(10,k)/(rik5) * { vibcoeff(ij,0)*vibcoeff(jk,0)*[1 - Cos2(θ)] - vibcoeff(ik,0)2*Cos(θ) } / bencoeff(1)2
            Uangle(θ) = benprefact(ijk) * [1 + bencoeff(2) * Cos(bencoeff(1)*θ)]
        • 10: Three Term Fourier Expansion
            This was originally implemented for use with the UFF forcefield. The prefactor is computed in a much more elaborate manner than is usual for angle force constants and only works in combination with the UFF 12-6 classical_potential. This prefactor is computed for each angle during the structure check at the start of a simulation. Note: this prefactor is nearly identical to the one for anglestyle 9, with the notable exception of not dividing by the bencoeff(1)2 term.
            θ0 = bencoeff(0)
            rik = [vibcoeff(ij,0)2*vibcoeff(jk,0)2 - 2*vibcoeff(ij,0)*vibcoeff(jk,0)*Cos(θ0) ]1/2
            benprefact(ijk) = kcaltok*664.12*nbcoeff(10,i)*nbcoeff(10,k)/(rik5) * { vibcoeff(ij,0)*vibcoeff(jk,0)*[1 - Cos2(θ0)] - vibcoeff(ik,0)2*Cos(θ0) }
            Uangle(θ) = benprefact(ijk) * [bencoeff(1) + bencoeff(2)*Cos(θ) + bencoeff(3)*Cos(2θ)]
        • 11: No interaction
            The angle energy is zero regardless of the distance or angle between the atoms.
        • 12: Sixth Power Angle with Autodetection of Angle-Bond Terms
            The Bond-Angle term is incorporated into this angle and that makes the potential assymetric with respect to the atom order. This version uses the first letter of the angle atom type definitions to decide the order of the terms. Then, it finds the appropriate bond terms for use in the bond-angle term.
            differ = θ - bencoeff(0)
            Uangle(θ) = bencoeff(1)*differ2 * (1 + bencoeff(2)*differ4)
            if there are Bond-Angle cross terms then
            Ubond-angle(θ) = bencoeff(3) * [ θijk - bencoeff(0) ] * [ ( rij - Equilibrium bond lengthij ) + ( rjk - Equilibrium bond lengthjk ) ]
        • 13: Infinite Square Well Angle
            The angle terms is a 1,3 distance based check with bounds.
            Uangle(θ) = infinity if distance1-3 < bencoeff(1)
            Uangle(θ) = bencoeff(3) if bencoeff(1) ≤ distance1-3 ≤ bencoeff(2)
            Uangle(θ) = infinity if bencoeff(2) < distance1-3
        • 14: Multiple Allowed Rigid Angles
            Uangle(θ) = Infinity, if |(θ - bencoeff(0))| > bencoeff(2)
              and |(θ - bencoeff(1))| > bencoeff(2).
            Uangle = 0 otherwise.
        • 15: MMFF Cubic with bond-angle cross terms
            differ = &theta - bencoeff(0)
            Uangle(θ) = bencoeff(1)*differ2 + bencoeff(2)*differ3 + Ubond-angle(rij,rjk,θ)
            if Bond-Angle Logical is true then
              Ubond-angle(rij,rjk,θ) = differ*[ bencoeff(3)*(rij - vibcoeff(0)ij) +bencoeff(4)*(rjk - vibcoeff(0)jk) ]
            otherwise
              Ubond-angle(rij,rjk,θ) = 0
        • 16: Harmonic Cosine plus 1-3 nonbonded interactions
            Uangle(θ) = bencoeff(1)*[ Cos(θ) - Cos(bencoeff(0)) ]2 + Unonbond(rik) + Uelectrostatic(rik)
      'Bond-Angle Logical' (a18) [logical]
      • This entry is only present for the Compass angle styles (4 and 8). It is true if there are bond-angle terms and false otherwise.
      'Bond-Angle Coefficients' (a23) [double precision]
      • This entry is only present for the Compass angle styles (4 and 8), and only if the Bond-Angle Logical is true. This lists the Bond-Angle coefficients (those appearing in Ubond-angle) one per line, in order starting from the lowest index coefficient.
      'Bond-Bond Logical' (a17) [logical]
      • This entry is only present for the Compass angle styles (4 and 8). This is true if there are bond-bond terms with this angle type.
      'Bond-Bond Coefficients' (a22) [double precision]
      • This entry is only present for the Compass angle styles (4 and 8), and only if the Bond-Bond Logical is true. This lists the Bond-Bond coefficients (those appearing in Ubond-bond) one per line, in order starting from the lowest index coefficient.
      'Angle Coefficients' (a18) [double precision]
      • The values of the bencoeffs listed in order (starting with 0) with a single coefficient per line. Note, that those bencoeffs that appear in the bond-angle or bond-bond cross terms are not listed here.
      'Angle Order' (a11) [character string a15]
      • This is used to distinguish between angle types that contain the same sets of angles, but in different environments. The environment is automatically detected according to the appropriate rules for the given forcefield. Required for all forcefields, but currently only used for MMFF94
      'Force Field Name' (a16) [character string a10]
      • The name of the Force Field
      'Number of Atoms with Same Parameters' (a38) [integer]
      • The number of triplets of atoms in this force field which use the same angle parameters.
      'Atom Names' (a10) [character string triplets a10,1x,a10,1x,a10]
      • All of the triplets of Atom Names for this Angle Type. These are listed with one triplet of atom names per line.
    End of the section repeated based on the Number of Angle Terms
    'Number of Torsion Terms' (a23) [integer]
    • The number of different regular torsion terms in the force field file. Regular torsions are defined for atoms that are connected through exactly three consecutive bonds.
    This is repeated based on the Number of Torsion Terms
      'Torsion Type Number' (a19) [integer]
      • The number coresponding to the torsion type. These numbers must run from 1 to the Number of Torsion Terms.
      'Torsion Style' (a13) [integer]
      • The number coresponding to the torsion style. Below is a list of Torsion Styles,their potential form, and the torcoeffs. The torsional potentials have quite a variety of forms (as shown below). Any energy terms need to use the standard code units of Kelvins. All torsion angles are computed in radians within Towhee, so any additive terms with the torsion angles should also be listed in radians. The torsional angle is listed as φ in the equations below. Please note that different authors define the torsional zero differently. Towhee uses the convention that a perfect cis bond has a torsional angle of 0.
        • 1: Simple Harmonic Potential
            Utorsion(φ) = torcoeff(0) * [φ - torcoeff(1)]^2
        • 2: OPLS style Cosine Series
            Utorsion(φ) = torcoeff(1) * [1 + Cos(φ)] + torcoeff(2) * [1 - Cos(2φ)] + torcoeff(3) * [1 + Cos(3φ)]
        • 3: Gromos/Charmm/Amber style Cosine Series
            Utorsion(φ) = ∑i=1,ntorloop torcoeff(3*(i-1)+1) * [ 1 + Cos(torcoeff(3*(i-1)+2)φ - torcoeff(3*(i-1)+3))]
        • 4: Gromos/Charmm/Amber style Cosine Series plus a harmonic term that Charmm traditionally calls an improper term despite the fact that the bonding pattern is that of a regular torsion.
            Utorsion(φ) = [∑i=1,ntorloop torcoeff(3*(i-1)+1) * [ 1 + Cos(torcoeff(3*(i-1)+2)φ - torcoeff(3*(i-1)+3))]]
            + torcoeff(3*ntorloop+1) * [ φ - torcoeff(3*ntorloop+2)]^2
        • 5: Compass Cosine series with cross terms and with autodetection.
            The bond and angle terms that appear in the cross terms use the appropriate vibcoeff and bencoeff parameters for the bond or angle in question. These are determined using the Compass naming conventions.
            Utorsion(φ) = torcoeff(0) * [1 - Cos(φ)] + torcoeff(1) * [1 - Cos(2φ)] + torcoeff(2) * [1 - Cos(3φ)]
            Ubond-torsion = (bond term 1-2) * [ torcoeff(3) * Cos(φ) + torcoeff(4) * Cos(2φ) + torcoeff(5) * Cos(3φ)]
            + (bond term 2-3) * [ torcoeff(6) * Cos(φ) + torcoeff(7) * Cos(2φ) + torcoeff(8) * Cos(3φ)]
            + (bond term 3-4) * [ torcoeff(9) * Cos(φ) + torcoeff(10) * Cos(2φ) + torcoeff(11) * Cos(3φ)]
            Uangle-torsion = (angle term 1-2-3) * [ torcoeff(12) * Cos(φ) + torcoeff(13) * Cos(2φ) + torcoeff(14) * Cos(3φ)]
            + (angle term 2-3-4) * [ torcoeff(15) * Cos(φ) + torcoeff(16) * Cos(2φ) + torcoeff(17) * Cos(3φ)]
            Uangle-angle-torsion = torcoeff(18) * (angle term 1-2-3) * (angle term 2-3-4)
            Ubond-bond = torcoeff(19) * (bond term 1-2) * (bond term 3-4)
        • 6: Compass Cosine series without cross terms and with autodetection.
            Utorsion(φ) = torcoeff(0) * [1 - Cos(φ)] + torcoeff(1) * [1 - Cos(2φ)] + torcoeff(2) * [1 - Cos(3φ)]
        • 7: TraPPE simple Cosine function.
            Note that unlike most force fields, TraPPE uses the convention that a trans bond is 0 φ. That is the reason for the shift of π in this potential as Towhee uses the zero cis convention.
            Utorsion(φ) = torcoeff(0) * [1 - Cos(2*[φ - &pi] + torcoeff(1)) ]
        • 8: Nonbonded Potential only.
            The nonbonded potential for van der Waals and coulombic terms is the only thing used.
        • 9: Compass cosine series with cross terms and explicit parameter declaration.
            Utorsion(φ) = torcoeff(0) * [1 - Cos(φ - torcoeff(3))] + torcoeff(1) * [1 - Cos(2*(φ - torcoeff(4)))] + torcoeff(2) * [1 - Cos(3*(φ - torcoeff(5)))]
            Ubond-torsion(φ) = [length(1-2) - torcoeff(9)] * [ torcoeff(6) * Cos(φ) + torcoeff(7) * Cos(2φ) + torcoeff(8) * Cos(3φ)]
            + [length(2-3) - torcoeff(13)] * [ torcoeff(10) * Cos(φ) + torcoeff(11) * Cos(2φ) + torcoeff(12) * Cos(3φ)]
            + [length(3-4) - torcoeff(17)] * [ torcoeff(14) * Cos(φ) + torcoeff(15) * Cos(2φ) + torcoeff(16) * Cos(3φ)]
            Uangle-torsion(φ) = [ angle(1-2-3) - torcoeff(24)] * [ torcoeff(18) * Cos(φ) + torcoeff(19) * Cos(2φ) + torcoeff(20) * Cos(3φ)]
            + [angle(2-3-4) - torcoeff(25)] * [ torcoeff(21) * Cos(φ) + torcoeff(22) * Cos(2φ) + torcoeff(23) * Cos(3φ)]
            Uangle-angle-torsion(φ) = torcoeff(26) * [angle(1-2-3) - torcoeff(27)] * [angle(2-3-4) - torcoeff(28)] * [Cos(φ - torcoeff(3))]
            Ubond-bond = torcoeff(29) * [length(1-2) - torcoeff(30)] * [length(3-4) - torcoeff(31)]
        • 10: Cosine Power Series.
            Note that there are actually ntorloop+1 parameters required as this sum runs from a power of 0 to a power of ntorloop.
            Utorsion(φ) = ∑i=0,ntorloop torcoeff(i)*Cosi(φ)
        • 11: Old style OPLS Cosine Series.
            Utorsion(φ) = torcoeff(0) + torcoeff(1)*[1 + Cos(φ)] + torcoeff(2)*[1 - Cos(2φ)] + torcoeff(3)*[1 + Cos(3φ)]
        • 12: Sum2003 Cosine Series (Sum et al. 2003)
            Note: this was not functioning properly prior to version 8.1.0.
            Utorsion(φ) = ∑i=0, ntorloop-1 torcoeff(i+1) * [ 1 - Cos{i*(φ - torcoeff(0))}]
        • 13: Old OPLS Two term Cosine Series.
            Utorsion(φ) = torcoeff(1)*[1 - Cos(φ)] + torcoeff(2)*[1 - Cos(2φ)]
        • 14: UFF 1 - Cosine divided by Total Number of Torsions
            This potential has the initial force constant divided by the total number of torsions that share the central two atoms. This value is computed on the fly based upon the number of atoms bonded to each of the two central atoms.
            totbond = [invib(atom2)-1]*[invib(atom3)-1]
            Utorsion(φ) = torcoeff(1) * [1 - torcoeff(2)*Cos(torcoeff(3)φ)]/totbond
        • 15: DREIDING 1 - Cos(n (φ - φ0)) divided by Total Number of Torsions
            This potential has the initial force constant divided by the total number of torsions that share the central two atoms. This value is computed on the fly based upon the number of atoms bonded to each of the two central atoms.
            totbond = [invib(atom2)-1]*[invib(atom3)-1]
            Utorsion(φ) = torcoeff(1) * [1 - Cos(torcoeff(2)*(φ - torcoeff(3)))]/totbond
        • 16: 2-fold Cosine
            Utorsion(φ) = torcoeff(1) * [1 - Cos( 2φ ) ]
        • 17: TraPPE planer cosine series
            Utorsion(φ) = torcoeff(1) * [1 + Cos(φ + torcoeff(3))] + torcoeff(2) * [1 - Cos2(φ)]
        • 18: Square Well
            Utorsion = torcoeff(3): if torcoeff(1) < rik < torcoeff(2)
            Utorsion = infinity : otherwise
        • 19: Amber / Number of Torsions
            tottor(jk) = ∑ of all torsions that share j and k as the central atoms
            Utorsion(φ) = ∑i=1,ntorloop torcoeff(3*(i-1)+1) * [ 1 + Cos(torcoeff(3*(i-1)+2)φ - torcoeff(3*(i-1)+3))] / tottor(jk)
        • 20: OPLS Fluorocarbon Four Term
            Utorsion(φ) = torcoeff(1)*(1 + Cos(φ)) + torcoeff(2)*(1 - Cos(2φ)) + torcoeff(3)*(1 + Cos(3φ)) + torcoeff(4)*(1 - Cos(4φ))
        • 21: Multiple Rigid Dihedrals
            Utorsion(φ) = Infinity if for each value of x as x goes from 1 to ntorloop
              | φ - torcoeff(x) | > torcoeff(0)
            Utorsion(φ) = 0 otherwise.
        • 22: Fluoroalkane series from Cui et al. 1998
            Utorsion(φ) = torcoeff(0) + torcoeff(1)*(1 - Cos(φ)) + torcoeff(2)*(1 - Cos(3 φ)) + torcoeff(3)*(1 - Cos(φ))5 + torcoeff(4)*exp[-torcoeff(5) φ2]
      'One-Four Nonbond Logical' (a24) [logical]
      • True if you wish to add nonbonded van der Waals and coulombic terms into the torsional potential.
      'One-Four Coulombic Scaling' (a26) [double precision]
      • This entry is only present if the One-Four Nonbond Logical is true. This is the scaling prefactor for the coulombic one-four terms in the torsional potential. A value of 1.0 uses the same coulombic energy as for all of the other nonbonded terms.
      'Number of Torsion Loops' (a23) [integer]
      • This entry is only present for those potentials that have the ntorloop variable (see above). In those cases, this is the number of terms in the Cosine series.
      'Torsion Coefficients' (a20) [double precision]
      • The torcoeff variables, listed one per line, starting with the lowest index variable and working upwards.
      'Torsion Order' (a13) [character string a15]
      • This is used to distinguish between torsion types that contain the same sets of atoms, but in different environments. The environment is automatically detected according to the appropriate rules for the given forcefield. Required for all forcefields, but currently only used for MMFF94
      'Force Field Name' (a16) [character string a10]
      • The name of the Force Field
      'Number of Atoms with Same Parameters' (a38) [integer]
      • The number of quartets of atoms in this force field which use the same torsion parameters.
      'Atom Names' (a10) [character string quartets a10,1x,a10,1x,a10,1x,a10]
      • All of the quartets of Atom Names for this Torsion Type. These are listed with one quartet of atom names per line.
    End of the section repeated based on the Number of Torsion Terms
    'Number of Improper Terms' (a24) [integer]
    • The number of different improper torsion terms in the force field file. Improper torsions are defined for three atoms that are all bonded to a single central atom. They are most often used to enforce planarity of sp2 centers.
    This section is repeated based on the Number of Improper Terms
      'Improper Type Number' (a20) [integer]
      • The number coresponding to the improper type. These numbers must run from 1 to the Number of Improper Terms.
      'Improper Form' (a13) [integer]
      • The number coresponding to the improper form.
        • 1: Amber / Stereocenter Form.
        • 2: Charmm / Simple out-of-plane form.
        • 3: Average over the three possible Wilson out-of-plane angles form.
        • 4: UFF angle between the IL vector and the IJK plane where I is the central atom bonded to atoms J, K, and L.
        • 5: Vector of the three possible Wilson out-of-plane angles
      'Improper Style' (a13) [integer]
      • The number coresponding to the improper style. Below is a list of Improper Styles and their impcoeffs. The improper potentials have a variety of forms (as shown below). Any energy terms need to use the standard code units of Kelvins. All improper angles are computed in radians within Towhee, so any additive terms with the improper angles should also be listed in radians. The improper angle is listed as φ in the equations below. Please note that different authors define the torsional zero differently.
        • 1: Simple Harmonic Potential
            Uimproper(φ) = impcoeff(0) * [φ - impcoeff(1)]^2
        • 2: OPLS style cosine series
            Uimproper(φ) = impcoeff(1) * [ 1 + Cos(φ)] + impcoeff(2) * [1 - Cos(2φ)] + impcoeff(3) * [ 1 + Cos(3φ)]
        • 3: Potentials for enforcing stereochemistry
            This is to be used only with the stereochemistry improper form. It is used to enforce R/S or D/L stereochemistry centers. The same potential is used in both cases, but the atom ordering is different.
            Uimproper(φ) = 0, if 0 < φ ;
            = impcoeff(0) + impcoeff(1)*[-φ], if -0.5 &pi < φ < 0
            = impcoeff(0) + impcoeff(1)*[&pi + φ], if φ < -0.5 &pi
        • 4: Amber cosine
            Uimproper(φ) = impcoeff(1) * [ 1 + Cos(impcoeff(2)φ - impcoeff(3))]
        • 5: Wilson Harmonic
            Uimproper(φ) = impcoeff(0) * [ φ - impcoeff(1) ]2
        • 6: UFF Cosine Series
            Uimproper(φ) = impcoeff(3)*[impcoeff(0) + imcoeff(1)*Cos(φ) + impcoeff(2)*Cos(2φ)]
        • 7: 1 - Cos(φ)
            Uimproper(φ) = impcoeff(1)*[1 - Cos(φ)]
        • 8: MMFF Wilson planar enforcer (only works with Improper Form 5)
            Uimproper(φ) = impcoeff(1)*[ φj2 + φk2 + φl2]
      'Improper Coefficients' (a21) [double precision]
      • The impcoeff variables, listed one per line, starting with the lowest index variable and working upwards.
      'Force Field Name' (a16) [character string a10]
      • The name of the Force Field
      'Number of Atoms with Same Parameters' (a38) [integer]
      • The number of quartets of atoms in this force field which use the same improper parameters
      'Atom Names' (a10) [character string quartets a10,1x,a10,1x,a10,1x,a10]
      • All of the quartets of Atom Names for this Improper Type. These are listed with one quartet of atom names per line. The naming conventions for the different force fields are not consistent. If you wish to use this feature you should look in the code for examples.
    End of the section repeated based on the Number of Improper Terms
    'Number of Angle-Angle Terms' (a27) [integer]
    • The number of different angle-angle terms in the force field file. Angle-Angle terms are defined for 3 atoms that are all bonded to a single central atom. The two angles are inaa2-current-inaa3 and inaa2-current-inaa4.
    This section is repeated based on the Number of Angle-Angle Terms
      'Angle-Angle Type Number' (a23) [integer]
      • The number coresponding to the angle-angle type. These numbers must run from 1 to the Number of Angle-Angle Terms.
      'Angle-Angle Style' (a17) [integer]
      • The number coresponding to the angle-angle style. Below is a list of Angle-Angle Styles, their potential form, and the aacoeffs. The angle-angle potentials have a variety of forms (as shown below). Any energy terms need to use the standard code units of Kelvins. All angles are computed in radians within Towhee, so any additive terms with the angle-angles should also be listed in radians.
        • 1: Compass Angle-Angle with autodetect
            Uangle-angle = aacoeff(0) * [ angle(inaa0-current-inaa1) - bencoeff(inaa0-current-inaa1) ] * [ angle(inaa0-current-inaa2) - bencoeff(inaa0-current-inaa2) ]
        • 2: Compass Angle-Angle with explicit parameters
            Uangle-angle = aacoeff(0) * [ angle(inaa0-current-inaa1) - aacoeff(1) ] * [ angle(inaa0-current-inaa2) - aacoeff(2) ]
      'Angle-Angle Coefficients' (a24) [double precision]
      • The aacoeff variables, listed one per line, starting with the lowest index variable and working upwards.
      'Force Field Name' (a16) [character string a10]
      • The name of the Force Field
      'Number of Atoms with Same Parameters' (a38) [integer]
      • The number of quartets of atoms in this force field which use the same improper parameters
      'Atom Names' (a10) [character string quartets a10,1x,a10,1x,a10,1x,a10]
      • All of the quartets of Atom Names for this Angle-Angle Type. These are listed with one quartet of atom names per line. The Compass naming convention is followed to add ease of transfer from their publications. The order for atom output in this convention is inaa1, current, inaa0, inaa2.
    End of the section repeated based on the Number of Angle-Angle Terms
    'Number of One-Five Types' (a23) [integer]
    • The number of different special One-Five interaction types the force field file.
    This section is repeated once for each of the Number of One-Five Types
      'One-Five Type Number' (a20) [integer]
      • The number coresponding to the special one-five interaction type. These numbers must run from 1 to the Number of One-Five Types.
      'One-Five Style' (a14) [integer]
      • The number coresponding to the special one-five interaction style. Below is a list of One-Five Styles, their potential form, and the ofcoeffs. Currently there is only one kind of one-five potential. Any energy terms need to use the standard code units of Kelvins.
        • 1: Lennard-Jones 12-6 interactions with special parameters
            • Unonbond = 4 * ofcoeff(2) * [ (ofcoeff(1)/r)12 - (ofcoeff(1)/r)6 - Shift]
            Where the parameters are as follows.
            • ofcoeff(1): sigma in units of Å.
            • ofcoeff(2): epsilon in units of Kelvin.
            • Shift: computed automatically by the code if we are using a shifted potential.
        • 2: 12th power repulsion
            • Unonbond = ofcoeff(1) / r12
            Where the parameter is as follows.
            • ofcoeff(1): energy term in units of Kelvin * Å12.
      'One-Five Coefficients' (a22) [double precision]
      • The ofcoeff variables, listed one per line, starting with the lowest index variable and working upwards.
      'Force Field Name' (a16) [character string a10]
      • The name of the Force Field
      'Atom Names' (a10) [character string quintet a10,1x,a10,1x,a10,1x,a10,1x,a10]
      • The five atom names that, when bonded successively, create a special one-five interaction. These names are compared with the nbnames to find a valid match. If no matches are found then it is assumed that no special one-five interaction are needed for these atoms and the normal method of computing nonbonded pair interactions is used instead.
    End of the section repeated based on the Number of One-Five Types
    'Number of Bond Increments' (a25) [integer]
    • The number of different Bond Increment terms in the force field file.
    This section is repeated once for each of the Number of Bond Increments
      'Bond Increment Type Number' (a26) [integer]
      • The number coresponding to the bond increment type. These numbers must run from 1 to the Number of Bond Increments.
      'Bond Increment Value' (a20) [double precision]
      • The value of the bond increment. This value is added to the charge of the first atom listed in the Atom Names and subtracted from the charge of the second atom listed in the Atom Names.
      'Bond Increment Order' (a20) [character string a10]
      • The bond increment order. This is used to distinguish between different types of bonds that can be made between the same pairs of atoms in the same manner as the 'Vibration Order'. Commonly used options are listed in the advanced connectivity map manual.
      'Force Field Name' (a16) [character string a10]
      • The name of the Force Field
      'Atom Names' (a10) [character string pair a10,1x,a10]
      • The pair of atoms that share this bond increment. The bond increment is added to the base charge of the first atom listed and subtracted from the base charge of the second atom listed when these atoms are connected by a bond.
    End of the section repeated based on the Number of Bond Increments
Return to the main towhee web page
Last updated: July 22, 2021 Send comments to: Marcus Martin