diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index ed4b73b82..046c0f454 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -17,6 +17,10 @@ organisation on `GitHub `__. * 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 `__ - October 2025 -------------------------------------------------------------------------------------------- diff --git a/src/sire/mol/_dynamics.py b/src/sire/mol/_dynamics.py index 7d0ae6592..b18e186c3 100644 --- a/src/sire/mol/_dynamics.py +++ b/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 diff --git a/tests/convert/test_openmm.py b/tests/convert/test_openmm.py index 90d4dff1b..4a6b58d69 100644 --- a/tests/convert/test_openmm.py +++ b/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) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 11d881663..85d79657f 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/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("clj", system); - std::vector 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("clj", system); - std::vector 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> _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> 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> all_coords; + + // Get the coordinates from all of the molecules. + for (const auto &mol : mols) + { + const auto coords = mol.property(map["coordinates"]).asA().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(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?