\documentclass[12pt]{article}
\usepackage{url,setspace,amsmath}
\usepackage{color}
%\setlength{\oddsidemargin}{-8mm}
%\setlength{\evensidemargin}{0mm}
%\setlength{\textwidth}{175mm}
%\setlength{\topmargin}{-5mm}
%\setlength{\textheight}{225mm}
%\setlength{\headheight}{0cm}
\setstretch{1}

\begin{document}
\title{\vspace{-2cm}Documentation for shearableWLC code. Simulations using the dssWLC model.}
\author{E.~F.~Koslover, A.~J.~Spakowitz}
\date{Last updated \today}
\maketitle

This code can be used to run Brownian Dynamics or Monte Carlo simulations of the dssWLC model.

\tableofcontents
\newpage

\section{Compilation Instructions}
To compile and run the program, you will need the following:
\begin{itemize}
\item a compiler capable of handling Fortran90.
The code has been tested with the gfortran compiler. 
\item BLAS and LAPACK libraries installed in a place where the compiler knows to look for them
\item Python (version 2.5 or higher) to run various auxiliary scripts
  (e.g., for visualization). The scripts have been tested with Python
  2.6.4 only. You will also need the NumPy extension package.
\item Recommended: PyMOL to visualize pdb files.
\end{itemize}

The code has been tested on Ubuntu Linux only.

To compile with gfortran, go into the \path=source= directory. Type \verb=make=.
To compile with any other compiler that can handle Fortran90, type
\begin{verbatim}
make FC=compiler
\end{verbatim}
substituting in the command you usually use to call the compiler. 

If the compilation works properly, the executable \path=shearableWLC.exe= will appear in the main directory.

%% To test that the code is running properly, go into the \path=testing= directory and type 
%% \begin{verbatim}
%% ./runalltests.py
%% \end{verbatim}
%% This will run a number of test calculations to make sure everything works properly. The tests will take a few minutes to complete.

\section{Usage Instructions}
To run the program in the main directory, type:
\begin{verbatim}
./shearableWLC.exe suffix
\end{verbatim}

Here, \verb=suffix= can be any string up to 100 characters in length. 
The program reads in all input information from a file named
\path=param.suffix= where, again, \verb=suffix= is the command-line
argument. If no argument is supplied, it will look for a file named
\path=param=. If the desired parameter file does not exist, the
program will exit with an error.

The parameters in the input file are given in the format "{\em KEYWORD} value" where the possible keywords and values are described
in Section \ref{sec:keywords}. Each keyword goes on a separate
line. Any line that starts with "\#" is treated as a comment and
ignored. Any blank line is also ignored. The keywords in the parameter
file are not case sensitive. For the most part, the order in which the
keywords are given does not matter. All parameters have default
values, so you need only specify keywords and values when you want to
change something from the default.

%Instructions for running specific calculations are given in more
detail in Section \ref{sec:tasks}. 
Example parameter files for the
different calculations are provided in the \path=examples= directory.

\section{Examples for a Quick Start}
Example parameter files are located in the \path=examples= directory

\subsection{Example 1: equilibrium sampling}
Sample a number of chain configurations from the equilibrium distribution. Run this example by typing
\begin{verbatim}
../shearableWLC.exe ex.equildistrib
\end{verbatim}
This example uses a chain where the 2 edge segments on either side have a segment length of $0.1$ and corresponding parameters, whereas all other segments have a segment length of $0.2$.
This will sample 10000 chains in total, and output end-to-end vectors in the first three columns of the output file ex.equildistrib.out. The next 3 columns give the first orientation vector.  Snapshots of every $1000^\text{th}$ chain will be output in the \path=ex.equildistrib.snap.out= file, which can then be converted to pdb format as follows:
\begin{verbatim}
../snapshot2pdb.py ex.equildistrib.snap.out
\end{verbatim}
The resulting pdb file (\path=ex.equildistrib.snap.pdb=) can be loaded into PyMOL and colored using the following commands within PyMOL itself.
\begin{verbatim}
load ex.equildistrib.snap.pdb, multiplex=0
../scripts/viewsnapshots.pml
\end{verbatim}

\subsection{Example 2: Brownian Dynamics}
Run a Brownian Dynamics simulation for 2000 chains in parallel, keeping track of shear stress correlations in the file \path=ex.browndyn.stress.out=.
The simulation starts from an equilibrium distribution and runs for 10000 steps. 
\subsection{Example 3: Looping first-passage time}
Run a Brownian Dynamics simulation to track the first passage time to looping for 1000 chains in parallel. This example uses a chain where the 2 edge segments on either side have a segment length of $0.1$ and corresponding parameters, whereas all other segments have a segment length of $0.2$. The friction coefficients are supplied per unit length, which means the edge beads have half the friction of the 2-nd beads from the edge. Those in turn have half the friction of the inner beads (which correspond to longer segment lengths). The time to looping for each chain is output in the 3rd column of the file \path=ex.loop.loop.out= as each chain loops.

\subsection{Example 4: Looping first-passage time, with dynamic rediscretization}
Run a Brownian Dynamics simulation to track the first passage time to looping for 1000 chains in parallel. Rediscretize the whole chain dynamically by a factor of 4 (to segment lengths of 0.05) whenever the end-to-end distance goes below 0.4. Coarsen the chain whenever the end-to-end distance goes above 0.5. Parameters are interpolated from the file \path=dssWLCparams.txt=.

\subsection{Example 5: Monte Carlo}
Run a Monte Carlo simulation of $10^8$ steps for a dssWLC. Output the end-to-end distance every 1000 steps into columns 5-7 of the file \path=ex.montecarlo.out=.

\section{Auxiliary Scripts}

\subsection{Visualizing structures}

The script \path=scripts/snapshot2pdb.py= will convert a snapshot file containing many chain configurations into a concatenated pdb file that can be loaded into pymol. Run this script without any arguments to get usage information. 
To load the resulting snapshots into different states to make a movie do ``load snap.pdb, multiplex=0'' in pymol. Color an visualize in pymol using  \path=scripts/viewsnapshots.pml=.

\section{Parameters for the dssWLC model}
\label{sec:getparams}
The Matlab code for calculating the energetic and dynamic parameters for the dssWLC model is located in \path=getparams/=. The relevant function for getting the parameters is \path=dssWLCcalcparams=

There are three approaches to getting the parameters. Firstly, as described in Koslover and Spakowitz, Soft Matter, 2013, one can optimize of the parameter $\alpha = \eta^2\epsilon_b/\epsilon_\perp$ to get the minimal length scale of accuracy. This can be done for a chain with Nseg segments of length $\Delta = \text{del}$ by running the following in Matlab:

\begin{verbatim}
[eb,gam,epar,eperp,eta,alpha,zetau,delt] = dssWLCcalcparams(del,Nseg)
\end{verbatim}

These calculations can be very slow (taking on the order of 15 min). The first time this is run will be even slower as it needs to generate a tensor of coupling coefficients for spherical harmonics which will then by saved in a .mat file for future use. The size of the matrix used for calculating the structure factor is set with the optional parameter LMAX. Default value is LMAX=10. Smaller values of $\Delta$ require higher values of LMAX. 

The output parameters are, in order, $\epsilon_b, \gamma, \epsilon_\parallel, \epsilon_\perp, \eta, \alpha, \zeta_{ub}, \delta t$. They should be input into the parameter files for the Fortran simulation code as follows:
\begin{eqnarray*}
\text{EB} & \quad & \epsilon_b \\
\text{GAM} & \quad & \gamma \\
\text{EPERP} & \quad & \epsilon_\perp + \eta^2\epsilon_b \\
\text{EPAR} & \quad & \epsilon_\parallel \\
\text{EC} & \quad & - \eta \epsilon_b \\
\text{FRICT} & \quad & 1/(\text{Nseg}+1) \quad \zeta_{ub} \\
\text{DELTSCL} & \quad & \delta t/\zeta_{ub}
\end{eqnarray*}
(assuming that  $\zeta_{ub} < 1/(\text{Nseg}+1)$). 

Alternately, one can input the friction coefficients per unit length, which needs to be done in the case where not all segment lengths are equal. If using this format, then even if all segment lengths are equal the edge beads will have have the translational friction of the central beads.
\begin{equation*}
\text{FRICT} \quad 1D0 \quad \zeta_{ub}(\text{Nseg}+1)/(\Delta \text{Nseg}) \quad T
\end{equation*}

Generally, if running a chain with different segment lengths, one should use the smallest segment length to determine the appropriate values of $\delta t$ and $\zeta_{ub}$. 

Another approach to selecting the parameters is to use a specific value of $\alpha$, which can be done as follows in Matlab for a value $\alpha=a$:
\begin{verbatim}
[eb,gam,epar,eperp,eta,alpha,zetau,delt] = dssWLCcalcparams(del,Nseg,`alpha',a)
\end{verbatim}
This will run much faster and is generally more stable. 

Finally, the recommended approach is to simply interpolate from the pretabulated values available in \path=getparams/dssWLCparams.txt=. This file lists on each line the values of $\Delta, \epsilon_b, \gamma, \epsilon_\parallel, \epsilon_\perp, \eta, \zeta_u, \delta t/\zeta_u$, respectively. The interpolation can be done via:
\begin{verbatim}
[eb,gam,epar,eperp,eta,alpha,zetau,delt] = dssWLCcalcparams(del,Nseg,`intepfile',`dssWLCparams.txt')
\end{verbatim}
The tabulated file was generated using optimization over $\alpha$ for $\Delta = 0.01 - 1$, and a constant value of $\alpha$ thereafter for $\Delta = 1 - 4$. The interpolation file was generated using \path=getparams/tabulateparams.m=
%The matrix size used was $\text{LMAX}=14$ for $\Delta < 0.045$ and $\text{LMAX}=10$ otherwise. 

\section{Description of Specific Calculations}
\label{sec:tasks}

The {\em ACTION} keyword specifies what type of calculation will be
done. The possible actions are EQUILDISTRIB, BROWNDYN, and MONTECARLO, as described below.

\subsection{EQUILDISTRIB action}

This generates a bunch of dssWLC chain configurations sampled from an equilibrium distributions. The sampling method is set by the {\em STARTEQUIL} keyword (2 types of rejection sampling or monte carlo). Number of chains sampled is set by {\em MCSTEPS} keyword. Dumping of snapshots is set by the {\em SNAPSHOTS} keyword. End-to-end vectors for all configurations are output into the file set by {\em OUTFILE}. Will also work with a Gaussian chain (see {\em GAUSSIANCHAIN} keyword) and a bead-rod chain (if {\em STRETCHABLE} and {\em SHEARABLE} are set to false).

\subsection{BROWNDYN action}

This runs a Brownian Dynamics simulation for a set of chains, keeping track of the shear stress correlation over time. Use {\em NCHAIN} to set the number of chains being run in parallel, {\em BDSTEPS} to set the total number of steps and the printing / output frequency. Use {\em LOOPING} keyword to tabulate first passage times for looping of chain ends. {\em FRICT} keyword sets friction coefficients. {\em DELTSCL} keyword sets the timestep as a multiple of the friction coefficients. Can use {\em STARTEQUIL} keyword to start from an equilibrated set of configurations. Can periodically dump out snapshots of chain configurations.

\subsection{MONTECARLO action}

Run a Monte Carlo simulation for a dssWLC chain. Only Monte Carlo of a single chain at a time has been recently tested. Set total number of steps and number of initialization steps with the {\em MCSTEPS} keyword. Set printint / output frequency with {\em MCPRINTFREQ}. Can output snapshots with {\em SNAPSHOTS} keyword. Used {\em ADJUSTRANGE} keyword to set how often step sizes are adjusted. For the most part, this has been supplanted by the equilibrium configuration sampling (EQUILDISTRIB action). However this general procedure can be implemented with additional complications (ie: chain meshes, non-local interactions) which EQUILDISTRIB cannot.
% ---------------------------------------------------------

%----------------------------------------------------------
\section{Keyword Index}
\label{sec:keywords}
The code will attempt to read parameters out of a file named \path=param.suffix= where ``suffix'' is the command line argument. If no command line arguments are supplied, it will look for a file named \path=param=. If multiple arguments are supplied, it will read multiple parameter files in sequence.

The parameter file should have one keyword per line and must end with a blank line. All blank lines and all lines beginning with \# are ignored. For the most part, the order of the lines and the capitalization of the keywords does not matter. All keywords except {\em ACTION} are optional. The default values for each parameter are listed below. If a keyword is supplied, then values may or may not be needed as well. Again, the required and optional value types are listed below. 

Keywords and multiple values are separated by spaces. 

When reading the parameter file, lines longer than 500 characters will be truncated. To continue onto the next line, add ``+++'' at the end of the line to be continued.
No individual keyword or  value should be longer than 100 characters.

Floating point numbers can be formated as $1.0$, $1.1D0$, $10e-1$, $-1.0E+01$, etc., where the exponential notation specifier must be D or E (case insensitive). Integer numbers can also be specified in exponential notation without decimal points (eg: 1000 or 1E3). Logical values can be specified as T, F, TRUE, FALSE, 1, or 0 (with 1 corresponding to true and 0 to false).

By default, all energy units are in kT. 

\begin{itemize}
%
\item {\it ACTION}
  \begin{itemize}
    \item  value: 1 string of at most 20 characters; no default
    \item This keyword sets the overall calculation performed by the program
 (see Sec.\ref{sec:tasks})
    \item Possible values are: MONTECARLO, BROWNDYN, EQUILDISTRIB
  \end{itemize}
%
\item {\it ADJUSTRANGE}
  \begin{itemize}
    \item  value: 1 required integer (ADJUSTEVERY), 3 optional floats (FACCTARGET, FACCTOL, ADJUSTSCL)
    \item When doing a Monte Carlo simulation, how to adjust the step size.
    \item The accepted fraction is checked every ADJUSTEVERY steps. If it is outside the range of FACCTARGET $\pm$ FACCTOL, then the step sizes are multiplied or divided by ADJUSTSCL
    \item defaults are ADJUSTEVERY=1000, FACCTARGET=0.5, FACCTOL=0.1, ADJUSTSCL=2
  \end{itemize}
%
\item {\it BDSTEPS}
  \begin{itemize}
    \item  value: 1 required integer (BDSTEPS), 1 optional float (BDPRINTEVERY), 1 optional logical (BDPRINTLOG)
    \item Sets the total number of Brownian Dynamics steps (BDSTEPS) and how often to print output.
    \item If BDPRINTLOG is true then print at logarithmically spaced step numbers, where BDPRINTEVERY sets the multiplicative factor for the spacing. Otherwise, print every BDPRINTEVERY steps.
    \item dafaults: BDSTEPS=1000, BDPRINTEVERY=1, BDPRINTLOG=false
  \end{itemize}
%
\item {\it BRCRELAX}
  \begin{itemize}
    \item  value: 1 float
    \item parameter for extra force to keep segment lengths fixed when running Brownian Dynamics with a bead-rod model
    \item {\color{red} Bead-rod BD are not debugged! do not use.}
  \end{itemize}
%
\item {\it CONNECT}
  \begin{itemize}
    \item  value: 4 integers
    \item Connection point in a mesh of many chains
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
  \end{itemize}
%
\item {\it CONNECTMOD}
  \begin{itemize}
    \item  value: 2 floats
    \item something about connecting chains in a mesh...
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
  \end{itemize}
%
\item {\it CONNECTTYPE}
  \begin{itemize}
    \item  value: 2 logicals
    \item whether to connect together bead positions and/or orientations
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
  \end{itemize}
%
\item {\it CONSTMOD}
  \begin{itemize}
    \item  value: 1 float
    \item not currently used
  \end{itemize}
%
\item {\it COUPLED}
  \begin{itemize}
    \item  value: 1 logical
    \item not currently used
  \end{itemize}
%
\item {\it DELTSCL}
  \begin{itemize}
    \item  value: 1 float
    \item scaling constant used to set the timestep in Brownian Dynamics simulations
    \item Ultimately, if the chain is designated as {\it SHEARABLE} then the timestep is $\delta t = \text{DELTSCL} * \text{min}(\zeta_{ub},\zeta_{rb})$, where the friction coefficients $\zeta_{ub},\zeta_{rb}$ are set using keyword {\it FRICT}. If the chain is not shearable, then $\delta t = \text{DELTSCL} * \zeta_{rb}$
  \end{itemize}
%
\item {\it DIAMONDLATTICE}
  \begin{itemize}
    \item  value: 3 integers; 1 optional float
    \item connect up a mesh of chains in a diamond lattice
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
  \end{itemize}
%
\item {\it DOLOCALMOVES}
  \begin{itemize}
    \item  value: no values
    \item When running a Monte Carlo calculation for a single chain, move individual beads rather than the default crank-shaft type moves
  \end{itemize}
%
\item {\it DYNAMICREDISC}
  \begin{itemize}
    \item  values: 2 floats (REDISCCUTFINE, REDISCCUTCOARSE), 1 integer (REDISCFACT), 1 optional integer (REDISCEDGE), 1 optional float (REDISCEQTIME)
    \item When running a first-passage looping time calculation, rediscretize the chain dynamically depending on how close the two ends are
    \item Switch to a finer discretization when the end-to-end distance is smaller than REDISCCUTFINE. Switch back to a coarse discretization when the end-to-end distance is larger than REDISCCUTCOARSE. 
    \item Rediscretization is by an integer factor of REDISCFACT. Each segment turns into REDISCFACT shorter segments
    \item REDISCEDGE sets how many of the edge segments are rediscretized (make it bigger than half chain length to rediscretize the whole chain). By default, REDISCEDGE=1, so only the first and last segment are rediscretized.
    \item REDISCEQTIME gives the total time for running a mini brownian dynamics simulation to equilibrate newly inserted beads. Default value is $\zeta_r\Delta^4/(3\pi/2)^4$, the timescale for bending equilibration of a single coarse segment.
    \item 
  \end{itemize}
%
\item {\it EC}
  \begin{itemize}
    \item  value: 1 float; default: 0
    \item The bend-shear coupling energetic parameter for the dssWLC
    \item EC = $-\eta\epsilon_b$ using the notation in the Soft Matter paper
  \end{itemize}
%
\item {\it EDGESEGS}
  \begin{itemize}
    \item  value: 1 integer (N); 6 floats
    \item set separate parameters for the first and last N segments
    \item parameters, in order, are the ones usually set with the {\em LS, LP, GAM, EPAR, EPERP, EC} keywords
  \end{itemize}
%
\item {\it EPAR}
  \begin{itemize}
    \item  value: 1 float; default: 1D3
    \item Stretch modulus for the dssWLC energetics
    \item EPAR = $\epsilon_\parallel$ using the notation in the Soft Matter paper
  \end{itemize}
%
\item {\it EPERP}
  \begin{itemize}
    \item  value: 1 float; default: 1D3
    \item Modified shear modulus for the dssWLC energetics
    \item EPERP = $\hat{\epsilon}_\perp = \epsilon_\perp + \eta^2\epsilon_b$ using the notation in the Soft Matter paper
  \end{itemize}
%
\item {\it EPAR}
  \begin{itemize}
    \item  value: 1 float; default: 1D3
    \item Stretch modulus for the dssWLC energetics
    \item EPAR = $\epsilon_\parallel$ using the notation in the Soft Matter paper
  \end{itemize}
%
\item {\it FINITEXT}
  \begin{itemize}
    \item  value: 1 optional float; default: 1D-3
    \item For Monte Carlo simulations, prevent individual segments from stretching beyond the contour length.
    \item optional float is a scaling factor F where the stretch energy is a fraction (1-F) of the usual gaussian  and a fraction F of a logarithmic term that prevents the overextension of the segment (as in the FENE model)
    \item {\color{red} finite extension not tested for a while. Use at own risk.}
  \end{itemize}
%
\item {\it FIXBEAD}
  \begin{itemize}
    \item  value: 1 integer; 1 optional integer; 2 optional logical
    \item first integer: which bead to hold fix
    \item second integer: on which chain? (for multi-chain mesh runs); default=1
    \item logicals: fix position and/or orientation
  \end{itemize}
%
\item {\it FIXBEAD1}
  \begin{itemize}
    \item  value: no value
    \item hold the first bead fixed, when running Brownian Dynamics
    \item only set up to work {\bf without} Runge-Kutta
  \end{itemize}
%
\item {\it FIXBEADMID}
  \begin{itemize}
    \item  value: no value
    \item hold the middle bead fixed, when running Brownian Dynamics
    \item only set up to work {\bf without} Runge-Kutta
  \end{itemize}
%
\item {\it FIXBOUNDARY}
  \begin{itemize}
    \item  value: 2 or 4 integers
    \item hold the boundary of a multi-chain mesh fixed
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
  \end{itemize}
%
\item {\it FIXBOUNDARY}
  \begin{itemize}
    \item  no value
    \item something to do with multichain mesh simulations...
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
  \end{itemize}
%
\item {\it FORCE}
  \begin{itemize}
    \item values: 2 integers; 3 optional floats
    \item Add a force on a bead in a particular chain
    \item first integer: which bead; second integer: which chain; floats: force vector
    \item Used in Monte Carlo calculations only
  \end{itemize}
%
\item {\it FRICT}
  \begin{itemize}
    \item values: 2 floats; one optional logical (FRICTPERLEN); default FRICTPERLEN is false
    \item if FRICTPERLEN is true, then the position and u-vector friction per length of chain ($\zeta_r, \zeta_b$). For a chain with identical segment sizes, the 2 edge beads have half the friction coefficients of the other beads
    \item if FRICTPERLEN is false, friction coefficients $\zeta_{rb}, \zeta_{ub}$ for the bead positions and orientation vectors explicitly      
    \item Used for Brownian Dynamics calculations only
  \end{itemize}
%
\item {\it GAM}
  \begin{itemize}
    \item values: 1 float; default=1
    \item Fractional preferred segment extension $\gamma$ for the dssWLC energetics.
  \end{itemize}
%
\item {\it GAUSSIANCHAIN}
  \begin{itemize}
    \item no values
    \item For Brownian Dynamics calculations, treat the chain as a plain bead-spring chain with modulus EPAR/2/LS for each spring.
    \item For EQUILDISTRIB calculations, treat it as a bead-spring chain with modulus EPAR/2/LS along $\vec{u}$ and modulus EPERP/2/LS perpendicular to $\vec{u}$.
  \end{itemize}
%
\item {\it INITRANGE}
  \begin{itemize}
    \item values: 4 floats; defaults: 1D0 1D0 1D0 1D0
    \item For Monte Carlo simulations, initial step size ranges
    \item angle range for type 1 moves; shift range for type 1 moves; angle range for type 2 moves; shift range for type 2 moves    
  \end{itemize}
%
\item {\it INTERPPARAMS}
  \begin{itemize}
    \item values: 1 logical; 1 string
    \item if logical is true then at the start of the calculation, extract the energetic and dynamic parameters of the chain by interpolating from a data file. 
    \item string supplies the file containing data for the parameters (see \ref{sec:getparams})
  \end{itemize}
%
\item {\it LOGRTERM}
  \begin{itemize}
    \item no values
    \item For an non-shearable chain, this includes the additional logarithmic terms to make it behave as a chain with $\epsilon_\perp\rightarrow \infty$ rather than just a plain stretchable chain
    \item See appendix in Soft Matter paper for details
    \item implemented in Brownian Dynamics and EQUILDISTRIB calculations only
  \end{itemize}
%
\item {\it LOOPING}
  \begin{itemize}
    \item optional float LOOPRAD; optional string LOOPFILE
    \item track the first-passage looping time in a Brownian Dynamics simulation
    \item LOOPRAD is the radius such that chain is looped if ends approach within this distance
    \item LOOPFILE is the filename in which to save the looping time of each chain 
    \item output columns in LOOPFILE: chain, step when first looped, loop time, end-to-end vector when first looped    
  \end{itemize}
%
\item {\it LP}
  \begin{itemize}
    \item 1 float; default 1
    \item bending modulus ($\epsilon_b$) for the dssWLC model
  \end{itemize}
%
\item {\it LS}
  \begin{itemize}
    \item 1 float; default 1
    \item Segment length ($\Delta$) for the dssWLC model
  \end{itemize}
%
\item {\it MCPRINTFREQ}
  \begin{itemize}
    \item 1 integer; 1 optional integer
    \item first number is how often to print output to screen during MC simulation (in terms of number of steps)
    \item second number is how often to output to the file set by {\em OUTFILE}
    \item When doing EQUILDISTRIB action, second number is how often to print to screen (in terms of number of chains)
    \item By default, second number is same as the first
  \end{itemize}
%
\item {\it MCSTEPS}
  \begin{itemize}
    \item 1 integer; 2 optional integers; defaults: 1000,100,100
    \item First number is total number of MC steps to run if doing the MONTECARLO action
    \item Second number is how often to update statistics, such as average end-to-end distance
    \item Third number is the number of initialization steps in the Monte Carlo before you start calculating statistics or outputting to file
    \item When running the EQUILDISTRIB action, the total number of sampled chains is given by the first parameter here
  \end{itemize}
%
\item {\it NCHAIN}
  \begin{itemize}
    \item one integer
    \item For Brownian Dynamics or Monte Carlo simulations, the number of chains to run in parallel 
    \item {\color{red} Warning: multi-chain Monte Carlo has not been tested in a long while. Use at own risk}
  \end{itemize}
%
\item {\it NOBROWN}
  \begin{itemize}
    \item no values
    \item Do not include the random Brownian forces in the Brownian Dynamics simulations
  \end{itemize}
%
\item {\it NPT}
  \begin{itemize}
    \item 1 integer; 1 optional integer (MAXNPT)
    \item Number of beads in each chain
    \item If second integer provided, can also set the maximum allowed number of beads for the case of dynamic rediscretization. By default MAXNPT is the same as NPT
    \item MAXNPT sets the sizes of all the different arrays in the chain object, so things will break massively if more beads than this ever show up
  \end{itemize}
%
\item {\it OBSTACLE}
  \begin{itemize}
    \item 3 floats
    \item radius, steric modulus, friction coefficient
    \item sets up an obstacle to interact with the chain
    \item {\color{red} Never fully implemented. Do not use!}
  \end{itemize}
%
\item {\it OUTFILE}
  \begin{itemize}
    \item string; default *.out
    \item General output file; * is replaced with suffix for the job
    \item for EQUILDISTRIB action: end-to-end vector and first orientation vector for each chain
    \item for BROWNDYN action: step, chain, energy, end-to-end vector, center of mass, first $\vec{u}$ vector; output frequency set by {\em BDPRINTFREQ}
    \item for MONTECARLO action: step, type 1 move acceptance frequency, type 2 move acceptance frequency, average $R^2$, end-to-end vector, correlation between first and last $\vec{u}$, radius of gyration, correlation between end-to-end vector and first $\vec{u}$, some other stuff; output frequency set by {\em MCPRINTFREQ}    
  \end{itemize}
%
\item {\it OUTPUTBEADWEIGHT}
  \begin{itemize}
    \item 2 optional integers; defaults 500 50
    \item For MC simulations, when dumping snapshots, output an extra two columns which contain, for each mobile bead, the partition function integrated over $\vec{u}$ and integrated over $\vec{r}$ respectively
    \item integers are number of integration points in each dimension; 2-dimensional integration over $\vec{u}$, 3-dimensional over $\vec{r}$
  \end{itemize}
%
\item {\it PARAMFROMSNAPSHOT}
  \begin{itemize}
    \item 1 optional logical; if not supplied then value set to true; if keyword is missing then default is false
    \item if true, then energetic parameters LS, LP, GAM, EPAR, EPERP, EC will be extracted from the input snapshot file (set with {\em RESTART}) rather than using the ones in the parameter file
  \end{itemize}
%
\item {\it REDISCRETIZE}
  \begin{itemize}
    \item 2 floats
    \item dynamic rediscretization based on segment lengths
    \item {\color{red} Not fully implemented. Do not use}
  \end{itemize}
%
\item {\it REDISCRETIZE}
  \begin{itemize}
    \item 2 floats
    \item dynamic rediscretization based on segment lengths
    \item {\color{red} Not fully implemented. Do not use}
  \end{itemize}
%
\item {\it RESTART}
  \begin{itemize}
    \item 1 optional string, 1 optional integer; defaults: start.out, 0
    \item restart calculation from a previously output chain snapshot
    \item first parameter is the snapshot file
    \item second parameters allows for skipping first few configurations in thefile
    \item Will attempt to reach NCHAIN chains from the snapshot file. If there are not enough will start cycling through the configs in the file
    \item For Monte Carlo simulations, the snapshot file contains a number on the first line that indicates what step to start from
  \end{itemize}
%
\item {\it RNGSEED}
  \begin{itemize}
    \item 1 integer; default: false
    \item seed for random number generator
    \item value of 0 will seed with system time in milliseconds
    \item value of -1 will use the last 5 characters in the suffix
    \item value of -2 will use the last 4 charactes in the suffix and the millisecond time
    \item otherwise: the seed is used directly (should be positive)
  \end{itemize}
%
\item {\it RUNGEKUTTA}
  \begin{itemize}
    \item 1 integer; default: 4
    \item what order of runge-kutta method to use with Brownian Dynamics simulations
    \item So far only 1st order (direct Euler's method) and 4th order Runge-Kutta are implemented; only 4-th order has been extensively tested
  \end{itemize}
%
\item {\it SETSHEAR}
  \begin{itemize}
    \item 1 float
    \item set a shear displasement for a chain mesh
    \item  \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
   \end{itemize}
%
\item {\it SHEARABLE}
  \begin{itemize}
    \item 1 logical; default true
    \item set whether chain is shearable
   \end{itemize}
%
\item {\it SNAPSHOTS}
  \begin{itemize}
    \item 1 optional integer, 1 optional string, 1 optional logical; defaults: 1, *.snap.out, false 
    \item Dump snapshots over the course of the calculation (for BROWNDYN, MONTECARLO, or EQUILDISTRIB actions)
    \item integer: how often to dump snapshots; string: snapshot file (* is replaced with suffix); logical: append rather than rewriting the snapshot file
   \end{itemize}
%
\item {\it STARTEQUIL}
  \begin{itemize}
    \item 1 optional integer (EQUILSAMPLETYPE), 1 optional float (STARTEQUILLP)
    \item For Brownian Dynamics simulations only, start from a set of equilibrium chain configurations
    \item EQUILSAMPLETYPE must be 1,2, or 3 and indicates how to generate the equilibrium sampling
    \item EQUILSAMPLETYPE=1, use rejection sampling with a Lorentz envelope. Becomes very inefficient for short segments, high shear modulus
    \item EQUILSAMPLETYPE=2, use rejection sampling with multivariate normal envelope; this is the preferred method for chains with short stiff segments; less efficient for highly flexible segments
    \item EQUILSAMPLETYPE=3, use Monte Carlo sampling; really inefficient
    \item If optional float is supplied, then sample from a bead-rod distribution with the given bend modulus (even if the actual chain for the BD simulations is a proper shearable chain)
   \end{itemize}
%
\item {\it SQUARELATTICE}
  \begin{itemize}
    \item no values
    \item set up a chain mesh with a square lattice
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
   \end{itemize}
%
\item {\it STARTCOLLAPSE}
  \begin{itemize}
    \item no values
    \item start with a collapsed linear chain mesh
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
   \end{itemize}
%
\item {\it STERICS}
  \begin{itemize}
    \item 1 float, 1 optional integer, 1 optional float
    \item  for steric inter-bead interactions
    \item First number is the steric radius, second is how many neighboring beads to skip in steric calculations, third is steric modulus
    \item {\color{red} Steric calculations are not tested! do not use.}
   \end{itemize}
%
\item {\it STRESSFILE}
  \begin{itemize}
    \item 1 string; default: *.stress.out
    \item For Brownian Dynamics simulations, file in which to output shear stress correlation ($C_\text{shear}$) over time
    \item How often to output is set by BDPRINTFREQ
   \end{itemize}
%
\item {\it STRETCHABLE}
  \begin{itemize}
    \item 1 logical; default true
    \item Chain is stretchable
    \item Chains that are shearable but not stretchable are not implemented
    \item Stretchable but not shearable is implemented
   \end{itemize}
%
\item {\it TRACKDIST}
  \begin{itemize}
    \item four integers
    \item For a multi-chain Monte Carlo calculation, track the average distance between two specific beads
    \item numbers are: bead 1, chain 1, bead 2, chain 2
    \item {\color{red} Many-chain mesh calculations are not tested! do not use.}
   \end{itemize}
%
\item {\it USEBDENERGY}
  \begin{itemize}
    \item no values
    \item When running a Monte Carlo simulation, use the energy calculation as it was defined for Brownian Dynamics as opposed to just calculating the local change in energy at each step.
    \item Makes for very inefficient simulations!
   \end{itemize}
%
\item {\it USEPSEUDOFORCE}
  \begin{itemize}
    \item no values
    \item Use pseudo-potential force when running bead-rod Brownian Dynamics simulations (SHEARABLE and STRETCHABLE set to false)
  \end{itemize}
%
\item {\it VERBOSE}
  \begin{itemize}
    \item 1 logical; default: false
    \item print extra output
    \item not really implemented
  \end{itemize}
%
% --------------------------

\end{itemize}

\bibliographystyle{aip} 
\bibliography{fiberModel}

\end{document}