Fixed 2 serious problems in all of the water examples: 1) replaced "f…

…ix rattle" with "fix shake" (rattle was not working again for some reason) 2) moved the "fix shake" command into the "In Settiongs" section (for compatibility with the README_remove__irrelevant_info.sh script). Hopefully now these examples are working again

README.md

@@ -23,41 +23,38 @@ 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
+[**ATB**](https://atb.uq.edu.au) molecule database, the
 [**OPLSAA**(2023)](./examples/all_atom/force_field_OPLSAA),
 [**OPLSUA**(2008)](./examples/all_atom/legacy_force_field_examples/force_field_OPLSUA_united_atom),
 [**LOPLS**(2015)](./examples/all_atom/force_field_OPLSAA/hexadecane),
 [**AMBER**(GAFF,GAFF2)](./examples/all_atom/force_field_AMBER),
 [**DREIDING**](./examples/all_atom/force_field_DREIDING),
 [**COMPASS**](./examples/all_atom/force_field_COMPASS),
 [**TraPPE**(1998)](./examples/coarse_grained/solvent_models/manybodywaterMW+hydrocarbonsTraPPE)
+[**EFF**](./examples/misc_examples/explicit_electrons/eff_CH4),
 force fields,
-the
-[**ATB**](https://atb.uq.edu.au) molecule database,
 and the
 [**MOLC**](https://pubs.rsc.org/en/content/articlelanding/2019/cp/c9cp04120f),
 [**mW**](https://doi.org/10.1021/jp805227c),
 [**ELBA**(water)](./examples/coarse_grained/solvent_models/ELBAwater%2Bmethanol),
-[**oxDNA2**](https://dna.physics.ox.ac.uk/index.php/DNA_model_introduction),
 and
-[**EFF**](./examples/misc_examples/explicit_electrons/eff_CH4)
+[**oxDNA2**](https://dna.physics.ox.ac.uk/index.php/DNA_model_introduction),
 molecular models (and others).
 (New force fields and examples are added continually by users.)
 Moltemplate is inter-operable with
 [**ATB**](https://atb.uq.edu.au),
-[**VMD/topotools**](https://www.ks.uiuc.edu/Research/vmd),
-[**PACKMOL**](http://m3g.iqm.unicamp.br/packmol/home.shtml),
-[**RED-server**](https://upjv.q4md-forcefieldtools.org),
-[**LigParGen**](https://moltemplate.org/doc/moltemplate_talk_2019-8-15.pdf#page=190),
-[**AmberTools**](https://ambermd.org/AmberTools.php),
 [**Open Babel**](https://open-babel.readthedocs.io/en/latest/FileFormats/The_LAMMPS_data_format.html),
+[**AmberTools**](https://ambermd.org/AmberTools.php),
+[**LigParGen**](https://moltemplate.org/doc/moltemplate_talk_2019-8-15.pdf#page=191),
+[**RED-server**](https://upjv.q4md-forcefieldtools.org),
+[**VMD**](https://www.ks.uiuc.edu/Research/vmd),
+[**topotools**](https://sites.google.com/site/akohlmey/software/topotools/tutorial-introduction),
+[**PACKMOL**](http://m3g.iqm.unicamp.br/packmol/home.shtml),
 [**EMC**](http://montecarlo.sourceforge.net/),
 [**CellPACK**](http://www.cellpack.org),
 [**Vipster**](https://sgsaenger.github.io/vipster),
 [**struc2lammpsdf**](https://nanohub.org/resources/struc2lammpsdf),
-and any other program that generates
-[**LAMMPS DATA**](https://docs.lammps.org/2001/data_format.html) (.lmpdat) files
-or
-[**MOL2**](http://chemyang.ccnu.edu.cn/ccb/server/AIMMS/mol2.pdf) files.
+and any other program that generates [**MOL2**](https://github.com/UnixJunkie/mol2-file-format-spec/blob/master/mol2.pdf) or [**LAMMPS DATA**](https://docs.lammps.org/2001/data_format.html) (.lmpdat) files *(by using the [mol22lt.py](https://github.com/jewettaij/moltemplate/blob/master/doc/doc_mol22lt.md) and [ltemplify.py](https://github.com/jewettaij/moltemplate/blob/master/doc/doc_ltemplify.md") file converters)*.
 
 
 This repository contains 3 folders:

examples/all_atom/force_field_AMBER/benzene/README_setup.sh

@@ -15,9 +15,6 @@ cd moltemplate_files
   # Move them to the directory where you plan to run LAMMPS (in this case "../")
   mv -f system.data system.in* ../
 
-  # Optional: Move the "log.cite" file (containing links to relevant papers).
-  mv -f log.cite* ../
-
   # Optional:
   # The "./output_ttree/" directory is full of temporary files generated by
   # moltemplate.  They can be useful for debugging, but are usually thrown away.

examples/all_atom/force_field_AMBER/ethylene+benzene/README_setup.sh

@@ -15,10 +15,6 @@ cd moltemplate_files
   mv -f system.data system.in* ../
 
 
-  # Optional: Move the "log.cite" file (containing links to relevant papers).
-  mv -f log.cite* ../
-
-
   # Optional:
   # The "./output_ttree/" directory is full of temporary files generated by
   # moltemplate.  They can be useful for debugging, but are usually thrown away.

examples/all_atom/force_field_AMBER/hexadecane/README_setup.sh

@@ -16,10 +16,6 @@ cd moltemplate_files
   mv -f system.data system.in* ../
 
 
-  # Optional: Move the "log.cite" file (containing links to relevant papers).
-  mv -f log.cite* ../
-
-
   # Optional:
   # The "./output_ttree/" directory is full of temporary files generated by
   # moltemplate.  They can be useful for debugging, but are usually thrown away.

examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/README_setup.sh

@@ -10,9 +10,6 @@ cd moltemplate_files
   # Move them to the directory where you plan to run LAMMPS (in this case "../")
   mv -f system.data system.in* ../
 
-  # Optional: Move the "log.cite" file (containing links to relevant papers).
-  mv -f log.cite* ../
-
   # Optional:
   # The "./output_ttree/" directory is full of temporary files generated by
   # moltemplate.  They can be useful for debugging, but are usually thrown away.

examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/run.in.npt

@@ -9,6 +9,15 @@
 
 include     "system.in.init"         # specify the style of force field used
 
+# Weird LAMMPS issue:
+# The "system.in.init" file contains definitions for bond, angle, dihedral,
+# and improper interactions that were defined in the AMBER/GAFF force-field.
+# But since the molecules in this example do not have any improper interactions,
+# LAMMPS will object and not allow us to run the simulation.  
+# To get around this, tell LAMMPS to disable improper interactions:
+
+improper_style none
+
 # ------------------------------- Atom Definition Section -------------------
 
 read_data   "system.data"            # load atom coordinates and topology
@@ -22,10 +31,9 @@ include      "system.in.settings"    # load the force field parameters
 
 # -- minimization protocol --
 
-minimize 1.0e-4 1.0e-6 100000 400000
+unfix fShakeTIP3P    # Disable SHAKE during minimization and pressure equilibr
 
-# (Apply constraints after minimization and after other fixes which exert force)
-include         system.in.constraints
+minimize 1.0e-4 1.0e-6 100000 400000
 
 
 # -- simulation protocol --
@@ -42,9 +50,6 @@ thermo_style custom step temp ke pe etotal epair ebond eangle edihed press vol
 
 fix          fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0
 
-include      "system.in.constraints" # apply constraints (after minimization
-                                     # and after all integration fixes)
-
-run          200000
+run          50000
 
 write_data   "system_after_npt.data"

examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/run.in.nvt

@@ -12,6 +12,15 @@
 
 include      "system.in.init"        # specify the style of force field used
 
+# Weird LAMMPS issue:
+# The "system.in.init" file contains definitions for bond, angle, dihedral,
+# and improper interactions that were defined in the AMBER/GAFF force-field.
+# But since the molecules in this example do not have any improper interactions,
+# LAMMPS will object and not allow us to run the simulation.  
+# To get around this, tell LAMMPS to disable improper interactions:
+
+improper_style none
+
 # ------------------------------- Atom Definition Section -------------------
 
 # Read the coordinates generated by an earlier NPT simulation
@@ -30,7 +39,7 @@ include      "system.in.settings"    # load the force field parameters
 # -- simulation protocol --
 
 
-timestep     1.0
+timestep     2.0
 
 thermo       500
 #thermo_modify  flush yes
@@ -39,9 +48,6 @@ dump         1 all custom 2000 traj_nvt.lammpstrj id mol type x y z ix iy iz
 
 fix          fxnvt all nvt temp 300.0 300.0 500.0 tchain 1
 
-include      "system.in.constraints" # apply constraints (after minimization
-                                     # and after all integration fixes)
-
-run          50000
+run          1000000
 
 write_data   "system_after_nvt.data"

examples/all_atom/force_field_OPLSAA/waterSPCE+methane/moltemplate_files/spce_oplsaa.lt

@@ -36,18 +36,13 @@ SPCE {
     angle_coeff  @angle:HOH       75.0    109.47
     pair_coeff   @atom:O @atom:O  0.1553  3.166 
     pair_coeff   @atom:H @atom:H  0.0     0.0
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
-    # parameters ("600.0", "75.0") do not matter. But the equilibrium bond
-    # length ("1.0") and equilibrium angle ("109.47") does matter. LAMMPS
-    # obtains these numbers from the bond_coeff and angle_coeff commands above.
-  }
 
-  write_once("In Constraints") {
+    # You can also define fixes and groups here as well
     group spce type  @atom:O  @atom:H
-    fix fRattleSPCE spce rattle 0.0001 10 100 b @bond:OH a @angle:HOH
-    # Remember to put this command in your LAMMPS input script:
-    #  include system.in.constraints
-    # ...after minimization and after all integration fixes
+    fix fShakeSPCE spce shake 0.0001 10 100 b @bond:OH a @angle:HOH
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeSPCE   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
   write_once("In Init") {
@@ -57,12 +52,6 @@ SPCE {
     bond_style   harmonic
     angle_style  harmonic
     kspace_style pppm 0.0001
-    # Don't specify the pair_style.  Why?
-    #   oplsaa2008.lt uses this pair_style:
-    #     pair_style   lj/cut/coul/long 11.0 11.0
-    #   oplsaa.lt (2023) uses this pair_style:
-    #     pair_style lj/charmm/coul/long 9.0 11.0
-    # pair_modify  mix arithmetic  # LEAVE THIS UNSPECIFIED!
   }
 
 } # end of definition of "SPCE" water molecule type

examples/all_atom/force_field_OPLSAA/waterSPCE+methane/run.in.npt

@@ -23,7 +23,8 @@ include      "system.in.charges"     # load the charge of each atom
 
 # -- minimization protocol --
 
-thermo        50
+unfix fShakeSPCE    # Disable SHAKE during minimization and pressure equilibr
+
 minimize 1.0e-4 1.0e-6 100000 400000
 
 # Optional: write the coordinates after minimization
@@ -36,11 +37,6 @@ write_data   "system_after_min.data"
 timestep     1.0
 dump         1 all custom 5000 traj_npt.lammpstrj id mol type x y z ix iy iz
 
-# COMMENTING OUT SHAKE/RATTLE CONSTRAINTS DURING PRESSURE EQUILIBRATION
-# (Pressure and SHAKE/RATTLE are currently not working well together (2024/11).
-#  Once the pressure has equilibrated, we will turn them on.  See "run.in.nvt")
-#include      "system.in.constraints" # apply constraints (after minimization but
-#                                     # before fixes that effect the box size)
 
 # Turn on the barostat and thermostat
 fix          fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0

examples/all_atom/force_field_OPLSAA/waterSPCE+methane/run.in.nvt

@@ -38,11 +38,8 @@ dump         1 all custom 2000 traj_nvt.lammpstrj id mol type x y z ix iy iz
 
 fix          fxnvt all nvt temp 300.0 300.0 500.0 tchain 1
 
-include      "system.in.constraints" # apply constraints (after minimization
-                                     # and after all integration fixes)
-
 #thermo_modify  flush yes
 
-run          500000
+run          1000000
 
 write_data   "system_after_nvt.data"

examples/all_atom/force_field_explicit_parameters/nanotube+water/moltemplate_files/spce.lt

@@ -33,19 +33,12 @@ SPCE {
     angle_coeff  @angle:HOH       harmonic          75.0    109.47
     pair_coeff   @atom:O @atom:O  lj/cut/coul/long  0.1553  3.166 
     pair_coeff   @atom:H @atom:H  lj/cut/coul/long  0.0     0.0
-    group spce type  @atom:O  @atom:H
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
-    # parameters ("600.0", "75.0") do not matter. But the equilibrium bond
-    # length ("1.0") and equilibrium angle ("109.47") does matter. LAMMPS
-    # obtains these numbers from the bond_coeff and angle_coeff commands above.
-  }
 
-  write_once("In Constraints") {
     group spce type  @atom:O  @atom:H
-    fix fRattleSPCE spce rattle 0.0001 10 100 b @bond:OH a @angle:HOH
-    # Remember to put this command in your LAMMPS input script:
-    #  include system.in.constraints
-    # ...after minimization and after all integration fixes
+    fix fShakeSPCE spce shake 0.0001 10 100 b @bond:OH a @angle:HOH
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeSPCE   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
   write_once("In Init") {

examples/all_atom/force_field_explicit_parameters/nanotube+water/run.in.npt

@@ -63,6 +63,8 @@ print "--------- beginning simulation (using fix nvt) ---------"
 
 dump            1 all custom 1000 traj_npt.lammpstrj id mol type x y z ix iy iz
 
+unfix fShakeSPCE    # Disable SHAKE during minimization and pressure equilibr
+
 thermo_style    custom step temp pe etotal press vol epair ebond eangle edihed
 thermo          200  # time interval for printing out "thermo" data
 
@@ -86,16 +88,11 @@ fix Ffreezestuff Cgraphene rigid single force * off off off torque * off off off
 #  so this may no longer be true.  Please use this example with caution.)
 
 
-
 # Thermostat+Barostat
 # Set temp=300K, pressure=200bar, and equilibrate volume only in the z direction
 
 fix fxMoveStuff mobile npt temp 300 300 100 z 200 200 1000.0 dilate mobile drag 2.0
 
-include      "system.in.constraints" # apply constraints (after minimization
-                                     # and after all integration fixes)
-
-
 #  ----------------------------------------
 
 # The next two lines recalculate the temperature using
@@ -112,31 +109,28 @@ reset_timestep 0
 
 timestep 0.25
 
-run    100000
+run    20000
 
 timestep 0.5
 
-run    200000
+run    50000
 
 
-# Hopefully the barostat is no longer oscillating. Increase the timestep and
-# also get get rid of "drag 2.0". (A non-zero drag parameter will result in
-# unrealistic fluctuations of volume under NPT conditions.)
-# drag 2.0  <-- commenting out
-#
-# Set temp=300K, pressure=0bar, and equilibrate volume only in the z direction
+# Hopefully the barostat is no longer oscillating.  Then we can try using fix
+# npt again and omit the "drag 2.0" argument. (A non-zero drag parameter will
+# result in unrealistic fluctuations of volume under NPT conditions.)
+# First, undo the previous "fix npt"
 
 unfix fxMoveStuff
-unfix fRattleSPCE  # (must undo all the fixes from "system.in.constraints")
+
+# Then re-apply "fix npt"
+# Set temp=300K, pressure=0bar, and equilibrate volume only in the z direction
 fix fxMoveStuff mobile npt temp 300 300 100 z 0 0 1000.0 dilate mobile
 fix_modify fxMoveStuff temp tempMobile
-include      "system.in.constraints" # apply constraints AFTER integration fixes
-                                     # (such as "fix npt" used above)
-                                     # https://lammps.sandia.gov/doc/fix_shake.html
 
-timestep 0.5
+timestep 0.5  # (Perhaps we can increase this timestep to 1.0?  I'm not sure.)
 
-run    1000000
+run    5000000
 
 write_data  system_after_npt.data
 

examples/all_atom/force_field_explicit_parameters/nanotube+water/run.in.nvt

@@ -47,9 +47,6 @@ thermo          500  # time interval for printing out "thermo" data
 # Integrate the equations of motion:
 fix fxMoveStuff mobile nvt temp 300.0 300.0 100.0
 
-include      "system.in.constraints" # apply constraints (after minimization
-                                     # and after all integration fixes)
-
 # The next two lines recalculate the temperature 
 # using only the mobile degrees of freedom:
 

examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/moltemplate_files/spce.lt

@@ -6,7 +6,7 @@
 
 SPCE {
 
-  # AtomID MoleculeID AtomType charge    X       Y       Z
+  # AtomID MoleculeID AtomType Charge    X       Y       Z
 
   write("Data Atoms") {
     $atom:o   $mol:w  @atom:O  -0.8476  0.0000000 0.00000 0.000000
@@ -33,26 +33,20 @@ SPCE {
     angle_coeff  @angle:HOH       75.0    109.47
     pair_coeff   @atom:O @atom:O  0.1553  3.166 
     pair_coeff   @atom:H @atom:H  0.0     0.0
-    group spce type  @atom:O  @atom:H
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
-    # parameters ("600.0", "75.0") do not matter. But the equilibrium bond
-    # length ("1.0") and equilibrium angle ("109.47") does matter. LAMMPS
-    # obtains these numbers from the bond_coeff and angle_coeff commands above.
-  }
 
-  write_once("In Constraints") {
     group spce type  @atom:O  @atom:H
-    fix fRattleSPCE spce rattle 0.0001 10 100 b @bond:OH a @angle:HOH
-    # Remember to put this command in your LAMMPS input script:
-    #  include system.in.constraints
-    # ...after minimization and after all integration fixes
+    fix fShakeSPCE spce shake 0.0001 10 100 b @bond:OH a @angle:HOH
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeSPCE   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
   write_once("In Init") {
     # -- Default styles (for solo "SPCE" water) --
     units        real
     atom_style   full
-    pair_style   lj/charmm/coul/long 9.0 10.0
+    # (Hybrid force fields were not necessary but are used for portability.)
+    pair_style   lj/cut/coul/long 9.0 10.0
     bond_style   harmonic
     angle_style  harmonic
     kspace_style pppm 0.0001