Fix issues #366 and #367

doc/source/changelog.rst

@@ -17,6 +17,10 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR.
 
+* Reset GCMC water state prior to minimisation following a dynamics crash.
+
+* Set appropriate default OpenMM box vectors when a system has no periodic space.
+
 `2025.2.0 <https://github.com/openbiosim/sire/compare/2025.1.0...2025.2.0>`__ - October 2025
 --------------------------------------------------------------------------------------------
 

src/sire/mol/_dynamics.py

@@ -966,6 +966,12 @@ def _rebuild_and_minimise(self):
 
         self._omm_mols = to(self._sire_mols, "openmm", map=self._map)
 
+        # reset the water state
+        if self._gcmc_sampler is not None:
+            self._gcmc_sampler.push()
+            self._gcmc_sampler._set_water_state(self._omm_mols)
+            self._gcmc_sampler.pop()
+
         self.run_minimisation()
 
     def run(
@@ -1419,12 +1425,6 @@ class NeedsMinimiseError(Exception):
                 self._prev_elapsed_time.to("picosecond") * picosecond
             )
 
-            # reset the water state
-            if self._gcmc_sampler is not None:
-                self._gcmc_sampler.push()
-                self._gcmc_sampler._set_water_state(self._omm_mols)
-                self._gcmc_sampler.pop()
-
             self.run(**orig_args)
             return
 

tests/convert/test_openmm.py

@@ -411,3 +411,39 @@ def test_openmm_skipped_constrained_bonds(zero_lj_mols, openmm_platform):
             line1 = lines1[i]
             assert "Bond" in line1
             i += 1
+
+
+@pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+def test_openmm_default_box_vectors(ala_mols, openmm_platform):
+    mols = ala_mols.clone()
+
+    # Remove the shared space property.
+    mols.remove_shared_property("space")
+
+    # Create a "vacuum" system.
+    omm = sr.convert.to(
+        mols[0],
+        "openmm",
+        map={
+            "cutoff": "12 A",
+            "platform": openmm_platform,
+        },
+    )
+
+    # Get the box vectors.
+    box_vectors = omm.getState().getPeriodicBoxVectors()
+
+    # Get the AABox of the first molecule.
+    aabox = sr.legacy.Vol.AABox(mols[0].property("coordinates").to_vector())
+
+    # Work out the box vectors from the AABox, i.e. the box size plus twice the cutoff.
+    box = 2.0 * aabox.half_extents() + sr.maths.Vector(24, 24, 24)
+
+    from openmm.unit import angstrom
+
+    # Check that the box vectors match the expected values.
+    for i, vec in enumerate(box_vectors):
+        assert vec[i].value_in_unit(angstrom) == pytest.approx(box[i].value(), abs=1e-3)

wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp

@@ -230,8 +230,6 @@ void _add_bond_restraints(const SireMM::BondRestraints &restraints,
     const double internal_to_nm = (1 * SireUnits::angstrom).to(SireUnits::nanometer);
     const double internal_to_k = (1 * SireUnits::kcal_per_mol / (SireUnits::angstrom2)).to(SireUnits::kJ_per_mol / (SireUnits::nanometer2));
 
-    auto cljff = lambda_lever.getForce<OpenMM::NonbondedForce>("clj", system);
-
     std::vector<double> custom_params = {1.0, 0.0, 0.0};
 
     for (const auto &restraint : atom_restraints)
@@ -308,8 +306,6 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res
     const double internal_to_k = (1 * SireUnits::kcal_per_mol / (SireUnits::angstrom2)).to(SireUnits::kJ_per_mol / (SireUnits::nanometer2));
     const double internal_to_de = (1 * SireUnits::kcal_per_mol).to(SireUnits::kJ_per_mol);
 
-    auto cljff = lambda_lever.getForce<OpenMM::NonbondedForce>("clj", system);
-
     std::vector<double> custom_params = {1.0, 0.0, 0.0, 0.0};
 
     for (const auto &restraint : atom_restraints)
@@ -797,7 +793,9 @@ void _set_clj_cutoff(OpenMM::NonbondedForce &cljff,
  */
 std::shared_ptr<std::vector<OpenMM::Vec3>>
 _set_box_vectors(OpenMM::System &system,
-                 const ForceFieldInfo &ffinfo)
+                 const SelectorMol &mols,
+                 const ForceFieldInfo &ffinfo,
+                 const PropertyMap &map)
 {
     // create the periodic box vectors
     std::shared_ptr<std::vector<OpenMM::Vec3>> boxvecs;
@@ -848,6 +846,45 @@ _set_box_vectors(OpenMM::System &system,
                                             boxvecs_data[1],
                                             boxvecs_data[2]);
     }
+    else
+    {
+        // Set the box vectors based on the AABox of the system and nonbonded
+        // cutoff distance.
+
+        QVector<QVector<SireMaths::Vector>> all_coords;
+
+        // Get the coordinates from all of the molecules.
+        for (const auto &mol : mols)
+        {
+            const auto coords = mol.property(map["coordinates"]).asA<SireMol::AtomCoords>().toVector();
+            all_coords.append(coords);
+        }
+
+        // Create an AABox.
+        const auto aabox = SireVol::AABox(all_coords);
+
+        // Work out the minimum box size and convert from Angstroms to nm.
+        auto min_box = 2.0 * aabox.halfExtents() * 0.1;
+
+        // Adjust the box size based on the cutoff distance.
+        if (ffinfo.hasCutoff())
+        {
+            // Get the nonbonded cutoff.
+            const auto double_cutoff = 2.0 * ffinfo.cutoff().to(SireUnits::nanometers);
+            min_box += SireMaths::Vector(double_cutoff, double_cutoff, double_cutoff);
+        }
+
+        boxvecs.reset(new std::vector<OpenMM::Vec3>(3));
+        auto boxvecs_data = boxvecs->data();
+
+        boxvecs_data[0] = OpenMM::Vec3(min_box.x(), 0, 0);
+        boxvecs_data[1] = OpenMM::Vec3(0, min_box.y(), 0);
+        boxvecs_data[2] = OpenMM::Vec3(0, 0, min_box.z());
+
+        system.setDefaultPeriodicBoxVectors(boxvecs_data[0],
+                                            boxvecs_data[1],
+                                            boxvecs_data[2]);
+    }
 
     return boxvecs;
 }
@@ -948,7 +985,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system,
     }
 
     // set the box vectors for periodic spaces
-    auto boxvecs = _set_box_vectors(system, ffinfo);
+    auto boxvecs = _set_box_vectors(system, mols, ffinfo, map);
 
     // Are any of the molecules perturbable, and if they are,
     // should we just ignore perturbations?