+Documentation, examples, and supporting code can be downloaded at:
+
+http://www.moltemplate.org
+
+## Requirements
+
+Moltemplate requires the Bourne-shell, and a recent version of python
+(2.7, 3.0 or higher), and can run on OS X, linux, or windows. (...if a
+suitable shell environment has been installed. See below.)
+
+
+## INSTALLATION INSTRUCTIONS
+
+This directory should contain 3 folders:
+
+ moltemplate/ <-- source code and force fields
+ doc/ <-- the moltemplate reference manual
+ examples/ <-- examples built with moltemplate
+
+There are two ways to install moltemplate:
+
+## Installation using pip
+If you are familiar with pip, then run the following command from within the directory where this README file is located:
+
+ pip install .
+
+If you receive an error regarding permissions, then run pip with the "--user" argument:
+
+ pip install . --user
+
+Make sure that your default pip install bin directory is in your PATH. (This is usually something like ~/.local/bin/ or ~/anaconda3/bin/. If you have installed anaconda, this will be done for you automatically.) Later, you can uninstall moltemplate using:
+
+ pip uninstall moltemplate
+
+If you continue to run into difficulty, try installing moltemplate into a temporary virtual environment by installing "virtualenv", downloading moltemplate (to "~/moltemplate" in the example below), and running these commands:
+
+ cd ~/moltemplate
+ virtualenv venv
+ source venv/bin/activate
+ pip install .
+ #(now do something useful with moltemplate...)
+
+(You will have to "run source ~/moltemplate/venv/bin/activate" beforehand every time you want to run moltemplate.) If all this fails, then try installing moltemplate by manually updating your \$PATH environment variable. Instructions for doing that are included below.
+
+## Manual installation:
+
+Alternatively, you can edit your $PATH environment variable manually to
+include the subdirectory where the "moltemplate.sh" script is located,
+as well as the subdirectory where most of the python scripts are located.
+Suppose the directory with this README file is named "moltemplate"
+and is located in your home directory:
+
+If you use the bash shell, typically you would edit your
+`~/.profile`, `~/.bash_profile` or `~/.bashrc` files
+After making these changes, you may need to start a new terminal (shell) for the changes to take effect. If you do not know what a `PATH` environment variable is and are curious, read:
+ http://www.linfo.org/path_env_var.html
+(I receive this question often.)
+
+
+### WINDOWS installation suggestions
+
+You can install both moltemplate and LAMMPS in windows, but you will first need to install the BASH shell environment on your computer. If you are using Windows 10 or later, try installing the "Windows Subsystem for Linux (WSL)"
+To use LAMMPS and moltemplate, You will also need to install (and learn how to use) a text editor. (Word, Wordpad, and Notepad will not work.) Popular free text editors which you can safely install and run from within the WSL terminal include: **nano**, **ne**, **emacs**, **vim**, and **jove**. (Unfortunately, as of 2017-5-17, [graphical unix-friendly text editors such as Atom, VSCode, Notepad++, and sublime won't work with WSL, and may cause file system corruption. Avoid these editors for now.](https://www.reddit.com/r/bashonubuntuonwindows/comments/6bu1d1/since_we_shouldnt_edit_files_stored_in_wsl_with/))
+
+## License
+
+Moltemplate is available under the terms of the open-source 3-clause BSD
+ %This manual (like moltemplate) is under development.
+
+\tableofcontents
+
+ %Additionally, several working examples of molecules created
+ %with moltemplate can be found in the ``examples/'' subdirectory
+ %(which is distributed with moltemplate).
+ %These were created to supplement the moltemplate documentation.
+
+\section{Introduction}
+
+Moltemplate is a general molecule builder and force-field database system for LAMMPS. A simple file format has been created to store molecule definitions and force-fields (the LAMMPS-template format, “LT”).
+LT files are templates containing \textit{all} of the text relevant to a particular molecule (including coordinates, bond-topology, angles, force-field parameters, constraints, groups and fixes). Moltemplate can then duplicate the molecule, customize it, and use it as a building-block for constructing larger, more complex molecules. (These molecules can be used to build even larger molecules.) Once built, individual molecules and subunits can be customized (atoms and bonds, and subunits can be inserted, moved, deleted and/or replaced).
+
+Popular force-fields such as AMBER GAFF and OPLS-AA have been converted into LT format, allowing users to quickly create molecules using moltemplate. (With help, more popular force-fields can be converted.) This way moltemplate users can build a molecule by specifying only a list of atoms in the molecule and the bonds connecting them. End-users are not required to manually specify all of the force-field parameters. However they still have the freedom to easily customize individual interactions when needed.
+ %or generate all of its angle, dihedral, improper interactions manually.
+
+Moltemplate is extremely flexible. It supports all LAMMPS force-field styles and nearly all atom-styles (now and in the future).
+
+ % OLD VERSION
+ %Moltemplate is a cross-platform text-based molecule builder for LAMMPS. It is typically used for building coarse-grained toy molecular models. Moltemplate users have access to (nearly) all of the standard and non-standard (custom, user-created) force-field and features available in LAMMPS.
+ %
+ %\textit{(Although optimized for LAMMPS, moltemplate is a general text manipulation tool which, in principle, could be used to generate topology and force-field files for other simulation programs. Please email \includegraphics[height=0.3cm]{author_email.png} if you want to attempt this.)}
+ %
+ %A file format has been created to store molecule definitions (the LAMMPS-template format, ``LT''). Typical ``.LT'' files contain atom coordinates, topology data (bonds), LAMMPS force-field data, and other LAMMPS settings (such as group definitions, fixes, and user-defined input files) for a type of molecule (or a molecular subunit). Molecules can be copied, combined, and linked together to define new molecules. (These can be used to define larger molecules.)
+ %%Unlimited levels of object composition, nesting, and inheritance are supported.)
+ %Once built, individual molecules and subunits can be customized (atoms and bonds, and subunits can be moved, deleted and replaced).
+
+
+Moltemplate requires the Bourne-shell, and a recent version of python (2.7 or 3.0 or higher), and can run on OS X, linux, or windows (if a suitable shell environment has been installed).
+\textbf{A substantial amount of memory is needed} to run moltemplate.
+For example, building a system of 1000000 atoms typically requires
+between 3 and 12 GB of \textit{available} memory.
+(This depends on the number of bonds, molecules, and angular interactions.
+ See section \ref{sec:limitations} for details.)
+%Memory requirements are discussed in section \ref{sec:limitations}.
+
+ %Moltemplate is a text-manipulation tool for generating
+ %input files for molecular dynamics simulation programs.
+ %Moltemplate has been optimized for constructing input files for LAMMPS.
+ %from constituent parts.
+ %Molecules are stored in a hierarchical,
+ %object-oriented,
+ %template-based file format (``.LT'').
+ % %using an object-oriented style
+ % %which can mimic many popular molecular file formats.
+ % %such as PDB, amber TOP, Gromacs TOP,
+ % %PSF files, and some limited xplor parameter files.
+% \subsubsection*{Using ``lttree.py'' instead of ``moltemplate.sh''}
+% The format of an ``.LT'' file closely mimics the syntax in
+% current LAMMPS data and input script files (as of early 2013).
+% However LAMMPS file formats are constantly changing
+% as users add their own custom features to LAMMPS.
+% (In addition, there are some currently known limitations of
+% moltemplate.sh, which are discussed in section \ref{sec:limitations}.)
+% %However this file format must be flexible enough
+% %to handle potentially radical syntax changes in the future.
+% %End users who add new features to LAMMPS may also modify the syntax
+% %of these input files, and will likely introduce new file formats.
+% Consequently, we also provide several simple python scripts:
+% %(which should remain useful when/if moltemplate.sh breaks)
+% ``ttree.py'', ``lttree.py'', and ``nbody\_by\_type.py''.
+% %\begin{list}
+% %\item
+% %``ttree.py'', is a general text manipulation
+% %tool which prints out the text contained in the
+% %``write()'' and ``write\_once()'', commands in an LT file,
+% %and substitutes numerical values into the \$ and \@ variables
+% %contained inside.
+% %\item
+% %``lttree.py'' is a variant of ``ttree.py''
+% % %understand LAMMPS atom\_style syntax and
+% %which also generates atomic coordinates.
+% %(It process the ``.move()'' and ``.rot()'' commands.)
+% %\item
+% %``nbody\_by\_type.py'' is a utility which generates
+% %many-body bonded interactions between atoms automatically,
+% %according to the atom and bond type.
+% %(It processes the ``Data Angles By Type'',
+% %``Data Dihedrals By Type'', and ``Data Impropers By Type'' sections.)
+% %\end{list}
+% The ``ttree.py'' program is a general text manipulation tool which
+% should continue to work in the distant future,
+% even if the LAMMPS syntax changes radically, and ``moltemplate.sh'' breaks.
+% (``ttree.py'' is nearly identical to and supports all the
+% command line options used by ``moltemplate.sh'',
+% with the exception of ``-pdb'', ``-xyz'', and ``-raw''.)
+% %However this tool is intentionally simple and ignorant about LAMMPS.
+% %This allows programmers to add features to LAMMPS without ever
+% %breaking the ``.LT'' file format.
+% %(although you may have to sacrifice some convenience
+% %that using moltemplate.sh provides).
+% A tutorial for using these programs
+% is provided in appendix \ref{sec:ttree}.
+
+
+
+
+\section{Installation}
+
+\subsection*{Obtaining Moltemplate}
+\textit{If you don't already have moltemplate},
+the most up-to-date version can be downloaded at
+\url{http://www.moltemplate.org}
+If you obtained moltemplate as a .tar.gz file,
+(as opposed to github or pip), you can unpack it using:
+\begin{verbatim}
+tar -xzvf moltemplate_2017-8-22.tar.gz
+\end{verbatim}
+(The date will vary from version to version.)
+Alternately, if you obtained moltemplate bundled with LAMMPS,
+then the \textit{``moltemplate''} directory will probably be located
+in the \textit{``tools''} subdirectory of your lammps installation.
+
+There are two ways to install moltemplate:
+
+\subsubsection*{Installation Method 1 (pip)}
+
+\textit{If you are familiar with pip}, then run the following command from within outermost directory:
+\begin{verbatim}
+pip install .
+\end{verbatim}
+\textit{In order for this to work, this directory should contain a file named ``\textbf{setup.py}''.} (If no such file exists, then either proceed to ``Installation Method 2'' below, or download a newer version of moltemplate.) If you receive an error regarding permissions, then run pip this way instead:
+\begin{verbatim}
+pip install . --user
+\end{verbatim}
+Make sure that your default pip install bin directory is in your PATH. (This is usually something like \textapprox/.local/bin/ or \textapprox/anaconda3/bin/. If you have installed anaconda, your PATH should have been updated for you automatically.) Later, you can uninstall moltemplate using:
+\begin{verbatim}
+pip uninstall moltemplate
+\end{verbatim}
+
+\textit{If you continue to run into difficulty}, try installing moltemplate into a temporary virtual environment by installing ``virtualenv'', downloading moltemplate (to ``\textapprox/moltemplate'' in the example below), and running these commands:
+\begin{verbatim}
+cd ~/moltemplate
+virtualenv venv
+source venv/bin/activate
+pip install .
+ #(now do something useful with moltemplate...)
+\end{verbatim}
+You will have to ``run source \textapprox/moltemplate/venv/bin/activate'' beforehand whenver you want to run moltemplate. If all this fails, then try installing moltemplate by manually updating your \$PATH environment variable. Instructions for doing that are included below.
+
+
+\subsubsection*{Installation Method 2}
+
+Alternatively, you can edit your PATH variable manually to include
+the subdirectory where the moltemplate.sh script is located
+(typically ``\textapprox/moltemplate/moltemplate/scripts/''), as well as
+the directory containing the most of the python scripts (``\textapprox/moltemplate/moltemplate/'').
+Suppose the directory where with the README file is named ``moltemplate''
+and it is located in your home directory:
+
+If you use the \textbf{bash} shell, typically you would edit your
+You can install both moltemplate and LAMMPS in windows, but you will first need to install the BASH shell environment on your computer. If you are using Windows 10 or later, try installing the "Windows Subsystem for Linux (WSL)"
+To use LAMMPS and moltemplate, You will also need to install (and learn how to use) a text editor. (Word, Wordpad, and Notepad will not work.) Popular free text editors which you can safely install and run from within the WSL terminal include: nano, ne, emacs, vim, and jove. (Unfortunately, as of 2017-5-17, graphical unix-friendly text editors such as Atom, VSCode, Notepad++, and sublime won't work with WSL, and may cause file system corruption. Avoid these editors for now. (\url{https://www.reddit.com/r/bashonubuntuonwindows/comments/6bu1d1/since_we_shouldnt_edit_files_stored_in_wsl_with/})
+
+
+%\pagebreak
+\section{Quick reference \textit{(skip on first reading)}}
+
+\section*{
+\textit{Note: New users should skip to section \ref{sec:tutorial}}
+}
+
+
+\subsection{Moltemplate commands}
+
+%\begin{table}
+\begin{longtable}[h]{l|p{9cm}}
+\textbf{command} & \textbf{meaning}
+\\
+\hline
+\hline
+\begin{tabular}[t]{l}
+\\
+\textit{MolType} \textbf{\{} \\
+\\
+\hspace{0.35cm} \textit{content} ... \\
+\\
+\textbf{\}} \\
+\end{tabular}
+&
+Define a new type of molecule (or namespace) named \textit{MolType}.
+The text enclosed in curly brackets (\textit{content})
+Allow alternate names for the same variable. This replaces all instances of \textit{oldvariable} with \textit{newvariable}. Both variable names must have a ``@'' prefix. This is typically used to reduce the length of long variables, for example to allow the shorthand ``@atom:C2'' to refer to ``@atom:C2\_bC2\_aC\_dC\_iC''
+\\
+\hline
+ \textbf{\#}\textit{commented text} &
+All text following a ``\#'' character is treated as a comment and ignored.
+\end{longtable}
+
+%\caption{List of moltemplate commands}
+%\label{tab:commands}
+%\end{table}
+
+
+
+%\pagebreak
+\subsection{Common \$ and @ variables}
+
+(See section \ref{sec:variables} for details.) \\
+\begin{tabular}[h]{l|p{11cm}}
+\textbf{variable type} & \textbf{meaning}
+\\
+\hline
+\hline
+\$atom:\textit{name} &
+A unique ID number assigned to atom \textit{name} in this molecule.
+(Note: The \textit{:name} suffix can be omitted if the molecule
+in which this variable appears only contains a single atom.)
+%(This number is unique even if there are multiple copies of this molecule.)
+\\
+\hline
+@atom:\textit{type} &
+A number which indicates an atom's \textit{type}
+ (typically used to lookup pair interactions.)
+\\
+\hline
+\$bond:\textit{name} &
+A unique ID number assigned to bond \textit{name}
+(Note: The \textit{:name} suffix can be omitted if the molecule
+in which this variable appears only contains a single bond.)
+\\
+\hline
+@bond:\textit{type} &
+A number which indicates a bond's \textit{type}
+\\
+\hline
+\$angle:\textit{name} &
+A unique ID number assigned to angle \textit{name}
+(Note: The \textit{:name} suffix can be omitted if the molecule
+in which this variable appears only contains a single angle interaction.)
+\\
+\hline
+@angle:\textit{type} &
+A number which indicates an angle's \textit{type}
+\\
+\hline
+\$dihedral:\textit{name} &
+A unique ID number assigned to dihedral \textit{name}
+(Note: The \textit{:name} suffix can be omitted if the molecule in which
+this variable appears only contains a single dihedral-angle interaction.)
+\\
+\hline
+@dihedral:\textit{type} &
+A number which indicates a dihedral's \textit{type}
+\\
+\hline
+\$improper:\textit{name} &
+A unique ID number assigned to improper \textit{name}
+(Note: The \textit{:name} suffix can be omitted if the molecule in which
+this variable appears only contains a single improper interaction.)
+\\
+\hline
+@improper:\textit{type} &
+A number which indicates an improper's \textit{type}
+\\
+\hline
+\$\textit{mol} \hspace{0.2cm} or \hspace{0.2cm} \$\textit{mol:.} &
+This variable refers to the ID number of \textit{this} molecule object.
+(See section \ref{sec:spce_example}.
+Note: \mbox{\textit{``\$mol''}} is shorthand for \mbox{\textit{``\$mol:.''}})
+\\
+\hline
+\$\textit{mol:}... &
+The ID number assigned to the molecule to which this object belongs
+(if applicable).
+See sections \ref{sec:2beadPeptide},
+\ref{sec:ellipsis_mol},
+%\ref{sec:paths},
+and appendix \ref{sec:adv_variable_syntax}.
+\\
+\hline
+\hline
+\multicolumn{2}{p{16.5cm}} {
+%Variable operations
+\textit{The numbers assigned to each variable are saved in the \textbf{output\_ttree/ttree\_assignments.txt} file}
+%See section \ref{sec:output_ttree}.
+}
+\\
+\hline
+\hline
+\multicolumn{2}{l} {
+%Variable operations
+\quad \textit{\textbf{Advanced variable usage}}
+}
+\\
+\hline
+\textit{\$category}:\textbf{query}()
+&
+Query the current value of the counter in this \textit{\$category}
+without incrementing it.
+(The ``\textit{\$category}'' is usually either \textit{\$atom}, \textit{\$bond}, \textit{\$angle}, \textit{\$dihedral}, \textit{\$improper}, or \textit{\$mol}.)
+This is useful for counting the number of
+atoms, bonds, angles, molecules, etc... created so far.
+\\
+\hline
+\textit{@category}:\textbf{query}()
+&
+Query the current value of the counter in this \textit{@category}
+without incrementing it.
+(The ``\textit{@category}'' is usually either \textit{@atom}, \textit{@bond}, \textit{@angle}, \textit{@dihedral}, or \textit{@improper}.)
+This is useful for counting the number of
+atom types, bond types, angle types, etc... declared so far.)
+\\
+\hline
+\begin{tabular}[t]{l}
+\textit{@\textbf{\{}category:variable\textbf{\}}} \ or \\
+Coordinate transformations introduced using the \textit{push()} command are applied to molecules instantiated later (using the \textit{new}) command, and remain in effect until they are removed using the \textit{pop()} command. (And transformations appearing in arrays accumulate as well, but do not need to be removed with \textit{pop()}.)
+%The \textit{push()} and \textit{pop()} commands allow the user to control exactly how coordinate transformations accumulate. The \textit{pop()} command undoes the transformations introduced in the most recent \textit{push()} command.
+In this example, the first transformation, ``rot()'', is applied to both ``monomer1'' and ``monomer2''. The last transformation, ``move()'', is applied after ``rot()'' and only acts on ``monomer2''.
+\\
+\hline
+\end{tabular}
+%\caption{Coordinate Transformation Commands}
+%\label{tab:transformation_commands}
+%\end{table}
+
+
+
+
+
+\subsection{moltemplate.sh command line arguments:}
+\label{sec:args_table}
+%\begin{table}
+\begin{tabular}[h]{l|p{10cm}}
+\textbf{argument} & \textbf{meaning}
+\\
+\hline
+\hline
+-atomstyle \textit{style}
+&
+Inform moltemplate which atom\_style you are using.
+(\textit{style} is "full" by default).
+Other styles like "molecular" or "hybrid full dipole" are supported.
+For custom atom styles, you can also specify the list of column
+names manually. For example:
+\textbf{-atomstyle "molid x y z atomid atomtype mux muy muz"}
+Atom styles should be enclosed in quotes (").
+\\
+\hline
+-raw coords.raw
+&
+Read all of the atomic coordinates from an external RAW file.
+(RAW files are simple 3-column ASCII files contain X Y Z coordinates
+ for every atom, separated by spaces.)
+\\
+\hline
+-xyz coords.xyz
+&
+Read all of the atomic coordinates from an external XYZ file
+(XYZ files are 4-column ascii files in ATOMTYPE X Y Z format.
+ The first column, ATOMTYPE, is skipped.
+ The first line should contain the number of atoms.
+ The second line is skipped. See section \ref{sec:coords_intro}.)
+\\
+\hline
+-pdb coords.pdb
+&
+Read all of the atomic coordinates from an external PDB file
+(Periodic boundary conditions are also read, if present.
+ Atoms are sorted by the chainID, resID, insertCode, and atomID
+ fields on every line beginning with ``ATOM'' or ``HETATM''.
+ This order must match the order that the atoms appear in the data file.
+ See section \ref{sec:coords_intro}.)
+\\
+\hline
+-a '\textit{variable} \textit{value}'
+&
+Assign \textit{variable} to \textit{value}.
+(The \textit{variable} should begin with either a @ character
+ or a \$ character.
+ Single-quotes and a space separator are required.
+ See appendix \ref{sec:manual_assignment}.)
+\\
+\hline
+-a bindings\_file'
+&
+The variables in column 1 of
+\textit{bindings\_file}
+(which is a text file)
+will be assigned to
+the values in column 2 of that file.
+(This is useful when there are many variable assignments to make.
+See appendix \ref{sec:manual_assignment}.)
+% \$-variables should \textit{not} be preceded by \textbackslash\ in this case.)
+\\
+\hline
+\begin{tabular}[t]{l}
+-b '\textit{variable} \textit{value}'
+\\
+\hspace{0.35cm} \textit{or} \\
+-b \textit{bindings\_file}
+\\
+\end{tabular}
+&
+Assign variables to values.
+Unlike assignments made with ``-a'',
+assignments made using ``-b''
+are non-exclusive.
+(They may overlap with other variables in the same category.
+ See appendix \ref{sec:manual_assignment}.)
+\\
+\hline
+
+\begin{tabular}[t]{l}
+-overlay-bonds
+\\
+-overlay-angles
+\\
+-overlay-dihedrals
+\\
+-overlay-impropers
+\\
+\end{tabular}
+&
+By default moltemplate overwrites
+duplicate bonded interactions which
+involve the same set of atoms.
+These flags disable that behavior.
+This can be useful when you want to superimpose
+multiple angular or dihedral forces on the same set of atoms
+(eg. to enable more complex force fields).
+\\
+\hline
+-nocheck &
+Do \textit{not} check for common LAMMPS/moltemplate syntax errors.
+(This might be useful when using moltemplate
+ with simulation software other than LAMMPS,
+ \textit{or} to build systems which need new non-standard LAMMPS features.)
+\\
+\hline
+-checkff &
+This forces moltemplate.sh to check that there
+are valid angle and dihedral interactions defined for every
+3 or 4 consecutively bonded atoms in the system
+(defined in "Data Angles By Type'' and ``Data Dihedrals By Type" sections).
+\\
+\hline
+-vmd &
+Invoke VMD after running moltemplate to view the system you have just created.
+(VMD must be installed.
+ %This feature uses Axel Kohlmeyer's topotools plugin.
+ See sections \ref{sec:vmd_topotools}, \ref{sec:vmd_advanced} for details.)
+\\
+\hline
+%-import-path LOCATION
+%&
+%When a user imports an .LT file, moltemplate first looks in the directory
+%where it was run, and then in the ``force\_fields'' subdirectory in the
+%moltemplate installation. Additional directories can be appended using
+%this command. (Multiple directories must be separated by ':' characters)
+%This allows moltemplate to look for .LT files
+%in other directories when using ``import''.
+%(Multiple directories must be separated by ':' characters.)
+%\\
+%\hline
+\end{tabular}
+
+\begin{tabular}[h]{l|p{10cm}}
+\hline
+\begin{tabular}[t]{l}
+-dihedral-sym file.py
+\\
+-improper-sym file.py
+\\
+-bond-symmetry file.py
+\\
+-angle-symmetry file.py
+\\
+\end{tabular}
+&
+Normally moltemplate.sh reorders the atoms in each bond, angle, dihedral, and improper interaction before writing them to the DATA file in order to help avoid duplicate interactions between the same atoms if listed in different but equivalent orders. Sometimes this is undesirable. \textit{\textbf{To disable this behavior, set ``file.py'' to ``None''.}} You can also manually choose alternate symmetry rules for unusual force fields. (Such as class2 force fields, dihedral\_style spherical, etc... For an example of the file format for ``file.py'', see the ``nbody\_Impropers.py'' file.)
+\\
+\hline
+\end{tabular}
+
+
+%\pagebreak
+
+
+
+
+\section{Introductory tutorial}
+\label{sec:tutorial}
+\subsection*{\textit{Summary}}
+\textit{Moltemplate is based on a very simple text generator (wrapper) which
+repetitively copies short text fragments into one (or more) files
+and keeps track of various kinds of counters.}
+
+LAMMPS is a powerful but complex program with many contributors.
+Moltemplate is a front-end for LAMMPS.
+Moltemplate users will have to tackle the same steep learning-curve
+(and occasional bugs) that other LAMMPS users must face.
+ %Moltemplate is (intentionally) ignorant about LAMMPS
+ %and molecular dynamics in general.
+ %Gradually other features have been added to moltemplate.sh which make
+ %it somewhat more convenient for generating LAMMPS simulation files.
+Moltemplate files (LT files) share the same file format and
+syntax structure as LAMMPS DATA files and INPUT scripts.
+%Moltemplate can understand some simple LAMMPS commands,
+%and it will attempt to correct user mistakes.
+Moltemplate will attempt to correct user mistakes,
+however users must still learn
+LAMMPS syntax and write LT files which obey it.
+For users who are new to LAMMPS, the easiest way
+to do this is to modify an existing example
+(such as the water box example in this section).
+(The official LAMMPS documentation
+\url{http://lammps.sandia.gov/doc/Manual.html}
+is an excellent reference to look up LAMMPS commands
+you see in these examples that you are not familiar with.)
+
+%In addition to the examples here, there are more complex examples
+%distributed with the moltemplate source code.}
+
+
+\subsection{Simulating a box of water using moltemplate and LAMMPS}
+%Later on when we build the final LAMMPS input and data files,
+%data from the different files (``Data Atoms'', ``Data Bonds'',
+%``In Init'', ``In Settings'', etc...)
+%must be pasted together in the correct order according to LAMMPS conventions.
+%This is a simple task that can be performed by the ``moltemplate'' script
+%(included with ttree),
+%or manually by the user, depending on their familiarity with LAMMPS.
+%LAMMPS has a complex and diverse syntax because
+%it supports a wide variety of force-field types.
+%The commands above are \textit{raw} LAMMPS commands,
+%augmented by ttree variables
+%(like ``@atom:O'' and ``\$atom:o''), which are explained below.
+
+
+ %and/or CHARMM27 parameter files.
+ %We hope ttree is flexible enough that it should remain useful in the future,
+ %even if the LAMMPS input syntax changes radically.
+
+ %\subsection{Before ttree}
+ %The ability to load and combine data from multiple different types of
+ %molecules together is missing from LAMMPS.
+ %\textit{Normally} LAMMPS users are required to manually assign unique
+ %id numbers to \textit{every} atom, bonded, 3-body, and 4-body interaction
+ %in the entire simulation.
+ %Each molecule is assigned a unique id number as well.
+ %For a system with 6000 water molecules, a user would be required to
+ %specify 18000 atom ids, 12000 bond ids, 6000 3-body angle ids
+ %and 6000 molecule ids.
+ %(This should not be done by hand.)
+
+
+\subsection{Visualizing Trajectories}
+\label{sec:vmd_trajectory}
+After you have run a simulation in LAMMPS, there are several programs which
+can visualize the system.
+If you have saved your trajectory in LAMMPS ``dump'' format,
+later you can view it in VMD \cite{VMD}.
+For the purpose of viewing trajectories in LAMMPS,
+I recommend using the following style of ``dump'' commands in the LAMMPS
+input-script that you use when you run LAMMPS:
+\begin{verbatim}
+dump 1 all custom 1000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
+\end{verbatim}
+(The ``all'' and ``1000'', refer to the atom selection and save interval, which may differ depending on the kind of simulation you are running. See \url{http://lammps.sandia.gov/doc/dump.html} for details.)
+
+
+Once you have a dump file, you can view it in VMD using:
+ltemplify.py has \textit{limited} support for ``fix'' and ``group'' commands,
+including ``fix shake'', ``fix rigid'', and ``fix poems''.
+Other fixes must be added manually to the file generated by ltemplify.py.
+(Such as fix ``restrain'', ``bond/create'', ``bond/break'', ``ttm'', etc...)
+
+ltemplify.py can understand simple (static) ``group'' commands, and will include them in the output file, if it can determine that they contain any relevant atoms. (Fixes depending on irrelevant groups are also deleted.)
+
+
+\textit{Note: This feature has not been tested carefully. So please review all of the group and fix commands generated by ltemplify.py to make sure they refer to the correct atoms. And please report any bugs you find. (-Andrew 2014-10-29)}
+
+
+
+
+
+\subsection{Known bugs and limitations (ltemplify.py):}
+\label{sec:ltemplify_limitations}
+%\subsubsection*{Wildcard characters ``*''}
+%Support for wildcards is not consistent throughout an LT file.
+%
+%Wildcard characters like ``*'' currently mean different things
+%in different places.
+%In the \textit{write\_once(``Data Angles By Type'') \{...\}} section,
+%for example, the ``*'' and ``?'' wildcard characters are interpreted
+%as \textit{string wildcards}.
+%This means that ``@atom:C?'' will match ``@atom:C1'', ``@atom:C2'', and
+%``@atom:CA'', but \textit{not} ``@atom:CA2''.
+%However ``@atom:CH*'' will match all of these examples.
+%(See appendix \ref{sec:nbody_by_type}.)
+%Moltemplate ignores ``*'' characters elsewhere in an LT file,
+%and leaves it up to LAMMPS.
+%
+%
+%This means that a ``*'' character appearing in a
+%\textit{pair\_coeff},
+%\textit{bond\_coeff},
+%\textit{angle\_coeff},
+%\textit{dihedral\_coeff},
+%\textit{improper\_coeff},
+%or
+%\textit{group}
+%command, for example,
+%is interpreted (by LAMMPS) as a \textit{numeric wildcard}.
+This means: find the closest ancestor of the \textit{\textbf{cpath}} object containing a category named ``\textit{\textbf{catname}}''. This ancestor determines the category's scope. Counter variables in this category are local to ancestors of that object. In this usage example, \textit{\textbf{lpath}} identifies the location of the variable's corresponding ``leaf'' object
+relative to the category scope object (\textit{\textbf{cpath}}).
+On the other hand, if the the category's scope (\textit{\textbf{cpath}})
+was not explicitly stated by the user (which is typical),
+then the \textit{\textbf{lpath}} identifies the location of the leaf object relative to
+the object in which the variable was referenced
+(the current-context ``.'').
+
+\subsection{Variable shorthand equivalents}
+\label{sec:variables_shorthand}
+
+\subsubsection*{\$\textit{\textbf{catname}}:\textit{\textbf{lpath}} is equivalent to ``\$.../\textit{\textbf{catname}}:\textit{\textbf{lpath}}''}
+ %\label{sec:variables_shorthand_catname:lpath}
+This means: find the closest direct ancestor of the current object containing a category whose name matches \textit{\textbf{catname}}. If not found, create a new category (at the global level). \textit{This is the syntax used most frequently in LT files.}
+
+If the colon is omitted, as in \$\textit{\textbf{lpath}}/\textit{\textbf{catname}},
+then it is equivalent to: \$\textit{\textbf{catname}}:\textit{\textbf{lpath}}.
+Again, in these cases, \textit{\textbf{lpath}} is a path which is relative to the object
+in which the variable was referenced.
+
+If \$\textit{\textbf{lpath}} is omitted, then this is equivalent to \$\textit{\textbf{catname}}:. In other words, the the leaf node is the current node, ``.''. (This syntax is often used to count keep track of molecule ID numbers. You can use the counter variable ``\$mol'' to keep track of the current molecule id number, because it counts the molecular objects in which this variable was defined. In this case the name of the category is ``mol''. As in most examples, the category object, \textit{\textbf{cpath}}, is not specified. This means the category object is automatically global. A global category object means that every molecule object is given a unique ID number which is unique for the entire system, not just unique within some local molecule. As a counter-example, consider amino acid residue counters. Each amino acid in a protein can be assigned a residue ID number which identifies it within a single protein chain. However because their category was defined locally at the protein level, these residue ID numbers are not global, and are not uniquely defined if there are multiple protein chains present.)
+Find the category name and object corresponding to ``\$\textit{\textbf{cpath}}/\textit{\textbf{catname}}:''
+(see above)
+If \$\textit{\textbf{cpath}}/ is blank, then search for an ancestor with a category whose name matches \textit{\textbf{catname}}, as described above.
+To find the variable's corresponding ``leaf object'', start from the CURRENT object (not the category object). If \textit{\textbf{lpath}} is not empty, follow \textit{\textbf{lpath}} to a new position in the tree. Otherwise, start at the current object. (An empty \textit{\textbf{lpath}} corresponds to the current object.) From this position in the object tree search for a direct ancestor which happens to also be ``leaf object'' for some other variable which belongs to the desired category. If no such variable is found, then ttree creates a new variable whose leaf object is the object at the \textit{\textbf{lpath}} position, and put it in the desired category.
+
+\subsubsection*{\$\textit{\textbf{lpath}}/.../\textit{\textbf{catname}} is equivalent to \$\textit{\textbf{catname}}:\textit{\textbf{lpath}}/...}
+If \textit{\textbf{lpath}} is omitted, then start from the current node.
+(In the molecular examples, ``\$.../mol'' is a variable whose category name is ``mol''. The ``leaf object'' for the variable is either the current object in which this variable was defined, OR a direct ancestor of this object which has been assigned to a variable belonging to the category named ``mol''. In this way large objects (large molecules) can be comprised of smaller objects, without corrupting the ``mol'' counter which keeps track of which molecule we belong to. In other words, ``\$.../mol'' unambiguously refers to the ID\# of the large molecule to which this sub-molecule belongs (regardless of however many layers up that may be).)
+\textit{Variables in the output\_ttree/ttree\_assignments.txt file
+ use the this syntax.}
+
+If the user explicitly specifies the path leading up to the cat node, and avoids using ``...'', then \textit{\textbf{lpath}} is interpreted relative to the category object, not the current object (however \textit{\textbf{cpath}} is interpreted relative to the current object). This happens to be the format used in the ``ttree\_assignments.txt'' file (although you can use it anywhere else in an ``.LT'' file). In ``ttree\_assignments.txt'' file, \textit{\textbf{cpath}} is defined relative to the global object. The variables in that file always begin with ``\$/'' or ``@/''. The slash at the beginning takes us to the global environment object (to which all the other objects belong). (Since the variables in the ``ttree\_assignments.txt'' always begin with ``\$/'' or ``@/'', this distinction is usually not important because the category object for most variables usually is the ``global'' root object.)
+ TITLE="Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes",
+ JOURNAL="J. Phys. Chem. B",
+ VOLUME=102,
+ NUMBER=14,
+ PAGES={2569--2577},
+ YEAR=1998
+}
+
+
+@ARTICLE{Raviv++SafinyaBiophysJ2007,
+ AUTHOR="Uri Raviv and Toan Nguyen and Rouzbeh Ghafouri and Daniel J. Needleman and Youli Li and Herbert P. Miller and Leslie Wilson and Robijn F. Bruinsma and Cyrus R. Safinya",
+ TITLE="Microtubule Protofilament Number Is Modulated in a Stepwise Fashion by the Charge Density of an Enveloping Layer",
+ JOURNAL="Biophys. J.",
+ VOLUME=92,
+ NUMBER=1,
+ PAGES={278--287},
+ YEAR="2007",
+ MONTH={January},
+ DAY={1}
+}
+
+
+@ARTICLE{Raviv++SafinyaPNAS2005,
+ AUTHOR="Uri Raviv and Daniel J. Needleman and Youli Li and Herbert P. Miller and Leslie Wilson and Cyrus R. Safinya",
+ TITLE="Cationic liposome-microtubule complexes: Pathways to the formation of two-state lipid-protein nanotubes with open or closed ends",
+ JOURNAL="Proc. Natl. Acad. Sci. USA",
+ VOLUME=102,
+ NUMBER=32,
+ PAGES={11167--11172},
+ YEAR="2005",
+ MONTH={August},
+ DAY={9}
+}
+
+@ARTICLE{Berendsen++StraatsmaJPhysChem1987,
+ AUTHOR="H. J. C. Berendsen and J. R. Grigera and T. P. Straatsma",
+ TITLE="The Missing Term in Effective Pair Potentials",
+ (Rough-draft documentation for dump2data.py and raw2data.py)
+
+ ---- Description ----
+
+"dump2data.py" was originally designed to convert dump files into LAMMPS DATA format (for restarting a simulation from where it left off). However it also reads and writes .XYZ and .RAW (simple 3-column text format) files also.
+
+ Comparison with pizza.py:
+This script duplicates some of the tools in pizza.py, but you don't have to learn python to use it. If you are willing to learn a little python, pizza.py, can handle more general dump files which might cause dump2data.py to crash (eg "atom_style tri"). Unlike "dump2data.py", pizza.py is maintained by the lammps team:
+http://pizza.sandia.gov/doc/Manual.html
+
+ ----- General Usage -----
+
+General usage:
+
+dump2data.py [old_data_file -xyz -raw -last -t time -tstart ta -tstop tb -interval n -multi -center -scale x -atomstyle style] < DUMP_FILE > OUTPUT_FILE
+
+ ----- examples -----
+
+ If your LAMMPS dump file is named "traj.lammpstrj", you can
+extract the coordinates this way:
+
+dump2data.py -xyz < traj.lammpstrj > traj.xyz
+
+This generates a 3-column text file containing the xyz coordinates on each line of each atom (sorted by atomid). If there are multiple frames in the trajectory file, it will concatenate them together this way:
+
+8192
+LAMMPS data from timestep 50000
+1 -122.28 -19.2293 -7.93705
+2 -121.89 -19.2417 -8.85591
+3 -121.6 -19.2954 -7.20586
+: : : :
+8192
+LAMMPS data from timestep 100000
+1 -121.59 -20.3273 -2.0079
+2 -122.2 -19.8527 -2.64669
+3 -120.83 -19.7342 -2.2393
+
+(When using the "-raw" argument to create simple 3-column .RAW files, blank lines are used to delimit different frames in the trajectory.)
+
+---- optional command line arguments ---
+
+If you want to select a particular frame from the trajectory, use:
+...to indicate the desired interval between frames (it must be a multiple of
+the save interval). You can also use "-tstart 500000 and "-tstop 1000000" arguments to limit the output to a particular range of time. (500000-1000000 in this example).
+
+--- creating DATA files ---
+
+"dump2data.py" can also create lammps DATA files. You must supply it with an existing DATA file containing the correct number of atoms and topology information.
+
+If your coordinates are stored in a DUMP file (eg "traj.lammpstrj"), you can create a new data file this way:
+Again, in this example, "10000" is the timestep for the frame you have selected. You can use "-last" to select the last frame. If you do not specify the frame you want, multiple data files may be created...
+
+Creating multiple data files:
+The "-multi" command line argument tells "dump2data.py" to generate a new data file for each frame in the trajectory/dump-file. Those files will have names ending in ".1", ".2", ".3", ... (If you use the "-interval" argument, frames in the trajectory whose timestep is not a multiple of the interval will be discarded.) I can't remember if this behavior is switched on by default.
+
+Reading simple 3-column coordinate files:
+If you have a file containing only the coordinates of the atoms (in sorted order), you can use "raw2data.py" to create a data file with those atoms coordinates.
+(where ATOMSTYLE is a quoted string, such as "full" or "hybrid sphere dipole" discussed earlier. Warning: "raw2data.py" is not a stand-alone script. Make sure raw2data.py is located in the same directory with dump2data.py.)
+
+--- scaling and centering coordinates ---
+
+-center
+ This will center the coordinates around the geometric center, so that the average position of the atoms in each frame is located at the origin. (This script attempts to pay attention to the periodic image flags. As such, I think this script works with triclinic cells, but I have not tested that feature carefully.)
+
+-scale 1.6
+ This will multiply the coordinates by a constant (eg "1.6") (Please email me if this fails with periodic image flags.)
+
+---- limitations ----
+
+Speed.
+The program is somewhat slow, although it should be able to handle big trajectories. If speed is important to you, you probably should write your own custom script or use pizza.py which might be faster.
+
+triclinic cells
+ Support for triclinic cells has been added, but not tested.
+
+exotic atom_styles
+
+ This script was designed to work with point-like atoms, and it extracts the x,y,z (and if present vx,vy,vz velocity) degrees of freedom and (by default) copies it to the new data being created by this script.
+
+By default, this script assumes you are using "atom_style full".
+If you are using some other atom style (eg "hybrid bond dipole"), then you can try to run it this way:
+
+dump2data.py -t 10000 \
+ -atomstyle "hybrid bond dipole" \
+ old_data_file < traj.lammpstrj > new_data_file
+
+In general, the -atomstyle argument can be any of the atom styles listed in the
+table at:
+http://lammps.sandia.gov/doc/atom_style.html
+...such as "angle", "bond", "charge", "full", "molecular", "dipole", "ellipsoid", or any hybrid combination of these styles. (When using hybrid atom styles, you must enclose the argument in quotes, for example: "hybrid sphere dipole")
+Warning: I have not tested using dump2data.py with exotic (non-point-like) atom
+styles. (I suspect that the script will not crash, but the dipole orientations
+will not be updated.)
+
+You can also customize the order columns you want to appear in that file using -atomstyle ”molid x y z atomid atomtype mux muy muz”, but again, I don't think the mux, muy, muz information in the new data file will be accurate.
+
+I also strongly suspect that "dump2data.py" does not currently work with the "tri", "ellipsoid", and new "body" styles.
+
+Again, try using pizza.py if you are simulating systems with exotic data types.
+This will create a new LAMMPS DATA file named "FILE_NEW.data" whose atom
+coordinates are copied from the COORDS.raw file, but is otherwise identical
+to the original DATA file (eg, "FILE_OLD.data"). The optional
+-atomstyle ATOMSTYLE argument tells raw2data.py about the format of the DATA
+file. If not specified, the atom style is "full" by default.
+
+
+Arguments:
+
+ATOMSTYLE is a quoted string, such as "full" or "hybrid sphere dipole" indicating the format of the data file. It can be any of the atom styles listed in the table at:
+http://lammps.sandia.gov/doc/atom_style.html
+...such as "angle", "bond", "charge", "full", "molecular", "dipole", "ellipsoid"
+or any hybrid combination of these styles.
+
+FILE_OLD.data
+The second argument to raw2data.py is the name of a DATA file you want to read.
+raw2data.py will replace the coordinates in the "Atoms" section of this file,
+while preserving the rest of the data file.
+
+COORDS.raw is a simple 3-column ASCII file containing the coordinates of the
+atoms in your system. It has a very simple format:
+-122.28 -19.2293 -7.93705
+-121.89 -19.2417 -8.85591
+-121.6 -19.2954 -7.20586
+-121.59 -20.3273 -2.0079
+-122.2 -19.8527 -2.64669
+-120.83 -19.7342 -2.2393
+ : : :
+
+The order of the atoms in this file should match the ATOM-ID number in the
+first column of the "Atoms" section of the FILE_OLD.data file.
+(...I THINK...
+ To be on the safe side, use a DATA file with the atoms in sorted order.)
+
+Exotic atom styles:
+ When using hybrid atom styles, you must enclose the argument in quotes,
+for example: "hybrid sphere dipole"
+
+ Warning 1: I have not tested using raw2data.py with exotic (non-point-like)
+atom styles. (I suspect that the script will not crash, but dipole orientations
+and other internal degrees of freedom will not be updated.)
+
+ Warning 2: "raw2data.py" is not a stand-alone script. Make sure dump2data.py is located in the same directory with raw2data.py.
+
+ Note: I have not tested it, but I suspect many of the other arguments that work with "dump2data.py", such as "-scale", and "-xyz" also work with raw2data.py.
+
+Try using pizza.py if you are simulating systems with exotic data types.
+These are examples for the "moltemplate" molecule builder for LAMMPS.
+http://www.moltemplate.org
+
+Each directory contains one or more examples.
+
+Each example directory contains:
+
+ images/ This folder has pictures of the molecules in the system
+ moltemplate_files/ This folder contains LT files and other auxiliary files
+ README_setup.sh Instructions for how to use moltemplate (executable)
+ README_visualize.txt Instructions for viewing in DATA/DUMP files in VMD
+
+ ...and one or more LAMMPS input scripts with names like
+
+ run.in.min
+ run.in.npt
+ run.in.nvt
+
+You can run these scripts using
+ lmp_linux -i run.in.npt
+(The name of your lammps binary, "lmp_linux" in this example, may vary.
+ Sometimes, these scripts must be run in a certain order. For example
+ it may be necessary to run run.in.min to minimize the system before you can use run.in.npt, and later run.in.nvt. The README_run.sh file in each subdirectory
+ specifies indicates the order. These files have not been optimized.)
+This is an implementation of the "two-stage" model used by Maxim Imakaev
+in the Naumova et Al 2013 Science paper on metaphase chromatin.
+(Download the supplemental materials section and scroll down to the section:
+ "Two-stage process: linear compaction - axial compression")
+
+---- SMALL MODIFICATION ----
+
+Unlike that study, I did not use "softened" Lennard-Jones potentials
+(which allow the chains to pass through each other).
+
+--- Why use moltemplate? ---
+
+Honestly, you don't need to use moltemplate to build this polymer.
+It is almost counter-productive to use moltemplate to build this kind of
+polymer because it is so simple. (The polymer has only 1 bead per atom.
+It just makes it more complicated to introduce all these extra
+files including monomer.lt, condensin.lt and system.lt, especially considering
+that system.lt is a complex file which is generated by a separate script.)
+
+However building the sytem using moltemplate may pay off if you
+replace each point-like monomer with a multi-atom molecule later on.
+(Right now, using moltemplate to build this system is sort of overkill.
+ I'll post an example of building more complex models of chromatin eventually.)
+
+Anyway, the two-stage model at the end of Naumova et al Science 2013 uses the "30nm-fiber" model, whose details are (somewhat vaguely) described in the supplemental materials section.
+
+---- 10-nm fiber model: ----
+
+For the 10nm model,
+ n=128000,
+ L=200,
+ U(alpha)=5*(1 - cos(alpha))
+ bond_length=1.0 (=10nm)
+ sigma=1.0 (particle radius = 10nm)
+
+---- 30-nm fiber model: ----
+
+"The 30nm-like fiber was modeled by increasing the volume of each monomer and the amount of DNA represented by each monomer by a factor of 4.25, while keeping other parameters the same at the monomer level."
+
+I interpret this to mean that, for the 30nm model,
+ n=128000/4.25~=30117 (however I rounded up to 32768=2^15)
+The average diameter of a mammalian cell nucleus is is 6 micrometers (μm),
+which occupies about 10% of the total cell volume.
+
+(See "Molecular Biology of the Cell, Chapter 4, pages 191–234 (4th ed.)", by
+ Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, Peter Walter, 2002)
+
+... of that, 25% of it is occupied by the nucleolus
+http://en.wikipedia.org/wiki/Nucleolus
+("citation needed")
+---------------
+
+From the supplemental material for the original HiC paper
+(Lieberman-Aiden et al., Science 2009)
+
+Appendix 1.
+ Estimate of the Volume Fraction of Chromatin in Human Cells
+In the simulations we sought to obtain an ensemble of structures that, in their statistical properties, resemble some of the features of chromatin arrangement in the cell. Below we demonstrate that chromatin occupies a significant fraction of cell volume, a property that we reproduced in simulations. Taking the nuclear diameter of a tissue culture cell to be 5-10um, and assuming close to a spherical shape we obtain the volume in the range 50-500 um^3, with a (geometric) mean of ~160 um^3. If we assume that the chromatin is built of DNA wrapped around nucleosomes, then we have 6x10^9bp/200bp=3x10^7 nucleosomes. Each may be approximated as a cylinder ~10nm in diameter and ~5nm in height, suggesting a volume of about 500nm3 each. The linker DNA after each nucleosome is about 50bps long, suggesting a volume of about 50*0.34nm*3.14*1nm^2=50nm^3. Thus the total volume of chromatin = 550x3x10^7 =16 um^3, or ~10% (3-23%) of the nuclear volume. This strikingly large volume fraction is itself a significant underestimate, since we ignored, among other things, all other DNA-bound proteins. Note that any further packing or localization of chromatin inside the nucleus will increase local density.
+---- This next section mostly only justifies why they ----
+---- they did not stop the simulation when the globules ----
+---- were fully crumpled (ie with uniform density) ----
+ In our simulations, the radius of the final crumpled globule was R≈12.5 and the volume V≈8000 cubic units. The total volume of the 4000 monomers, 1 unit in diameters each, is V≈2000. This implies a volume fraction of about 25%, which is consistent with the volume fraction estimated above.
+ ---- ----
+
+Appendix 2.
+ Monomer length in base pairs
+Each monomer of the chain corresponds to a fragment of chromatin that equals the Kuhn length of the chromatin fiber, i.e. approximately twice the persistence length of the fiber. Although the persistence length of the chromatin fiber is unknown it can be estimated using the following arguments. DNA is packed into nucleosomes, where 150 bps are wrapped around the histone core and do not contribute to flexibility of the fiber. The linker DNA of about 50 bps that connects consecutive nucleosomes is bendable, and is the source of flexibility in the fiber. Since the persistence length of double-stranded DNA is 150 bps, an equally flexible region of the nucleosomal DNA should contain 3 linkers, i.e. 3 consecutive nucleosomes packing about 600 bps of DNA. The excluded volume of the nucleosomes, nucleosome interactions, and other DNA-bound proteins can make the fiber less flexible or prohibit certain conformation and may tend to increase the persistence length of the fiber. Using this estimated lower bound estimate for the persistence length, we obtain the Kuhn length of the equivalent freely-jointed chain to be 6 nucleosomes, or ~ 1200bp. A simulated chain of 4000 monomers corresponds to 4.8Mb of packed DNA. The size of each monomer was chosen such that its volume is equal to (or slightly above) that of 6 nucleosomes (V=6 x 600 nm^3); thus the radius of the spherical monomer is R=10nm. The diameter of each globule shown in Figure 4 is about 200 nm.
+The two-stage model at the end of Naumova et al Science 2013 uses the "30nm-fiber" model, whose details are (somewhat vaguely) described in the supplemental materials section.
+
+For the 10nm model,
+ n=128000,
+ L=200,
+ U(alpha)=5*(1 - cos(alpha))
+ bond_length=1.0 (=10nm)
+ sigma=1.0 (particle radius = 10nm)
+
+30nm-fiber model details:
+"The 30nm-like fiber was modeled by increasing the volume of each monomer and the amount of DNA represented by each monomer by a factor of 4.25, while keeping other parameters the same at the monomer level."
+
+I interpret this to mean that, for the 30nm model,
+To increase the volume by a factor o 4.25, originally I thought I should increase the "sigma" parameter from 1.0 to 4.25^(1/3)~=1.6198. But I suspect that the bond-lengths between monomers should be fixed at 1.0. If that is the case, then, perhaps I should increase "sigma" from 1.0 to 4.25^(1/2)~=2.061552, and keep the bond-length fixed at 1.0 (which in the units used by thsi paper, corresponds to 10nm). (That would increase the volume of a cylinder of radius "sigma" and length="bond-length" by a factor of 4.25)
+ bond_length=1.0 (10nm again)
+ sigma=2.061552 (Yes, this is less than 3.0<-->30nm. See below.)
+
+
+
+
+--- Excerpts from the Supplemental section of Naumova et al Science 2013 ---
+
+From p. 18 of the supplemental materials section of Naumova et al Science 2013.
+
+ (This section was probably written by Maxim Imakaev.)
+
+ In vivo, the structure of the chromatin fiber can be complicated and many details remain unknown, particularly in metaphase. Given this uncertainty, we simulated chromatin as a homogeneous “beads-on-a-string” polymer fiber. We consider a 10nm fiber, as the pervasiveness of the 30nm fiber in vivo has become increasingly contested. In our simulations, 77Mb is represented by a densely-packed 10nm fiber of 128,000 monomers. Each monomer represents a 10nm-sized DNA-histone complex containing 3 nucleosomes (around 600bp). The fiber has a persistence length of 4 monomers (~2.4Kb), which is based on earlier estimates of 5-10 nucleosomes for interphase (14). Those estimates arise from the assumption that 5-10 linker DNA fragments, each of 20-40bp, can collectively provide flexibility equal to that of the 150bp persistence length of DNA. Binding of proteins to the linker DNA (e.g. histone H1) and interactions between neighboring nucleosomes can further constrain dynamics, requiring more linkers to provide the persistence length. Due to the tight packing of nucleosomes in metaphase, we use the upper limit of this range, i.e. 12 nucleosomes.
+ For the consecutive loops on a scaffold model (the final folded state model with the best agreement with Hi-C data), we also performed simulations with a more flexible 10nm fiber, or with a 30nm fiber, and found similar results. The more flexible 10nm fiber was modeled by decreasing persistence length to 1.8 monomers. The 30nm-like fiber was modeled by increasing the volume of each monomer and the amount of DNA represented by each monomer by a factor of 4.25, while keeping other parameters the same at the monomer level. We note that a classic model of a 30nm fiber is much less dense than a compact metaphase chromosome. A textbook model of a 30nm fiber assumes packing of about 6 nucleosomes per 10nm of fiber length. This model predicts that only 28% of the volume of the fiber (a 30nm-diameter cylinder) is occupied by nucleosomes, assuming a nucleosome shell volume of 328 nm^3. This is much less than the estimated 30-50% density of nucleosomes in a metaphase chromosome, assuming a diameter of 600nm, a packing density of 50-70 Mb/um, and the same nucleosome volume. (See also (15), which gives an estimate of 0.14-0.18 pg/μm for DNA only, and would give about twice the density if DNA is counted with nucleosomes). As follows, these fibers would have to interdigitate, and fill in gaps within each other. We account for this overlap by assuming the effective diameter of the fiber to be less than 30nm. The effective diameter was chosen to make the volume of the fiber equal the total volume of all the nucleosomes.
+
+ We accounted for topoisomerase II activity by allowing chromatin fibers to pass through each other, while still having excluded volume interactions. This was achieved by using a soft-core Lennard-Jones potential with 1kT energy cost for monomer overlap (see below). This allows for changes in the topological state of a chromosome that are known to occur during compaction in vivo.
+
+Our simulations of a two-step folding process show that Hi-C data for mitotic chromosomes is consistent with a linearly compressed array of consecutive chromatin loops. Whereas mechanisms for formation of consecutive chromatin loops have been proposed, the process of axial compression is less understood. Chromatid compression cannot be accomplished by increased chromatin-chromatin affinity alone, as this would lead to condensation into a globular geometry (14, 16, 17). However, mechanisms which locally compress the fiber of loop bases naturally allow for anisotropic compression into a shorter and thicker fiber, with the same width regardless of chromosome length (18). Differences in the duration or efficiency of the first and second stages of chromosomal condensation provide a natural mechanism for condensation-related proteins to separately affect mitotic chromosome length and width (19). We also note that the axis of loop-bases in our two-stage model does not necessarily form a continuous and rigid scaffold (Figure S26). As follows, we remain agnostic about the molecular details of the chromosomal scaffold, which might for example be formed by a network consisting of protein-protein and/or protein-DNA interactions (20).
+
+ 1. Polymer simulations
+
+ To perform Langevin dynamics polymer simulations we used OpenMM, a high-performance GPU-assisted molecular dynamics API (21, 22). To represent chromatin fibers as polymers, we used a sequence of spherical monomers of 1 unit of length in diameter. Here and below all distances are measured in monomer sizes, set to be 10nm unless specified otherwise. Neighboring monomers are connected by harmonic bonds, with a potential U = 100*(r - 1)^2 (here and below in units of kT). Polymer stiffness is modeled with a three point interaction term, with the potential U = 5*(1 - cos(alpha)), where alpha is the angle between neighboring bonds. All monomers interact via either a shifted Lennard-Jones (LJ) repulsive potential, or an attractive Lennard-Jones potential. At high densities in a confined volume, the details of the inter-monomer interactions become negligible due to screening (23), and we therefore used the computationally efficient shifted LJ potential. The shifted LJ potential allows for a short-range repulsion by truncating the LJ potential at its minimum and shifting the minimum to zero: U = 4 * (1/r^12 - 1/r^6) + 1, for r<2^(1/6); U=0 for r > 2^(1/6). The shifted LJ potential is one of the most computationally efficient repulsive potentials due to a very short cutoff radius.
+
+ To allow chain passing, which represents activity of topoisomerase II, we softened the shifted LJ potential by truncating the interaction energy at Ecutoff = 1 kT. At energies more than 0.5 Ecutoff, the LJ potential was softened via: Usoftened = 0.5 * Ecutoff * (1 + tanh(2*U/Ecutoff - 1)). To avoid numerical 19instabilities in the calculation of U at r ~ 0, the interaction radius r was truncated at r=0.3 via: rtruncated = (r^10 + (0.3)^10)^0.1, which introduced negligible shift in a final softened potential. For an attractive LJ potential, we used: U = 4 * e * (1/r^12 - 1/r^6), with e = 0.46 kT, slightly below the theta-temperature. The attractive potential was similarly softened at 2 kT and cut off at r=2.5. Unless noted, we used a softened shifted LJ repulsive potential.
+
+ Polymer models were visualized using Pymol and Rasmol. For images with loop bases highlighted, a base of each loop and 3 monomers surrounding it in each direction were labeled in red.
+
+ SECTIONS 2-5 SKIPPED
+
+6. Two-stage process: linear compaction - axial compression
+
+To simulate the two-stage process of metaphase chromosome folding, we used the 30nm fiber representation described above for its computational efficiency. Simulations were initialized from 30000 monomer fractal globule conformations; fractal globule is a model for interphase chromatin organization. First, random consecutive loops with L=100 monomers (see above) were introduced, and anchors of neighboring loops were brought together using harmonic springs with a potential U = k * (r – r0)2; r0=0.5. To avoid abrupt motion of the loop anchors, the force was gradually turned on over the first
+300000 timesteps, with k linearly increasing in time from 0 to 10 kT. We used softened shifted repulsive LJ potential for inter-monomer interaction.
+
+Upon completion of linear compaction, axial compression was initiated. This involves following changes: the repulsive LJ force is replaced with an attractive LJ force for all monomers, and the chromosomal core of loop anchors is homogeneously compressed. To achieve the latter, all anchor pairs separated by less than 30 anchors were attracted via a potential U = step(d-3) * abs(d-3) * 10 kT, which implements a constant attractive force between two anchors if they are separated by a distance larger than 3. The interactions between neighboring loop anchors were kept throughout this process.
+
+To obtain the contact map from this simulation, 50 independent runs of 1.5e7 timesteps were performed, and 250 conformations were collected from the second half of each run. The contact map was calculated from all conformations of all runs at a 30-monomer resolution, and was further averaged over three 10000-monomer blocks along the diagonal of the heatmap. The latter was done to show contact map at a relevant length scale (0 to 25 Mb), and to achieve a better averaging of the contact map.
+Msg {write() { On the ${day}th day of Christmas, my true love gave to me:
+}
+}
+Gifts1 {write(){@day partridge in a pear tree.
+
+}} Gifts2 {write(){@day turtle doves, and
+} gifts = new Gifts1 } Gifts3 {write(){@day french hens,
+} gifts = new Gifts2 } Gifts4 {write(){@day calling birds,
+} gifts = new Gifts3 } Gifts5 {write(){@day golden rings,
+} gifts = new Gifts4 } Gifts6 {write(){@day geese a-laying,
+} gifts = new Gifts5 } Gifts7 {write(){@day swans a-swimming,
+} gifts = new Gifts6 } Gifts8 {write(){@day maids a-milking,
+} gifts = new Gifts7 } Gifts9 {write(){@day ladies dancing,
+} gifts = new Gifts8 } Gifts10 {write(){@day lords a-leaping,
+} gifts = new Gifts9 } Gifts11 {write(){@day pipers piping,
+} gifts = new Gifts10 } Gifts12 {write(){@day drummers drumming,
+} gifts = new Gifts11 }
+Msg1 inherits Msg{g = new Gifts1} Msg2 inherits Msg{g = new Gifts2} Msg3 inherits Msg{g = new Gifts3} Msg4 inherits Msg{g = new Gifts4} Msg5 inherits Msg{g = new Gifts5} Msg6 inherits Msg{g = new Gifts6} Msg7 inherits Msg{g = new Gifts7} Msg8 inherits Msg{g = new Gifts8} Msg9 inherits Msg{g = new Gifts9} Msg10 inherits Msg{g = new Gifts10} Msg11 inherits Msg{g = new Gifts11} Msg12 inherits Msg{g = new Gifts12}
+msg1 = new Msg1 msg2 = new Msg2 msg3 = new Msg3 msg4 = new Msg4 msg5 = new Msg5 msg6 = new Msg6 msg7 = new Msg7 msg8 = new Msg8 msg9 = new Msg9 msg10 = new Msg10 msg11 = new Msg11 msg12 = new Msg12
+This folder contains an lt file for a Cooke-type 3-bead lipid, as described in:
+
+"Tunable generic model for fluid bilayer membranes"
+ Cooke IR, Kremer K, Deserno M, Physical Review E, 2005
+
+Due to the form of the forcefield, this requires a 'tabulated potential' style in lammps. This is easily generated using the included python script. Usage as follows:
+
+ gen_potential-cooke.py w_c
+
+where w_c is an optional parameter as described in the original paper (10.1103/PhysRevE.72.011506) with default value of 1.6
+
+This creates the 'tabulated_potential' file which is needed by lammps during the simulation.