Merge changes for 8.3.1 (#5004)

* Fixed bug in computing pressure (#4980)

* Add function parsing to custom multiparticle force (#4986)

* Merge pull request #4989 from epretti/fix-charmm-water-models

Update CHARMM36 2024 and add water models

---------
Co-authored-by: Ezra Greenberg <120955867+egreenberg7@users.noreply.github.com>
Co-authored-by: Evan Pretti <pretti@stanford.edu>

docs-source/usersguide/application/02_running_sims.rst

@@ -590,31 +590,31 @@ for full references).
 
 .. tabularcolumns:: |l|L|
 
-=================================  ============================================
-File                               Parameters
-=================================  ============================================
-:file:`charmm36_2024.xml`          Protein, DNA, RNA, lipids, carbohydrates, and small molecules
-:file:`charmm36/water.xml`         Default CHARMM water model (a modified version of TIP3P\ :cite:`Jorgensen1983`) and ions
-:file:`charmm36/spce.xml`          SPC/E water model\ :cite:`Berendsen1987` and ions
-:file:`charmm36/tip3p-pme-b.xml`   TIP3P-PME-B water model\ :cite:`Price2004` and ions
-:file:`charmm36/tip3p-pme-f.xml`   TIP3P-PME-F water model\ :cite:`Price2004` and ions
-:file:`charmm36/tip4pew.xml`       TIP4P-Ew water model\ :cite:`Horn2004` and ions
-:file:`charmm36/tip4p2005.xml`     TIP4P-2005 water model\ :cite:`Abascal2005` and ions
-:file:`charmm36/tip5p.xml`         TIP5P water model\ :cite:`Mahoney2000` and ions
-:file:`charmm36/tip5pew.xml`       TIP5P-Ew water model\ :cite:`Rick2004` and ions
-=================================  ============================================
+=====================================  =========================================
+File                                   Parameters
+=====================================  =========================================
+:file:`charmm36_2024.xml`              Protein, DNA, RNA, lipids, carbohydrates, and small molecules
+:file:`charmm36_2024/water.xml`        Default CHARMM water model (a modified version of TIP3P\ :cite:`Jorgensen1983`) and ions
+:file:`charmm36_2024/spce.xml`         SPC/E water model\ :cite:`Berendsen1987` and ions
+:file:`charmm36_2024/tip3p-pme-b.xml`  TIP3P-PME-B water model\ :cite:`Price2004` and ions
+:file:`charmm36_2024/tip3p-pme-f.xml`  TIP3P-PME-F water model\ :cite:`Price2004` and ions
+:file:`charmm36_2024/tip4pew.xml`      TIP4P-Ew water model\ :cite:`Horn2004` and ions
+:file:`charmm36_2024/tip4p2005.xml`    TIP4P-2005 water model\ :cite:`Abascal2005` and ions
+:file:`charmm36_2024/tip5p.xml`        TIP5P water model\ :cite:`Mahoney2000` and ions
+:file:`charmm36_2024/tip5pew.xml`      TIP5P-Ew water model\ :cite:`Rick2004` and ions
+=====================================  =========================================
 
 The file :file:`charmm36_2024.xml` bundles everything but the water and ions into a single
 file.  In most cases, you can simply include that file, plus one of the water models,
-such as :file:`charmm36/water.xml`, which specifies the default CHARMM water model
+such as :file:`charmm36_2024/water.xml`, which specifies the default CHARMM water model
 (a modified version of TIP3P\ :cite:`Jorgensen1983`) and ions:
 ::
 
-    forcefield = ForceField('charmm36_2024.xml', 'charmm36/water.xml')
+    forcefield = ForceField('charmm36_2024.xml', 'charmm36_2024/water.xml')
 
-.. tip:: The solvent model XML files included under the :file:`charmm36/` directory
+.. tip:: The solvent model XML files included under the :file:`charmm36_2024/` directory
          include both water *and* ions compatible with that water model, so if you
-         mistakenly specify :file:`tip3p.xml` instead of :file:`charmm36/water.xml`,
+         mistakenly specify :file:`tip3p.xml` instead of :file:`charmm36_2024/water.xml`,
          you run the risk of having :class:`ForceField` raise an exception due to
          missing parameters for ions in your system.
 
@@ -641,6 +641,13 @@ such as :file:`charmm36/water.xml`, which specifies the default CHARMM water mod
          `CHARMM-GUI <http://charmm-gui.org/>`_, it's easiest to load
          the PSF file directly, as discussed in Section :numref:`using-charmm-files`.
 
+.. warning:: Solvent model files in the :file:`charmm36_2024/` directory are
+             designed for use with :file:`charmm36_2024.xml`, and those in the
+             :file:`charmm36/` directory are designed for use with an older
+             release of the CHARMM force field found in :file:`charmm36.xml` and
+             described below.  Be sure to load the correct version of the
+             solvent model you wish to use along with the main force field file.
+
 .. warning:: Some residues and patches are not included in this port of
              CHARMM36, either due to a lack of support for certain CHARMM
              features in `ParmEd <https://github.com/parmed/parmed>`_, or
@@ -764,6 +771,7 @@ File                           Force Field
 :code:`amber03.xml`            Amber03\ :cite:`Duan2003`
 :code:`amber10.xml`            Amber10 (documented in the AmberTools_ manual as `ff10`)
 :code:`charmm36.xml`           July 2017 release of the CHARMM36 force field without CHARMM36m protein parameters
+:code:`charmm36/*.xml`         Parameters for water models and ions specifically for use with :code:`charmm36.xml`
 :code:`charmm_polar_2013.xml`  2013 version of the CHARMM polarizable force field\ :cite:`Lopes2013`
 :code:`charmm_polar_2019.xml`  2019 version of the CHARMM polarizable force field\ :cite:`Lopes2013`
 =============================  ==================================================================================

platforms/common/src/kernels/monteCarloBarostat.cc

@@ -79,7 +79,7 @@ KERNEL void computeMolecularKineticEnergy(int numMolecules, GLOBAL mixed4* RESTR
             molVel += mass*trimTo3(v);
             molMass += mass;
         }
-        molVel *= RECIP((mixed) numAtoms);
+        molVel *= RECIP((mixed) molMass);
 #if COMPONENTS == 1
         ke[0] += 0.5f*molMass*dot(molVel, molVel);
 #else

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp

@@ -133,7 +133,7 @@ void ReferenceMonteCarloBarostat::computeMolecularKineticEnergy(const vector<Vec
             molVel += masses[atom]*velocities[atom];
             molMass += masses[atom];
         }
-        molVel /= molecule.size();
+        molVel /= molMass;
         if (components == 1)
             ke[0] += 0.5*molMass*molVel.dot(molVel);
         else {

tests/TestMonteCarloBarostat.h

@@ -167,8 +167,8 @@ void testMolecularGas() {
     init_gen_rand(0, sfmt);
     for (int i = 0; i < numMolecules; ++i) {
         system.addParticle(1.0);
-        system.addParticle(1.0);
-        system.addParticle(1.0);
+        system.addParticle(2.0);
+        system.addParticle(3.0);
         Vec3 pos(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
         bonds->addBond(positions.size(), positions.size()+1, 0.1, 10.0);
         system.addConstraint(positions.size(), positions.size()+2, 0.1);

wrappers/python/openmm/app/data/charmm36.xml

@@ -103763,4 +103763,12 @@
     <NBFixPair epsilon="0.8631592000000001" sigma="0.25134925534893393" type1="MG" type2="SO"/>
     <NBFixPair epsilon="0.8414024" sigma="0.25219560913116723" type1="MG" type2="SS1"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>

wrappers/python/openmm/app/data/charmm36/spce.xml

@@ -121,4 +121,12 @@
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OC2D2" type2="SOD"/>
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OG2D2" type2="SOD"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>

wrappers/python/openmm/app/data/charmm36/tip3p-pme-b.xml

@@ -121,4 +121,12 @@
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OC2D2" type2="SOD"/>
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OG2D2" type2="SOD"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>

wrappers/python/openmm/app/data/charmm36/tip3p-pme-f.xml

@@ -121,4 +121,12 @@
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OC2D2" type2="SOD"/>
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OG2D2" type2="SOD"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>

wrappers/python/openmm/app/data/charmm36/tip4p2005.xml

@@ -127,4 +127,12 @@
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OC2D2" type2="SOD"/>
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OG2D2" type2="SOD"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>

wrappers/python/openmm/app/data/charmm36/tip4pew.xml

@@ -127,4 +127,12 @@
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OC2D2" type2="SOD"/>
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OG2D2" type2="SOD"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>

wrappers/python/openmm/app/data/charmm36/tip5p.xml

@@ -129,4 +129,12 @@
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OC2D2" type2="SOD"/>
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OG2D2" type2="SOD"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>

wrappers/python/openmm/app/data/charmm36/tip5pew.xml

@@ -129,4 +129,12 @@
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OC2D2" type2="SOD"/>
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28419669108676826" type1="OG2D2" type2="SOD"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>

wrappers/python/openmm/app/data/charmm36/water.xml

@@ -127,4 +127,12 @@
     <NBFixPair epsilon="0.31388368000000005" sigma="0.2877602859593296" type1="OG2D2" type2="SOD"/>
     <NBFixPair epsilon="0.31388368000000005" sigma="0.28152399493234725" type1="OG2P1" type2="SOD"/>
   </LennardJonesForce>
+  <InitializationScript>
+EXPECTED_VERSION = "charmm36"
+if hasattr(self, "_charmm_ff_version"):
+    if self._charmm_ff_version != EXPECTED_VERSION:
+        raise RuntimeError(f"CHARMM force field versions ({self._charmm_ff_version} and {EXPECTED_VERSION}) incompatible")
+else:
+    self._charmm_ff_version = EXPECTED_VERSION
+  </InitializationScript>
 </ForceField>