Add restricted angle bending for Boresch restraint and automated restraint search

corelib/src/libs/SireMM/boreschrestraints.cpp

@@ -74,9 +74,7 @@ QDataStream &operator>>(QDataStream &ds, BoreschRestraint &borrest)
     {
         SharedDataStream sds(ds);
 
-        sds >> borrest.receptor_atms >> borrest.ligand_atms >> borrest._r0
-            >> borrest._theta0 >> borrest._phi0 >> borrest._kr >> borrest._ktheta
-            >> borrest._kphi >> static_cast<Property &>(borrest);
+        sds >> borrest.receptor_atms >> borrest.ligand_atms >> borrest._r0 >> borrest._theta0 >> borrest._phi0 >> borrest._kr >> borrest._ktheta >> borrest._kphi >> static_cast<Property &>(borrest);
     }
     else
         throw version_error(v, "1", r_borrest, CODELOC);
@@ -285,7 +283,7 @@ QString BoreschRestraint::toString() const
             .arg(k.join(", "))
             .arg(_r0.toString())
             .arg(t.join(", "))
-            .arg(p.join(', '));
+            .arg(p.join(", "));
     }
 }
 
@@ -345,11 +343,11 @@ static const RegisterMetaType<BoreschRestraints> r_borrests;
 
 QDataStream &operator<<(QDataStream &ds, const BoreschRestraints &borrests)
 {
-    writeHeader(ds, r_borrests, 2);
+    writeHeader(ds, r_borrests, 3);
 
     SharedDataStream sds(ds);
 
-    sds << borrests.r << borrests.use_pbc
+    sds << borrests.r << borrests.use_pbc << borrests.angle_potential << borrests.restraint_lever
         << static_cast<const Restraints &>(borrests);
 
     return ds;
@@ -364,16 +362,29 @@ QDataStream &operator>>(QDataStream &ds, BoreschRestraints &borrests)
         SharedDataStream sds(ds);
 
         sds >> borrests.r >> static_cast<Restraints &>(borrests);
+
+        borrests.use_pbc = false;
+        borrests.angle_potential = "harmonic";
+        borrests.restraint_lever = "combined";
     }
     else if (v == 2)
     {
         SharedDataStream sds(ds);
 
-        sds >> borrests.r >> borrests.use_pbc
-            >> static_cast<Restraints &>(borrests);
+        sds >> borrests.r >> borrests.use_pbc >> static_cast<Restraints &>(borrests);
+
+        borrests.angle_potential = "harmonic";
+        borrests.restraint_lever = "combined";
+    }
+    else if (v == 3)
+    {
+        SharedDataStream sds(ds);
+
+        sds >> borrests.r >> borrests.use_pbc >> borrests.angle_potential >> borrests.restraint_lever >>
+            static_cast<Restraints &>(borrests);
     }
     else
-        throw version_error(v, "1,2", r_borrests, CODELOC);
+        throw version_error(v, "1,2,3", r_borrests, CODELOC);
 
     return ds;
 }
@@ -426,7 +437,8 @@ BoreschRestraints::BoreschRestraints(const QString &name,
 }
 
 BoreschRestraints::BoreschRestraints(const BoreschRestraints &other)
-    : ConcreteProperty<BoreschRestraints, Restraints>(other), r(other.r), use_pbc(other.use_pbc)
+    : ConcreteProperty<BoreschRestraints, Restraints>(other), r(other.r), use_pbc(other.use_pbc),
+      angle_potential(other.angle_potential), restraint_lever(other.restraint_lever)
 {
 }
 
@@ -438,14 +450,17 @@ BoreschRestraints &BoreschRestraints::operator=(const BoreschRestraints &other)
 {
     r = other.r;
     use_pbc = other.use_pbc;
+    angle_potential = other.angle_potential;
+    restraint_lever = other.restraint_lever;
     Restraints::operator=(other);
     return *this;
 }
 
 bool BoreschRestraints::operator==(const BoreschRestraints &other) const
 {
     return r == other.r and Restraints::operator==(other) and
-           use_pbc == other.use_pbc;
+           use_pbc == other.use_pbc and angle_potential == other.angle_potential and
+           restraint_lever == other.restraint_lever;
 }
 
 bool BoreschRestraints::operator!=(const BoreschRestraints &other) const
@@ -499,11 +514,14 @@ QString BoreschRestraints::toString() const
         }
     }
 
-    return QObject::tr("BoreschRestraints( name=%1, size=%2, use_pbc=%3\n%4\n)")
-            .arg(this->name())
-            .arg(n)
-            .arg(this->use_pbc ? "true" : "false")
-            .arg(parts.join("\n"));
+    return QObject::tr("BoreschRestraints( name=%1, size=%2, use_pbc=%3, angle_potential=%4, "
+                       "restraint_lever=%5\n%6\n)")
+        .arg(this->name())
+        .arg(n)
+        .arg(this->use_pbc ? "true" : "false")
+        .arg(this->angle_potential)
+        .arg(this->restraint_lever)
+        .arg(parts.join("\n"));
 }
 
 /** Return whether or not this is empty */
@@ -610,3 +628,50 @@ bool BoreschRestraints::usesPbc() const
 {
     return this->use_pbc;
 }
+
+/** Set the functional form used for the two Boresch angle restraint terms.
+ *  Must be either "harmonic" (the default) or "restricted_bending". */
+void BoreschRestraints::setAnglePotential(const QString &angle_potential)
+{
+    if (angle_potential != "harmonic" and angle_potential != "restricted_bending")
+    {
+        throw SireError::invalid_arg(QObject::tr(
+                                         "'angle_potential' must be either 'harmonic' or "
+                                         "'restricted_bending', got '%1'.")
+                                         .arg(angle_potential),
+                                     CODELOC);
+    }
+
+    this->angle_potential = angle_potential;
+}
+
+/** Return the functional form used for the two Boresch angle restraint terms,
+ *  either "harmonic" or "restricted_bending". */
+QString BoreschRestraints::anglePotential() const
+{
+    return this->angle_potential;
+}
+
+/** Set how the restraint's six degrees of freedom are grouped into
+ *  lambda-addressable OpenMM Forces. Must be either "combined" (the
+ *  default) or "split". */
+void BoreschRestraints::setRestraintLever(const QString &restraint_lever)
+{
+    if (restraint_lever != "combined" and restraint_lever != "split")
+    {
+        throw SireError::invalid_arg(QObject::tr(
+                                         "'restraint_lever' must be either 'combined' or "
+                                         "'split', got '%1'.")
+                                         .arg(restraint_lever),
+                                     CODELOC);
+    }
+
+    this->restraint_lever = restraint_lever;
+}
+
+/** Return how the restraint's six degrees of freedom are grouped into
+ *  lambda-addressable OpenMM Forces, either "combined" or "split". */
+QString BoreschRestraints::restraintLever() const
+{
+    return this->restraint_lever;
+}

corelib/src/libs/SireMM/boreschrestraints.h

@@ -197,12 +197,35 @@ namespace SireMM
         void setUsesPbc(bool use_pbc);
         bool usesPbc() const;
 
+        void setAnglePotential(const QString &angle_potential);
+        QString anglePotential() const;
+
+        void setRestraintLever(const QString &restraint_lever);
+        QString restraintLever() const;
+
     private:
         /** The actual list of restraints*/
         QList<BoreschRestraint> r;
 
         /** Whether the restraints use periodic boundary conditions */
         bool use_pbc = false;
+
+        /** The functional form used for the two Boresch angle restraint
+         *  terms (thetaA, thetaB). Either "harmonic" (default) or
+         *  "restricted_bending" - the latter uses
+         *  k*(cos(theta)-cos(theta0))^2/sin(theta)^2, which diverges as
+         *  theta approaches 0 or pi, preventing the restraint angles from
+         *  ever reaching the Boresch collinearity singularity. */
+        QString angle_potential = "harmonic";
+
+        /** How the restraint's six degrees of freedom are grouped into
+         *  lambda-addressable OpenMM Forces. Either "combined" (default,
+         *  all six terms share a single scale factor / lever) or "split"
+         *  (the distance and two angle terms share one scale factor, the
+         *  three dihedral terms share a second, independent scale factor -
+         *  allowing them to be turned on according to different lambda
+         *  schedules). */
+        QString restraint_lever = "combined";
     };
 
 }

doc/source/changelog.rst

@@ -26,6 +26,21 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
   name '_fix_siremm' from 'sire.mm'``) if it was the first ``sire`` submodule touched in
   a process, due to ``use_new_api()`` being called before ``_fix_siremm`` was defined.
 
+* Added an optional ``angle_potential="restricted_bending"`` mode to ``BoreschRestraints``
+  (default remains ``"harmonic"``), which uses a ``sin(theta)^2``-weighted potential for the
+  two angle restraint terms to prevent the restraint angle from ever reaching the Boresch
+  collinearity singularity at 0/180 degrees.
+
+* Added ``sire.restraints.boresch_search()``, which automatically generates a
+  ``BoreschRestraints`` object and its standard state correction from a trajectory of
+  a protein-ligand complex, using either a hydrogen-bond-driven anchor search or a
+  reference distance-variance-driven protocol.
+
+* Added an optional ``restraint_lever="split"`` mode to ``BoreschRestraints`` (default
+  remains ``"combined"``), which converts the distance/angle and dihedral restraint
+  terms into two independently lambda-addressable OpenMM Forces, allowing them to be
+  turned on according to different lambda schedule equations.
+
 `2026.1.0 <https://github.com/openbiosim/sire/compare/2025.4.0...2026.1.0>`__ - June 2026
 -----------------------------------------------------------------------------------------
 

doc/source/tutorial/part06/03_restraints.rst

@@ -428,7 +428,7 @@ BoreschRestraint( [1574, 1554, 1576] => [4, 3, 5],
                   k=[5 kcal mol-1 Å-2, 0.0152309 kcal mol-1 °-2, 0.0152309 kcal mol-1 °-2,
                   0.0152309 kcal mol-1 °-2, 0.0152309 kcal mol-1 °-2, 0.0152309 kcal mol-1 °-2]
                   r0=15.1197 Å, theta0=[80.5212°, 59.818°],
-                  phi0=[170.562°Ⱐ128.435°Ⱐ192.21°] )
+                  phi0=[170.562°, 128.435°, 192.21°] )
 
 creates a Boresch restraint where the receptor anchor atoms are r1 = 1574, r2 = 1554, and r3 = 1576,
 and the ligand anchor atoms are l1 = 4, l2 = 3, and l3 = 5. The default half force constants have been set
@@ -461,14 +461,17 @@ Alternatively, we could have explicitly set the half force constants and equilib
 BoreschRestraint( [1574, 1554, 1576] => [4, 3, 5],
                   k=[6.2012 kcal mol-1 Å-2, 0.00876339 kcal mol-1 °-2, 0.00756073 kcal mol-1 °-2, 0.0182352 kcal mol-1 °-2, 0.000241348 kcal mol-1 °-2, 0.016808 kcal mol-1 °-2]
                   r0=16 Å, theta0=[68.7549°, 74.4845°],
-                  phi0=[126.051°Ⱐ143.239°Ⱐ85.9437°] )
+                  phi0=[126.051°, 143.239°, 85.9437°] )
 
 .. note::
 
    :func:`sire.restraints.boresch` returns a list of Boresch restraints. If you are only
    interested in a single Boresch restraint, you can extract it with the index, e.g.
    ``boresch_restraint = boresch_restraints[0]``.
 
+   :func:`sire.restraints.boresch` always creates just a single restraint, but the
+   returned :class:`~sire.mm.BoreschRestraints` container can hold several.
+
 When performing an alchemical absolute binding free energy calculation, it is necessary to
 calculate the free energy of releasing the decoupled ligand to the standard state volume.
 The analytical Boresch correction is almost always accurate if stable restraints have been
@@ -481,6 +484,42 @@ selected (see 10.26434/chemrxiv-2023-8s9dz-v2). This can be calculated with
 >>> print(correction)
 -6.2399 kcal mol-1
 
+Automatic Boresch Restraint Search
+-----------------------------------
+
+Choosing suitable Boresch anchor atoms and equilibrium values by hand is
+tedious and error-prone. The :func:`sire.restraints.boresch_search` function
+automates this by analysing a trajectory of the protein-ligand complex
+(e.g. from a short run at :math:`\lambda=0`) and returning a
+:class:`~sire.mm.BoreschRestraints` object together with its standard state
+correction, ready to pass to ``restraints``.
+
+>>> mols = sr.load_test_files("boresch_restraints.prm7", "boresch_restraints.dcd")
+>>> mols.update(sr.morph.decouple(mols.molecule(1), as_new_molecule=False))  # molecule 1 is the ligand
+>>> restraints, correction = sr.restraints.boresch_search(mols, temperature="298 K")
+>>> print(restraints[0])
+BoreschRestraint( [692, 702, 704] => [1496, 1498, 1499],
+                  k=[1 kcal mol-1 Å-2, 0.0243694 kcal mol-1 °-2, 0.0243694 kcal mol-1 °-2,
+                  0.0243694 kcal mol-1 °-2, 0.0243694 kcal mol-1 °-2, 0.0243694 kcal mol-1 °-2]
+                  r0=4.56908 Å, theta0=[82.5581°, 94.9595°],
+                  phi0=[27.9429°, 125.68°, -107.008°] )
+>>> print(correction)
+-10.5706 kcal mol-1
+
+Two search protocols are available via the ``protocol`` argument. The
+default, ``"rxrx"``, seeds candidate anchor atoms from protein-ligand
+hydrogen bonds and scores them to avoid angles close to the Boresch
+collinearity singularity. The alternative, ``"aldeghi"``, is a reference
+implementation of the distance-variance-driven approach used by
+MDRestraintsGenerator/BioSimSpace, kept for comparison.
+
+.. note::
+
+   :func:`sire.restraints.boresch_search` defaults to
+   ``angle_potential="restricted_bending"`` for the restraints it generates
+   (unlike :func:`sire.restraints.boresch`, which defaults to
+   ``"harmonic"``), since this directly avoids the collinearity instability
+   that the search protocols are otherwise designed to steer away from.
 
 Using restraints in minimisation or dynamics
 --------------------------------------------

src/sire/restraints/CMakeLists.txt

@@ -7,6 +7,7 @@
 # Add your script to this list
 set ( SCRIPTS
         __init__.py
+        _boresch_search.py
         _restraints.py
         _standard_state_correction.py
     )

src/sire/restraints/__init__.py

@@ -5,6 +5,7 @@
     "dihedral",
     "distance",
     "boresch",
+    "boresch_search",
     "inverse_bond",
     "inverse_distance",
     "rmsd",
@@ -25,3 +26,4 @@
     rmsd,
 )
 from ._standard_state_correction import get_standard_state_correction
+from ._boresch_search import boresch_search