diff --git a/README.md b/README.md
index ef1af3e6..0caa520b 100644
--- a/README.md
+++ b/README.md
@@ -23,6 +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
+[**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),
@@ -30,34 +31,30 @@ It currently supports the
 [**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:
diff --git a/examples/all_atom/force_field_AMBER/benzene/README_setup.sh b/examples/all_atom/force_field_AMBER/benzene/README_setup.sh
index cd04e94d..5cd2142a 100755
--- a/examples/all_atom/force_field_AMBER/benzene/README_setup.sh
+++ b/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.
diff --git a/examples/all_atom/force_field_AMBER/ethylene+benzene/README_setup.sh b/examples/all_atom/force_field_AMBER/ethylene+benzene/README_setup.sh
index 94804451..fdfb19eb 100755
--- a/examples/all_atom/force_field_AMBER/ethylene+benzene/README_setup.sh
+++ b/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.
diff --git a/examples/all_atom/force_field_AMBER/hexadecane/README_setup.sh b/examples/all_atom/force_field_AMBER/hexadecane/README_setup.sh
index 5ff0bf3f..e8d8be3f 100755
--- a/examples/all_atom/force_field_AMBER/hexadecane/README_setup.sh
+++ b/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.
diff --git a/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/README_setup.sh b/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/README_setup.sh
index 0dd730b1..4fc8dda6 100755
--- a/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/README_setup.sh
+++ b/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.
diff --git a/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/run.in.npt b/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/run.in.npt
index 8b43671f..05661412 100644
--- a/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/run.in.npt
+++ b/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"
diff --git a/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/run.in.nvt b/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/run.in.nvt
index c36cf8db..48d25d54 100644
--- a/examples/all_atom/force_field_AMBER/waterTIP3P+isobutane/run.in.nvt
+++ b/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"
diff --git a/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/moltemplate_files/spce_oplsaa.lt b/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/moltemplate_files/spce_oplsaa.lt
index ad66463d..38effff0 100644
--- a/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/moltemplate_files/spce_oplsaa.lt
+++ b/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
diff --git a/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/run.in.npt b/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/run.in.npt
index 08703ee3..f53ba541 100644
--- a/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/run.in.npt
+++ b/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
diff --git a/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/run.in.nvt b/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/run.in.nvt
index 5f82a063..781f28b0 100644
--- a/examples/all_atom/force_field_OPLSAA/waterSPCE+methane/run.in.nvt
+++ b/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"
diff --git a/examples/all_atom/force_field_explicit_parameters/nanotube+water/moltemplate_files/spce.lt b/examples/all_atom/force_field_explicit_parameters/nanotube+water/moltemplate_files/spce.lt
index 234e0f40..9e61b3a5 100644
--- a/examples/all_atom/force_field_explicit_parameters/nanotube+water/moltemplate_files/spce.lt
+++ b/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") {
diff --git a/examples/all_atom/force_field_explicit_parameters/nanotube+water/run.in.npt b/examples/all_atom/force_field_explicit_parameters/nanotube+water/run.in.npt
index 9e8389bc..94abde9d 100644
--- a/examples/all_atom/force_field_explicit_parameters/nanotube+water/run.in.npt
+++ b/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
 
diff --git a/examples/all_atom/force_field_explicit_parameters/nanotube+water/run.in.nvt b/examples/all_atom/force_field_explicit_parameters/nanotube+water/run.in.nvt
index 6cdb1a68..851e868d 100644
--- a/examples/all_atom/force_field_explicit_parameters/nanotube+water/run.in.nvt
+++ b/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:
 
diff --git a/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/moltemplate_files/spce.lt b/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/moltemplate_files/spce.lt
index fa4426b7..b5c4fa83 100644
--- a/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/moltemplate_files/spce.lt
+++ b/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
diff --git a/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/run.in.npt b/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/run.in.npt
index b50bd64f..f189f5df 100644
--- a/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/run.in.npt
+++ b/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/run.in.npt
@@ -21,20 +21,20 @@ include      "system.in.settings"    # load the force field parameters
 
 # -- minimization protocol --
 
+unfix fShakeSPCE    # Disable SHAKE during minimization and pressure equilibr
 
 # Note: The minimization step is not necessary in this example.  However
 #       in general, it's always a good idea to minimize the system beforehand.
 minimize 1.0e-3 1.0e-5 100000 400000
 
+
 # -- simulation protocol --
 
 
-timestep     1.0
+timestep     0.8
 dump         1 all custom 200 traj_npt.lammpstrj id mol type x y z ix iy iz
-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)
+fix          fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0
 
 thermo       100
 
@@ -43,7 +43,7 @@ run          10000
 # Now that the system's temperature has become more equilibrated,
 # we can increase the timestep:
 
-timestep     2.0
-run          50000
+timestep     1.0   # (Note: You can increase this to 2.0 if SHAKE/RATTLE is on.)
+#run          50000
 
 write_data   system_after_npt.data
diff --git a/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/run.in.nvt b/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/run.in.nvt
index c43221a9..120b3284 100644
--- a/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/run.in.nvt
+++ b/examples/all_atom/force_field_explicit_parameters/waterSPCE+Na+Cl/run.in.nvt
@@ -35,9 +35,6 @@ timestep     1.0
 dump         1 all custom 500 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       500
 
 run          10000
diff --git a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/alkane_chain_single/moltemplate_files/run.in.EXAMPLE b/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/alkane_chain_single/moltemplate_files/run.in.EXAMPLE
deleted file mode 100644
index 28e373bc..00000000
--- a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/alkane_chain_single/moltemplate_files/run.in.EXAMPLE
+++ /dev/null
@@ -1,19 +0,0 @@
-
-
-# ----------------- Init Section -----------------
-
-include "system.in.init"
-
-
-# ----------------- Atom Definition Section -----------------
-
-read_data "system.data"
-
-# ----------------- Settings Section -----------------
-
-include "system.in.settings"
-
-
-# ----------------- Charges Section -----------------
-
-include "system.in.charges"
diff --git a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/butane/moltemplate_files/run.in.EXAMPLE b/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/butane/moltemplate_files/run.in.EXAMPLE
deleted file mode 100644
index 28e373bc..00000000
--- a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/butane/moltemplate_files/run.in.EXAMPLE
+++ /dev/null
@@ -1,19 +0,0 @@
-
-
-# ----------------- Init Section -----------------
-
-include "system.in.init"
-
-
-# ----------------- Atom Definition Section -----------------
-
-read_data "system.data"
-
-# ----------------- Settings Section -----------------
-
-include "system.in.settings"
-
-
-# ----------------- Charges Section -----------------
-
-include "system.in.charges"
diff --git a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/ethylene+benzene/moltemplate_files/run.in.EXAMPLE b/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/ethylene+benzene/moltemplate_files/run.in.EXAMPLE
deleted file mode 100644
index 28e373bc..00000000
--- a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/ethylene+benzene/moltemplate_files/run.in.EXAMPLE
+++ /dev/null
@@ -1,19 +0,0 @@
-
-
-# ----------------- Init Section -----------------
-
-include "system.in.init"
-
-
-# ----------------- Atom Definition Section -----------------
-
-read_data "system.data"
-
-# ----------------- Settings Section -----------------
-
-include "system.in.settings"
-
-
-# ----------------- Charges Section -----------------
-
-include "system.in.charges"
diff --git a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/functionalized_nanotubes_NH2/moltemplate_files/run.in.EXAMPLE b/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/functionalized_nanotubes_NH2/moltemplate_files/run.in.EXAMPLE
deleted file mode 100644
index 28e373bc..00000000
--- a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/functionalized_nanotubes_NH2/moltemplate_files/run.in.EXAMPLE
+++ /dev/null
@@ -1,19 +0,0 @@
-
-
-# ----------------- Init Section -----------------
-
-include "system.in.init"
-
-
-# ----------------- Atom Definition Section -----------------
-
-read_data "system.data"
-
-# ----------------- Settings Section -----------------
-
-include "system.in.settings"
-
-
-# ----------------- Charges Section -----------------
-
-include "system.in.charges"
diff --git a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/methane/moltemplate_files/run.in.EXAMPLE b/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/methane/moltemplate_files/run.in.EXAMPLE
deleted file mode 100644
index 28e373bc..00000000
--- a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/methane/moltemplate_files/run.in.EXAMPLE
+++ /dev/null
@@ -1,19 +0,0 @@
-
-
-# ----------------- Init Section -----------------
-
-include "system.in.init"
-
-
-# ----------------- Atom Definition Section -----------------
-
-read_data "system.data"
-
-# ----------------- Settings Section -----------------
-
-include "system.in.settings"
-
-
-# ----------------- Charges Section -----------------
-
-include "system.in.charges"
diff --git a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/waterSPCE+methane/moltemplate_files/run.in.EXAMPLE b/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/waterSPCE+methane/moltemplate_files/run.in.EXAMPLE
deleted file mode 100644
index af494be0..00000000
--- a/examples/all_atom/legacy_force_field_examples/force_field_OPLSAA_2008/waterSPCE+methane/moltemplate_files/run.in.EXAMPLE
+++ /dev/null
@@ -1,23 +0,0 @@
-
-
-# ----------------- Init Section -----------------
-
-include "system.in.init"
-
-
-# ----------------- Atom Definition Section -----------------
-
-read_data "system.data"
-
-# ----------------- Settings Section -----------------
-
-include "system.in.settings"
-
-
-# ----------------- Charges Section -----------------
-
-include "system.in.charges"
-
-# ----------------- Constraints Section -----------------
-
-include "system.in.constraints"
diff --git a/moltemplate/force_fields/spc_oplsaa2008.lt b/moltemplate/force_fields/spc_oplsaa2008.lt
index 2f458aff..26d11014 100644
--- a/moltemplate/force_fields/spc_oplsaa2008.lt
+++ b/moltemplate/force_fields/spc_oplsaa2008.lt
@@ -40,21 +40,24 @@ SPC_oplsaa inherits OPLSAA {
     $bond:OH2  $atom:O $atom:H2
   }
 
-  write_once("In Constraints") {
+  write_once("In Settings") {
 
     # Define a group for the spc water molecules. (We will use "spc" below.)
     group spc type  @atom:76  @atom:77
 
     # Constrain the angles and distances in spc water:
-    # (See "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    # (See "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
 
-    fix fRattleSPC spc rattle 0.0001 10 100 b @bond:042_043 a @angle:043_042_043
+    fix fShakeSPC spc shake 0.0001 10 100 b @bond:042_043 a @angle:043_042_043
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeSPC   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
 
     # Note:    Bonds of type "@bond:042_043"
     #  and    Angles of type "@angle:043_042_043"
     # ...are defined in the "oplsaa.lt" file. These codes are used to lookup the
     # OH bond length (1.0 Angstroms) and HOH angle (109.47°) used in SPC water.
-    # (Also see "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    # (Also see "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
 
     # ------- Remember to put this command in your LAMMPS input script! -------
     #
diff --git a/moltemplate/force_fields/spce.lt b/moltemplate/force_fields/spce.lt
index 586ffb87..0bba4f7d 100644
--- a/moltemplate/force_fields/spce.lt
+++ b/moltemplate/force_fields/spce.lt
@@ -33,18 +33,12 @@ 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") {
diff --git a/moltemplate/force_fields/spce_amber.lt b/moltemplate/force_fields/spce_amber.lt
index 60d7274a..563bf9aa 100644
--- a/moltemplate/force_fields/spce_amber.lt
+++ b/moltemplate/force_fields/spce_amber.lt
@@ -35,18 +35,15 @@ 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.
+    # You can also define fixes and groups here as well:
+    group spce type  @atom:O  @atom:H
+    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 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
+  write_once("In Settings") {
   }
 
   write_once("In Init") {
diff --git a/moltemplate/force_fields/spce_dreiding.lt b/moltemplate/force_fields/spce_dreiding.lt
index 936ae7f0..abe240b6 100644
--- a/moltemplate/force_fields/spce_dreiding.lt
+++ b/moltemplate/force_fields/spce_dreiding.lt
@@ -62,19 +62,12 @@ SPCE {
     angle_coeff  @angle:HOH                        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
+    # You can also define fixes and groups here as well:
     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") {
diff --git a/moltemplate/force_fields/spce_hybrid.lt b/moltemplate/force_fields/spce_hybrid.lt
index 83b5adcb..1006b015 100644
--- a/moltemplate/force_fields/spce_hybrid.lt
+++ b/moltemplate/force_fields/spce_hybrid.lt
@@ -43,18 +43,12 @@ SPCE {
     angle_coeff  @angle:HOH   harmonic                 75.0    109.47
     pair_coeff   @atom:O @atom:O  lj/charmm/coul/long  0.1553  3.166 
     pair_coeff   @atom:H @atom:H  lj/charmm/coul/long  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") {
diff --git a/moltemplate/force_fields/spce_more_comments.lt b/moltemplate/force_fields/spce_more_comments.lt
index 2f855d0c..1cbfa2d3 100644
--- a/moltemplate/force_fields/spce_more_comments.lt
+++ b/moltemplate/force_fields/spce_more_comments.lt
@@ -55,7 +55,7 @@ SPCE {
   # Sometimes it is convenient to nest them within a molecule's definition.
   # The "In Settings" section contains LAMMPS commands which typically define 
   # force-field parameters (eg "pair_coeff", "bond_coeff", "angle_coeff") and
-  # constraints (eg "fix rattle") which are specific to this type of molecule.
+  # constraints (eg "fix shake") which are specific to this type of molecule.
   # (For more complicated simulations, it's convenient to define a "ForceField"
   # object containing force field parameters shared by many different types of 
   # molecules. Since this is a simple example, there is no need to do that.)
@@ -70,13 +70,11 @@ SPCE {
     pair_coeff   @atom:O @atom:O    0.1553  3.166 
     pair_coeff   @atom:H @atom:H    0.0     2.058
     # explanation:  http://lammps.sandia.gov/doc/pair_charmm.html
-  }
 
-  write_once("In Constraints") {
     # OPTIONAL: We want the water molecules to behave like rigid objects, so 
-    # we apply the "RATTLE" constraint to the group of atoms in water molecules.
-    # (RATTLE is used to make bonds and angles rigid.)
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
+    # we apply the "SHAKE" constraint to the group of atoms in water molecules.
+    # (SHAKE is used to make bonds and angles rigid.)
+    # Note: Since we are using SHAKE 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.
@@ -84,11 +82,12 @@ SPCE {
     group spce type  @atom:O  @atom:H
     # Explanation:  http://lammps.sandia.gov/doc/group.html
 
-    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.
-    # Explanation:  http://lammps.sandia.gov/doc/fix_shake.html
+    # You can also define fixes and groups here as well:
+    group spce type  @atom:O  @atom:H
+    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)
   }
 
   # LAMMPS requires that some commands must be issued in a certain order.
diff --git a/moltemplate/force_fields/spce_oplsaa.lt b/moltemplate/force_fields/spce_oplsaa.lt
index ad66463d..38effff0 100644
--- a/moltemplate/force_fields/spce_oplsaa.lt
+++ b/moltemplate/force_fields/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
diff --git a/moltemplate/force_fields/tip3p_1983.lt b/moltemplate/force_fields/tip3p_1983.lt
index 95c01098..51bed4bc 100644
--- a/moltemplate/force_fields/tip3p_1983.lt
+++ b/moltemplate/force_fields/tip3p_1983.lt
@@ -47,7 +47,7 @@ TIP3P_1983 {
   write_once("In Settings") {
     bond_coeff   @bond:OH        450.0  0.9572
     angle_coeff  @angle:HOH      55.0   104.52
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
+    # Note: Since we are using SHAKE constraints, the bond and angle strength
     # parameters ("450.0", "55.0") do not matter.  But the equilibrium bond
     # length ("0.9572") and equilibrium angle ("104.52") does matter.  LAMMPS
     # obtains these numbers from the bond_coeff and angle_coeff commands.
@@ -83,18 +83,18 @@ TIP3P_1983 {
   }
 
 
-  write_once("In Constraints") {
+  write_once("In Settings") {
     # Optional: Constrain the angles and distances.
     #           (Most implementations use this, but it is optional.)
     # First, define a group for the tip3p water molecules:
     group tip3p type  @atom:O  @atom:H
 
-    # Then use "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
-    fix fRattleTIP3P tip3p rattle 0.0001 10 100 b @bond:OH a @angle:HOH
+    # Then use "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    fix fShakeTIP3P tip3p shake 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.
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeTIP3P   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
 
diff --git a/moltemplate/force_fields/tip3p_1983_charmm.lt b/moltemplate/force_fields/tip3p_1983_charmm.lt
index eca59185..af7308a8 100644
--- a/moltemplate/force_fields/tip3p_1983_charmm.lt
+++ b/moltemplate/force_fields/tip3p_1983_charmm.lt
@@ -48,7 +48,7 @@ TIP3P_1983_CHARMM {
     bond_coeff   @bond:OH         450.0  0.9572
     angle_coeff  @angle:HOH       55.0   104.52
 
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
+    # Note: Since we are using SHAKE constraints, the bond and angle strength
     # parameters ("450.0", "55.0") do not matter.  But the equilibrium bond
     # length ("0.9572") and equilibrium angle ("104.52") does matter.  LAMMPS
     # obtains these numbers from the bond_coeff and angle_coeff commands.
@@ -84,18 +84,18 @@ TIP3P_1983_CHARMM {
   }
 
 
-  write_once("In Constraints") {
+  write_once("In Settings") {
     # Optional: Constrain the angles and distances.
     #           (Most implementations use this, but it is optional.)
     # First, define a group for the tip3p water molecules:
     group tip3p type  @atom:O  @atom:H
 
-    # Then use "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
-    fix fRattleTIP3P tip3p rattle 0.0001 10 100 b @bond:OH a @angle:HOH
+    # Then use "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    fix fShakeTIP3P tip3p shake 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.
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeTIP3P   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
 
diff --git a/moltemplate/force_fields/tip3p_1983_charmm_hybrid.lt b/moltemplate/force_fields/tip3p_1983_charmm_hybrid.lt
index fa221019..b9a6090b 100644
--- a/moltemplate/force_fields/tip3p_1983_charmm_hybrid.lt
+++ b/moltemplate/force_fields/tip3p_1983_charmm_hybrid.lt
@@ -93,18 +93,18 @@ TIP3P_1983_CHARMM {
   }
 
 
-  write_once("In Constraints") {
+  write_once("In Settings") {
     # Optional: Constrain the angles and distances.
     #           (Most implementations use this, but it is optional.)
     # First, define a group for the tip3p water molecules:
     group tip3p type  @atom:O  @atom:H
 
-    # Then use "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
-    fix fRattleTIP3P tip3p rattle 0.0001 10 100 b @bond:OH a @angle:HOH
+    # Then use "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    fix fShakeTIP3P tip3p shake 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.
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeTIP3P   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
 
diff --git a/moltemplate/force_fields/tip3p_1983_hybrid.lt b/moltemplate/force_fields/tip3p_1983_hybrid.lt
index fd1b6980..91cfffa1 100644
--- a/moltemplate/force_fields/tip3p_1983_hybrid.lt
+++ b/moltemplate/force_fields/tip3p_1983_hybrid.lt
@@ -57,7 +57,7 @@ TIP3P_1983 {
     bond_coeff   @bond:OH         harmonic       450.0  0.9572
     angle_coeff  @angle:HOH       harmonic       55.0   104.52
 
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
+    # Note: Since we are using SHAKE constraints, the bond and angle strength
     # parameters ("450.0", "55.0") do not matter.  But the equilibrium bond
     # length ("0.9572") and equilibrium angle ("104.52") does matter.  LAMMPS
     # obtains these numbers from the bond_coeff and angle_coeff commands.
@@ -93,18 +93,18 @@ TIP3P_1983 {
   }
 
 
-  write_once("In Constraints") {
+  write_once("In Settings") {
     # Optional: Constrain the angles and distances.
     #           (Most implementations use this, but it is optional.)
     # First, define a group for the tip3p water molecules:
     group tip3p type  @atom:O  @atom:H
 
-    # Then use "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
-    fix fRattleTIP3P tip3p rattle 0.0001 10 100 b @bond:OH a @angle:HOH
+    # Then use "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    fix fShakeTIP3P tip3p shake 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.
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeTIP3P   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
 
diff --git a/moltemplate/force_fields/tip3p_1983_oplsaa.lt b/moltemplate/force_fields/tip3p_1983_oplsaa.lt
index 5f815589..3f1e93e8 100644
--- a/moltemplate/force_fields/tip3p_1983_oplsaa.lt
+++ b/moltemplate/force_fields/tip3p_1983_oplsaa.lt
@@ -49,27 +49,26 @@ TIP3P_1983_oplsaa inherits OPLSAA {
     $bond:oh2  $atom:o $atom:h2
   }
 
-  write_once("In Constraints") {
+  write_once("In Settings") {
 
     # Define a group for the tip3p water molecules. (We will use "tip3p" below.)
     group tip3p type  @atom:63  @atom:64
 
     # Optional: Constrain the angles and distances.
     #           (Most implementations use this, but it is optional.)
-    # (See "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
-    fix fRattleTIP3P tip3p rattle 0.0001 10 100 b @bond:031_032 a @angle:032_031_032
+    # (See "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    fix fShakeTIP3P tip3p shake 0.0001 10 100 b @bond:tipO_tipH a @angle:tipH_tipO_tipH
 
-    # Note:    Bonds of type "@bond:031_032"
-    #  and    Angles of type "@angle:032_031_032"
+    # Note:    Bonds of type "@bond:tipO_tipH"
+    #  and    Angles of type "@angle:tipH_tipO_tipH"
     # ...are defined in the "oplsaa.lt" file. These codes are used to lookup the
     # OH bond length (0.9572 Angstroms) and HOH angle (104.52°) used in
     # the TIP3P water model.
 
 
-
-    # Remember to put this command in your LAMMPS input script:
-    #  include system.in.constraints
-    # ...AFTER minimization and after all integration fixes.
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeTIP3P   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
 } # "TIP3P_1983_oplsaa" water molecule type
diff --git a/moltemplate/force_fields/tip3p_2004.lt b/moltemplate/force_fields/tip3p_2004.lt
index 7e3f1046..89810903 100644
--- a/moltemplate/force_fields/tip3p_2004.lt
+++ b/moltemplate/force_fields/tip3p_2004.lt
@@ -48,7 +48,7 @@ TIP3P_2004 {
     bond_coeff   @bond:OH         450.0  0.9572
     angle_coeff  @angle:HOH       55.0   104.52
 
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
+    # Note: Since we are using SHAKE constraints, the bond and angle strength
     # parameters ("450.0", "55.0") do not matter.  But the equilibrium bond
     # length ("0.9572") and equilibrium angle ("104.52") does matter.  LAMMPS
     # obtains these numbers from the bond_coeff and angle_coeff commands.
@@ -84,18 +84,18 @@ TIP3P_2004 {
   }
 
 
-  write_once("In Constraints") {
+  write_once("In Settings") {
     # Optional: Constrain the angles and distances.
     #           (Most implementations use this, but it is optional.)
     # First, define a group for the tip3p water molecules:
     group tip3p type  @atom:O  @atom:H
 
-    # Then use "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
-    fix fRattleTIP3P tip3p rattle 0.0001 10 100 b @bond:OH a @angle:HOH
+    # Then use "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    fix fShakeTIP3P tip3p shake 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.
+    # It's a good idea to put this command in your LAMMPS input scripts:
+    #  unfix fShakeTIP3P   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }
 
 
diff --git a/moltemplate/force_fields/tip3p_2004_hybrid.lt b/moltemplate/force_fields/tip3p_2004_hybrid.lt
index c0fc6816..cc862bbd 100644
--- a/moltemplate/force_fields/tip3p_2004_hybrid.lt
+++ b/moltemplate/force_fields/tip3p_2004_hybrid.lt
@@ -59,11 +59,6 @@ TIP3P_2004 {
     bond_coeff   @bond:OH         harmonic       450.0  0.9572
     angle_coeff  @angle:HOH       harmonic       55.0   104.52
 
-    # Note: Since we are using RATTLE constraints, the bond and angle strength
-    # parameters ("450.0", "55.0") do not matter.  But the equilibrium bond
-    # length ("0.9572") and equilibrium angle ("104.52") does matter.  LAMMPS
-    # obtains these numbers from the bond_coeff and angle_coeff commands.
-
     #########################################################################
     ####  There are two choices for for the O-O interactions
     #########################################################################
@@ -95,18 +90,17 @@ TIP3P_2004 {
   }
 
 
-  write_once("In Constraints") {
+  write_once("In Settings") {
     # Optional: Constrain the angles and distances.
     #           (Most implementations use this, but it is optional.)
     # First, define a group for the tip3p water molecules:
     group tip3p type  @atom:O  @atom:H
 
-    # Then use "fix rattle" (See https://lammps.sandia.gov/doc/fix_shake.html)
-    fix fRattleTIP3P tip3p 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.
+    # Then use "fix shake" (See https://lammps.sandia.gov/doc/fix_shake.html)
+    fix fShakeTIP3P tip3p 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 fShakeTIP3P   # <-- this disables fix shake
+    # ...before minimization and pressure equilibration.  (compatibility issues)
   }