Merge pull request #163 from OpenBioSim/feature_minimise

Feature minimise

Author: chryswoods

corelib/src/libs/SireBase/progressbar.cpp

@@ -1158,6 +1158,11 @@ static void check_raise_interrupt()
     }
 }
 
+void ProgressBar::silentTick()
+{
+    check_raise_interrupt();
+}
+
 void ProgressBar::tick()
 {
     check_raise_interrupt();

corelib/src/libs/SireBase/progressbar.h

@@ -86,6 +86,7 @@ namespace SireBase
 
         void tick();
         void tick(const QString &text);
+        void silentTick();
 
         void setProgress(quint32 value);
         void setProgress(quint32 value, const QString &text);

corelib/src/libs/SireCAS/lambdaschedule.cpp

@@ -238,10 +238,12 @@ QString LambdaSchedule::toString() const
         auto keys = this->stage_equations[i].keys();
         std::sort(keys.begin(), keys.end());
 
-        for (auto lever : keys)
+        for (const auto &lever : keys)
         {
+            auto output_name = lever;
+            output_name.replace("*::", "");
             lines.append(QString("    %1: %2")
-                             .arg(lever.replace("*::", ""))
+                             .arg(output_name)
                              .arg(this->stage_equations[i][lever].toOpenMMString()));
         }
     }
@@ -956,20 +958,35 @@ SireCAS::Expression LambdaSchedule::_getEquation(int stage,
                                      CODELOC);
 
     const auto default_lever = _get_lever_name("*", lever);
+    const auto default_force = _get_lever_name(force, "*");
+    const auto lever_name = _get_lever_name(force, lever);
 
-    if (force == "*")
+    const auto equations = this->stage_equations[stage];
+
+    // search from most specific to least specific
+    auto it = equations.find(lever_name);
+
+    if (it != equations.end())
     {
-        return this->stage_equations[stage].value(
-            default_lever, this->default_equations[stage]);
+        return it.value();
     }
-    else
+
+    it = equations.find(default_force);
+
+    if (it != equations.end())
+    {
+        return it.value();
+    }
+
+    it = equations.find(default_lever);
+
+    if (it != equations.end())
     {
-        return this->stage_equations[stage].value(
-            _get_lever_name(force, lever),
-            this->stage_equations[stage].value(
-                default_lever,
-                this->default_equations[stage]));
+        return it.value();
     }
+
+    // we don't have any match, so return the default equation for this stage
+    return this->default_equations[stage];
 }
 
 /** Return the equation used to control the specified 'lever'
@@ -1142,15 +1159,33 @@ QHash<QString, QVector<double>> LambdaSchedule::getLeverValues(
     QVector<double> values(lambda_values.count(), NAN);
 
     QHash<QString, QVector<double>> lever_values;
-    lever_values.reserve(this->lever_names.count() + 1);
+
+    // get all of the lever / force combinations in use
+    QSet<QString> all_levers;
+
+    for (const auto &equations : this->stage_equations)
+    {
+        for (const auto &lever : equations.keys())
+        {
+            all_levers.insert(lever);
+        }
+    }
+
+    QStringList levers = all_levers.values();
+    std::sort(levers.begin(), levers.end());
+
+    lever_values.reserve(levers.count() + 2);
 
     lever_values.insert("λ", values);
 
     lever_values.insert("default", values);
 
-    for (const auto &lever_name : this->lever_names)
+    for (const auto &lever : levers)
     {
-        lever_values.insert(lever_name, values);
+        if (lever.startsWith("*::"))
+            lever_values.insert(lever.mid(3), values);
+        else
+            lever_values.insert(lever, values);
     }
 
     if (this->nStages() == 0)
@@ -1174,12 +1209,16 @@ QHash<QString, QVector<double>> LambdaSchedule::getLeverValues(
         const auto equation = this->default_equations[stage];
         lever_values["default"][i] = equation(input_values);
 
-        for (const auto &lever_name : lever_names)
+        for (const auto &lever : levers)
         {
-            const auto equation = this->stage_equations[stage].value(
-                lever_name, this->default_equations[stage]);
+            auto parts = lever.split("::");
+
+            const auto equation = this->_getEquation(stage, parts[0], parts[1]);
 
-            lever_values[lever_name][i] = equation(input_values);
+            if (lever.startsWith("*::"))
+                lever_values[lever.mid(3)][i] = equation(input_values);
+            else
+                lever_values[lever][i] = equation(input_values);
         }
     }
 

doc/source/changelog.rst

@@ -110,6 +110,20 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Fixed compile error using Python 3.12. This fixes issue #147.
 
+* Optimised the OpenMM minimisation code and making it more robust.
+  This includes vectorising for Apple Silicon and adding more tests for
+  convergence so that we can have more confidence that the structures
+  output are sensible. Also made sure that optimised compilation (-O3) is
+  used for all of the plugins (SireOpenMM, SireGemmi and SireRDKit).
+  They were previously compiled with wrapper options (e.g. -Os).
+  Minimisation now gives better progress updates, using a progress
+  bar to show progress towards the maximum number of iterations.
+  This has been reduced to 1500 by default. Also, if the minimisation
+  fails to create a structure that obeys constraints on the first pass,
+  then the minimisation is repeated, with the maximum number of
+  iterations reset. If it fails again, then this structure, with
+  constraints re-applied, is returned.
+
 * Please add an item to this changelog when you create your PR
 
 `2023.5.1 <https://github.com/openbiosim/sire/compare/2023.5.0...2023.5.1>`__ - January 2024

src/sire/mol/_dynamics.py

@@ -564,20 +564,78 @@ def step(self, num_steps: int = 1):
 
         self._omm_mols.getIntegrator().step(num_steps)
 
-    def run_minimisation(self, max_iterations: int):
+    def get_minimisation_log(self):
+        """
+        Return the log from the last minimisation
+        """
+        if self.is_null():
+            return None
+        else:
+            try:
+                return self._minimisation_log
+            except Exception:
+                return None
+
+    def run_minimisation(
+        self,
+        max_iterations: int = 10000,
+        tolerance: float = 10.0,
+        max_restarts: int = 10,
+        max_ratchets: int = 20,
+        ratchet_frequency: int = 500,
+        starting_k: float = 100.0,
+        ratchet_scale: float = 2.0,
+        max_constraint_error: float = 0.001,
+    ):
         """
         Internal method that runs minimisation on the molecules.
 
+        If the system is constrained, then a ratcheting algorithm is used.
+        The constraints are replaced by harmonic restraints with an
+        force constant based on `tolerance` and `starting_k`. Minimisation
+        is performed, with the actual constrained bond lengths checked
+        whenever minimisation converges, or when ratchet_frequency steps
+        have completed (whichever is sooner). The force constant of
+        the restraints is ratcheted up by `ratchet_scale`, and minimisation
+        continues until there is no large change in energy or the maximum
+        number of ratchets has been reached. In addition, at each ratchet,
+        the actual bond lengths of constrained bonds are compared against
+        the constrained values. If these have drifted too far away from
+        the constrained values, then the minimisation is restarted,
+        going back to the starting conformation and starting minimisation
+        at one higher ratchet level. This will repeat a maximum of
+        `max_restarts` times.
+
+        If a stable structure cannot be reached, then an exception
+        will be raised.
+
         Parameters:
 
         - max_iterations (int): The maximum number of iterations to run
+        - tolerance (float): The tolerance to use for the minimisation
+        - max_restarts (int): The maximum number of restarts before giving up
+        - max_ratchets (int): The maximum number of ratchets before giving up
+        - ratchet_frequency (int): The maximum number of steps between ratchets
+        - starting_k (float): The starting value of k for the minimisation
+        - ratchet_scale (float): The amount to scale k at each ratchet
+        - max_constraint_error (float): The maximum error in the constraint in nm
         """
         from ..legacy.Convert import minimise_openmm_context
 
         if max_iterations <= 0:
             max_iterations = 0
 
-        minimise_openmm_context(self._omm_mols, max_iterations=max_iterations)
+        self._minimisation_log = minimise_openmm_context(
+            self._omm_mols,
+            tolerance=tolerance,
+            max_iterations=max_iterations,
+            max_restarts=max_restarts,
+            max_ratchets=max_ratchets,
+            ratchet_frequency=ratchet_frequency,
+            starting_k=starting_k,
+            ratchet_scale=ratchet_scale,
+            max_constraint_error=max_constraint_error,
+        )
 
     def _rebuild_and_minimise(self):
         if self.is_null():
@@ -600,7 +658,7 @@ def _rebuild_and_minimise(self):
 
         self._omm_mols = to(self._sire_mols, "openmm", map=self._map)
 
-        self.run_minimisation(max_iterations=10000)
+        self.run_minimisation()
 
     def run(
         self,
@@ -1032,18 +1090,57 @@ def __repr__(self):
     def minimise(
         self,
         max_iterations: int = 10000,
+        tolerance: float = 10.0,
+        max_restarts: int = 10,
+        max_ratchets: int = 20,
+        ratchet_frequency: int = 500,
+        starting_k: float = 100.0,
+        ratchet_scale: float = 2.0,
+        max_constraint_error: float = 0.001,
     ):
         """
-        Perform minimisation on the molecules, running a maximum
-        of max_iterations iterations.
+        Internal method that runs minimisation on the molecules.
+
+        If the system is constrained, then a ratcheting algorithm is used.
+        The constraints are replaced by harmonic restraints with an
+        force constant based on `tolerance` and `starting_k`. Minimisation
+        is performed, with the actual constrained bond lengths checked
+        whenever minimisation converges, or when ratchet_frequency steps
+        have completed (whichever is sooner). The force constant of
+        the restraints is ratcheted up by `ratchet_scale`, and minimisation
+        continues until there is no large change in energy or the maximum
+        number of ratchets has been reached. In addition, at each ratchet,
+        the actual bond lengths of constrained bonds are compared against
+        the constrained values. If these have drifted too far away from
+        the constrained values, then the minimisation is restarted,
+        going back to the starting conformation and starting minimisation
+        at one higher ratchet level. This will repeat a maximum of
+        `max_restarts` times.
+
+        If a stable structure cannot be reached, then an exception
+        will be raised.
 
         Parameters:
 
         - max_iterations (int): The maximum number of iterations to run
+        - tolerance (float): The tolerance to use for the minimisation
+        - max_restarts (int): The maximum number of restarts before giving up
+        - max_ratchets (int): The maximum number of ratchets before giving up
+        - ratchet_frequency (int): The maximum number of steps between ratchets
+        - starting_k (float): The starting value of k for the minimisation
+        - ratchet_scale (float): The amount to scale k at each ratchet
+        - max_constraint_error (float): The maximum error in the constraints in nm
         """
         if not self._d.is_null():
             self._d.run_minimisation(
                 max_iterations=max_iterations,
+                tolerance=tolerance,
+                max_restarts=max_restarts,
+                max_ratchets=max_ratchets,
+                ratchet_frequency=ratchet_frequency,
+                starting_k=starting_k,
+                ratchet_scale=ratchet_scale,
+                max_constraint_error=max_constraint_error,
             )
 
         return self