Fixed regression caused by the removal of constraints from water mole…

…cules

with no internal bonding. Added a better algorithm to detect and fix
this.

Also cleaned up the code to get a shared constraint length.

Added the changelog entries ready for the PR

doc/source/changelog.rst

@@ -15,6 +15,23 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 `2023.4.0 <https://github.com/openbiosim/sire/compare/2023.4.0...2023.5.0>`__ - December 2023
 ---------------------------------------------------------------------------------------------
 
+* Fixed regression introduced in 2023.4.0 that meant that removed the constraints
+  from water molecules that had no internal bonds. These waters would blow up
+  as there was nothing holding them together. The need for these constraints is
+  now better detected and explicitly added.
+
+* Significantly sped up the OpenMM layer by checking for similar constraint lengths
+  and matching them all to be the same (within 0.05 A for calculated constraints,
+  e.g. unbonded atoms or angle constraints) or to R0 for bonds where the bond
+  length is within 0.1 A of R0 and the molecule isn't perturbable.
+
+* Added a custom minimiser that is based on OpenMM's LocalEnergyMinimizer,
+  but that copes better with exclusion errors, and that has deep integration
+  with the progress bar / interuption system.
+
+* Fixed a bug where the exclusions and exceptions were mismatched for the
+  OpenMM CPU platform, leading to exclusion errors.
+
 * Please add an item to this changelog when you create your PR
 
 `2023.4.0 <https://github.com/openbiosim/sire/compare/2023.3.0...2023.4.0>`__ - October 2023

tests/convert/test_openmm_minimise.py

@@ -6,8 +6,8 @@
     "openmm" not in sr.convert.supported_formats(),
     reason="openmm support is not available",
 )
-def test_openmm_simple_minimise(kigaki_mols):
-    mols = kigaki_mols
+def test_openmm_simple_minimise(ala_mols):
+    mols = ala_mols
 
     nrg = mols.energy()
 
@@ -40,3 +40,44 @@ def test_openmm_minimise_lambda(merged_ethane_methanol):
     mols = (
         mols[0].minimisation(platform="cpu", lambda_value=0.5).run(5).commit()
     )
+
+
+@pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+def test_openmm_minimise_unbonded_water(kigaki_mols):
+    mols = kigaki_mols
+
+    atoms = mols[100].atoms()
+
+    # the water molecules have no internal bonds, so this tests
+    # whether or not constraints have been added correctly
+    mols = mols.minimisation(platform="cpu").run(10).commit()
+
+    new_atoms = mols[100].atoms()
+
+    # make sure the water has not blown up... (low tolerance as
+    # the water bond lengths will be constrained to a common shared
+    # value)
+    assert (
+        sr.measure(atoms[0], atoms[1]).value()
+        == pytest.approx(
+            sr.measure(new_atoms[0], new_atoms[1]).value(), abs=1e-2
+        )
+    ) or (
+        sr.measure(new_atoms[0], new_atoms[2]).value() == pytest.approx(0.9572)
+    )
+
+    assert (
+        sr.measure(atoms[0], atoms[2]).value()
+        == pytest.approx(
+            sr.measure(new_atoms[0], new_atoms[2]).value(), abs=1e-2
+        )
+    ) or (
+        sr.measure(new_atoms[0], new_atoms[2]).value() == pytest.approx(0.9572)
+    )
+
+    assert sr.measure(atoms[1], atoms[2]).value() == pytest.approx(
+        sr.measure(new_atoms[1], new_atoms[2]).value(), abs=1e-2
+    )

wrapper/Convert/SireOpenMM/openmmmolecule.cpp

@@ -312,15 +312,21 @@ std::tuple<int, int, double, double, double> OpenMMMolecule::getException(
                            charge, sigma, epsilon);
 }
 
-// distance below which constraints are considered to be equal (e.g. to
-// r0 or to another angle - units in nanometers)
-const double constraint_length_tolerance = 0.01;
-
 /** Return closest constraint length to 'length' based on what
  *  we have seen before and the constraint_length_tolerance
  */
 double getSharedConstraintLength(double length)
 {
+    // distance below which constraints are considered to be equal (e.g. to
+    // r0 or to another angle - units in nanometers)
+    const double constraint_length_tolerance = 0.005;
+
+    // is this close to a O-H bond length of water?
+    if (std::abs(length - 0.09572) < constraint_length_tolerance)
+    {
+        return 0.09572;
+    }
+
     static QVector<double> angle_constraint_lengths;
     static QReadWriteLock l;
 
@@ -330,7 +336,7 @@ double getSharedConstraintLength(double length)
     // is this close to any of the existing angle constraints?
     for (const auto &cl : angle_constraint_lengths)
     {
-        if (std::abs(cl - length) < 0.01)
+        if (std::abs(cl - length) < constraint_length_tolerance)
         {
             // this is close enough to an existing constraint
             // so we will use that instead
@@ -346,7 +352,7 @@ double getSharedConstraintLength(double length)
 
     for (const auto &cl : angle_constraint_lengths)
     {
-        if (std::abs(cl - length) < 0.01)
+        if (std::abs(cl - length) < constraint_length_tolerance)
         {
             // this is close enough to an existing constraint
             // so we will use that instead
@@ -494,6 +500,14 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol,
     this->bond_params.clear();
     this->constraints.clear();
 
+    // initialise all atoms as being unbonded
+    this->unbonded_atoms.reserve(nats);
+
+    for (int i = 0; i < nats; ++i)
+    {
+        this->unbonded_atoms.insert(i);
+    }
+
     // now the bonds
     const double bond_k_to_openmm = 2.0 * (SireUnits::kcal_per_mol / (SireUnits::angstrom * SireUnits::angstrom)).to(SireUnits::kJ_per_mol / (SireUnits::nanometer * SireUnits::nanometer));
     const double bond_r0_to_openmm = SireUnits::angstrom.to(SireUnits::nanometer);
@@ -523,6 +537,12 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol,
         const double k = bondparam.k() * bond_k_to_openmm;
         const double r0 = bondparam.r0() * bond_r0_to_openmm;
 
+        if (k != 0)
+        {
+            this->unbonded_atoms.remove(atom0);
+            this->unbonded_atoms.remove(atom1);
+        }
+
         const bool has_light_atom = (masses_data[atom0] < 2.5 or masses_data[atom1] < 2.5);
         const bool has_massless_atom = masses_data[atom0] < 0.5 or masses_data[atom1] < 0.5;
 
@@ -545,7 +565,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol,
 
             // use the r0 for the bond if this is close to the measured length and this
             // is not a perturbable molecule
-            if (not is_perturbable and std::abs(constraint_length - r0) < constraint_length_tolerance)
+            if (not is_perturbable and std::abs(constraint_length - r0) < 0.01)
             {
                 constraint_length = r0;
             }
@@ -832,6 +852,9 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map)
             bond_params_1.append(bond1);
 
             found_index_1[j] = true;
+
+            unbonded_atoms.remove(atom0);
+            unbonded_atoms.remove(atom1);
         }
     }
 
@@ -1081,10 +1104,13 @@ void OpenMMMolecule::buildExceptions(const Molecule &mol,
     // we need to allow this set to morph
     exception_params.clear();
 
+    // save memory
+    if (unbonded_atoms.isEmpty())
+        unbonded_atoms = QSet<qint32>();
+
     const int nats = this->cljs.count();
 
     const auto &nbpairs = mol.property(map["intrascale"]).asA<CLJNBPairs>();
-    const auto &connectivity = mol.property(map["connectivity"]).asA<Connectivity>();
 
     const int ncgs = mol.nCutGroups();
 
@@ -1103,16 +1129,15 @@ void OpenMMMolecule::buildExceptions(const Molecule &mol,
                 // are any of these atoms unbonded? If so, then
                 // we will need to add a constraint to hold them
                 // in place
-                if (connectivity.nConnections(AtomIdx(i)) == 0 or
-                    connectivity.nConnections(AtomIdx(j)) == 0)
+                if (unbonded_atoms.contains(i) or unbonded_atoms.contains(j))
                 {
                     if (not constrained_pairs.contains(to_pair(i, j)))
                     {
                         const auto delta = coords[j] - coords[i];
                         const auto length = std::sqrt((delta[0] * delta[0]) +
                                                       (delta[1] * delta[1]) +
                                                       (delta[2] * delta[2]));
-                        constraints.append(std::make_tuple(i, j, length));
+                        constraints.append(std::make_tuple(i, j, getSharedConstraintLength(length)));
                         constrained_pairs.insert(to_pair(i, j));
                     }
                 }
@@ -1131,13 +1156,51 @@ void OpenMMMolecule::buildExceptions(const Molecule &mol,
         {
             // two atoms must be excluded
             exception_params.append(std::make_tuple(0, 1, 0.0, 0.0));
+
+            if (unbonded_atoms.contains(0) or unbonded_atoms.contains(1))
+            {
+                // these atoms are not bonded - we need to add a constraint
+                // to keep them in place
+                const auto delta = coords[1] - coords[0];
+                const auto length = std::sqrt((delta[0] * delta[0]) +
+                                              (delta[1] * delta[1]) +
+                                              (delta[2] * delta[2]));
+                constraints.append(std::make_tuple(0, 1, getSharedConstraintLength(length)));
+                constrained_pairs.insert(to_pair(0, 1));
+            }
         }
-        else if (nats <= 3)
+        else if (nats == 3)
         {
             // three atoms must be excluded
             exception_params.append(std::make_tuple(0, 1, 0.0, 0.0));
             exception_params.append(std::make_tuple(1, 2, 0.0, 0.0));
             exception_params.append(std::make_tuple(0, 2, 0.0, 0.0));
+
+            if (unbonded_atoms.contains(0) or unbonded_atoms.contains(1) or unbonded_atoms.contains(2))
+            {
+                // these atoms are not bonded - we need to add a constraint
+                // to keep them in place
+                auto delta = coords[1] - coords[0];
+                auto length = std::sqrt((delta[0] * delta[0]) +
+                                        (delta[1] * delta[1]) +
+                                        (delta[2] * delta[2]));
+                constraints.append(std::make_tuple(0, 1, getSharedConstraintLength(length)));
+                constrained_pairs.insert(to_pair(0, 1));
+
+                delta = coords[2] - coords[0];
+                length = std::sqrt((delta[0] * delta[0]) +
+                                   (delta[1] * delta[1]) +
+                                   (delta[2] * delta[2]));
+                constraints.append(std::make_tuple(0, 2, getSharedConstraintLength(length)));
+                constrained_pairs.insert(to_pair(0, 2));
+
+                delta = coords[2] - coords[1];
+                length = std::sqrt((delta[0] * delta[0]) +
+                                   (delta[1] * delta[1]) +
+                                   (delta[2] * delta[2]));
+                constraints.append(std::make_tuple(1, 2, getSharedConstraintLength(length)));
+                constrained_pairs.insert(to_pair(1, 2));
+            }
         }
         else
         {

wrapper/Convert/SireOpenMM/openmmmolecule.h

@@ -144,6 +144,9 @@ namespace SireOpenMM
          */
         SireBase::PropertyMap perturtable_map;
 
+        /** The atoms that are missing any internal parameters */
+        QSet<qint32> unbonded_atoms;
+
         /** Alpha values for all of the atoms. This is equal to zero for
          *  non-ghost atoms, and one for ghost atoms
          */