+----- DETAILS: Editing the "oplsaa_subset.prm" file -------
+
+Again, before you run "oplsaa_moltemplate.py", you must edit the "oplsaa.prm"
+file (or "oplsaa_subset.prm file) and eliminate atom types which do not
+correspond to any of the atoms in your simulation. This means you must
+look for lines near the beginning of this file which begin with the word "atom"
+and refer to atom types which appear in the simulation you plan to run. All
+other lines (beginning with the word "atom") must be deleted or commented out.
+(Leave the rest of the file alone.)
+
+For example:
+If you were working with ethylene and benzene you would delete every line
+beginning with the word "atom", except for these four lines:
+
+# for ethylene:
+atom 88 47 CM "Alkene H2-C=" 6 12.011 3
+atom 89 46 HC "Alkene H-C=" 1 1.008 1
+# for benzene:
+atom 90 48 CA "Aromatic C" 6 12.011 3
+atom 91 49 HA "Aromatic H-C" 1 1.008 1
+
+Then you are ready to run oplsaa_moltemplate.py on this file.
+
+(You can try to delete more irrelevant information, but be careful.
+ It is not necessary, and it is easy to make mistakes.)
+
+
+----- Using the "oplsaa.lt" file -----
+
+Once you have created the "oplsaa.lt" file, you can create files (like
+ethylene.lt) which define molecules that refer to these atom types.
+Here is an excerpt from "ethyelene.lt":
+
+Ethylene inherits OPLSAA {
+ write('Data Atoms') {
+ list of atoms goes here ...
+ }
+ write('Data Bond List') {
+ list of bonds goes here ...
+ }
+}
+
+And then run moltemplate.
+
+
+----------- CHARGE: -----------
+
+By default, the OPLSAA force-field assigns atom charge according to atom type.
+When you run moltemplate, it will create a file named "system.in.charges",
+containing commands like:
+
+set type 2 charge -0.42
+set type 3 charge 0.21
+
+(This assumes your main moltemplate file is named "system.lt". If it was
+named something else, eg "polymer.lt", then the file created by moltemplate
+will be named "polymer.in.charges".)
+
+Include these commands somewhere in your LAMMPS input script
+(or use the LAMMPS "include" command to load the commands in system.in.charges)
+
+Note that the atom numbers (eg "2", "3") in this file will not match the
+OPLS atom numbers. (Check the output_ttree/ttree_assignments.txt file,
+created by moltemplate, to see a table of "@atom" type numbers translated
+from OPLSAA into LAMMPS.)
+
+----------- CREDIT -----------
+
+If you use these tools and you publish a paper using OPLSAA, please also cite
+the TINKER program. (Because the original "oplsaa.prm" file which we depend on
+is distributed with TINKER.) I think these are the relevant citations:
+
+1) Ponder, J. W., & Richards, F. M. (1987). "An efficient newtonâlike method for molecular mechanics energy minimization of large molecules. Journal of Computational Chemistry", 8(7), 1016-1024.
+
+2) Ponder, J. W, (2004) "TINKER: Software tools for molecular design", http://dasher.wustl.edu/tinker/
+
+-------------------------------
+
+Jason Lambert and Andrew Jewett
+April, 2014
+
+Please email bugs to jewett.aij@gmail.com and jlamber9@gmail.com
+The original oplsaa.txt from moltemplate version 1.20 file is located here.
+Feel free to pick whichever one works well for your application, but
+if you publish a paper using these parameters please cite TINKER,
+since that appears to be the source of all of these files.
+I think these are the relevant citations:
+
+Ponder, J. W., & Richards, F. M. (1987). "An efficient newtonâlike method for molecular mechanics energy minimization of large molecules. Journal of Computational Chemistry", 8(7), 1016-1024.
+
+Ponder, J. W, (2004) "TINKER: Software tools for molecular design", http://dasher.wustl.edu/tinker/
+
+
+Here is a discussion of a difference between these two files:
+
+>> On Fri, Apr 11, 2014 at 6:44 PM, Andrew Jewett <jewett.aij@gmail.com> wrote:
+>> Jason
+>> Lambert graciously wrote the "oplsaa_moltemplate.py" script and
+>> provided the "oplsaa.txt" file which it reads as an input file (after
+>> some minor manual editing on the part of the user). I did some
+>> googling, and I noticed that this file is very similar to the
+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.)
+
+---- 30-nm fiber model: ----
+
+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.
+
+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,
+ 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.
+and save this file as "oplsaa_subset.prm". Then you must EDIT THIS FILE
+so that it only contains atom types you plan to have in your simulation
+(see below for details). Then run the opls_moltemplate.py script this way:
+
+
+python oplsaa_moltemplate.py oplsaa_subset.prm
+
+
+This will create a file named "oplsaa.lt"
+Look over the newly created "oplsaa.lt" file.
+Then move this file to wherever you plan to run moltemplate. For example:
+
+mv -f oplsaa.lt ..
+
+----- DETAILS: Editing the "oplsaa_subset.prm" file -------
+
+Again, before you run "oplsaa_moltemplate.py", you must edit the "oplsaa.prm"
+file (or "oplsaa_subset.prm file) and eliminate atom types which do not
+correspond to any of the atoms in your simulation. This means you must
+look for lines near the beginning of this file which begin with the word "atom"
+and refer to atom types which appear in the simulation you plan to run. All
+other lines (beginning with the word "atom") must be deleted or commented out.
+(Leave the rest of the file alone.)
+
+For example:
+If you were working with ethylene and benzene you would delete every line
+beginning with the word "atom", except for these four lines:
+
+# for ethylene:
+atom 88 47 CM "Alkene H2-C=" 6 12.011 3
+atom 89 46 HC "Alkene H-C=" 1 1.008 1
+# for benzene:
+atom 90 48 CA "Aromatic C" 6 12.011 3
+atom 91 49 HA "Aromatic H-C" 1 1.008 1
+
+Then you are ready to run oplsaa_moltemplate.py on this file.
+
+(You can try to delete more irrelevant information, but be careful.
+ It is not necessary, and it is easy to make mistakes.)
+
+
+----- Using the "oplsaa.lt" file -----
+
+Once you have created the "oplsaa.lt" file, you can create files (like
+ethylene.lt) which define molecules that refer to these atom types.
+Here is an excerpt from "ethyelene.lt":
+
+Ethylene inherits OPLSAA {
+ write('Data Atoms') {
+ list of atoms goes here ...
+ }
+ write('Data Bond List') {
+ list of bonds goes here ...
+ }
+}
+
+And then run moltemplate.
+
+
+----------- CHARGE: -----------
+
+By default, the OPLSAA force-field assigns atom charge according to atom type.
+When you run moltemplate, it will create a file named "system.in.charges",
+containing commands like:
+
+set type 2 charge -0.42
+set type 3 charge 0.21
+
+(This assumes your main moltemplate file is named "system.lt". If it was
+named something else, eg "polymer.lt", then the file created by moltemplate
+will be named "polymer.in.charges".)
+
+Include these commands somewhere in your LAMMPS input script
+(or use the LAMMPS "include" command to load the commands in system.in.charges)
+
+Note that the atom numbers (eg "2", "3") in this file will not match the
+OPLS atom numbers. (Check the output_ttree/ttree_assignments.txt file,
+created by moltemplate, to see a table of "@atom" type numbers translated
+from OPLSAA into LAMMPS.)
+
+----------- CREDIT -----------
+
+If you use these tools and you publish a paper using OPLSAA, please also cite
+the TINKER program. (Because these examples use the "oplsaa.prm" file which
+is distributed with TINKER.) I think these are the relevant citations:
+
+1) Ponder, J. W., & Richards, F. M. (1987). "An efficient newtonâlike method for molecular mechanics energy minimization of large molecules. Journal of Computational Chemistry", 8(7), 1016-1024.
+
+2) Ponder, J. W, (2004) "TINKER: Software tools for molecular design", http://dasher.wustl.edu/tinker/
+
+-------------------------------
+
+Jason Lambert and Andrew Jewett
+April, 2014
+
+Please email bugs to jewett.aij@gmail.com and jlamber9@gmail.com
+and save this file as "oplsaa_subset.prm". Then you must EDIT THIS FILE
+so that it only contains atom types you plan to have in your simulation
+(see below for details). Then run the opls_moltemplate.py script this way:
+
+
+python oplsaa_moltemplate.py oplsaa_subset.prm
+
+
+This will create a file named "oplsaa.lt"
+Look over the newly created "oplsaa.lt" file.
+Then move this file to wherever you plan to run moltemplate. For example:
+
+mv -f oplsaa.lt ..
+
+----- DETAILS: Editing the "oplsaa_subset.prm" file -------
+
+Again, before you run "oplsaa_moltemplate.py", you must edit the "oplsaa.prm"
+file (or "oplsaa_subset.prm file) and eliminate atom types which do not
+correspond to any of the atoms in your simulation. This means you must
+look for lines near the beginning of this file which begin with the word "atom"
+and refer to atom types which appear in the simulation you plan to run. All
+other lines (beginning with the word "atom") must be deleted or commented out.
+(Leave the rest of the file alone.)
+
+For example:
+If you were working with ethylene, you would delete every line
+beginning with the word "atom", except for these two lines:
+
+
+atom 88 47 CM "Alkene H2-C=" 6 12.011 3
+atom 89 46 HC "Alkene H-C=" 1 1.008 1
+
+
+Then you are ready to run oplsaa_moltemplate.py on this file.
+
+(You can try to delete more irrelevant information, but be careful.
+ It is not necessary, and it is easy to make mistakes.)
+
+
+----- Using the "oplsaa.lt" file -----
+
+Once you have created the "oplsaa.lt" file, you can create files (like
+ethylene.lt) which define molecules that refer to these atom types.
+Here is an excerpt from "ethyelene.lt":
+
+Ethylene inherits OPLSAA {
+ write('Data Atoms') {
+ list of atoms goes here ...
+ }
+ write('Data Bond List') {
+ list of bonds goes here ...
+ }
+}
+
+And then run moltemplate.
+
+
+----------- CHARGE: -----------
+
+By default, the OPLSAA force-field assigns atom charge according to atom type.
+When you run moltemplate, it will create a file named "system.in.charges",
+containing commands like:
+
+set type 2 charge -0.42
+set type 3 charge 0.21
+
+(This assumes your main moltemplate file is named "system.lt". If it was
+named something else, eg "polymer.lt", then the file created by moltemplate
+will be named "polymer.in.charges".)
+
+Include these commands somewhere in your LAMMPS input script
+(or use the LAMMPS "include" command to load the commands in system.in.charges)
+
+Note that the atom numbers (eg "2", "3") in this file will not match the
+OPLS atom numbers. (Check the output_ttree/ttree_assignments.txt file,
+created by moltemplate, to see a table of "@atom" type numbers translated
+from OPLSAA into LAMMPS.)
+
+----------- CREDIT -----------
+
+If you use these tools and you publish a paper using OPLSAA, please also cite
+the TINKER program. (Because these examples use the "oplsaa.prm" file which
+is distributed with TINKER.) I think these are the relevant citations:
+
+1) Ponder, J. W., & Richards, F. M. (1987). "An efficient newtonâlike method for molecular mechanics energy minimization of large molecules. Journal of Computational Chemistry", 8(7), 1016-1024.
+
+2) Ponder, J. W, (2004) "TINKER: Software tools for molecular design", http://dasher.wustl.edu/tinker/
+
+-------------------------------
+
+Jason Lambert and Andrew Jewett
+April, 2014
+
+Please email bugs to jewett.aij@gmail.com and jlamber9@gmail.com
+ (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. 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"). Furthermore, 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 -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 will be created.)
+
+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.
+(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
+
+
+I have not tested this, and have no idea what will happen. (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.