
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 7.1.0.
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 BondAngle 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 OneFive 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 OneFive 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 'LorentzBerthelot' 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 3part': The pair potential is a complicated 3 part function as described in
Mendelev et al. 2003.
If r < nbcoeff(1)
 x = r/nbcoeff(4)
 u_{pair} = ( 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) * r^{2} + nbcoeff(8) * r^{3}
 u_{pair} = nbcoeff(9) * exp[x]
 else
 u_{pair} = ∑_{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
 u_{pair} = ( 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) * r^{2} + nbcoeff(15) * r^{3}
 u_{pair} = nbcoeff(16) * exp[x]
 else
 u_{pair} = ∑_{i = 1,table_npair} table_pair(2,i) *
[table_pair(1,i)  r]^{i+3} * Heaviside[table_pair(1,i)r]
 'Belashchenko 58': Two part pair potential as described in Belashchenko 2013.
 if r < r_{1} (where r_{1} is nbcoeff(1)) then
 u_{pair} = nbcoeff(7) + nbcoeff(8) * (r_{1}  r) + nbcoeff(9)*(exp[nbcoeff(10)*(r_{1}r)]1)
 else
 u_{pair} = ∑_{i = 1,5} ∑_{m=0,8} a_{i,m} (r  r_{i+1})^{m} Heaviside[r_{i},r_{i+1}]
 where the r_{i} are stored in nbcoeff(i) for i=1,6
 and a_{i,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.
 u_{pair} = nbcoeff(1) e^{[nbcoeff(2) rij]}.
 'morse': The pair potential is described by the Morse function.
 u_{pair} = nbcoeff(1) * { exp[2 nbcoeff(2) (r  nbcoeff(3))]  2 exp[nbcoeff(2) (rnbcoeff(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 3part'
 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 nontabular 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 selfinteraction 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'
 ρ(r_{ij}) = ∑_{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.
 ρ(r_{ij}) = eam_dens(1,2) e^{eam_dens(2,2) rij} for r_{ij} < 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': u_{embed}(ρ(r_{ij}))

= eam_embed_data(4,1) * √ ρ (r_{ij})
+ eam_embed_data(4,2) * ρ(r_{ij})
for 0 ≤ ρ(r_{ij}) < eam_embed_data(5,1)

= eam_embed_data(1,2) + eam_embed_data(2,2) * [ρ(r_{ij})  eam_embed_data(5,2)]^{2} + eam_embed_data(3,2)*[ρ(r_{ij})  eam_embed_data(5,2)]^{3}
for eam_embed_data(5,1) ≤ ρ(r_{ij})
 '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)
u_{embed}(ρ(r_{ij}))

= [eam_embed_data(8,1) + eam_embed_data(8,2) * (ρ(r_{ij})  belrho(7)) + eam_embed_data(8,3) * (ρ(r_{ij})  belrho(7))^{2}]
* [2*ρ(r_{ij})/belrho(7)  {ρ(r_{ij})/belrho(7)}^{2}]
for 0 < ρ(r_{ij}) < belrho(7)

= eam_embed_data(i,1) + (ρ(r_{ij})  belrho(i1)) * [eam_embed_data(i,2) + eam_embed_data(i,3)*(ρ(r_{ij})belrho(i1))]
for belrho(i) ≤ ρ(r_{ij}) < belrho(i1) where 2 ≤ i ≤ 7

= eam_embed_data(1,1) + eam_embed_data(3,1) * (ρ(r_{ij})  belrho(0))^{2}
for belrho(1) ≤ ρ(r_{ij}) < belrho(8).

= eam_embed_data(9,1) + eam_embed_data(9,2) * [ρ(r_{ij})  belrho(8)] + eam_embed_data(9,3)*[ρ(r_{ij})  belrho(8)]^{eam_embed_data(14,2)}
for belrho(8) ≤ ρ(r_{ij}) < belrho(9)

= eam_embed_data(10,1) + eam_embed_data(10,2) * [ρ(r_{ij})  belrho(9)] + eam_embed_data(10,3)*[ρ(r_{ij})  belrho(9)]^{eam_embed_data(14,3)}
for belrho(9) ≤ ρ(r_{ij})
 'logarithmic': u_{embed}(ρ(r_{ij})) = ρ(r_{ij}) * (eam_embed_data(1,2) * ln[ρ(r_{ij})] + eam_embed_data(2,2))
 'power 0.5 and 2': u_{embed}(ρ(r_{ij})) =
eam_embed_data(1,2) √ ρ (r_{ij})
+ eam_embed_data(2,2) *[ρ(r_{ij})]^{2}
 'power 0.5, 2, and 4': u_{embed}(ρ(r_{ij})) =
eam_embed_data(1,2) √ ρ (r_{ij})
+ eam_embed_data(2,2) *[ρ(r_{ij})]^{2}
+ eam_embed_data(3,2) *[ρ(r_{ij})]^{4}
 'squareroot': u_{ebmed}(ρ(r_{ij})) = eam_embed_data(2,2) * ρ(r_{ij})
 eam_embed_data(1,2) √ ρ (r_{ij})
 '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/k_{b}).
 'Ackland 2term': 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 configurationalbias moves that utilize a nonuniform generation of
the bond angles and dihedrals. See the configurationalbias 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 configurationalbias 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
 U_{bond} = Infinity, if ( length  vibcoeff(0) ) / vibcoeff(0) > 0.01
 U_{bond} = 0 otherwise.
 2: Standard Harmonic
 U_{bond} = vibcoeff(1)*[ length  vibcoeff(0) ]^{2}
 3: Gromos Quartic
 U_{bond} = vibcoeff(1)*[ length^{2}  vibcoeff(0)^{2} ]^{2}
 4: Nonlinear
 distance = [ length  vibcoeff(0) ]^{2}
 U_{bond} = [ vibcoeff(1)*distance ] / [ vibcoeff(2)  distance ]
 5: MM2 Quadradic and Triplet
 distance = [ length  vibcoeff(0) ]
 U_{bond} = vibcoeff(1)*distance^{2} * [ 1.0  2.0*distance ], if distance ≤ 1/3
 U_{bond} = Infinity, if distance > 1/3
 6: Compass Quartic
 distance = [ length  vibcoeff(0) ]
 U_{bond} = vibcoeff(1)*distance^{2} + vibcoeff(2)*distance^{3} + vibcoeff(3)*distance^{4}
 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) ]
 U_{bond} = vibcoeff(1)*[ e^{vibcoeff(2)*distance}  1 ]^{2}
 10: Infinite Square Well
 vibcoeff(0) is no longer listed for this potential starting with version 7.1.0
 U_{bond} = infinity if length < vibcoeff(1)
 U_{bond} = vibcoeff(3) if vibcoeff(1) <= length <= vibcoeff(2)
 U_{bond} = infinity if vibcoeff(2) < length
 11: Standard Harmonic plus nonbonded interactions
 U_{bond} = vibcoeff(1)*[ length  vibcoeff(0) ]^{2}
+ U_{nonbond}(r_{ij}) + U_{electrostatic}(r_{ij})
 12: FENE plus nonbonded interactions (Reference: Kremer and Grest, 1990
 U_{bond} =  0.5 * vibcoeff(2) * vibcoeff(1)^{2} * ln[ 1.0  (length/(vibcoeff(1)))^{2} ]
+ U_{nonbond}(r_{ij}) + U_{electrostatic}(r_{ij}) if length < vibcoeff(1).
 U_{bond} = 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, sp^{2}
carbons bond to each other with either single or double bonds). Commonly used options are listed in the
inpstyle 6 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
 U_{angle}(θ) = Infinity, if (θ  bencoeff(0)) > bencoeff(1)
 U_{angle}(θ) = 0 otherwise.
 1: Standard Harmonic
 U_{angle}(θ) = bencoeff(1)*[ θ  bencoeff(0) ]^{2}
 2: DREIDING 1 + Cos(θ)
 U_{angle}(θ) = bencoeff(1)*[ 1.0 + Cos(θ) ]
 3: Harmonic Cosine
 U_{angle}(θ) = bencoeff(1)*[ Cos(θ)  Cos(bencoeff(0)) ]^{2}
 4: Compass Quartic Angle with autodetection
The BondAngle and BondBond (12) 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 bondangle and bondbond cross terms. Two logicals control whether there are any BondAngle or
BondBond cross terms.
Definitions: r_{ij} is the distance between the first and second atoms in the angle, r_{jk} is the distance between the second
and third atoms of the angle, and the appropriate vibcoeff(0) values are determined automatically.
 differ = θ  bencoeff(0)
 U_{angle}(θ) = bencoeff(1)*differ^{2} + bencoeff(2)*differ^{3} + bencoeff(3)*differ^{4}
+ U_{bondangle}(r_{ij},r_{jk},θ)
+ U_{bondbond}(r_{ij},r_{jk})
 if BondAngle Logical is true then
U_{bondangle}(θ,r_{ij},r_{jk}) = differ*[ bencoeff(4)*(r_{ij}  vibcoeff(0)_{ij})
+bencoeff(5)*(r_{jk}  vibcoeff(1)_{jk}) ]
otherwise
U_{bondangle}(r_{ij},r_{jk},θ) = 0
 if BondBond Logical is true then
U_{bondbond}(r_{ij},r_{jk}) = bencoeff(6)*[
(r_{ij}  vibcoeff(0)_{ij}) *(r_{kj}  vibcoeff(0)_{jk}) ]
otherwise
U_{bondbond}(r_{ij},r_{jk}) = 0
 5: Charmm Harmonic with UreyBradley
U_{angle}(θ,d_{ik}) = bencoeff(1)*[ θ  bencoeff(0) ]^{2}
+ bencoeff(3)*[d_{ik}  bencoeff(2) ]^{2}
 where d_{ik} is the distance between the first and third atoms of the angle.
 6: Nonbonded Terms only between the 13 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 13 atoms plus a harmonic angle
The angle energies are the nonbonded van der Waals, the nonbonded coulomb terms, and the following.
 U_{angle}(θ) = bencoeff(1)*[ θ  bencoeff(0)]^{2}
 8: Compass Quartic Angle with explicit ordering of terms
The BondAngle and BondBond (12) 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 BondAngle or BondBond cross terms.
 differ = anglebencoeff(0)
 U_{angle}(θ) = bencoeff(1)*differ^{2} + bencoeff(2)*differ^{3} + bencoeff(3)*differ^{4}
 if there are BondAngle cross terms then
 U_{bondangle}(θ,d_{ij},d_{jk}) = [ θ  bencoeff(0) ] * {bencoeff(4)*[ d_{ij}  bencoeff(5) ]
+ bencoeff(6)*[ d_{jk}  bencoeff(7) ]}
 U_{bondbond}(d_{ij},d_{jk}) = bencoeff(8) * [ d_{ij}  bencoeff(9) ]
* [ d_{jk}  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 126 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)
 r_{ik} = [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)/(r_{ik}^{5})
* { vibcoeff(ij,0)*vibcoeff(jk,0)*[1  Cos^{2}(θ)]  vibcoeff(ik,0)^{2}*Cos(θ) } / bencoeff(1)^{2}
 U_{angle}(θ) = 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 126 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)
 r_{ik} = [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)/(r_{ik}^{5})
* { vibcoeff(ij,0)*vibcoeff(jk,0)*[1  Cos^{2}(θ_{0})]  vibcoeff(ik,0)^{2}*Cos(θ_{0}) }
 U_{angle}(θ) = 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 AngleBond Terms
The BondAngle 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 bondangle term.
 differ = θ  bencoeff(0)
 U_{angle}(θ) = bencoeff(1)*differ^{2} * (1 + bencoeff(2)*differ^{4})
 if there are BondAngle cross terms then
 U_{bondangle}(θ) = bencoeff(3) * [ θ_{ijk}  bencoeff(0) ]
* [ ( r_{ij}  Equilibrium bond length_{ij} ) + ( r_{jk}  Equilibrium bond length_{jk} ) ]
 13: Infinite Square Well Angle
The angle terms is a 1,3 distance based check with bounds.
 U_{angle}(θ) = infinity if distance_{13} < bencoeff(1)
 U_{angle}(θ) = bencoeff(3) if bencoeff(1) <= distance_{13} <= bencoeff(2)
 U_{angle}(θ) = infinity if bencoeff(2) < distance_{13}
 14: Multiple Allowed Rigid Angles
U_{angle}(θ) = Infinity, if (θ  bencoeff(0)) > bencoeff(2)
and (θ  bencoeff(1)) > bencoeff(2).
 U_{angle} = 0 otherwise.
 15: MMFF Cubic with bondangle cross terms
differ = &theta  bencoeff(0)
 U_{angle}(θ) = bencoeff(1)*differ^{2} + bencoeff(2)*differ^{3}
+ U_{bondangle}(r_{ij},r_{jk},θ)
 if BondAngle Logical is true then
U_{bondangle}(r_{ij},r_{jk},θ) = differ*[ bencoeff(3)*(r_{ij}  vibcoeff(0)_{ij})
+bencoeff(4)*(r_{jk}  vibcoeff(0)_{jk}) ]
otherwise
U_{bondangle}(r_{ij},r_{jk},θ) = 0
 16: Harmonic Cosine plus 13 nonbonded interactions
 U_{angle}(θ) = bencoeff(1)*[ Cos(θ)  Cos(bencoeff(0)) ]^{2}
+ U_{nonbond}(r_{ik}) + U_{electrostatic}(r_{ik})
 'BondAngle Logical' (a18) [logical]
 This entry is only present for the Compass angle styles (4 and 8). It is true if there are bondangle terms and false otherwise.
 'BondAngle Coefficients' (a23) [double precision]
 This entry is only present for the Compass angle styles (4 and 8), and only if the BondAngle Logical is true. This lists the BondAngle
coefficients (those appearing in U_{bondangle}) one per line, in order starting from the lowest index coefficient.
 'BondBond Logical' (a17) [logical]
 This entry is only present for the Compass angle styles (4 and 8). This is true if there are bondbond terms with this angle type.
 'BondBond Coefficients' (a22) [double precision]
 This entry is only present for the Compass angle styles (4 and 8), and only if the BondBond Logical is true. This lists the BondBond
coefficients (those appearing in U_{bondbond}) 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 bondangle or bondbond 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
 U_{torsion}(φ) = torcoeff(0) * [φ  torcoeff(1)]^2
 2: OPLS style Cosine Series
 U_{torsion}(φ) = torcoeff(1) * [1 + Cos(φ)] + torcoeff(2) * [1  Cos(2φ)] + torcoeff(3) * [1 + Cos(3φ)]
 3: Gromos/Charmm/Amber style Cosine Series
 U_{torsion}(φ) = ∑^{i=1,ntorloop} torcoeff(3*(i1)+1) * [ 1 + Cos(torcoeff(3*(i1)+2)φ  torcoeff(3*(i1)+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.
 U_{torsion}(φ) = [∑^{i=1,ntorloop} torcoeff(3*(i1)+1)
* [ 1 + Cos(torcoeff(3*(i1)+2)φ  torcoeff(3*(i1)+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.
 U_{torsion}(φ) = torcoeff(0) * [1  Cos(φ)] + torcoeff(1) * [1  Cos(2φ)] + torcoeff(2) * [1  Cos(3φ)]
 U_{bondtorsion} = (bond term 12) * [ torcoeff(3) * Cos(φ) + torcoeff(4) * Cos(2φ) + torcoeff(5) * Cos(3φ)]
 + (bond term 23) * [ torcoeff(6) * Cos(φ) + torcoeff(7) * Cos(2φ) + torcoeff(8) * Cos(3φ)]
 + (bond term 34) * [ torcoeff(9) * Cos(φ) + torcoeff(10) * Cos(2φ) + torcoeff(11) * Cos(3φ)]
 U_{angletorsion} = (angle term 123) * [ torcoeff(12) * Cos(φ) + torcoeff(13) * Cos(2φ) + torcoeff(14) * Cos(3φ)]
 + (angle term 234) * [ torcoeff(15) * Cos(φ) + torcoeff(16) * Cos(2φ) + torcoeff(17) * Cos(3φ)]
 U_{angleangletorsion} = torcoeff(18) * (angle term 123) * (angle term 234)
 U_{bondbond} = torcoeff(19) * (bond term 12) * (bond term 34)
 6: Compass Cosine series without cross terms and with autodetection.
 U_{torsion}(φ) = 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.
 U_{torsion}(φ) = 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.
 U_{torsion}(φ) = torcoeff(0) * [1  Cos(φ  torcoeff(3))] + torcoeff(1) * [1  Cos(2*(φ  torcoeff(4)))]
+ torcoeff(2) * [1  Cos(3*(φ  torcoeff(5)))]
 U_{bondtorsion}(φ) = [length(12)  torcoeff(9)] * [ torcoeff(6) * Cos(φ) + torcoeff(7) * Cos(2φ)
+ torcoeff(8) * Cos(3φ)]
 + [length(23)  torcoeff(13)] * [ torcoeff(10) * Cos(φ) + torcoeff(11) * Cos(2φ) + torcoeff(12) * Cos(3φ)]
 + [length(34)  torcoeff(17)] * [ torcoeff(14) * Cos(φ) + torcoeff(15) * Cos(2φ) + torcoeff(16) * Cos(3φ)]
 U_{angletorsion}(φ) = [ angle(123)  torcoeff(24)] * [ torcoeff(18) * Cos(φ)
+ torcoeff(19) * Cos(2φ) + torcoeff(20) * Cos(3φ)]
 + [angle(234)  torcoeff(25)] * [ torcoeff(21) * Cos(φ) + torcoeff(22) * Cos(2φ) + torcoeff(23) * Cos(3φ)]
 U_{angleangletorsion}(φ) = torcoeff(26) * [angle(123)  torcoeff(27)] * [angle(234)  torcoeff(28)]
* [Cos(φ  torcoeff(3))]
 U_{bondbond} = torcoeff(29) * [length(12)  torcoeff(30)] * [length(34)  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.
 U_{torsion}(φ) = ∑^{i=0,ntorloop} torcoeff(i)*Cos^{i}(φ)
 11: Old style OPLS Cosine Series.
 U_{torsion}(φ) = torcoeff(0) + torcoeff(1)*[1 + Cos(φ)] + torcoeff(2)*[1  Cos(2φ)] + torcoeff(3)*[1 + Cos(3φ)]
 12: Sun2003 Cosine Series.
 U_{torsion}(φ) = ∑^{i=0,ntorloop} torcoeff(i+1) * [ 1  Cos{i*(φ  torcoeff(0))}]
 13: Old OPLS Two term Cosine Series.
 U_{torsion}(φ) = 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(atom_{2})1]*[invib(atom_{3})1]
 U_{torsion}(φ) = 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(atom_{2})1]*[invib(atom_{3})1]
 U_{torsion}(φ) = torcoeff(1) * [1  Cos(torcoeff(2)*(φ  torcoeff(3)))]/totbond
 16: 2fold Cosine
 U_{torsion}(φ) = torcoeff(1) * [1  Cos( 2φ ) ]
 17: TraPPE planer cosine series
 U_{torsion}(φ) = torcoeff(1) * [1 + Cos(φ + torcoeff(3))] + torcoeff(2) * [1  Cos^{2}(φ)]
 18: Square Well
 U_{torsion} = torcoeff(3): if torcoeff(1) < r_{ik} < torcoeff(2)
 U_{torsion} = infinity : otherwise
 19: Amber / Number of Torsions
 tottor(jk) = ∑ of all torsions that share j and k as the central atoms
 U_{torsion}(φ) = ∑^{i=1,ntorloop} torcoeff(3*(i1)+1) *
[ 1 + Cos(torcoeff(3*(i1)+2)φ  torcoeff(3*(i1)+3))] / tottor(jk)
 20: OPLS Fluorocarbon Four Term
 U_{torsion}(φ) = torcoeff(1)*(1 + Cos(φ)) + torcoeff(2)*(1  Cos(2φ)) + torcoeff(3)*(1 + Cos(3φ)) +
torcoeff(4)*(1  Cos(4φ))
 21: Multiple Rigid Dihedrals
 U_{torsion}(φ) = Infinity if for each value of x as x goes from 1 to ntorloop
 φ  torcoeff(x)  > torcoeff(0)
 U_{torsion}(φ) = 0 otherwise.
 22: Fluoroalkane series from Cui et al. 1998
 U_{torsion}(φ) = torcoeff(0) + torcoeff(1)*(1  Cos(φ)) + torcoeff(2)*(1  Cos(3 φ))
+ torcoeff(3)*(1  Cos(φ))^{5} + torcoeff(4)*exp[torcoeff(5) φ^{2}]
 'OneFour Nonbond Logical' (a24) [logical]
 True if you wish to add nonbonded van der Waals and coulombic terms into the torsional potential.
 'OneFour Coulombic Scaling' (a26) [double precision]
 This entry is only present if the OneFour Nonbond Logical is true. This is the scaling prefactor for the coulombic onefour 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 sp^{2} 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 outofplane form.
 3: Average over the three possible Wilson outofplane 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 outofplane 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. Towhee uses the
convention that a perfect cis bond has a torsional angle of 0.
 1: Simple Harmonic Potential
 U_{improper}(φ) = impcoeff(0) * [φ  impcoeff(1)]^2
 2: OPLS style cosine series
 U_{improper}(φ) = 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.
 U_{improper}(φ) = 0, if 0 < φ ;
 = impcoeff(0) + impcoeff(1)*[φ], if 0.5 &pi < φ < 0
 = impcoeff(0) + impcoeff(1)*[&pi + φ], if φ < 0.5 &pi
 4: Amber cosine
 U_{improper}(φ) = impcoeff(1) * [ 1 + Cos(impcoeff(2)φ  impcoeff(3))]
 5: Wilson Harmonic
 U_{improper}(φ) = impcoeff(1) * [ 1 + Cos(impcoeff(2)φ  impcoeff(3))]
 6: UFF Cosine Series
 U_{improper}(φ) = impcoeff(3)*[impcoeff(0) + imcoeff(1)*Cos(φ) + impcoeff(2)*Cos(2φ)]
 7: 1  Cos(φ)
 U_{improper}(φ) = impcoeff(1)*[1  Cos(φ)]
 8: MMFF Wilson planar enforcer (only works with Improper Form 5)
 U_{improper}(φ) = impcoeff(1)*[ φ_{j}^{2} + φ_{k}^{2} + φ_{l}^{2}]
 '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 AngleAngle Terms' (a27) [integer]
 The number of different angleangle terms in the force field file. AngleAngle terms are defined for 3 atoms that are all bonded to a
single central atom. The two angles are inaa2currentinaa3 and inaa2currentinaa4.
 This section is repeated based on the Number of AngleAngle Terms
 'AngleAngle Type Number' (a23) [integer]
 The number coresponding to the angleangle type. These numbers must run from 1 to the Number of AngleAngle Terms.
 'AngleAngle Style' (a17) [integer]
 The number coresponding to the angleangle style. Below is a list of AngleAngle Styles, their potential form, and the aacoeffs.
The angleangle 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 angleangles should also be listed in radians.
 1: Compass AngleAngle with autodetect
 U_{angleangle} = aacoeff(0) * [ angle(inaa0currentinaa1)  bencoeff(inaa0currentinaa1) ]
* [ angle(inaa0currentinaa2)  bencoeff(inaa0currentinaa2) ]
 2: Compass AngleAngle with explicit parameters
 U_{angleangle} = aacoeff(0) * [ angle(inaa0currentinaa1)  aacoeff(1) ] * [ angle(inaa0currentinaa2)  aacoeff(2) ]
 'AngleAngle 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 AngleAngle 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 AngleAngle Terms
 'Number of OneFive Types' (a23) [integer]
 The number of different special OneFive interaction types the force field file.
 This section is repeated once for each of the Number of OneFive Types
 'OneFive Type Number' (a20) [integer]
 The number coresponding to the special onefive interaction type. These numbers must run from 1 to the Number of OneFive Types.
 'OneFive Style' (a14) [integer]
 The number coresponding to the special onefive interaction style. Below is a list of OneFive Styles, their potential form, and the
ofcoeffs. Currently there is only one kind of onefive potential. Any energy terms need to use the standard code units of Kelvins.
 1: LennardJones 126 interactions with special parameters
 U_{nonbond} = 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
 U_{nonbond} = ofcoeff(1) / r^{12}
 Where the parameter is as follows.
 ofcoeff(1): energy term in units of Kelvin * Å^{12}.
 'OneFive 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 onefive 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 onefive 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 OneFive 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 inpstyle 6 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
