updated (most of) the OPLS and LOPLS examples for compatibility with …

…the 2023 version of OPLSAA

README.md

@@ -23,7 +23,7 @@ is a *general* cross-platform text-based molecule builder for
 Moltemplate was intended for building custom coarse-grained molecular models,
 but it can be used to prepare realistic all-atom simulations as well.
 It currently supports the
-[**OPLSAA**(2008)](./examples/all_atom/force_field_OPLSAA),
+[**OPLSAA**(2023)](./examples/all_atom/force_field_OPLSAA),
 [**OPLSUA**](./examples/all_atom/force_field_OPLSUA_united_atom),
 [**LOPLS**(2015)](./examples/all_atom/force_field_OPLSAA/hexadecane),
 [**AMBER**(GAFF,GAFF2)](./examples/all_atom/force_field_AMBER),

examples/all_atom/force_field_OPLSAA/README.md

@@ -0,0 +1,93 @@
+This directory contains some examples of all-atom simulations using the OPLSAA force field.
+
+### WARNING
+
+There is no gaurantee that simulations prepared using moltemplate will reproduce the behavior of other MD codes.  If you notice a problem with these examples, please report it. Peer-review is the only way to improve this software (or any software).  (jewett.aij @ gmail.com)
+
+
+### Duplicate dihedrals, angles, and bonds
+
+Sometimes OPLSAA has multiple plausible choices for the dihedral forces
+between atoms in your molecules.
+(This also sometimes occurs for angles and bonds.)
+When this happens, a "warning_duplicate_dihedrals.txt" file will be
+generated.  *You can ignore these files*.  But if you want to maximize
+the accuracy of your simulation, you should read the warning messages
+in those files and modify your .lt files accordingly.
+
+Several example .lt files demonstrate how to do that:
+- butane/moltemplate_files/butane.lt
+- benzene+benzoic_acid/moltemplate_files/benzoic_acid.lt
+
+
+### Atomic charges
+
+In most of the OPLSAA examples,
+the atomic charges are determined by their @atom types
+*(...according to a lookup table located at the beginning of the
+["oplsaa.lt"](../../../moltemplate/force_fields/oplsaa.lt) file)*.
+*(Any atomic charges listed in the "Data Atoms" section of your molecule's
+LT files will be ignored.)*
+**These charges can be overridden.**
+
+
+### Customizing atomic charges in OPLSAA molecules
+
+#### Background information
+
+LAMMPS provides two different methods to specify atomic charges:
+1) **Specify charges in a DATA file** (eg "system.data").
+*(This is the most popular way to specify atomic charges.
+After running moltemplate.sh, the information in the "Data Atoms" section
+is copied into the "Atoms" section of the "system.data" file created by
+moltemplate.sh, and later read by LAMMPS.)*
+2) **Specify charges using "set" commands.**
+*(This is how the OPLSAA atom charges are specified.
+After running moltemplate, atom charge information in the
+["oplsaa.lt" file](../../../moltemplate/force_fields/oplsaa.lt)
+is copied into the "system.in.charges" file created by moltemplate.sh.
+A LAMMPS input script file (eg. "run.in.nvt" or "run.in.npt")
+is included with all of the OPLSAA examples.  It tells LAMMPS to read
+this "system.in.charges" file after reading the "system.data" file,
+This overrides the atom charges from the "system.data" file.)*
+
+
+#### How to customize atomic charge
+*(without modifying "oplsaa.lt")*
+
+1) If you use use a 3rd-party program to calculate each atom's charge, you can
+copy those charges into the "Data Atoms" section of your molecule's LT files.
+To prevent LAMMPS from ignoring these charges, delete or comment-out the line
+containing: **"include system.in.charges"** from your LAMMPS input script
+(such as "run.in.min", "run.in.nvt", and "run.in.npt").
+2) Alternatively, if you only want to override the charges of *some* of the
+atoms in your molecules (and use default "oplsaa.lt" charges for the remaining
+atoms), then you can do this by adding an "In Charges" section to your LT file
+and providing a list of custom charges for the \$atoms you want to modify.
+This is demonstrated in the ["graphene_nh2.lt"](functionalized_nanotubes_NH2/moltemplate_files/graphene_nh2.lt)
+file located in [this example](functionalized_nanotubes_NH2).
+*(In that example, the charge of one of the carbon atoms in the "Graphene_NH2"
+object was modified.  If you use this strategy, do not comment out
+"include system.in.charges" from your "run.in\*" script files.)*
+3) The discussion so far only applies to molecules that use the OPLSAA force
+field *(i.e. molecules whose definition begins with "inherits OPLSAA")*.
+You can also mix molecules that use OPLSAA with other molecules
+that don't.  In the [waterSPCE+methane](waterSPCE+methane) example,
+the SPC/E water molecules do not use OPLSAA.
+Hence, their atomic charges are located in the "Data Atoms" section
+of the [spce.lt](waterSPCE+methane/moltemplate_files/spce.lt) file.
+*(The same is true of most of the carbon atoms in the
+[carbon nanotube](functionalized_nanotubes_NH2) example.)*
+
+
+
+### Improper angles
+
+The style of improper interaction used by OPLS force fields depends on an angle which depends on the order of the atoms surrounding the central atom. When multiple atoms have the same type, this creates ambiguity in atom order. Since there is no guarantee that moltemplate will choose the same atom order as other molecule builders (such as VMD), this can lead to small unavoidable discrepancies in energies and forces computed by LAMMPS and NAMD.  But their effect should be neglegible.
+*(Please let us know if this is not the case.)*
+
+### Bloated lammps input scripts
+
+By default, LAMMPS input scripts prepared using moltemplate contain the entire contents of the OPLS force-field, even when simulating small systems with just a few atom types.
+
+This is harmless, but if you want to get rid of this extra information, follow the instructions in the "README_remove_irrelevant_info.sh" files.

examples/all_atom/force_field_OPLSAA/alkane_chain_single/README.md

@@ -0,0 +1,65 @@
+Long alkane chain example
+==============
+This example is a simple simulation of a long alkane chain, in a vacuum at room temperature using the OPLSAA force field.  The molecule in this example was constructed from monomeric subunits (named "CH2", and "CH3").  In this example, the polymer was initially straight.  However you can use the ["genpoly_lt.py" script](../../../../doc/doc_genpoly_lt.md) to generate .lt files describing long polymers that can wrap around any curve.  (An example of "genpoly_lt.py" usage can be found [here](../../../coarse_grained/DNA_models/dsDNA_only/2strands/3bp_2particles/confined_viral_DNA).)  *Note: This particular example uses the a variant of the OPLSAA force-field suitable for long alkane chains (sometimes called the "LOPLSAA" force field).*
+
+
+#### Images
+
+<img src="images/ch2_ry90.jpg" width=110> <img src="images/plus.svg" height=80> <img src="images/ch3_ry60.jpg" width=110>
+<img src="images/rightarrow.svg" height=80> <img src="images/alkane50_t=0_straight.jpg" width=250> <img src="images/rightarrow.svg" height=80> <img src="images/alkane50_t=1ns_equilibrated.jpg" width=150>
+
+The length of the polymer can be controlled by editing the [alkane50.lt file](moltemplate_files/alkane50.lt).  The simulation contitions can be controlled by editing the [run.in.nvt file](run.in.nvt).
+
+
+### *Suggestion: Start with the "butane" example*
+
+If this is your first time learning how to build a polymer in moltemplate,
+I suggest starting with the [butane](../butane) example instead.
+
+
+### Instructions
+
+1) To build the files which LAMMPS needs, follow the instructions in:
+[README_setup.sh](README_setup.sh)
+
+2) To run LAMMPS with these files, follow these instructions:
+[README_run.sh](README_run.sh)
+
+(The instructions in "README_remove_irrelevant_info.sh" are optional.  *(If you notice a problem with this example, please [report it](../README.md).*)
+
+
+### Details
+
+The "Alkane50" molecule, as well as the "CH2", and "CH3" monomers it contains, use the LOPLSAA force-field.  As with all of the OPLSAA examples, when we define these molecules, we only specify the atom names, bond list, and coordinates.  We do not have to list the atom charges, angles, dihedrals, or impropers.  The rules for creating atomic charge and angle topology are contained in the ["loplsaa.lt"](../../../../moltemplate/force_fields/loplsaa.lt) and  ["oplsaa.lt"](../../../../moltemplate/force_fields/oplsaa.lt) files.  To let moltemplate know that you want to use these rules, define your molecules (and molecular subunits) this way:
+
+
+```
+import "loplsaa.lt"
+CH2 inherits OPLSAA { ... }         # (see "ch2group.lt")
+CH3 inherits OPLSAA { ... }         # (see "ch3group.lt")
+Alkane50 inherits OPLSAA { ... }    # (see "alkane50.lt")
+```
+
+
+#### OPLSAA or LOPLSAA"?
+
+There are only a few differences between LOPLSAA and OPLSAA.  The LOPLSAA force field contains a few extra atom types and dihedral interactions which improve the accuracy of long alkane chains.  The ["loplsaa.lt"](../../../../moltemplate/force_fields/loplsaa.lt) file (referenced above) incorporates these extra atom and dihedral types in the existing OPLSAA force field *instead* of creating a new force field named "LOPLSAA".  *(There is no separate "LOPLSAA" force field object.  I apologize if this is confusing.)*  You can mix LOPLSAA and OPLSAA atoms in the same molecule.
+
+
+
+### Customizing atomic charges
+
+In this example, atomic charge for OPLSAA atoms is determined by @atom type
+*(...according to a lookup table located at the beginning of the
+["oplsaa.lt"](../../../moltemplate/force_fields/oplsaa.lt) file)*.
+*(Any atomic charges listed in the "Data Atoms" section of your molecules'
+LT files will be ignored.)*
+**These charges can be overridden.**
+See [here](../README.md#Customizing-atomic-charges-in-OPLSAA-molecules)
+for instructions explaining how to customize atomic charge.
+
+
+### Manual control of bond and angle interactions
+
+If necessary, you can customize existing bonds, angles, dihedrals etc. in your molecule (eg. *Alkane50*), or add new ones (if the force field does not define them).  To do this, edit the corresponding LT file (eg. ["alkane50.lt"](./moltemplate_files/alkane50.lt)), and add extra sections to that file (eg. *write("Data Bonds")* or *write("Data Angles")*).  Then add a list of bonded interactions to these sections (containing lines similar to *"\$bond:c7h5 @bond:CustomType \$atom:c7 \$atom:h5"*).  By default, this will override the bond and bonded angular interactions created by the force field.  For more details, read the chapter in the moltemplate manual named "Customizing molecule position and topology".)
+

examples/all_atom/force_field_OPLSAA/alkane_chain_single/README_remove_irrelevant_info.sh

@@ -0,0 +1,11 @@
+
+# Note: By default, the system.data and system.in.settings files contain
+# extra information for atoms defined in OPLSAA which you are not using
+# in this simulation.
+# This is harmless, but if you to delete this information from your
+# system.in.settings and system.in.data files, run this script:
+
+cleanup_moltemplate.sh
+
+# (Note: Removing unecessary atom types will make it easier to visualize the
+#        simulation in VMD.)

examples/all_atom/force_field_OPLSAA/alkane_chain_single/README_run.sh

@@ -0,0 +1,34 @@
+# --- Running LAMMPS ---
+#
+# The 2 files "run.in.npt", and "run.in.nvt" are LAMMPS
+# input scripts which link to the input scripts and data files
+# you hopefully have created earlier with moltemplate.sh:
+#   system.in.init, system.in.settings, system.data
+# If not, carry out the instructions in "README_setup.sh".
+#
+#  -- Instructions: --
+# If "lmp_mpi" is the name of the command you use to invoke lammps,
+# then you would run lammps on these files this way:
+
+
+lmp_mpi -i run.in.min  # minimization
+lmp_mpi -i run.in.nvt  # simulation at constant volume
+
+#(Note: The constant volume simulation lacks pressure equilibration. These are
+#       completely separate simulations. The results of the constant pressure
+#       simulation might be ignored when beginning the simulation at constant
+#       volume.  (This is because restart files in LAMMPS don't always work,
+#       and I was spending a lot of time trying to convince people it was a
+#       LAMMPS bug, instead of a moltemplate bug, so I disabled restart files.)
+#       Read the "run.in.nvt" file to find out how to use the "read_restart"
+#       command to load the results of the pressure-equilibration simulation,
+#       before beginning a constant-volume run.
+
+
+
+
+
+# If you have compiled the MPI version of lammps, you can run lammps in parallel
+#mpirun -np 4 lmp_mpi -i run.in.npt
+#mpirun -np 4 lmp_mpi -i run.in.nvt
+# (assuming you have 4 processors available)

examples/all_atom/force_field_OPLSAA/alkane_chain_single/README_setup.sh

@@ -0,0 +1,51 @@
+
+# Create LAMMPS input files this way:
+cd moltemplate_files
+
+  # run moltemplate
+
+  # This was the original (simple) way to run moltemplate:
+  # moltemplate.sh  system.lt  <-- COMMENTING OUT
+  # Instead, this is the recommended way to run moltemplate with OPLSAA:
+
+  moltemplate.sh  system.lt  -report-duplicates bytype __
+
+  # (The optional "-report-duplicates bytype __" arguments check to make
+  #  sure that there was no ambiguity in the dihedrals that were generated.
+  #  This is an issue with OPLSAA.  If there was, then moltemplate will create
+  #  a file named "warning_duplicate_dihedrals.txt".)
+  #
+  # (Note: You can also check for missing angle,dihedral params this way:)
+  # moltemplate.sh -checkff system.lt 
+
+
+  # Moltemplate generates various files with names ending in *.in* and *.data.
+  # Move them to the directory where you plan to run LAMMPS (in this case "../")
+  mv -f system.data system.in* ../
+
+  # Optional:
+  # The "./output_ttree/" directory is full of temporary files generated by
+  # moltemplate.  They can be useful for debugging, but are usually thrown away.
+  rm -rf output_ttree/
+
+  # Optional:
+  rm -f run.in.EXAMPLE  # not needed.  We have several run.in... files already.
+
+  # Optional:
+  # If any warnings or log files were generated, move them to the parent folder
+  # (so they get noticed).
+  mv -f warning*.txt ../ 2> /dev/null
+  mv -f log.* ../ 2> /dev/null
+
+cd ../
+
+
+
+
+# Optional:
+# Note: The system.data and system.in.settings files contain extra information
+# for atoms defined in OPLSAA which you are not using in this simulation.
+# This is harmless, but if you to delete this information from your
+# system.in.settings and system.in.data files, run this script:
+#
+# cleanup_moltemplate.sh

examples/all_atom/force_field_OPLSAA/alkane_chain_single/README_visualize.txt

@@ -0,0 +1,87 @@
+
+ ------- To view a lammps trajectory in VMD --------
+
+
+1) Build a PSF file for use in viewing with VMD.
+
+This step works with VMD 1.9 and topotools 1.2.
+(Older versions, like VMD 1.8.6, don't support this.)
+
+
+a) Start VMD
+b) Menu  Extensions->Tk Console
+c) Enter:
+
+(I assume that the the DATA file is called "system.data")
+
+   topo readlammpsdata system.data full
+   animate write psf system.psf
+
+2)
+
+Later, to Load a trajectory in VMD:
+
+  Start VMD
+  Select menu: File->New Molecule
+ -Browse to select the PSF file you created above, and load it.
+  (Don't close the window yet.)
+ -Browse to select the trajectory file.
+  If necessary, for "file type" select: "LAMMPS Trajectory"
+  Load it.
+
+   ----  A note on trajectory format: -----
+If the trajectory is a DUMP file, then make sure the it contains the
+information you need for pbctools (see below.  I've been using this
+command in my LAMMPS scripts to create the trajectories:
+
+  dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz
+
+It's a good idea to use an atom_style which supports molecule-ID numbers
+so that you can assign a molecule-ID number to each atom.  (I think this
+is needed to wrap atom coordinates without breaking molecules in half.)
+
+Of course, you don't have to save your trajectories in DUMP format,
+(other formats like DCD work fine)  I just mention dump files
+because these are the files I'm familiar with.
+
+3) -----  Wrap the coordinates to the unit cell
+          (without cutting the molecules in half)
+
+a) Start VMD
+b) Load the trajectory in VMD (see above)
+c) Menu  Extensions->Tk Console
+d) Try entering these commands:
+
+    pbc wrap -compound res -all
+    pbc box
+
+    ----- Optional ----
+    Sometimes the solvent or membrane obscures the view of the solute.
+    It can help to shift the location of the periodic boundary box
+    To shift the box in the y direction (for example) do this:
+
+    pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
+    pbc box -shiftcenterrel {0.0 0.15 0.0}
+
+    Distances are measured in units of box-length fractions, not Angstroms.
+
+    Alternately if you have a solute whose atoms are all of type 1,
+    then you can also try this to center the box around it:
+
+    pbc wrap -sel type=1 -all -centersel type=2 -center com
+
+4)
+    You should check if your periodic boundary conditions are too small.
+    To do that:
+       select Graphics->Representations menu option
+       click on the "Periodic" tab, and
+       click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.
+
+5) Optional: If you like, change the atom types in the PSF file so
+   that VMD recognizes the atom types, use something like:
+
+sed -e 's/   1    1      /   C    C      /g' < system.psf > temp1.psf
+sed -e 's/   2    2      /   H    H      /g' < temp1.psf  > temp2.psf
+sed -e 's/   3    3      /   P    P      /g' < temp2.psf  > system.psf
+
+(If you do this, it might effect step 2 above.)