<html>
<head>
<title>MCCCS Towhee (towhee_input Version 3.6.x)</title>
</head>
<body bgcolor="#FFFFFF" text="#000000">
<table width="675" border="0" cellspacing="0" cellpadding="0" height="590">
<tr>
<td colspan="2" height="29">
<div align="center"> <font size="5"> <b><font face="Arial, Helvetica, sans-serif"><a name="top"></a>MCCCS
Towhee (towhee_input Version 3.6.x)</font></b> </font> </div>
</td>
</tr>
<A href="http://sourceforge.net">
<IMG src="http://sourceforge.net/sflogo.php?group_id=87039&type=5" width="210" height="62" border="0"
alt="SourceForge.net Logo" align="right"/>
</A>
<tr>
<td width="18" height="371" valign="top">
<p> </p>
<p> </p>
</td>
<td width="697" valign="top"> <b>Overview</b>
<ul>
This section covers the variables that are set in the towhee_input file
Version 3.6.x. Each variable is listed along with its type (logical, character,
integer, or double precision). towhee_input is the main input file for
Towhee and is generally the only file that needs to be edited on a regular
basis. It has a regimented style to the input. The variables are described
here in the order they appear in this file. Please look at one of the
example files (available with the code package) for the precise file format.
</ul>
<b>Bug reports / feature enhancements for 3.6.x versions</b>
<ul>
<li>3.6.9: Moved all function definitions into a new include file called functions.inc</li>
<li>3.6.8: Fixed some bugs in the readquest functionallity to work better for fitting multiple
atom types.</li>
<li>3.6.7: Modified some of the example input file and readquest to make the code compatable with
a DEC machine. Need to have single quotes around any unformatted text strings.</li>
<li>3.6.6: Increased the precision of the writelammps.F output. Modified the eetokcal constant
to give agreement with the current value in LAMMPS (which was modified to agree with Charmm).
<li>3.6.5: Made some changes to readquest.F to accomodate the new quest_input format.</li>
<li>3.6.4: Added some safety checking to the parallel functionallity. Fixed a parallel bug in
readquest.F. Added two new types of "aminio acids" to the Charmm22 force field. These are the
commonly seen N-methylamide and acetyl capping groups.</li>
<li>3.6.3: Added new parallel functionallity to the readquest routine that allows the force field
parameter search to be automatically divided up on a large number of processors.</li>
<li>3.6.2: Fixed an initialization bug in readquest that led to some strange errors when running the
SW_Silicon test case in parallel with some of the coulombic test cases.</li>
<li>3.6.1: Cleaned up the minimizer routines. Fixed a small bug in the readquest routine that resulted in
an inaccurate atom count.</li>
<li>3.6.0: Added a minimizer feature (coded by Aidan Thompson) to the Towhee code. Replaced the <b>lcycles</b>
variable with the <b>stepstyle</b> variable. If <b>stepstyle</b> is set to 'minimize' then there are also
some new variables to control the minimization.</li>
</ul>
<b>towhee_input file differences from version 3.5.x</b>
<ul>
<li> Replaced the <b>lcycles</b> variable with <b>stepstyle</b>. Added a minimizer.</li>
</ul>
<b>Variable explanations for towhee_input</b>
<ul>
<dt><a name="inputformat"><b>inputformat (character string)</b></a>
<ul>
<li>'Towhee' : reads in the input variables following the format for Towhee. This is the format
that is described in this file.</li>
<li>'LAMMPS' : reads in the input variables from the lammps_input and lammps_data files. Outputs files
suitable for use with Towhee.</li>
<li>'Quest' : reads in the input variables from the quest_input file. Runs energy calculations for a database
of conformations. This feature is not documented as it is still under development.</li>
</ul>
<dt><a name="temperature"><b>temperature (double precision)</b></a>
<ul>
<li> The temperature in Kelvin.</li>
</ul>
<dt><a name="pressure"><b>pressure (double precision)</b></a>
<ul>
<li> The external pressure in kPa. This in only used if lnpt = .true.,
but you have to list a value regardless.</li>
</ul>
<dt><a name="nmolty"><b>nmolty (integer)</b></a>
<ul>
<li> The total number of molecule types in the simulation. This must
be less then or equal to NTMAX (see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="nmolectyp"><b>nmolectyp (integer)</b></a>
<ul>
<li> The number of molecules of each type (listed sequentially on
a single line). This variable was formerly known as <b>moltyp</b>.</li>
</ul>
<dt><a name="numboxes"><b>numboxes (integer)</b></a>
<ul>
<li> The number of simulation boxes in the system. This value must
be less than or equal to MAXBOXES (set in <a href="../code/code_manual.html#preproc">preproc.h</a>). Note that
many of the variables below depend upon numboxes as you will have
to provide information for each box (such as box lengths).</li>
</ul>
<dt><a name="stepstyle"><b>stepstyle (character string of length 10)</b></a>
<ul>
The different settings for stepstyle require a different set of variables afterwards. For each
<b>stepstyle</b> I list a description of the resutling step style and the set of variables that must follow.
<li> 'cycles': Run a Monte Carlo simulation for <b>nstep</b> Monte Carlo cycles.
A cycle is equal to <b>N</b> Monte Carlo moves, where <b>N</b> is the number of molecules in the system.</li>
<ul>
<dt><a name="stepstyle_cycles_nstep"><b>nstep (integer)</b></a>
<ul>
<li> The number of Monte Carlo steps to perform where each step is a cycle.</li>
</ul>
</dt>
</ul>
<li> 'moves': Run a Monte Carlo simulation for <b>nstep</b> Monte Carlo moves.</li>
<ul>
<dt><a name="stepstyle_moves_nstep"><b>nstep (integer)</b></a>
<ul>
<li> The number of Monte Carlo steps to perform where each step is a single move.</li>
</ul>
</dt>
</ul>
<li> 'minimize': Perform a minimization.</li>
<ul>
<dt><a name="stepstyle_minimize_optstyle"><b>optstyle (integer)</b></a>
<ul>
<li> 1: Use the Broyden-Fletcher-Goldfarb-Shanno variant of the variable-metric
or quasi-newton method for minimization. The suggested reference in Numirical
Recipes was <a href="../references.html#polak_1971">Polak 1971</a>.</li>
</ul>
</dt>
<dt><a name="stepstyle_minimize_mintol"><b>mintol (double precision)</b></a>
<ul>
<li> The convergance tolerance for the minimization.</li>
</ul>
</dt>
</ul>
</ul>
<dt><a name="printfreq"><b>printfreq (integer)</b></a>
<ul>
<li> The step frequency for outputting information about the system to stdout (fort.6).
The information is the number of Monte
Carlo steps performed thus far during the run, the total energy
in each box, the x-box length of each box, the pressure of each
box, and the number of molecules of each type in each box. This variable was formerly
known as <b>iprint</b></li>
</ul>
<dt><a name="blocksize"><b>blocksize (integer)</b></a>
<ul>
<li> The size of the blocks for computing block averages. If you want
this to be meaningful then blocksize should divide cleanly into nstep.
The quantities that are averaged (in each simulation box) are the
specific density, the pressure, all of the energy terms, the chemical
potential of each molecule type, number density of each molecule
type, and the mole fractions. This variable was formerly known as <b>iblock</b></li>
</ul>
<dt><a name="moviefreq"><b>moviefreq (integer)</b></a></a>
<ul>
<li> The step frequency for outputting the system conformations
to the towhee_movie file. This file is analyzed after the run using
the analyze_movie.F routine to compute a variety of distribution
functions. This file can get pretty big if you output frequently
so be careful if you have a limited amount of hard disk space available. This variable was
formerly known as <b>imovie</b></li>
</ul>
<dt><a name="backupfreq"><b>backupfreq (integer)</b></a>
<ul>
<li> The step frequency for writing a file named
towhee_backup that is suitable for use as a restart file. It overwrites
the previous version of towhee_backup each time so it does not take
up much space. Typically I set backupfreq so that I get around 10 backups
during a run. For more information about restart files look at the
manual entries for towhee_initial, towhee_backup, and towhee_final. This variable was formerly
known as <b>ibackup</b></li>
</ul>
<dt><a name="loutpdb"><b>loutpdb (logical)</b></a>
<ul>
<li> .true. if you wish to output Protein Data Bank (pdb) files for
each simulation box at the end of the run. These files are named
box_xx.pdb where xx is the simulation box number.</li>
<li> .false. if you do not want to output pdb files.</li>
</ul>
<dt><a name="loutdft"><b>loutdft (logical)</b></a>
<ul>
<li> .true. if you wish to output files for use with the Tramonto
classical density functional theory code. This outputs dft_surfaces.dat
and dft_decode.dat. See the Tramonto code for information on what
these files mean.</li>
<li> .false. if you do not want to output dft files.</li>
</ul>
<dt><a name="loutlammps"><b>loutlammps (logical)</b></a>
<ul>
<li> .true. if you wish to output files for use with the LAMMPS massively
parallel molecular dynamics code. This outputs lammps_input and
lammps_data# where the number is each of the simulation box numbers.
See the LAMMPS documentation for more information on how to read
in these files.</li>
<li> .false. if you do not want to output LAMMPS files.</li>
</ul>
<dt><a name="pressurefreq"><b>pressurefreq (integer)</b></a>
<ul>
<li> The step frequency for computing the pressure in each simulation box.
Be aware that computing the pressure is a fairly expensive
task (especially for large systems) so if you don't really care
about the computed pressure then it will pay to set <b>pressurefreq</b> to a high
value. This variable formerly known as <b>iratp</b></li>
</ul>
<dt><a name="trmaxdispfreq"><b>trmaxdispfreq (integer)</b></a>
<ul>
<li> The step frequency for updating the maximum translational (atom and center-of-mass) and rotational
displacements. They are
adjusted to try and achieve the target acceptance rates (see <b>tatraa</b>,
<b>tatrac</b>, and <b>tarot</b>). It is a good idea to do this fairly frequently
at the start of the simulation (every step or every 10 steps) in
order to get good values for the maximum displacements. Once the
acceptance rates are near their desired values I typically set <b>trmaxdispfreq</b>
to do 10 updates during a run. This variable formerly known as <b>iratio</b></li>
</ul>
<dt><a name="volmaxdispfreq"><b>volmaxdispfreq (integer)</b></a>
<ul>
<li> The step frequency for updating the maximum volume displacements.
They are adjusted to try and achieve the target acceptance
rates (see <b>tavol</b>). It is a good idea to do this fairly frequently
at the start of the simulation (every few steps) in order to get
good values for the maximum displacements. Once the acceptance rates
are near their desired values I typically set <b>volmaxdisp</b> to do 10 updates
during a run. This variable formerly known as <b>iratv</b></li>
</ul>
<dt><a name="ffnumber"><b>ffnumber (integer)</b></a>
<ul>
<li> 1 or more: reads the force field information from this number of file(s) listed in the <b>ff_filename</b>.</li>
</ul>
<dt><a name="ff_filename"><b>ff_filename (formatted character*70)</b></a>
<ul>
<li>A list of the filenames (one per line) that contain the force field information.</li>
</ul>
<dt><a name="potentyp"><b>potentyp (integer)</b></a>
<ul>
The different settings for <b>potentyp</b> require a different set of variables afterwards. For each
<b>potentyp</b> I list a description of the nonbonded potential and the set of variables that must be
specified below the <b>potentyp</b>
<li> potentyp = 0: 12-6 Lennard-Jones van der Waals potential.</li>
<ul>
<dt>If the two atoms are separated by more than 3 bonds, or are on different molecules then</dt>
<dt>U<sub>nonbond</sub> = 4 * nbcoeff(2) * [ (nbcoeff(1)/r)^12 - (nbcoeff(1)/r)^6 ]</dt>
<dt>else if the two atoms are separated by exactly 3 bonds then</dt>
<dt>U<sub>nonbond</sub> = 4 * nbcoeff(4) * [ (nbcoeff(3)/r)^12 - (nbcoeff(3)/r)^6 ]</dt>
<dt><a name="p0_mixrule"><b>mixrule (integer)</b></a>
<ul>
<li> mixrule = 0: Lorentz-Berthelot (arithmetic mean of sigma, geometric
mean of epsilon) mixing rules.</li>
<li> mixrule = 1: Geometric (geometric mean of sigma and epsilon)
mixing rules.</li>
<li> mixrule = 3: Gromos (geometric mean of sigma and epsilon with
some special cases) mixing rules.</li>
<li> mixrule = 4: Explicit (defined in towhee_ff files) mixing rules.</li>
</ul>
<dt><a name="lshift"><b>lshift (logical)</b></a>
<ul>
<li> .true. if you want the nonbonded van der Waals potential to be
shifted so that it is zero at the cutoff.</li>
<li> .false. if you do not want to shift the nonbonded van der Waals
potential.</li>
</ul>
<dt><a name="ltailc"><b>ltailc (logical)</b></a>
<ul>
<li>.true. if you want to apply analytical tail corrections for the
portion of the van der Waals potential that is past the cutoff.
Note that you cannot have a shifted potential and tail corrections
at the same time.</li>
<li>.false. if you do not want analytical tail corrections for van
der Waals.</li>
</ul>
<dt><a name="rmin"><b>rmin (double precision)</b></a>
<ul>
<li> A hard inner cutoff that can speed computation for Lennard-Jones
systems, and is required to avoid the potential hitting infinity
for exponential repulsion systems which also contain point charges.
This should be set smaller than the smallest radius of any atom.
Generally I set this to 0.5 or 1.0 Angstroms.</li>
</ul>
<dt><a name="rcut"><b>rcut (double precision)</b></a>
<ul>
<li> The van der Waals potential cutoff in Angstroms.</li>
</ul>
<dt><a name="rcutin"><b>rcutin (double precision)</b></a>
<ul>
<li> The inner nonbonded cutoff used in configurational-bias Monte
Carlo moves. This dual-cutoff method can speed configurational-bias
computations by at least a factor of 2, without affecting the acceptance
rate. The inner cutoff is used during the growth procedure, and
the full potential is calculated at the end of the move and everything
is fixed up in the acceptance criteria. I typically set this to
5 Angstroms for noncoulombic simulations, and to 10 Angstroms for
coulombic simulations.</li>
</ul>
</ul>
<li> potentyp = 1: 9-6 Lennard-Jones van der Waals potential</li>
<ul>
<dt>If the two atoms are separated 3 or more bonds, or are on different molecules then</dt>
<dt>U<sub>nonbond</sub> = nbcoeff(2) * [ 2*(nbcoeff(1)/r)^9 - 3*(nbcoeff(1)/r)^6 ]</dt>
<dt><a name="p0_mixrule"><b>mixrule (integer)</b></a>
<ul>
<li> mixrule = 2: Compass (sixth order combination of sigma and epsilon)
mixing rules.</li>
</ul>
<dt><a name="lshift"><b>lshift (logical)</b></a>
<ul>
<li> .true. if you want the nonbonded van der Waals potential to be
shifted so that it is zero at the cutoff.</li>
<li> .false. if you do not want to shift the nonbonded van der Waals
potential.</li>
</ul>
<dt><a name="ltailc"><b>ltailc (logical)</b></a>
<ul>
<li>.true. if you want to apply analytical tail corrections for the
portion of the van der Waals potential that is past the cutoff.
Note that you cannot have a shifted potential and tail corrections
at the same time.</li>
<li>.false. if you do not want analytical tail corrections for van
der Waals.</li>
</ul>
<dt><a name="rmin"><b>rmin (double precision)</b></a>
<ul>
<li> A hard inner cutoff that can speed computation for Lennard-Jones
systems, and is required to avoid the potential hitting infinity
for exponential repulsion systems which also contain point charges.
This should be set smaller than the smallest radius of any atom.
Generally I set this to 0.5 or 1.0 Angstroms.</li>
</ul>
<dt><a name="rcut"><b>rcut (double precision)</b></a>
<ul>
<li> The van der Waals potential cutoff in Angstroms.</li>
</ul>
<dt><a name="rcutin"><b>rcutin (double precision)</b></a>
<ul>
<li> The inner nonbonded cutoff used in configurational-bias Monte
Carlo moves. This dual-cutoff method can speed configurational-bias
computations by at least a factor of 2, without affecting the acceptance
rate. The inner cutoff is used during the growth procedure, and
the full potential is calculated at the end of the move and everything
is fixed up in the acceptance criteria. I typically set this to
5 Angstroms for noncoulombic simulations, and to 10 Angstroms for
coulombic simulations.</li>
</ul>
</ul>
<li> potentyp = 2: Exponential-6 van der Waals potential</li>
<ul>
<dt>If the two atoms are separated by more than 3 bonds, or are on different molecules then</dt>
<dt>U<sub>nonbond</sub> = nbcoeff(1)/r^6 + nbcoeff(2) * exp[nbcoeff(3)*r]</dt>
<dt><a name="p0_mixrule"><b>mixrule (integer)</b></a>
<ul>
<li> mixrule = 4: Explicit (defined in towhee_ff files) mixing rules.</li>
</ul>
<dt><a name="lshift"><b>lshift (logical)</b></a>
<ul>
<li> .true. if you want the nonbonded van der Waals potential to be
shifted so that it is zero at the cutoff.</li>
<li> .false. if you do not want to shift the nonbonded van der Waals
potential.</li>
</ul>
<dt><a name="ltailc"><b>ltailc (logical)</b></a>
<ul>
<li>.true. if you want to apply analytical tail corrections for the
portion of the van der Waals potential that is past the cutoff.
Note that you cannot have a shifted potential and tail corrections
at the same time.</li>
<li>.false. if you do not want analytical tail corrections for van
der Waals.</li>
</ul>
<dt><a name="rcut"><b>rcut (double precision)</b></a>
<ul>
<li> The van der Waals potential cutoff in Angstroms.</li>
</ul>
<dt><a name="rcutin"><b>rcutin (double precision)</b></a>
<ul>
<li> The inner nonbonded cutoff used in configurational-bias Monte
Carlo moves. This dual-cutoff method can speed configurational-bias
computations by at least a factor of 2, without affecting the acceptance
rate. The inner cutoff is used during the growth procedure, and
the full potential is calculated at the end of the move and everything
is fixed up in the acceptance criteria. I typically set this to
5 Angstroms for noncoulombic simulations, and to 10 Angstroms for
coulombic simulations.</li>
</ul>
</ul>
<li> potentyp = 3: Hard Sphere potential</li>
<ul>
<dt>If the two atoms are separated by more than 3 bonds, or are on different molecules then</dt>
<dt>U<sub>nonbond</sub> = Infinity if r <= nbcoeff(1), or 0 otherwise</dt>
<dt><a name="p0_mixrule"><b>mixrule (integer)</b></a>
<ul>
<li> mixrule = 5: Hard sphere (arithmetic mean of sigmas) mixing rules.</li>
</ul>
<dt><a name="rcutin"><b>rcutin (double precision)</b></a>
<ul>
<li> The inner nonbonded cutoff used in configurational-bias Monte
Carlo moves. This dual-cutoff method can speed configurational-bias
computations by at least a factor of 2, without affecting the acceptance
rate. The inner cutoff is used during the growth procedure, and
the full potential is calculated at the end of the move and everything
is fixed up in the acceptance criteria. If you are using the Hard Sphere potential
without coulombic interations then just set this to something large (like 100), if you are
using coulombic interactions then I would suggest a value of 5 sigma.</li>
</ul>
</ul>
<li> potentyp = -3: Repulsive Sphere potential. Added to towhee in version 1.4.6. This
is used to help setup and equilibrate a hard sphere system where it is sometimes challenging to create an
initial conformation with no overlaps. Use the -3 option to equilibrate until the nonbonded potential energy
is 0.0, and then switch back to the normal hard sphere potential.</li>
<ul>
<dt>If the two atoms are separated by more than 3 bonds, or are on different molecules then</dt>
<dt>U<sub>nonbond</sub> = 1d5 + 1d5 * (nbcoeff(1)^2 - r^2) if r <= nbcoeff(1), or 0 otherwise</dt>
<dt><a name="p0_mixrule"><b>mixrule (integer)</b></a>
<ul>
<li> mixrule = 5: Hard sphere (arithmetic mean of sigmas) mixing rules.</li>
</ul>
<dt><a name="rcutin"><b>rcutin (double precision)</b></a>
<ul>
<li> The inner nonbonded cutoff used in configurational-bias Monte
Carlo moves. This dual-cutoff method can speed configurational-bias
computations by at least a factor of 2, without affecting the acceptance
rate. The inner cutoff is used during the growth procedure, and
the full potential is calculated at the end of the move and everything
is fixed up in the acceptance criteria. If you are using the Hard Sphere potential
without coulombic interations then just set this to something large (like 100), if you are
using coulombic interactions then I would suggest a value of 5 sigma.</li>
</ul>
</ul>
<li> potentyp = 4: Exponential plus 12-6 Lennard-Jones van der Waals
potential</li>
<ul>
<dt>If the two atoms are separated by more than 3 bonds, or are on different molecules then</dt>
<dt>U<sub>nonbond</sub> = nbcoeff(1)/r^6 + nbcoeff(2)/r^12 + nbcoeff(3)*exp[nbcoeff(4)*r]</dt>
<dt><a name="p0_mixrule"><b>mixrule (integer)</b></a>
<ul>
<li> mixrule = 4: Explicit (defined in towhee_ff files) mixing rules.</li>
</ul>
<dt><a name="lshift"><b>lshift (logical)</b></a>
<ul>
<li> .true. if you want the nonbonded van der Waals potential to be
shifted so that it is zero at the cutoff.</li>
<li> .false. if you do not want to shift the nonbonded van der Waals
potential.</li>
</ul>
<dt><a name="ltailc"><b>ltailc (logical)</b></a>
<ul>
<li>.true. if you want to apply analytical tail corrections for the
portion of the van der Waals potential that is past the cutoff.
Note that you cannot have a shifted potential and tail corrections
at the same time.</li>
<li>.false. if you do not want analytical tail corrections for van
der Waals.</li>
</ul>
<dt><a name="rmin"><b>rmin (double precision)</b></a>
<ul>
<li> A hard inner cutoff that can speed computation for Lennard-Jones
systems, and is required to avoid the potential hitting infinity
for exponential repulsion systems which also contain point charges.
This should be set smaller than the smallest radius of any atom.
Generally I set this to 0.5 or 1.0 Angstroms.</li>
</ul>
<dt><a name="rcut"><b>rcut (double precision)</b></a>
<ul>
<li> The van der Waals potential cutoff in Angstroms.</li>
</ul>
<dt><a name="rcutin"><b>rcutin (double precision)</b></a>
<ul>
<li> The inner nonbonded cutoff used in configurational-bias Monte
Carlo moves. This dual-cutoff method can speed configurational-bias
computations by at least a factor of 2, without affecting the acceptance
rate. The inner cutoff is used during the growth procedure, and
the full potential is calculated at the end of the move and everything
is fixed up in the acceptance criteria. I typically set this to
5 Angstroms for noncoulombic simulations, and to 10 Angstroms for
coulombic simulations.</li>
</ul>
</ul>
<li> potentyp = 5: Stillinger-Weber potential (see <a href="../references.html#stillinger_weber_1985">
Stillinger and Weber 1985</a>) </li>
<ul>
<dt>This is an atomic potential and can only be used with monatomic molecules in Towhee</dt>
<dt>U = nbcoeff(1)*[nbcoeff(2)*Sum u2(rij) + nbcoeff(7)*Sum u3(rij,rjk)]</dt>
<dt>u2(r_ij) = [nbcoeff(3)*(r_ij/nbcoeff(4)^-nbcoeff(5) - 1] * exp{1/[(r_ij/nbcoeff(4) - nbcoeff(6))]}
* Heaviside(nbcoeff(6) - [r_ij/nbcoeff(4)])</dt>
<dt>u3(r_ij,r_jk) = exp[ nbcoeff(8)/{r_ij/nbcoeff(4) - nbcoeff(6)} + nbcoeff(8)/{r_jk/nbcoeff(4) - nbcoeff(6)}]
* (cos(theta_ijk)-nbcoeff(9))^2 * Heaviside(nbcoeff(6)
- r_ij/nbcoeff(4)) * Heaviside(nbcoeff(6) - r_jk/nbcoeff(4))</dt>
<dt><a name="p0_mixrule"><b>mixrule (integer)</b></a>
<ul>
<li> mixrule = 4: Explicit (defined in towhee_ff files) mixing rules.</li>
</ul>
</ul>
</ul>
<dt><a name="lcoulomb"><b>lcoulomb (logical)</b></a>
<ul>
<li> .true. if you want to use point charges in the simulation. If you are using
point charges then the Ewald sum handles the long range corrections. Note that you can
essentially disable the Ewald sum by setting both <b>kalp</b> and <b>kmax</b> to zero.</li>
<li> .false. if you do want to use point charges.</li>
</ul>
<dt><a name="kalp"><b>kalp (double precision)</b></a>
<ul>
<li> Value used in the Ewald sum to compute alpha. If you set <b>kalp</b> and <b>kmax</b> both
to zero then you will effectively disable the Ewald sum. The actual Ewald sum alpha term is equal to <b>kalp</b>
divided by the shortest box length. The recommended value for <b>kalp</b> is 5.6.</li>
</ul>
<dt><a name="kmax"><b>kmax (integer)</b></a>
<ul>
<li> Maximum number of inverse space vectors to use in any dimension
for the Ewald sum. Recommended value of this parameter is 5. If
you want to set this to a larger value to may have to increase VECTORMAX
(see <a href="../code/code_manual.html#preproc">preproc.h</a>). Note that you can effectively
disable the Ewald sum by setting both <b>kalp</b> and <b>kmax</b> to zero.</li>
</ul>
<dt><a name="dielect"><b>dielect (double precision)</b></a>
<ul>
<li> The dielectric constant used when computing coulombic interactions. Generally this
should be set to 1.0 as the solvated system will act as the screening that the dielectric
constant is intended to represent. If you are performing a simulation without any solvent
(for example a protein without the water) you might want to set the dielectric constant to represent
the missing solvent.</li>
</ul>
<dt><a name="nhrdfld"><b>nhrdfld (integer)</b></a>
<ul>
<li> Number of hard walls you wish to include in the simulation. As
this is not typically used, it has a slightly different input style
from normal. If nhrdfld = 0 then none of the hrd* variables are
required in the input file. Otherwise for next nhrdfld lines each
line should contain the values for hrdbox(ifield), hrdxyz(ifield),
hrdcen(ifield), hrdrad(ifield) where those variables have the following
meanings.</li>
<li><b>hrdbox (integer)</b></a>
<ul>
<li> This is the number of the simulation box which contains this
hard wall. Must range from 1 to numboxes.</li>
</ul>
</li>
<li><b>hrdxyz (integer)</b></a>
<ul>
<li> 1: hard wall is perpendicular to the x-axis (in the yz plane)</li>
<li> 2: hard wall is perpendicular to the y-axis (in the xz plane)</li>
<li> 3: hard wall is perpendicular to the z-axis (in the xy plane)</li>
</ul>
</li>
<li><b>hrdcen (double precision)</b></a>
<ul>
<li> Position of the center of the hard wall. Must be between
0.0 and the box length of the axis that is perpendicular to
the wall.</li>
</ul>
</li>
<li><b>hrdrad (double precision)</b></a>
<ul>
<li> Radius of the hard wall. The wall will exclude all atoms
whose centers are within this radius regardless of the potentyp
or any of the atom parameters. Yes I know this is a strange
cd way to run a hard wall. The wall is felt through the periodic
boundaries.</li>
</ul>
</li>
</ul>
<dt><a name="nljfld"><b>nljfld (integer)</b></a>
<ul>
<li> Number of 9-3 Lennard-Jones walls you wish to include in the
simulation. As this is not typically used, it has a slightly different
input style from normal. If nljfld = 0 then none of the ljf* variables
are required in the input file. Otherwise for next nljfld lines
each line should contain the values for ljfbox(ifield), ljfxyz(ifield),
ljfcen(ifield), ljfsig(ifield), ljfeps(ifield) ,ljfcut(ifield),
and ljfdir(ifield) where those variables have the following meanings.</li>
<li><b>ljfbox (integer)</b></a>
<ul>
<li> This is the number of the simulation box which contains this
Lennard-Jones wall. Must range from 1 to numboxes.</li>
</ul>
</li>
<li><b>ljfxyz (integer)</b></a>
<ul>
<li> 1: Lennard-Jones wall is perpendicular to the x-axis (in
the yz plane)</li>
<li> 2: Lennard-Jones wall is perpendicular to the y-axis (in
the xz plane)</li>
<li> 3: Lennard-Jones wall is perpendicular to the z-axis (in
the xy plane)</li>
</ul>
</li>
<li><b>ljfcen (double precision)</b></a>
<ul>
<li> Position of the center of the Lennard-Jones wall. Must be
between 0.0 and the box length of the axis that is perpendicular
to the wall.</li>
</ul>
</li>
<li><b>ljfsig (double precision)</b></a>
<ul>
<li> Sigma parameter for the 9-3 Lennard-Jones wall.</li>
</ul>
</li>
<li><b>ljfeps (double precision)</b></a>
<ul>
<li> Epsilon parameter for the 9-3 Lennard-Jones wall. All atoms
in the system interact with the wall via the potential U = sqrt(2/5)
ljfeps [ (1/5) (ljfsig/r)^9 - (3/2) (ljfsig/r)^3 ] regardless
of the potentyp and parameters of each atom. Yes I know this
is a strange way to run a 9-3 Lennard-Jones wall, but I don't
use it much.</li>
</ul>
</li>
<li><b>ljfcut (double precision)</b></a>
<ul>
<li> The distance beyond which the wall-atom interactions are
not computed. This potential is cut, not shifted, and never
has tail corrections no matter how lshift and ltailc are set.</li>
</ul>
</li>
<li><b>ljfdir (integer)</b></a>
<ul>
<li> -1: Atoms only interact with the "left" face of this wall.
This extends through the periodic boundary.</li>
<li> 1: Atoms only interact with the "right" face of this wall.
This extends through the periodic boundary.</li>
</ul>
</li>
</ul>
<dt><a name="nragfld"><b>nlragfld (integer)</b></a>
<ul>
<li> Number of <a href="../references.html#raghavan_et_al_1991">Raghavan <i>et al.</i></a>
style 111 metal walls you wish to include in the
simulation. Note that all of the field variables have a slightly different
input style from normal. If nragfld = 0 then none of the rag* variables
are required in the input file. Otherwise for next nragfld lines
each line should contain the values for ragbox(ifield), ragxyz(ifield),
ragcen(ifield), ragcut(ifield), ragdir(ifield) , and raglat(ifield)
where those variables have the following meanings.</li>
<li><b>ragbox (integer)</b></a>
<ul>
<li> This is the number of the simulation box which contains this
Raghavan wall. Must range from 1 to numboxes.</li>
</ul>
</li>
<li><b>ragxyz (integer)</b></a>
<ul>
<li> 1: Raghavan wall is perpendicular to the x-axis (in
the yz plane)</li>
<li> 2: Raghavan wall is perpendicular to the y-axis (in
the xz plane)</li>
<li> 3: Raghavan wall is perpendicular to the z-axis (in
the xy plane)</li>
</ul>
</li>
<li><b>ragcen (double precision)</b></a>
<ul>
<li> Position of the center of the Raghavan wall. Must be
between 0.0 and the box length of the axis that is perpendicular
to the wall.</li>
</ul>
</li>
<li><b>ragcut (double precision)</b></a>
<ul>
<li> The distance beyond which the Raghavan wall-atom interactions are
not computed. This potential is cut, not shifted, and never
has tail corrections no matter how lshift and ltailc are set.</li>
</ul>
</li>
<li><b>ragdir (integer)</b></a>
<ul>
<li> -1: Atoms only interact with the "left" face of this wall.
This extends through the periodic boundary.</li>
<li> 1: Atoms only interact with the "right" face of this wall.
This extends through the periodic boundary.</li>
</ul>
</li>
<li><b>raglat (double precision)</b></a>
<ul>
<li> The lattice spacing of the metal atoms in the 111 wall. This is the "a" parameter
in the Raghavan wall potential.</li>
</ul>
</li>
</ul>
<dt><a name="ldftfld"><b>ldftfld (logical)</b></a>
<ul>
<li> .true. if you are using the classical density functional theory code Tramonto to
compute a solvation free energy.</li>
<li> .false. otherwise</li>
</ul>
<dt><a name="linit"><b>linit (logical)</b></a>
<ul>
<li> .true. if you are starting the simulation and wish to generate
the positions of all of the atoms, assign initial box lengths and
maximum displacements.</li>
<li> .false. if you want to continue the simulation by reading in
box lengths, maximum displacements, and coordinates from towhee_initial.</li>
</ul>
<dt><a name="initstyle"><b>initstyle (integer)</b></a>
<ul>
One line for each simulation box in the system. Each line contains
a value for each molecule type.
<li> 0: A template for this molecule type is created using configurational-bias.
This template is then replicated throughout the simulation box to
generate an initial configuration.</li>
<li> 1: A template for this molecule type is read from towhee_template.
This template is then replicated throughout the simulation box to
generate an initial configuration.</li>
<li> 2: The coordinates for each atom are read from towhee_coords.
This is useful if you are starting from a different file format
(such as pdb), or have another code for building an initial configuration.</li>
</ul>
<dt><a name="hmatrix"><b>hmatrix (double precision)</b></a>
<ul>
<li> The initial box dimensions (Angstroms) for the three
vectors that describe the simulation box. There are nine
entries (3 for each of the 3 vectors) in total for each
simulation box. These are listed one vector at a time, with the three numbers
which make up each vector listed on the same line. Note that the coordinate system
you choose does not have to be orthogonal, but it must follow the right hand rule. The
three vectors must also all be at least 45 degrees apart. Note that if you wish to use
a rectangular box then only the diagonal elements of hmatrix will be non-zero, and these
will be equal to the boxlengths in the x, y, and z dimensions.</li>
</ul>
<dt><a name="initmol"><b>initmol (integer)</b></a>
<ul>
<li> The initial number of each type of molecule in each box (one
line per box).</li>
</ul>
<dt><a name="inix"><b>inix, iniy, iniz (integer)</b></a>
<ul>
<li> The initial number of molecules in each direction in each box.
The product of inix*iniy*iniz must be greater than or equal to the
initial number of molecules in that box (the sums of initmol). While these
are labeled x, y, and z they actually correspond to the three coordinate
vectors.</li>
</ul>
<dt><a name="inimix"><b>inimix (integer)</b></a>
<ul>
One line for each simulation box in the system.
<li> -1: molecules are initially placed in each box in alternating
order.</li>
<li> 0: molecules are initially placed in each box in random order.</li>
<li> 1: molecules are initially placed in each box order. Thus all
molecules of type 1 are placed in a box before any molecules of
type 2. If you are using initstyle = 2 then this is the only valid
option and the other options will be reset to this option by the
code.</li>
</ul>
<dt><a name="lnpt"><b>lnpt (logical)</b></a>
<ul>
<li> .true. if you want the volume moves to be performed in an exchange
with an external pressure bath (isobaric isothermal ensemble or
isobaric-isothermal Gibbs ensemble).</li>
<li> .false. if you do want the total volume of the system to be conserved.
In the case of a multi-box simulation this exchanges volume between
pairs of boxes (canonical Gibbs ensemble), in a single box case
no volume moves are allowed (canonical ensemble).</li>
</ul>
<li> Note: the pm* variables are used to determine which move type to
perform every time we want to do a Monte Carlo move. A move is selected
by choosing a random number between 0.0 and 1.0 and then going down
the list of pm* until you find one which has a value higher than the
random number. At least one of the variables must be set to 1.0.
A similar procedure is
performed when we want to determine which boxes or molecule types
to perform the selected move upon. These are done using the pm**pr
and pm**mt arrays.</li>
<hr>
<dt><a name="pmvol"><b>pmvol (double precision)</b></a>
<ul>
<li> Probability of performing a volume move. If (lnpt = .true.) then
a single box is selected and it exchanges volume with an external
pressure bath (see pressure). If (lnpt = .false. and numboxes >
1) a pair of boxes are selected and volume is exchanged between
them.</li>
</ul>
<dt><a name="pmvlpr"><b>pmvlpr (double precision)</b></a>
<ul>
<li> Probability of performing a volume move on a particular box,
or box pair. All of these variables are listed on a single line
If (lnpt = .true.) then a value of pmvlpr is listed for each box.
If (lnpt = .false.) then a value is listed for each pair of simulation
boxes where the pairs are ordered (1,2), (1,3), ... (1,numboxes),
(2,3), ... (numboxes-1,numboxes).</li>
</ul>
<dt><a name="rmvol"><b>rmvol (double precision)</b></a>
<ul>
<li> The initial volume maximum displacement. If this is an isobaric-isothermal
ensemble (lnpt = .true.) then this is the initial maximum volume
displacement (cubic Angstroms) in each box. If this is the canonical
Gibbs ensemble (lnpt = .false. and numboxes > 1 ) then this is the
maximum displacement (logarithmic space) for each pair of boxes.
As the simulation progresses, these values will be updated for each
box, or each pair of boxes (see iratv).</li>
</ul>
<dt><a name="tavol"><b>tavol (double precision)</b></a>
<ul>
<li> The target acceptance rate for the volume move. Must be a value
between 0.0 and 1.0. The volume displacement (rmvol) is periodically
adjusted (see iratv) to yield this acceptance rate. I typically use
a value of 0.5, though some researchers prefer smaller values.</li>
</ul>
<hr>
<dt><a name="pmcell"><b>pmcell (double precision)</b></a>
<ul>
<li> Probability of performing a unit cell adjustment move. If (lnpt = .true.) then
a single box is selected and a single hmatrix element is changed. This results in a volume
exchange with a fictional external pressure bath (see pressure). If (lnpt = .false. and numboxes >
1) a pair of boxes are selected. One of the boxes is then selected according to the pmcellpt
variable and a single hmatrix element is changed in that box. This results in a change of volume for
the first box which is countered by isotropically changing the volume in the second box.</li>
</ul>
<dt><a name="pmcellpr"><b>pmcellpr (double precision)</b></a>
<ul>
<li> Probability of performing a unit cell adjustment move on a particular box,
or box pair. All of these variables are listed on a single line
If (lnpt = .true.) then a value of pmvlpr is listed for each box.
If (lnpt = .false.) then a value is listed for each pair of simulation
boxes where the pairs are ordered (1,2), (1,3), ... (1,numboxes),
(2,3), ... (numboxes-1,numboxes).</li>
</ul>
<dt><a name="pmcellpt"><b>pmcellpt (double precision)</b></a>
<ul>
<li> Probability of selecting the first box of the pair as the box to perform the non-isotropic
volume move upon, while its partner undergoes an isotropic volume move. This variable is only
meaningful if (lnpt=.false.). Note that you can choose to perform the non-isotropic volume
move always on the same box and this might be useful if you are doing a solid-vapor equilibria
calculation.</li>
</ul>
<dt><a name="rmcell"><b>rmcell (double precision)</b></a>
<ul>
<li> The initial unit cell adjustment maximum displacement. In all cases, this is the maximum amount
(in Angstroms) that a single element of the hmatrix can change in a single unit cell move.
Note, the in the canonical Gibbs ensemble case it is possible for the isotropic box to undergo
an hmatrix change that is larger than this value as that box simply makes up for the volume change
caused by the non-isotropic adjustment in the first box. As the simulation progresses, these values are
updated for each box with a frequency controlled by <b>iratv</b>.</li>
</ul>
<dt><a name="tacell"><b>tacell (double precision)</b></a>
<ul>
<li> The target acceptance rate for the unit cell adjustment move. Must be a value
between 0.0 and 1.0. The unit cell displacement (rmcell) is periodically
adjusted (see iratv) to yield this acceptance rate. I typically use
a value of 0.5.</li>
</ul>
<hr>
<dt><a name="pm2boxrbswap"><b>pm2boxrbswap (double precision)</b></a>
<ul>
<li>Probability of performing a rotational-bias interbox molecule
transfer move. This move takes a molecule out of one box and tries
to place it in another box. The molecule is grown using <b>nchoi1</b> attempted different orientations
and position (of the center-of-mass) for the new molecule.</li>
</ul>
<dt><a name="pm2rbswmt"><b>pm2rbswmt (double precision)</b></a>
<ul>
<li>Probability of performing a rotational-bias interbox molecule
transfer move on each type of molecule in the system.</li>
</ul>
<dt><a name="pm2rbswpr"><b>pm2rbswpr (double precision)</b></a>
<ul>
<li> Probability of performing a rotational-bias interbox molecule transfer move
between each pair of boxes in the system. The box pairs are ordered
(1,2), (1,3), ... (1,numboxes), (2,3), ... (numboxes-1,numboxes).</li>
</ul>
<hr>
<dt><a name="pm2boxcbswap"><b>pm2boxcbswap (double precision)</b></a>
<ul>
<li>Probability of performing a configurational-bias interbox molecule
transfer move. This move takes a molecule out of one box and tries
to place it in another box. The molecule is grown using <a href="../algorithm/cbmc.html">coupled-decoupled
configurational-bias Monte Carlo</a>. This variable was formerly known as <b>pmswap</b></li>
</ul>
<dt><a name="pm2cbswmt"><b>pm2cbswmt (double precision)</b></a>
<ul>
<li>Probability of performing a configurational-bias interbox molecule
transfer move on each type of molecule in the system. This variable was formerly known as
<b>pmswmt</b>.</li>
</ul>
<dt><a name="pm2cbswpr"><b>pm2cbswpr (double precision)</b></a>
<ul>
<li> Probability of performing a configurational-bias interbox molecule transfer move
between each pair of boxes in the system. The box pairs are ordered
(1,2), (1,3), ... (1,numboxes), (2,3), ... (numboxes-1,numboxes).
This variable was formerly known as <b>pmswpr</b></li>
</ul>
<hr>
<dt><a name="pm1boxcbswap"><b>pm1boxcbswap (double precision)</b></a>
<ul>
<li> Probability of performing an intrabox configurational-bias molecule
transfer move. This move takes a molecule out of one box and tries
to place it back into the same box. The molecule is grown using
<a href="../algorithm/cbmc.html">coupled-decoupled configurational-bias Monte Carlo</a>
This variable was formerly known as <b>pmiswp</b></li>
</ul>
<dt><a name="pm1bcbswmt"><b>pm1bcbswmt (double precision)</b></a>
<ul>
<li> Probability of performing an intrabox configurational-bias molecule
transfer move on each type of molecule in the system.
This variable was formerly known as <b>pmismt</b></li>
</ul>
<hr>
<dt><a name="pmavb1"><b>pmavb1 (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 1, as described in
<a href="../references.html#chen_siepmann_2000">Chen and Siepmann 2000</a>.
This is useful for forming and destroying clusters in simulations with molecules that tend to aggregate together.
</ul>
<dt><a name="pmavb1in"><b>pmavb1in (double precision)</b></a>
<ul>
<li> Probability of trying to move a molecule into an inner region for aggregation volume bias move of type 1.</li>
</ul>
<dt><a name="pmavb1mt"><b>pmavb1mt (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 1 where a molecule of a certain type is moved.
This is an array with one element for each molecule type.</li>
</ul>
<dt><a name="pmavb1ct"><b>pmavb1ct (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 1 where the molecule target is of a certain type.
The molecule that is moved is chosen according to <b>pmavb1mt</b> and then the type of molecule that is used as a reference
for determining the inner and outer regions is found using this variable.
This is a two dimensional array and uses one line of text for each type of molecule in the system.</li>
</ul>
<dt><a name="avb1rad"><b>avb1rad (double precision)</b></a>
<ul>
<li> The radius used to define the inner and outer volumes in the aggregation volume bias move of type 1.
The distance is specified in Angstroms and must be greater than zero, but less than or equal to <b>rcut</b>.</li>
</ul>
<hr>
<dt><a name="pmavb2"><b>pmavb2 (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 2, as described in
<a href="../references.html#chen_siepmann_2001">Chen and Siepmann 2001</a>.
This is useful for forming and destroying clusters in simulations with molecules that tend to aggregate together.
</ul>
<dt><a name="pmavb2in"><b>pmavb2in (double precision)</b></a>
<ul>
<li> Probability of trying to move a molecule into an inner region for aggregation volume bias move of type 2.</li>
</ul>
<dt><a name="pmavb1mt"><b>pmavb2mt (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 2 where a molecule of a certain type is moved.
This is an array with one element for each molecule type.</li>
</ul>
<dt><a name="pmavb1ct"><b>pmavb2ct (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 2 where the molecule target is of a certain type.
The molecule that is moved is chosen according to <b>pmavb2mt</b> and then the type of molecule that is used as a reference
for determining the inner and outer regions is found using this variable.
This is a two dimensional array and uses one line of text for each type of molecule in the system.</li>
</ul>
<dt><a name="avb1rad"><b>avb2rad (double precision)</b></a>
<ul>
<li> The radius used to define the inner and outer volumes in the aggregation volume bias move of type 2.
The distance is specified in Angstroms and must be greater than zero, but less than or equal to <b>rcut</b>.</li>
</ul>
<hr>
<dt><a name="pmavb1"><b>pmavb3 (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 3, as described in
<a href="../references.html#chen_siepmann_2001">Chen and Siepmann 2001</a>.
This is useful for transfering molecules between clusters.
</ul>
<dt><a name="pmavb1mt"><b>pmavb3mt (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 3 where a molecule of a certain type is moved.
This is an array with one element for each molecule type.</li>
</ul>
<dt><a name="pmavb1ct"><b>pmavb3ct (double precision)</b></a>
<ul>
<li> Probability of performing an aggregation volume bias move of type 3 where the molecule target is of a certain type.
The molecule that is moved is chosen according to <b>pmavb1mt</b> and then the type of molecule that is used as a reference
for determining the inner and outer regions is found using this variable.
This is a two dimensional array and uses one line of text for each type of molecule in the system.</li>
</ul>
<dt><a name="avb1rad"><b>avb3rad (double precision)</b></a>
<ul>
<li> The radius used to define the inner and outer volumes in the aggregation volume bias move of type 3.
The distance is specified in Angstroms and must be greater than zero, but less than or equal to <b>rcut</b>.</li>
</ul>
<hr>
<dt><a name="pmcb"><b>pmcb (double precision)</b></a>
<ul>
<li> Probability of performing a molecule regrowth move
on a molecule without regard to which box the molecule is currently
located in. This move chooses a molecule of the appropriate type
at random, selects an atom of the molecule at random, and then regrows
the molecule either entirely (if a random number < pmall) or in
all directions except for one. The molecule is regrown using
<a href="../algorithm/cbmc.html">configurational-bias</a>.</li>
</ul>
<dt><a name="pmcbmt"><b>pmcbmt (double precision)</b></a>
<ul>
<li> Probability of performing a molecule regrowth on
each type of molecule in the system.</li>
</ul>
<dt><a name="pmall"><b>pmall (double precision)</b></a>
<ul>
<li> pmall is the probability that a molecule regrowth move will regrow
the entire molecule. This is listed for each molecule type in the
simulation.</li>
</ul>
<hr>
<dt><a name="pmback"><b>pmback (double precision)</b></a>
<ul>
<li> Probability of performing configurational-bias fixed-endpoint
regrowth of a portion of the protein backbone. This selects an atom along the
peptide backbone, chooses another backbone atom that is connected by three bonds
to the first atom, and then regrows all of the atoms inbetween these two atoms.</li>
</ul>
<dt><a name="pmbkmt"><b>pmbkmt (double precision)</b></a>
<ul>
<li> Probability of performing a backbone regrowth
move on each type of molecule in the system.</li>
</ul>
<hr>
<dt><a name="pmplane"><b>pmplane (double precision)</b></a>
<ul>
<li> Probability of performing a plane shift move. This move displaces all of the
molecules whose center of mass lies in a plane of width <b>planewidth</b>. A new trial
position for the center of the plane of atoms is generated uniformly across the available
plane.</li>
</ul>
<dt><a name="pmplanebox"><b>pmplanebox (double precision)</b></a>
<ul>
<li> Probability of performing a plane shift in each of the simulation boxes. List one
value for each simulation box. At least one of the boxes must have a value of 1.0d0.</li>
</ul>
<dt><a name="planewidth"><b>planewidth (double precision)</b></a>
<ul>
<li> The width of the plane for the plane shift move. Any molecule whose center of mass is
within a plane of this thickness (whose position is chosen uniformly along one axis) will
move during the plane shift move. The value of planewidth must be greater than 0.0d0 and less
than the shortest boxlength.</li>
</ul>
<hr>
<dt><a name="pmrow"><b>pmrow (double precision)</b></a>
<ul>
<li> Probability of performing a row shift move. This move displaces all of the
molecules whose center of mass lies in a row of diameter <b>rowwidth</b>. A new trial
position for the center of the row of atoms is generated uniformly across the available
row.</li>
</ul>
<dt><a name="pmrowbox"><b>pmrowbox (double precision)</b></a>
<ul>
<li> Probability of performing a row shift in each of the simulation boxes. List one
value for each simulation box. At least one of the boxes must have a value of 1.0d0.</li>
</ul>
<dt><a name="rowwidth"><b>rowwidth (double precision)</b></a>
<ul>
<li> The width of the plan for the row shift move. Any molecule whose center of mass is
within a row of this thickness (whose position is chosen uniformly along one axis) will
move during the row shift move. The value of rowwidth must be greater than 0.0d0 and less
than the shortest boxlength.</li>
</ul>
<hr>
<dt><a name="pmtraat"><b>pmtraat (double precision)</b></a>
<ul>
<li> Probability of performing a single-atom translation
move a molecule without regard to which box the molecule is currently
located in. This move chooses a molecule of the appropriate type
at random, selects an atom of the molecule at random, chooses the
x,y, or z direction at random, and then attempts to displace the
atom a random distance between -rmtraa and +rmtraa in that direction.</li>
</ul>
<dt><a name="pmtamt"><b>pmtamt (double precision)</b></a>
<ul>
<li> Probability of performing a single-atom translation
move on each type of molecule in the system.</li>
</ul>
<dt><a name="rmtraa"><b>rmtraa (double precision)</b></a>
<ul>
<li> The initial Atom-translation maximum displacement (Angstroms)
for all molecules types in all boxes. As the simulation progresses,
these values will be updated independently to give the desired acceptance
rate for each molecule type in each dimension of each box (see iratio).</li>
</ul>
<dt><a name="tatraa"><b>tatraa (double precision)</b></a>
<ul>
<li> The target acceptance rate for the atom translation move. Must
be a value between 0.0 and 1.0. The maximum atom translational displacement
(rmtraa) is periodically adjusted (see iratio) to yield this acceptance
rate. I typically use a value of 0.5, though some researchers prefer
smaller values.</li>
</ul>
<hr>
<dt><a name="pmtracm"><b>pmtracm (double precision)</b></a>
<ul>
<li> Probability of performing a center-of-mass translation
move a molecule without regard to which box the molecule is currently
located in. This move chooses a molecule of the appropriate type
at random, chooses the x,y, or z direction at random, and then attempts
to displace the entire molecule a random distance between -rmtrac
and +rmtrac in that direction.</li>
</ul>
<dt><a name="pmtcmt"><b>pmtcmt (double precision)</b></a>
<ul>
<li> Probability of performing a center-of-mass translation
move on each type of molecule in the system.</li>
</ul>
<dt><a name="rmtrac"><b>rmtrac (double precision)</b></a>
<ul>
<li> The initial Center-of-mass translation maximum displacement (Angstroms)
for all molecule types in all boxes. As the simulation progresses,
these values will be updated independently to give the desired acceptance
rate for each molecule type in each dimension of each box (see iratio).</li>
</ul>
<dt><a name="tatrac"><b>tatrac (double precision)</b></a>
<ul>
<li> The target acceptance rate for the center-of-mass translation
move. Must be a value between 0.0 and 1.0. The maximum center-of-mass
translational displacement (rmtrac) is periodically adjusted (see
iratio) to yield this acceptance rate. I typically use a value of
0.5, though some researchers prefer smaller values.</li>
</ul>
<hr>
<dt><a name="pmrotate"><b>pmrotate (double precision)</b></a>
<ul>
<li> Probability of performing a rotation about the center-of-mass move for
a molecule without regard to which box the molecule is currently
located in.
This move chooses a molecule of the
appropriate type at random, chooses the x,y, or z direction at random,
and then attempts to rotate the entire molecule about an x,y, or
z axis that runs through the center-of-mass a random number of radians
between -<b>rmrot</b> and +<b>rmrot</b>.</li>
</ul>
<dt><a name="pmromt"><b>pmromt (double precision)</b></a>
<ul>
<li> Probability of performing a rotation move on each
type of molecule in the system.</li>
</ul>
<dt><a name="rmrot"><b>rmrot (double precision)</b></a>
<ul>
<li> The initial molecular rotation maximum displacement (radians)
for all molecule types in all boxes. As the simulation progresses,
these values will be updated independently to give the desired acceptance
rate for each molecule type about each axis of each box (see iratio).</li>
</ul>
<dt><a name="tarot"><b>tarot (double precision)</b></a>
<ul>
<li> The target acceptance rate for the rotation move. Must be a value
between 0.0 and 1.0. The rotation displacement (rmrot) is periodically
adjusted (see iratio) to yield this acceptance rate. I typically
use a value of 0.5, though some researchers prefer smaller values.</li>
</ul>
<hr>
<dt><a name="tor_cbstyle"><b>tor_cbstyle (integer)</b></a>
<ul>
<li> 0: When performing a <a href="../algorithm/cbmc.html">configurational-bias</a> move, generate trial dihedral angles
according to the true, ideal probability distribution. This is the method described in
<a href="../references.html#martin_siepmann_1999">Martin and Siepmann 1999</a></li>
<li> 1: When performing a <a href="../algorithm/cbmc.html">configurational-bias</a> move, generate trial dihedral angles
according to a different probability density and then
fix this up in the acceptance rules. This work is still in progress and is not yet published.</li>
</ul>
<dt><a name="bend_cbstyle"><b>bend_cbstyle (integer)</b></a>
<ul>
<li> 0: When performing a <a href="../algorithm/cbmc.html">configurational-bias</a> move, generate trial bending angles
according to the true, ideal probability distribution.
This is the method described in <a href="../references.html#martin_siepmann_1999">Martin and Siepmann 1999</a></li>
<li> 1: When performing a <a href="../algorithm/cbmc.html">configurational-bias</a> move, generate trial bending angles
according to a different probability density and then
fix this up in the acceptance rules. This work is still in progress and is not yet published.</li>
</ul>
<dt><a name="vib_cbstyle"><b>vib_cbstyle (integer)</b></a>
<ul>
<li> 0: When performing a <a href="../algorithm/cbmc.html">configurational-bias</a> move, generate trial bond lengths
according to the true, ideal probability distribution within the ranges set by the <b>vibrang</b> variable.</li>
<li> 1: When performing a <a href="../algorithm/cbmc.html">configurational-bias</a> move, generate trial bond lengths
according to a different probability density and then
fix this up in the acceptance rules. This work is still in progress and is not yet published.</li>
</ul>
<dt><a name="sdevtor"><b>sdevtor (double precision)</b></a>
<ul>
<li> This is the standard deviation of a gaussian distribution that is used to
sample dihedral angles [on 0,360] during a <a href="../algorithm/cbmc.html">configurational-bias</a>
regrowth for <b>tor_cbstyle</b> 1. Specify
a value in degrees. Right now I am using a value of 20.0. If <b>tor_cbstyle</b> is not 1 then
this value is not used in the code, but must still be entered into the input file.</li>
</ul>
<dt><a name="sdevbena"><b>sdevbena (double precision)</b></a>
<ul>
<li> This is the standard deviation of a gaussian distribution that is used to
generate trials for the part A bending angles [on 0,180] during a <a href="../algorithm/cbmc.html">configurational-bias</a>
regrowth for <b>bend_cbstyle</b> 1. Specify a value in degrees. Right now I am using a value of 10.0.
If <b>bend_cbstyle</b>
is not 1 then this value is not used in the code, but must still be entered into the input file.
</ul>
<dt><a name="sdevbenb"><b>sdevbenb (double precision)</b></a>
<ul>
<li> This is the standard deviation of a gaussian distribution that is used to
generate trials for the part B bending angles [on 0,360] during a <a href="../algorithm/cbmc.html">configurational-bias</a>
regrowth for <b>bend_cbstyle</b> 1. Specify a value in degrees. Right now I am using a value of 20.0.
If <b>bend_cbstyle</b>
is not 1 then this value is not used in the code, but must still be entered into the input file.
</ul>
<dt><a name="sdevvib"><b>sdevvib (double precision)</b></a>
<ul>
<li> This is the standard deviation of a gaussian distribution that is used to
sample bond lengths during a <a href="../algorithm/cbmc.html">configurational-bias</a> regrowth for <b>vib_cbstyle</b> 1. Specify
a value in Angtroms. Right now I am using a value of 0.1. If <b>vib_cbstyle</b> is not 1 then
this value is not used in the code, but must still be entered into the input file.
</ul>
<dt><a name="vibrang"><b>vibrang (double precision, double precision)</b></a>
<ul>
<li> This is the range of bond lengths to sample via <a href="../algorithm/cbmc.html">configurational-bias</a>
Monte Carlo. The range is expressed as a fraction of the equilibrium
bond length for the lower bound and the upper bound. I usually values
of 0.85 and 1.15 for <b>vib_cbstyle</b> 0, and the values have no meaning for <b>vib_cbstyle</b> 1.</li>
</ul>
<dt><a name="cdform"><b>cdform (integer)</b></a>
<ul>
<li> cdform = 0: When performing a <a href="../algorithm/cbmc.html">configurational-bias</a> move use the coupled-decoupled
formulation presented in the appendix of <a href="../references.html#martin_siepmann_1999">Martin and Siepmann 1999</a>.</li>
<li> cdform = 1: Uses a new algorithm that is not yet published. Still in early stages of testing this algorithm
to make sure it is producing correct answers.</li>
</ul>
<dt><a name="nch_nb_one"><b>nch_nb_one (integer)</b></a>
<ul>
<li> This is the number of trial positions that are sampled for the
first atom inserted during a <a href="../algorithm/cbmc.html">configurational-bias</a> molecule exchange
move (see pmswp and pmiswp). I typically use a value of 10. The
value must be less than or equal to NCHMAX (see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="nch_nb"><b>nch_nb (integer)</b></a>
<ul>
<li> This is the number of trial positions that are sampled for all
atoms except for the first atom inserted during a <a href="../algorithm/cbmc.html">configurational-bias</a>
molecule exchange move (see pmswp and pmiswp). This is used for
all atoms in a <a href="../algorithm/cbmc.html">configurational-bias</a> regrowth move. I
typically use a value of 10. The value must be less than or equal
to NCHMAX (see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="nch_tor_out"><b>nch_tor_out (integer)</b></a>
<ul>
<li> This is the number of outer loops over the dihedral angles that are sampled
during <a href="../algorithm/cbmc.html">configurational-bias</a> moves with <b>cdform</b> = 1. This has no meaning for
<b>cdform</b> = 0. I typically use a value in the range 1 to 10 with <b>cdform</b> = 1.
The value must be less than or equal to NCHTOR_MAX
(see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="nch_tor_in"><b>nch_tor_in (integer)</b></a>
<ul>
<li> This is the number of trial dihedral angles that are sampled
during <a href="../algorithm/cbmc.html">configurational-bias</a> moves. I typically use a value in the
range 100 to 360 for <b>tor_cbstyle</b> 0 and in the range 10 to 100 for <b>tor_cbstyle</b> 1.
The value must be less than or equal to NCHTOR_MAX
(see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="nch_tor_in_con"><b>nch_tor_in_con (integer)</b></a>
<ul>
<li> This is the number of trial dihedral angles that are sampled
during <a href="../algorithm/cbmc.html">configurational-bias</a> moves when we have grown the molecule
such that we need to connect back up with atoms that already exist.
This is needed in order to regrow cyclic molecules, and also could
be used to regrow the interiors of large molecules.
I typically use a value in the range 100 to 360. The value must
be less than or equal to NCHTOR_MAX (see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="nch_bend_a"><b>nch_bend_a (integer)</b></a>
<ul>
<li> This is the number of trial angles that are sampled during <a href="../algorithm/cbmc.html">configurational-bias</a>
moves when we are selecting the iugrow-iufrom-iuprev angle. I typically use a value of
1000 for <b>bend_cbstyle</b> of 0, and 10 for <b>bend_cbstyle</b> of 1.</li>
</ul>
<dt><a name="nch_bend_b"><b>nch_bend_b (integer)</b></a>
<ul>
<li> This is the number of trial angles that are sampled during <a href="../algorithm/cbmc.html">configurational-bias</a>
moves when we are selecting the rotation about a cone of one of
the iugrow angles relative to the others.
I typically use a value of 1000 for <b>bend_cbstyle</b> of 0, and 10 for <b>bend_cbstyle</b> of 1.</li>
</ul>
<dt><a name="nch_vib"><b>nch_vib (integer)</b></a>
<ul>
<li> This is the number of trial bond lengths that are sampled during
<a href="../algorithm/cbmc.html">configurational-bias</a> moves when we are growing atoms. I typically
use a value of 1000 for <b>vib_cbstyle</b> 0 and 10 for <b>vib_cbstyle</b> 1, unless I am using a
fixed-bond length force
field, in which case you might as well just use 1.</li>
</ul>
<li> The final section of towhee_input contains the information that
is used to construct the forcefield for the molecule types in the
system. There are three different ways to set up this section that
depend on inpstyle. The first way requires an explicit declaration
of the entire force field and while it allows the most flexibility,
is quite a burden to set up for any molecule with more than 10 atoms.
The second option is an automated builder which works for any molecule,
so long as the force field you wish to use has been implemented properly
into the code. The final option is a builder for proteins, or other
amino acid sequences. The choice of inpstyle determines what happens
below it so all of the remaining variables are described in each case
of inpstyle.</li>
<dt><a name="inpstyle"><b>inpstyle (integer)</b></a>
<ul>
<li> 0: Explicit declaration of all terms of the force field that describes
this molecule type. This is the hardest way to set up the molecule
information, but it also allows for the most flexibility. I do not
recommend using this option unless it is not possible to use the more
sophisticated options for the molecule of interest.</li>
<dt><a name="0_nunit"><b>nunit (integer)</b></a>
<ul>
<li> The number of atoms (or united-atoms) in this molecule. Must
be less than or equal to NUMAX (see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="0_nmaxcbmc"><b>nmaxcbmc (integer)</b></a>
<ul>
<li> The maximum number of atoms to regrow during a <a href="../algorithm/cbmc.html">configurational-bias</a>
regrowth. Generally I set this to the same value as <b>nunit</b>, but
occasionally a molecule is so large that almost none of the moves
which regrow large portions are accepted. In this case I set <b>nmaxcbmc</b>
to be around 1/2 of <b>nunit</b> so that only smaller regrowths are performed.</li>
</ul>
<dt><a name="0_lpdb"><b>lpdb (logical)</b></a>
<ul>
<li> .true. if you want to input information about the pdb (protein
data bank) atom name, amino acid number, amino acid name. This
information is needed in order to use the cartoon feature of certain
pdb viewers (such as Rasmol).</li>
<li> .false. if you don't need that feature. A pdb file will still
be output using the atom types defined in ffnonbond so you can
still view the system, you just can not use the cartoon feature.</li>
</ul>
<p><font color="red"> The variables listed immediately below (unit
through improper) are listed as a group for each atom in
the molecule. Thus, you input all of the information about the first
atom before you list information for the subsequent atoms.</font>
<dt><a name="0_unit"><b>unit (integer), ntype (integer), qqatom (double precision)</b></a>
<ul>
<li> unit is the number of the atom in order starting from atom
number 1. This is only used to help the user keep track of the
molecule as they are building it in the input file. If the unit
number listed in towhee_input does not match the running total
of unit numbers in Towhee then the code will stop with an error
message.</li>
<li> ntype is the integer that ffnonbond uses to determine the force
field parameters. You will need to look in ffnonbond.F in order
to find the number that you wish to use.</li>
<li> qqatom is the charge on this atom. This value is not used (but
must still be input) if lcoulomb is false.</li>
</ul><dt><a name="0_pdbname"><b>pdbname (character), aminonum (integer), aminoshort (character)</b></a>
<ul><dt>
These variables only need to be listed if lpdb = .true.
<li> pdbname: A four letter/number string that is output in the pdb
file. The precise spacing is important if you want to get most pdb
viewers to recognize the atoms as the pdb file is extremely specific.</li>
<li> aminonum: The number of each amino acid starting from the N-terminus.</li>
<li> aminoshort: The three letter code for each amino acid, or other
group (such as caps on the C or N termini).</li>
</ul>
<dt><b>vibration</b>
<ul><dt><a name="0_invib"><b>invib (integer)</b></a>
<ul>
<li> The first line under the vibration heading is the number of
atoms that are bonded to the current atom. Must be a number between
0 and NNBOND (see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="0_ijvib"><b>ijvib (integer), itvib (integer)</b></a>
<ul>
<li> The next <b>invib</b> lines underneath the vibration heading
are a list of the bond partner and bond force field number (see
ffbond.F) for the <b>invib</b> atoms that are bonded to the
current atom. Thus if you have 4 vibrations the next 4 lines
will list the bond partner and bond force field number for each
bond.</li>
</ul>
</ul><dt>
<b>bending</b>
<ul>
<dt><a name="0_inben"><b>inben (integer)</b></a>
<ul>
<li> The first line under the bending heading is the number of bond
bending angles that terminate at the current atom. You must list
all bond bending angles which have the current atom at one of
the ends, but you do not list bond bending angles which contain
the current atom in the center. Must be a number between 0 and
NNBOND*(NNBOND-1) (see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="0_ijben2"><b>ijben2 (integer), ijben3(integer), itben (integer)</b></a>
<ul>
<li> The next <b>inben</b> lines underneath the bending heading
are a list of the other atoms in the bond bending angle and
the bond bending angle force field number (see ffangle.F) for
the <b>inben</b> angles that contain the current atom at one
of the ends. The format for listing the angle is to consider
the current atom in the first position of an angle between atoms
current-ijben2-ijben3 so you only need to list atoms ijben2
and ijben3 and the bending type on each line, where there is
one angle per line of towhee_input.</li>
</ul>
</ul><dt>
<b>torsion</b>
<ul>
<dt><a name="0_intor"><b>intor (integer)</b></a>
<ul>
<li> The first line under the torsion heading is the number of dihedral
angles (regular torsions) that terminate at the current atom.
You must list all regular torsion angles which have the current
atom at one of the ends, but you do not list regular torsion angles
which contain the current atom in the center. Must be a number
between 0 and NNBOND*(NNBOND-1)*(NNBOND-1) (see <a href="../code/code_manual.html#preproc">preproc.h</a>).</li>
</ul>
<dt><a name="0_ijtor2"><b>ijtor2 (integer), ijtor3(integer), ijtor4 (integer), ittor (integer)</b></a>
<ul>
<li> The next <b>intor</b> lines underneath the torsion heading
are a list of the other atoms in the regular torsion angle and
the torsion force field number (see fftorsion.F) for the <b>intor</b>
torsions that contain the current atom at one of the ends. The
format for listing the regular torsion angle is to consider
the current atom in the first position of a dihedral between
atoms current-ijtor2-ijtor3-ijtor4 so you only need to list
atoms ijtor2, ijtor3, and ijtor4 and the torsion type on each
line, where there is one regular torsion per line of towhee_input.</li>
</ul>
</ul><dt>
<b>angle-angle</b>
<ul><dt><a name="0_inaa"><b>inaa (integer)</b></a>
<ul>
<li> The first line under the angle-angle heading is the number
of angle-angle terms which have their shared central atom located
at the current atom. You must list all angle-angle terms which
have the current atom at the shared central position, but you
do not list angle-angle terms which contain the current atom at
one of the ends. Must be a number between 0 and (NNBOND*(NNBOND-1))/2
(see <a href="../code/code_manual.html#preproc">preproc.h</a>). </li>
</ul>
<dt><a name="0_ijaa0"><b>ijaa0 (integer), ijaa1(integer), ijaa2 (integer), itaa (integer)</b></a>
<ul>
<li> The next <b>inaa</b> lines underneath the angle-angle heading
are a list of the other atoms in the angle-angle term and and
the angle-angle force field number (see ffangang.F) for the
<b>inaa</b> angle-angle terms that contain the current atom
at the shared central atom. The format for listing the regular
torsion angle is to consider the current atom as the central
shared atom in the angle-angle term between the angles ijaa0-current-ijaa1
and ijaa0-current-ijaa2. Each angle-angle term is listed on
one line according to the format ijaa0, ijaa1, ijaa2 and itaa
(the angle-angle type). </li>
</ul>
</ul><dt>
<b>improper torsion</b>
<ul><dt><a name="0_inimp"><b>inimprop (integer)</b></a>
<ul>
<li> The first line under the improper heading is the number
of improper torsions (any form) which have the central atom
located at the current atom. You must list all improper torsions
which have the current atom at the central position,
but you do not list improper torsions which contain
the current atom at one of the ends. Must be a number between
0 and ((NNBOND-1)*(NNBOND-2))/2 (see <a href="../code/code_manual.html#preproc">preproc.h</a>). </li>
</ul>
<dt><a name="0_ijimprop2"><b>ijimprop2 (integer), ijimprop3(integer), ijimprop4 (integer), itimprop (integer)</b></a>
<ul>
<li> The next <b>inimprop</b> lines underneath the improper torsion
heading are a list of the other atoms in the
improper torsion and the improper force field number (see
ffimproper.F) for the <b>inimprop</b> improper torsions that
contain the current atom at the central atom. There are currently three different
forms of improper torsions, and these forms are specified in the force field (either
in ffimproper.F or in towhee_ff). In all cases, the three atoms must all be bound to the
current atom.
<ul>
<li>Form 1: Amber / Stereo Improper.
This type of improper
torsion works by computing the dihedral angle between the planes
formed by the atoms ijimprop2-current-ijimprop3 and ijimprop2-current-ijimprop4.
This was originally used as a means to enforce planarity, but
is also useful for enforcing stereochemistry. The format for
the stereo improper torsion is to list atoms ijimprop2, ijimprop3,
and ijimprop4 and the torsion type (itimprop) on each line, where
there is one improper torsion per line of towhee_input.</li>
<li>Form 2: Charmm / Simple Out-of-Plane.
The format for listing the simple out-of-plane improper torsion is to consider
the current atom in the first position of a dihedral between
atoms current-ijimprop2-ijimprop3-ijimprop4 so you only need to list
atoms ijimprop2, ijimprop3, and ijimprop4 and the torsion type on
each line, where there is one improper torsion per line
of towhee_input. </li>
<li>Form 3: Averaged Out-of-Plane.
The format for listing the averaged out-of-plane improper torsion is identical to the
simple out-of-plane. The only difference is that the dihedral used is not just the
current-ijimprop2-ijimprop3-ijimprop4, but is instead an average of the following three
angles: current-ijimprop2-ijimprop3-ijimprop4; current-ijimprop3-ijimprop4-ijimprop2;
current-ijimprop4-ijimprop2-ijimprop3.</li>
</ul>
</ul>
</ul><dt>
<dt><font color="red"> This is the end of the section that is repeated
for each atom for inpstyle=0 </font>
<p> </p>
<li> 1: Polypeptide builder which sets up everything using buildprot
and assemble routines. This sets up everything for proteins from
a list of the amino acids. This routine also could be used to build
polymers, but I have not yet coded in any non-peptide monomers.
I am working on some documentation of the currently available monomers,
but until that is complete you will have to look into buildprot.F
to see what is available. </li>
<dt><a name="1_nunit"><b>nunit (integer)</b></a>
<ul>
<li> The number of monomers in this molecule. Must be less than
or equal to NUMAX (see <a href="../code/code_manual.html#preproc">preproc.h</a>) and there is also a check in
the builder to make sure you do not exceed NUMAX with the number
of atoms that all of your monomers require. </li>
</ul>
<dt><a name="1_nmaxcbmc"><b>nmaxcbmc (integer)</b></a>
<ul>
<li> The nmaxcbmc variable determines the largest number of atoms
that are regrown during a <a href="../algorithm/cbmc.html">configurational-bias</a> regrowth. Since
we only know the number of amino acids at this stage, and not
the number of atoms, the final value of <b>nmaxcbmc</b> is computed as
the product of the final number of atoms in the molecule times
the input value of <b>nmaxcbmc</b> divided by the input value of <b>nunit</b>. </li>
</ul>
<dt><a name="1_forcefield"><b>forcefield (character)</b></a>
<ul>
<li> The forcefield that you want to use to build this molecule.
I'm working on better documentation of which force fields are
implemented, but until that is finished just take a look in ffnonbond.F </li>
</ul>
<dt><a name="1_protgeom"><b>protgeom (character)</b></a>
<ul>
<li> linear: The polypeptide will be generated starting from the
N-terminus and finishing at the C-terminus. </li>
<li> cyclic: The polypeptide is cyclic so the first amino acid is
bonded to the final amino acid. </li>
</ul>
<p><font color="red"> The variables listed immediately below are listed
as a group for each amino acid in the molecule. There is only one
line of input per building block in the molecule.</font>
<dt><a name="1_pepname"><b>pepname (character), stereochem (character), bondpartner (integer)</b></a>
<ul>
<li> pepname is a two letter code for the building blocks of the
molecule. Peptides use a code where the first letter is the one
letter code for each amino acid. The second character is generally
the charge state of the side chain, but has some different meanings
for histidine and cysteine. Until the documentation is in place
you need to look in buildprot.F to find a complete list of valid
pepnames. </li>
<li> stereochem is a single letter which determines the stereochemistry
of the peptide. The valid options are l (L-amino acid), d (D-amino
acid) and r (racemic). </li>
<li> bondpartner is the number of the amino acid which is bonded
to this amino acid through the side chain. This is only used for
cysteines which have a disulfide bridge (type cs), otherwise just
set this value to 0. The bondpartner integers refer to the amino
acids starting with the first amino acid listed in towhee_input
as amino acid #1. </li>
</ul>
<dt><font color="red"> This is the end of the section that is repeated
for each atom for inpstyle=1 </font>
<p> </p>
<li> 2: Atom-based connectivity map, with the details of the force
field parameters determined via the buildmolec and assemble routines.
This is the easiest way to set up any molecule that is not a series
of amino acids. You must choose a forcefield and then input the
appropriate parameter names for that force field. I'm working on
better documentation of the atom names used in Towhee, but until
that is complete you have to look in ffnonbond.F The program will
sort out all of the vibration types, bending angles, bending types,
regular torsion angles, regular torsion types, angle-angle terms,
angle-angle types, and improper torsion types that are implied by the bonding structure of the
molecule. I have not automated the atom assignment for the improper
torsions as some of the force fields do not have general rules (at
least not rules I can decipher) for determining where the improper
torsions belong. </li>
<dt><a name="2_nunit"><b>nunit (integer)</b></a>
<ul>
<li> The number of atoms (or united-atoms) in this molecule. Must
be less than or equal to NUMAX (see <a href="../code/code_manual.html#preproc">preproc.h</a>). </li>
</ul>
<dt><a name="2_nmaxcbmc"><b>nmaxcbmc (integer)</b></a>
<ul>
<li> The maximum number of atoms to regrow during a <a href="../algorithm/cbmc.html">configurational-bias</a>
regrowth. Generally I set this to the same value as nunit, but
occasionally a molecule is so large that almost none of the moves
which regrow large portions are accepted. In this case I set <b>nmaxcbmc</b>
to be around 1/2 of <b>nunit</b> so that only smaller regrowths are performed. </li>
</ul>
<dt><a name="2_forcefield"><b>forcefield (character)</b></a>
<ul>
<li> The forcefield that you want to use to build this molecule.
I'm working on better documentation of which force fields are
implemented, but until that is finished just take a look in ffnonbond.F </li>
</ul>
<p><font color="red"> The variables listed immediately below (unit
through improper) are listed as a group for each atom in
the molecule. Thus, you input all of the information about the first
atom before you list information for the subsequent atoms.</font>
<dt><a name="2_unit"><b>unit (integer), type (character), qqatom (double precision)</b></a>
<ul>
<li> unit is the number of the atom in order starting from atom
number 1. This is only used to help the user keep track of the
molecule as they are building it in the input file. If the unit
number listed in towhee_input does not match the running total
of unit numbers in Towhee then the code will stop with an error
message. </li>
<li> ntype is the character string that contains the atom type for
the forcefield. I'm working on better documentation for the available
atom types, but until that is ready you will need to look in ffnonbond.F.
I generally use the same atom types as the original authors of
the force field, but sometimes I have to add information to those
atom types in order to automate the assembly process. </li>
<li> qqatom is the charge on this atom. This value is not used (but
must still be input) if lcoulomb is false. </li>
</ul><dt>
<b>vibration</b>
<ul>
<dt><a name="2_invib"><b>invib (integer)</b></a>
<ul>
<li> The first line under the vibration heading is the number of
atoms that are bonded to the current atom. Must be a number between
0 and NNBOND (see <a href="../code/code_manual.html#preproc">preproc.h</a>). </li>
</ul>
<dt><a name="2_ijvib"><b>ijvib (integer)</b></a>
<ul>
<li> The second line contains the atom numbers for all <b>invib</b>
atoms that are bonded to the current atom. </li>
</ul>
</ul><dt>
<b>improper torsion</b>
<ul><dt><a name="2_inimp"><b>inimprop (integer)</b></a>
<ul>
<li> The first line under the improper heading is the number
of improper torsions (any form) which have the central atom
located at the current atom. You must list all improper torsions
which have the current atom at the central position,
but you do not list improper torsions which contain
the current atom at one of the ends. Must be a number between
0 and ((NNBOND-1)*(NNBOND-2))/2 (see <a href="../code/code_manual.html#preproc">preproc.h</a>). </li>
</ul>
<dt><a name="2_ijimprop2"><b>ijimprop2 (integer), ijimprop3(integer), ijimprop4 (integer), itimprop (integer)</b></a>
<ul>
<li> The next <b>inimprop</b> lines underneath the improper torsion
heading are a list of the other atoms in the
improper torsion and the improper force field number (see
ffimproper.F) for the <b>inimprop</b> improper torsions that
contain the current atom at the central atom. There are currently three different
forms of improper torsions, and these forms are specified in the force field (either
in ffimproper.F or in towhee_ff). In all cases, the three atoms must all be bound to the
current atom.
If you want the program to determine the improper
torsion type then just enter an itimprop value of 0.</li>
<ul>
<li>Form 1: Amber / Stereo Improper.
This type of improper
torsion works by computing the dihedral angle between the planes
formed by the atoms ijimprop2-current-ijimprop3 and ijimprop2-current-ijimprop4.
This was originally used as a means to enforce planarity, but
is also useful for enforcing stereochemistry. The format for
the stereo improper torsion is to list atoms ijimprop2, ijimprop3,
and ijimprop4 and the torsion type (itimprop) on each line, where
there is one improper torsion per line of towhee_input.</li>
<li>Form 2: Charmm / Simple Out-of-Plane.
The format for listing the simple out-of-plane improper torsion is to consider
the current atom in the first position of a dihedral between
atoms current-ijimprop2-ijimprop3-ijimprop4 so you only need to list
atoms ijimprop2, ijimprop3, and ijimprop4 and the torsion type on
each line, where there is one improper torsion per line
of towhee_input. </li>
<li>Form 3: Averaged Out-of-Plane.
The format for listing the averaged out-of-plane improper torsion is identical to the
simple out-of-plane. The only difference is that the dihedral used is not just the
current-ijimprop2-ijimprop3-ijimprop4, but is instead an average of the following three
angles: current-ijimprop2-ijimprop3-ijimprop4; current-ijimprop3-ijimprop4-ijimprop2;
current-ijimprop4-ijimprop2-ijimprop3.</li>
</ul>
</ul>
</ul><dt>
<dt><font color="red"> This is the end of the section that is repeated
for each atom for inpstyle=2 </font>
</ul>
</ul>
<a href="../index.html">Return to the main towhee web page</a>
<p> </p>
</td>
</tr>
</table>
<hr width="715" align="left">
<i><font size="2">Send comments to:</font></i> <font size="2"> <a href="mailto:marcus_martin@users.sourceforge.net">Marcus
G. Martin</a><br>
<i>Last updated:</i>
<!-- #BeginDate format:Am1 -->September 15, 2004<!-- #EndDate -->
</font> <br>
</body>
</html>
