From ebd0e09ac9193b7f2a259b7a965deadc4f4f2bec Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 1 Jul 2026 15:57:56 +0100 Subject: [PATCH 1/4] Add optional restricted-bending angle potential to Boresch restraints. --- corelib/src/libs/SireMM/boreschrestraints.cpp | 64 +++- corelib/src/libs/SireMM/boreschrestraints.h | 11 + doc/source/changelog.rst | 5 + src/sire/restraints/_restraints.py | 33 ++ tests/restraints/test_boresch.py | 64 ++++ .../SireOpenMM/sire_to_openmm_system.cpp | 82 +++-- wrapper/MM/BoreschRestraints.pypp.cpp | 340 ++++++++---------- 7 files changed, 366 insertions(+), 233 deletions(-) diff --git a/corelib/src/libs/SireMM/boreschrestraints.cpp b/corelib/src/libs/SireMM/boreschrestraints.cpp index 25be2ebf6..aa8e95599 100644 --- a/corelib/src/libs/SireMM/boreschrestraints.cpp +++ b/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(borrest); + sds >> borrest.receptor_atms >> borrest.ligand_atms >> borrest._r0 >> borrest._theta0 >> borrest._phi0 >> borrest._kr >> borrest._ktheta >> borrest._kphi >> static_cast(borrest); } else throw version_error(v, "1", r_borrest, CODELOC); @@ -345,11 +343,11 @@ static const RegisterMetaType 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 << static_cast(borrests); return ds; @@ -364,16 +362,26 @@ QDataStream &operator>>(QDataStream &ds, BoreschRestraints &borrests) SharedDataStream sds(ds); sds >> borrests.r >> static_cast(borrests); + + borrests.use_pbc = false; + borrests.angle_potential = "harmonic"; } else if (v == 2) { SharedDataStream sds(ds); - sds >> borrests.r >> borrests.use_pbc - >> static_cast(borrests); + sds >> borrests.r >> borrests.use_pbc >> static_cast(borrests); + + borrests.angle_potential = "harmonic"; + } + else if (v == 3) + { + SharedDataStream sds(ds); + + sds >> borrests.r >> borrests.use_pbc >> borrests.angle_potential >> static_cast(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 +434,8 @@ BoreschRestraints::BoreschRestraints(const QString &name, } BoreschRestraints::BoreschRestraints(const BoreschRestraints &other) - : ConcreteProperty(other), r(other.r), use_pbc(other.use_pbc) + : ConcreteProperty(other), r(other.r), use_pbc(other.use_pbc), + angle_potential(other.angle_potential) { } @@ -438,6 +447,7 @@ BoreschRestraints &BoreschRestraints::operator=(const BoreschRestraints &other) { r = other.r; use_pbc = other.use_pbc; + angle_potential = other.angle_potential; Restraints::operator=(other); return *this; } @@ -445,7 +455,7 @@ BoreschRestraints &BoreschRestraints::operator=(const BoreschRestraints &other) 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; } bool BoreschRestraints::operator!=(const BoreschRestraints &other) const @@ -499,11 +509,12 @@ 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\n%5\n)") + .arg(this->name()) + .arg(n) + .arg(this->use_pbc ? "true" : "false") + .arg(this->angle_potential) + .arg(parts.join("\n")); } /** Return whether or not this is empty */ @@ -610,3 +621,26 @@ 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; +} diff --git a/corelib/src/libs/SireMM/boreschrestraints.h b/corelib/src/libs/SireMM/boreschrestraints.h index 7592daf5f..4b54d8938 100644 --- a/corelib/src/libs/SireMM/boreschrestraints.h +++ b/corelib/src/libs/SireMM/boreschrestraints.h @@ -197,12 +197,23 @@ namespace SireMM void setUsesPbc(bool use_pbc); bool usesPbc() const; + void setAnglePotential(const QString &angle_potential); + QString anglePotential() const; + private: /** The actual list of restraints*/ QList 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"; }; } diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 098551274..2170d6565 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -26,6 +26,11 @@ organisation on `GitHub `__. 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. + `2026.1.0 `__ - June 2026 ----------------------------------------------------------------------------------------- diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 51a97b919..3c3db78eb 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -143,6 +143,7 @@ def boresch( name=None, map=None, temperature=u("298 K"), + angle_potential=None, ): """ Create a set of Boresch restraints that will restrain the 6 @@ -232,6 +233,16 @@ def boresch( The temperature to use when checking for unstable restraints. If None, then this will default to 298 K. Default is None. + angle_potential : str, optional + The functional form used for the two angle restraint terms + (thetaA, thetaB), either "harmonic" or "restricted_bending". The + "restricted_bending" form (see the GROMACS manual, "Restricted + Bending Potential") diverges as the angle approaches 0 or pi, + preventing the restraint angles from ever reaching the Boresch + collinearity singularity, at the cost of no longer being a simple + harmonic potential away from theta0. Default is None, which is + equivalent to "harmonic". + Returns ------- BoreschRestraints : SireMM::BoreschRestraints @@ -277,6 +288,16 @@ def boresch( temperature = ( temperature if temperature is not None else map_dict.get("temperature", None) ) + angle_potential = ( + angle_potential + if angle_potential is not None + else map_dict.get("angle_potential", None) + ) + # Values retrieved from the map are wrapped as PropertyName, not a plain + # str, which the strict BoreschRestraints.set_angle_potential(QString) + # signature doesn't accept directly. + if angle_potential is not None: + angle_potential = str(angle_potential) receptor = _to_atoms(mols, receptor) ligand = _to_atoms(mols, ligand) @@ -296,6 +317,15 @@ def boresch( else: use_pbc = False + if angle_potential is not None: + if angle_potential not in ("harmonic", "restricted_bending"): + raise ValueError( + "'angle_potential' must be either 'harmonic' or " + f"'restricted_bending', got {angle_potential!r}" + ) + else: + angle_potential = "harmonic" + from .. import measure default_distance_k = u("5 kcal mol-1 A-2") @@ -449,6 +479,9 @@ def boresch( # Set the use_pbc flag. b.set_uses_pbc(use_pbc) + # Set the functional form used for the two angle restraint terms. + b.set_angle_potential(angle_potential) + return b diff --git a/tests/restraints/test_boresch.py b/tests/restraints/test_boresch.py index 2ca8ca2cc..8365469ba 100644 --- a/tests/restraints/test_boresch.py +++ b/tests/restraints/test_boresch.py @@ -297,3 +297,67 @@ def test_boresch_creation_with_map(thrombin_complex): ) boresch_restraint = boresch_restraints[0] assert boresch_restraint.kr().value() == 44.0 + + +def _make_default_boresch(thrombin_complex, **kwargs): + return boresch( + thrombin_complex, + receptor=thrombin_complex["protein"][ + BORESCH_PARAMS_DEFAULT["receptor_selection"] + ], + ligand=thrombin_complex["resname LIG"][ + BORESCH_PARAMS_DEFAULT["ligand_selection"] + ], + kr=BORESCH_PARAMS_DEFAULT["kr"], + ktheta=BORESCH_PARAMS_DEFAULT["ktheta"], + kphi=BORESCH_PARAMS_DEFAULT["kphi"], + r0=BORESCH_PARAMS_DEFAULT["r0"], + theta0=BORESCH_PARAMS_DEFAULT["theta0"], + phi0=BORESCH_PARAMS_DEFAULT["phi0"], + name=BORESCH_PARAMS_DEFAULT["name"], + **kwargs, + ) + + +def test_boresch_angle_potential_defaults_to_harmonic(thrombin_complex): + """ + Check that the angle_potential defaults to "harmonic" when not specified, + matching the pre-existing behaviour (no change for existing callers). + """ + boresch_restraints = _make_default_boresch(thrombin_complex) + assert boresch_restraints.angle_potential() == "harmonic" + + +def test_boresch_angle_potential_restricted_bending(thrombin_complex): + """ + Check that angle_potential="restricted_bending" is set correctly, and + that it doesn't change any of the other restraint parameters. + """ + boresch_restraints = _make_default_boresch( + thrombin_complex, angle_potential="restricted_bending" + ) + assert boresch_restraints.angle_potential() == "restricted_bending" + + boresch_restraint = boresch_restraints[0] + assert boresch_restraint.kr().value() == 6.2012 + assert boresch_restraint.theta0()[0].value() == 1.3031 + assert boresch_restraint.theta0()[1].value() == 1.4777 + + +def test_boresch_angle_potential_invalid_raises(thrombin_complex): + """ + Check that an invalid angle_potential value raises a ValueError. + """ + with pytest.raises(ValueError, match="angle_potential"): + _make_default_boresch(thrombin_complex, angle_potential="not_a_real_potential") + + +def test_boresch_angle_potential_via_map(thrombin_complex): + """ + Check that angle_potential can also be set via the map, matching every + other boresch() option. + """ + boresch_restraints = _make_default_boresch( + thrombin_complex, map={"angle_potential": "restricted_bending"} + ) + assert boresch_restraints.angle_potential() == "restricted_bending" diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 7e1f498b5..387bfe5de 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -94,31 +94,71 @@ void _add_boresch_restraints(const SireMM::BoreschRestraints &restraints, // // e_restraint = rho * (e_bond + e_angle + e_torsion) // e_bond = kr (r - r0)^2 - // e_angle_i = ktheta_i (theta_i - theta0_i)^2 // e_torsion_i = k_phi_i (min(dphi_i, 2pi-dphi_i))^2 where // dphi_i = abs(phi_i - phi0_i) // - const auto energy_expression = QString( - "rho*(e_bond + e_angle_A + e_angle_B + e_torsion_A + e_torsion_B + e_torsion_C);" - "e_bond=kr*(r-r0)^2;" - "e_angle_B=ktheta_B*(theta_B-theta0_B)^2;" - "e_angle_A=ktheta_A*(theta_A-theta0_A)^2;" - "e_torsion_C=kphi_C*(min(dphi_C, two_pi-dphi_C))^2;" - "e_torsion_B=kphi_B*(min(dphi_B, two_pi-dphi_B))^2;" - "e_torsion_A=kphi_A*(min(dphi_A, two_pi-dphi_A))^2;" - "dphi_C=abs(phi_C-phi0_C);" - "dphi_B=abs(phi_B-phi0_B);" - "dphi_A=abs(phi_A-phi0_A);" - "two_pi=6.283185307179586;" - "phi_C=dihedral(p1, p4, p5, p6);" - "phi_B=dihedral(p2, p1, p4, p5);" - "phi_A=dihedral(p3, p2, p1, p4);" - "theta_B=angle(p1, p4, p5);" - "theta_A=angle(p2, p1, p4);" - "r=distance(p1, p4);") - .toStdString(); + // The two angle terms use one of two functional forms, selected by + // restraints.anglePotential(): + // + // "harmonic" (default): e_angle_i = ktheta_i (theta_i - theta0_i)^2 + // + // "restricted_bending": e_angle_i = ktheta_i * (cos(theta_i) - cos(theta0_i))^2 / sin(theta_i)^2 + // (see the GROMACS manual, "Restricted Bending Potential"). The + // sin(theta)^2 denominator diverges as theta approaches 0 or pi, + // preventing the restraint angle from ever reaching the Boresch + // collinearity singularity, unlike the harmonic form which only + // penalises deviation from theta0 with no protection near the poles. + // To leading order in a large-force-constant expansion around theta0 + // (away from the poles) this form has the same local curvature as + // the harmonic form with the same ktheta, so no change is needed to + // the analytic standard state correction formula. + // + QString energy_expression; + + if (restraints.anglePotential() == "restricted_bending") + { + energy_expression = QString( + "rho*(e_bond + e_angle_A + e_angle_B + e_torsion_A + e_torsion_B + e_torsion_C);" + "e_bond=kr*(r-r0)^2;" + "e_angle_A=ktheta_A*(cos(theta_A)-cos(theta0_A))^2/sin(theta_A)^2;" + "e_angle_B=ktheta_B*(cos(theta_B)-cos(theta0_B))^2/sin(theta_B)^2;" + "e_torsion_C=kphi_C*(min(dphi_C, two_pi-dphi_C))^2;" + "e_torsion_B=kphi_B*(min(dphi_B, two_pi-dphi_B))^2;" + "e_torsion_A=kphi_A*(min(dphi_A, two_pi-dphi_A))^2;" + "dphi_C=abs(phi_C-phi0_C);" + "dphi_B=abs(phi_B-phi0_B);" + "dphi_A=abs(phi_A-phi0_A);" + "two_pi=6.283185307179586;" + "phi_C=dihedral(p1, p4, p5, p6);" + "phi_B=dihedral(p2, p1, p4, p5);" + "phi_A=dihedral(p3, p2, p1, p4);" + "theta_B=angle(p1, p4, p5);" + "theta_A=angle(p2, p1, p4);" + "r=distance(p1, p4);"); + } + else + { + energy_expression = QString( + "rho*(e_bond + e_angle_A + e_angle_B + e_torsion_A + e_torsion_B + e_torsion_C);" + "e_bond=kr*(r-r0)^2;" + "e_angle_B=ktheta_B*(theta_B-theta0_B)^2;" + "e_angle_A=ktheta_A*(theta_A-theta0_A)^2;" + "e_torsion_C=kphi_C*(min(dphi_C, two_pi-dphi_C))^2;" + "e_torsion_B=kphi_B*(min(dphi_B, two_pi-dphi_B))^2;" + "e_torsion_A=kphi_A*(min(dphi_A, two_pi-dphi_A))^2;" + "dphi_C=abs(phi_C-phi0_C);" + "dphi_B=abs(phi_B-phi0_B);" + "dphi_A=abs(phi_A-phi0_A);" + "two_pi=6.283185307179586;" + "phi_C=dihedral(p1, p4, p5, p6);" + "phi_B=dihedral(p2, p1, p4, p5);" + "phi_A=dihedral(p3, p2, p1, p4);" + "theta_B=angle(p1, p4, p5);" + "theta_A=angle(p2, p1, p4);" + "r=distance(p1, p4);"); + } - auto *restraintff = new OpenMM::CustomCompoundBondForce(6, energy_expression); + auto *restraintff = new OpenMM::CustomCompoundBondForce(6, energy_expression.toStdString()); restraintff->setName("BoreschRestraintForce"); restraintff->addPerBondParameter("rho"); diff --git a/wrapper/MM/BoreschRestraints.pypp.cpp b/wrapper/MM/BoreschRestraints.pypp.cpp index da907dcd9..9f57733b8 100644 --- a/wrapper/MM/BoreschRestraints.pypp.cpp +++ b/wrapper/MM/BoreschRestraints.pypp.cpp @@ -2,9 +2,9 @@ // (C) Christopher Woods, GPL >= 3 License -#include "boost/python.hpp" -#include "Helpers/clone_const_reference.hpp" #include "BoreschRestraints.pypp.hpp" +#include "Helpers/clone_const_reference.hpp" +#include "boost/python.hpp" namespace bp = boost::python; @@ -24,7 +24,7 @@ namespace bp = boost::python; #include "boreschrestraints.h" -SireMM::BoreschRestraints __copy__(const SireMM::BoreschRestraints &other){ return SireMM::BoreschRestraints(other); } +SireMM::BoreschRestraints __copy__(const SireMM::BoreschRestraints &other) { return SireMM::BoreschRestraints(other); } #include "Helpers/copy.hpp" @@ -36,232 +36,178 @@ SireMM::BoreschRestraints __copy__(const SireMM::BoreschRestraints &other){ retu #include "Helpers/len.hpp" -void register_BoreschRestraints_class(){ +void register_BoreschRestraints_class() +{ { //::SireMM::BoreschRestraints - typedef bp::class_< SireMM::BoreschRestraints, bp::bases< SireMM::Restraints, SireBase::Property > > BoreschRestraints_exposer_t; - BoreschRestraints_exposer_t BoreschRestraints_exposer = BoreschRestraints_exposer_t( "BoreschRestraints", "This class provides the information for a collection of positional\nrestraints that can be added to a collection of molecues. Each\nrestraint can act on a particle or the centroid of a collection\nof particles. The restaints are spherically symmetric, and\nare either flat-bottom harmonics or harmonic potentials\n", bp::init< >("Null constructor") ); - bp::scope BoreschRestraints_scope( BoreschRestraints_exposer ); - BoreschRestraints_exposer.def( bp::init< QString const & >(( bp::arg("name") ), "") ); - BoreschRestraints_exposer.def( bp::init< SireMM::BoreschRestraint const & >(( bp::arg("restraint") ), "") ); - BoreschRestraints_exposer.def( bp::init< QList< SireMM::BoreschRestraint > const & >(( bp::arg("restraints") ), "") ); - BoreschRestraints_exposer.def( bp::init< QString const &, SireMM::BoreschRestraint const & >(( bp::arg("name"), bp::arg("restraint") ), "") ); - BoreschRestraints_exposer.def( bp::init< QString const &, QList< SireMM::BoreschRestraint > const & >(( bp::arg("name"), bp::arg("restraints") ), "") ); - BoreschRestraints_exposer.def( bp::init< SireMM::BoreschRestraints const & >(( bp::arg("other") ), "") ); + typedef bp::class_> BoreschRestraints_exposer_t; + BoreschRestraints_exposer_t BoreschRestraints_exposer = BoreschRestraints_exposer_t("BoreschRestraints", "This class provides the information for a collection of positional\nrestraints that can be added to a collection of molecues. Each\nrestraint can act on a particle or the centroid of a collection\nof particles. The restaints are spherically symmetric, and\nare either flat-bottom harmonics or harmonic potentials\n", bp::init<>("Null constructor")); + bp::scope BoreschRestraints_scope(BoreschRestraints_exposer); + BoreschRestraints_exposer.def(bp::init((bp::arg("name")), "")); + BoreschRestraints_exposer.def(bp::init((bp::arg("restraint")), "")); + BoreschRestraints_exposer.def(bp::init const &>((bp::arg("restraints")), "")); + BoreschRestraints_exposer.def(bp::init((bp::arg("name"), bp::arg("restraint")), "")); + BoreschRestraints_exposer.def(bp::init const &>((bp::arg("name"), bp::arg("restraints")), "")); + BoreschRestraints_exposer.def(bp::init((bp::arg("other")), "")); { //::SireMM::BoreschRestraints::add - - typedef void ( ::SireMM::BoreschRestraints::*add_function_type)( ::SireMM::BoreschRestraint const & ) ; - add_function_type add_function_value( &::SireMM::BoreschRestraints::add ); - - BoreschRestraints_exposer.def( - "add" - , add_function_value - , ( bp::arg("restraint") ) - , bp::release_gil_policy() - , "Add a restraint onto the list" ); - + + typedef void (::SireMM::BoreschRestraints::*add_function_type)(::SireMM::BoreschRestraint const &); + add_function_type add_function_value(&::SireMM::BoreschRestraints::add); + + BoreschRestraints_exposer.def( + "add", add_function_value, (bp::arg("restraint")), bp::release_gil_policy(), "Add a restraint onto the list"); } { //::SireMM::BoreschRestraints::add - - typedef void ( ::SireMM::BoreschRestraints::*add_function_type)( ::SireMM::BoreschRestraints const & ) ; - add_function_type add_function_value( &::SireMM::BoreschRestraints::add ); - - BoreschRestraints_exposer.def( - "add" - , add_function_value - , ( bp::arg("restraints") ) - , bp::release_gil_policy() - , "Add a restraint onto the list" ); - + + typedef void (::SireMM::BoreschRestraints::*add_function_type)(::SireMM::BoreschRestraints const &); + add_function_type add_function_value(&::SireMM::BoreschRestraints::add); + + BoreschRestraints_exposer.def( + "add", add_function_value, (bp::arg("restraints")), bp::release_gil_policy(), "Add a restraint onto the list"); } { //::SireMM::BoreschRestraints::at - - typedef ::SireMM::BoreschRestraint const & ( ::SireMM::BoreschRestraints::*at_function_type)( int ) const; - at_function_type at_function_value( &::SireMM::BoreschRestraints::at ); - - BoreschRestraints_exposer.def( - "at" - , at_function_value - , ( bp::arg("i") ) - , bp::return_value_policy() - , "Return the ith restraint" ); - + + typedef ::SireMM::BoreschRestraint const &(::SireMM::BoreschRestraints::*at_function_type)(int) const; + at_function_type at_function_value(&::SireMM::BoreschRestraints::at); + + BoreschRestraints_exposer.def( + "at", at_function_value, (bp::arg("i")), bp::return_value_policy(), "Return the ith restraint"); } { //::SireMM::BoreschRestraints::count - - typedef int ( ::SireMM::BoreschRestraints::*count_function_type)( ) const; - count_function_type count_function_value( &::SireMM::BoreschRestraints::count ); - - BoreschRestraints_exposer.def( - "count" - , count_function_value - , bp::release_gil_policy() - , "Return the number of restraints" ); - + + typedef int (::SireMM::BoreschRestraints::*count_function_type)() const; + count_function_type count_function_value(&::SireMM::BoreschRestraints::count); + + BoreschRestraints_exposer.def( + "count", count_function_value, bp::release_gil_policy(), "Return the number of restraints"); } { //::SireMM::BoreschRestraints::usesPbc - - typedef bool ( ::SireMM::BoreschRestraints::*usesPbc_function_type)( ) const; - usesPbc_function_type usesPbc_function_value( &::SireMM::BoreschRestraints::usesPbc ); - - BoreschRestraints_exposer.def( - "usesPbc" - , usesPbc_function_value - , bp::release_gil_policy() - , "Return whether or not periodic boundary conditions are to be used" ); - + + typedef bool (::SireMM::BoreschRestraints::*usesPbc_function_type)() const; + usesPbc_function_type usesPbc_function_value(&::SireMM::BoreschRestraints::usesPbc); + + BoreschRestraints_exposer.def( + "usesPbc", usesPbc_function_value, bp::release_gil_policy(), "Return whether or not periodic boundary conditions are to be used"); + } + { //::SireMM::BoreschRestraints::anglePotential + + typedef ::QString (::SireMM::BoreschRestraints::*anglePotential_function_type)() const; + anglePotential_function_type anglePotential_function_value(&::SireMM::BoreschRestraints::anglePotential); + + BoreschRestraints_exposer.def( + "anglePotential", anglePotential_function_value, bp::release_gil_policy(), "Return the functional form used for the two Boresch angle restraint terms,\neither \"harmonic\" or \"restricted_bending\"."); } { //::SireMM::BoreschRestraints::isEmpty - - typedef bool ( ::SireMM::BoreschRestraints::*isEmpty_function_type)( ) const; - isEmpty_function_type isEmpty_function_value( &::SireMM::BoreschRestraints::isEmpty ); - - BoreschRestraints_exposer.def( - "isEmpty" - , isEmpty_function_value - , bp::release_gil_policy() - , "Return whether or not this is empty" ); - + + typedef bool (::SireMM::BoreschRestraints::*isEmpty_function_type)() const; + isEmpty_function_type isEmpty_function_value(&::SireMM::BoreschRestraints::isEmpty); + + BoreschRestraints_exposer.def( + "isEmpty", isEmpty_function_value, bp::release_gil_policy(), "Return whether or not this is empty"); } { //::SireMM::BoreschRestraints::isNull - - typedef bool ( ::SireMM::BoreschRestraints::*isNull_function_type)( ) const; - isNull_function_type isNull_function_value( &::SireMM::BoreschRestraints::isNull ); - - BoreschRestraints_exposer.def( - "isNull" - , isNull_function_value - , bp::release_gil_policy() - , "Return whether or not this is empty" ); - + + typedef bool (::SireMM::BoreschRestraints::*isNull_function_type)() const; + isNull_function_type isNull_function_value(&::SireMM::BoreschRestraints::isNull); + + BoreschRestraints_exposer.def( + "isNull", isNull_function_value, bp::release_gil_policy(), "Return whether or not this is empty"); } { //::SireMM::BoreschRestraints::nRestraints - - typedef int ( ::SireMM::BoreschRestraints::*nRestraints_function_type)( ) const; - nRestraints_function_type nRestraints_function_value( &::SireMM::BoreschRestraints::nRestraints ); - - BoreschRestraints_exposer.def( - "nRestraints" - , nRestraints_function_value - , bp::release_gil_policy() - , "Return the number of restraints" ); - + + typedef int (::SireMM::BoreschRestraints::*nRestraints_function_type)() const; + nRestraints_function_type nRestraints_function_value(&::SireMM::BoreschRestraints::nRestraints); + + BoreschRestraints_exposer.def( + "nRestraints", nRestraints_function_value, bp::release_gil_policy(), "Return the number of restraints"); } - BoreschRestraints_exposer.def( bp::self != bp::self ); - BoreschRestraints_exposer.def( bp::self + bp::other< SireMM::BoreschRestraint >() ); - BoreschRestraints_exposer.def( bp::self + bp::self ); + BoreschRestraints_exposer.def(bp::self != bp::self); + BoreschRestraints_exposer.def(bp::self + bp::other()); + BoreschRestraints_exposer.def(bp::self + bp::self); { //::SireMM::BoreschRestraints::operator= - - typedef ::SireMM::BoreschRestraints & ( ::SireMM::BoreschRestraints::*assign_function_type)( ::SireMM::BoreschRestraints const & ) ; - assign_function_type assign_function_value( &::SireMM::BoreschRestraints::operator= ); - - BoreschRestraints_exposer.def( - "assign" - , assign_function_value - , ( bp::arg("other") ) - , bp::return_self< >() - , "" ); - + + typedef ::SireMM::BoreschRestraints &(::SireMM::BoreschRestraints::*assign_function_type)(::SireMM::BoreschRestraints const &); + assign_function_type assign_function_value(&::SireMM::BoreschRestraints::operator=); + + BoreschRestraints_exposer.def( + "assign", assign_function_value, (bp::arg("other")), bp::return_self<>(), ""); } - BoreschRestraints_exposer.def( bp::self == bp::self ); + BoreschRestraints_exposer.def(bp::self == bp::self); { //::SireMM::BoreschRestraints::operator[] - - typedef ::SireMM::BoreschRestraint const & ( ::SireMM::BoreschRestraints::*__getitem___function_type)( int ) const; - __getitem___function_type __getitem___function_value( &::SireMM::BoreschRestraints::operator[] ); - - BoreschRestraints_exposer.def( - "__getitem__" - , __getitem___function_value - , ( bp::arg("i") ) - , bp::return_value_policy() - , "" ); - + + typedef ::SireMM::BoreschRestraint const &(::SireMM::BoreschRestraints::*__getitem___function_type)(int) const; + __getitem___function_type __getitem___function_value(&::SireMM::BoreschRestraints::operator[]); + + BoreschRestraints_exposer.def( + "__getitem__", __getitem___function_value, (bp::arg("i")), bp::return_value_policy(), ""); } { //::SireMM::BoreschRestraints::restraints - - typedef ::QList< SireMM::BoreschRestraint > ( ::SireMM::BoreschRestraints::*restraints_function_type)( ) const; - restraints_function_type restraints_function_value( &::SireMM::BoreschRestraints::restraints ); - - BoreschRestraints_exposer.def( - "restraints" - , restraints_function_value - , bp::release_gil_policy() - , "Return all of the restraints" ); - + + typedef ::QList (::SireMM::BoreschRestraints::*restraints_function_type)() const; + restraints_function_type restraints_function_value(&::SireMM::BoreschRestraints::restraints); + + BoreschRestraints_exposer.def( + "restraints", restraints_function_value, bp::release_gil_policy(), "Return all of the restraints"); } { //::SireMM::BoreschRestraints::setUsesPbc - - typedef void ( ::SireMM::BoreschRestraints::*setUsesPbc_function_type)( bool ) ; - setUsesPbc_function_type setUsesPbc_function_value( &::SireMM::BoreschRestraints::setUsesPbc ); - - BoreschRestraints_exposer.def( - "setUsesPbc" - , setUsesPbc_function_value - , ( bp::arg("use_pbc") ) - , bp::release_gil_policy() - , "Set whether or not periodic boundary conditions are to be used" ); - + + typedef void (::SireMM::BoreschRestraints::*setUsesPbc_function_type)(bool); + setUsesPbc_function_type setUsesPbc_function_value(&::SireMM::BoreschRestraints::setUsesPbc); + + BoreschRestraints_exposer.def( + "setUsesPbc", setUsesPbc_function_value, (bp::arg("use_pbc")), bp::release_gil_policy(), "Set whether or not periodic boundary conditions are to be used"); + } + { //::SireMM::BoreschRestraints::setAnglePotential + + typedef void (::SireMM::BoreschRestraints::*setAnglePotential_function_type)(::QString const &); + setAnglePotential_function_type setAnglePotential_function_value(&::SireMM::BoreschRestraints::setAnglePotential); + + BoreschRestraints_exposer.def( + "setAnglePotential", setAnglePotential_function_value, (bp::arg("angle_potential")), bp::release_gil_policy(), "Set the functional form used for the two Boresch angle restraint terms.\nMust be either \"harmonic\" or \"restricted_bending\"."); } { //::SireMM::BoreschRestraints::size - - typedef int ( ::SireMM::BoreschRestraints::*size_function_type)( ) const; - size_function_type size_function_value( &::SireMM::BoreschRestraints::size ); - - BoreschRestraints_exposer.def( - "size" - , size_function_value - , bp::release_gil_policy() - , "Return the number of restraints" ); - + + typedef int (::SireMM::BoreschRestraints::*size_function_type)() const; + size_function_type size_function_value(&::SireMM::BoreschRestraints::size); + + BoreschRestraints_exposer.def( + "size", size_function_value, bp::release_gil_policy(), "Return the number of restraints"); } { //::SireMM::BoreschRestraints::toString - - typedef ::QString ( ::SireMM::BoreschRestraints::*toString_function_type)( ) const; - toString_function_type toString_function_value( &::SireMM::BoreschRestraints::toString ); - - BoreschRestraints_exposer.def( - "toString" - , toString_function_value - , bp::release_gil_policy() - , "" ); - + + typedef ::QString (::SireMM::BoreschRestraints::*toString_function_type)() const; + toString_function_type toString_function_value(&::SireMM::BoreschRestraints::toString); + + BoreschRestraints_exposer.def( + "toString", toString_function_value, bp::release_gil_policy(), ""); } { //::SireMM::BoreschRestraints::typeName - - typedef char const * ( *typeName_function_type )( ); - typeName_function_type typeName_function_value( &::SireMM::BoreschRestraints::typeName ); - - BoreschRestraints_exposer.def( - "typeName" - , typeName_function_value - , bp::release_gil_policy() - , "" ); - + + typedef char const *(*typeName_function_type)(); + typeName_function_type typeName_function_value(&::SireMM::BoreschRestraints::typeName); + + BoreschRestraints_exposer.def( + "typeName", typeName_function_value, bp::release_gil_policy(), ""); } { //::SireMM::BoreschRestraints::what - - typedef char const * ( ::SireMM::BoreschRestraints::*what_function_type)( ) const; - what_function_type what_function_value( &::SireMM::BoreschRestraints::what ); - - BoreschRestraints_exposer.def( - "what" - , what_function_value - , bp::release_gil_policy() - , "" ); - - } - BoreschRestraints_exposer.staticmethod( "typeName" ); - BoreschRestraints_exposer.def( "__copy__", &__copy__); - BoreschRestraints_exposer.def( "__deepcopy__", &__copy__); - BoreschRestraints_exposer.def( "clone", &__copy__); - BoreschRestraints_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::BoreschRestraints >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); - BoreschRestraints_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::BoreschRestraints >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); - BoreschRestraints_exposer.def_pickle(sire_pickle_suite< ::SireMM::BoreschRestraints >()); - BoreschRestraints_exposer.def( "__str__", &__str__< ::SireMM::BoreschRestraints > ); - BoreschRestraints_exposer.def( "__repr__", &__str__< ::SireMM::BoreschRestraints > ); - BoreschRestraints_exposer.def( "__len__", &__len_size< ::SireMM::BoreschRestraints > ); - } + typedef char const *(::SireMM::BoreschRestraints::*what_function_type)() const; + what_function_type what_function_value(&::SireMM::BoreschRestraints::what); + + BoreschRestraints_exposer.def( + "what", what_function_value, bp::release_gil_policy(), ""); + } + BoreschRestraints_exposer.staticmethod("typeName"); + BoreschRestraints_exposer.def("__copy__", &__copy__); + BoreschRestraints_exposer.def("__deepcopy__", &__copy__); + BoreschRestraints_exposer.def("clone", &__copy__); + BoreschRestraints_exposer.def("__rlshift__", &__rlshift__QDataStream<::SireMM::BoreschRestraints>, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1, 2>>()); + BoreschRestraints_exposer.def("__rrshift__", &__rrshift__QDataStream<::SireMM::BoreschRestraints>, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1, 2>>()); + BoreschRestraints_exposer.def_pickle(sire_pickle_suite<::SireMM::BoreschRestraints>()); + BoreschRestraints_exposer.def("__str__", &__str__<::SireMM::BoreschRestraints>); + BoreschRestraints_exposer.def("__repr__", &__str__<::SireMM::BoreschRestraints>); + BoreschRestraints_exposer.def("__len__", &__len_size<::SireMM::BoreschRestraints>); + } } From 7560b0b2d11c78dfca5dbd8a796c4fdb1bf48095 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 2 Jul 2026 09:38:17 +0100 Subject: [PATCH 2/4] Add sire.restraints.boresch_search() for automatic Boresch restraint search --- doc/source/changelog.rst | 5 + doc/source/tutorial/part06/03_restraints.rst | 43 +- src/sire/restraints/CMakeLists.txt | 1 + src/sire/restraints/__init__.py | 2 + src/sire/restraints/_boresch_search.py | 1381 ++++++++++++++++++ tests/restraints/test_boresch_search.py | 160 ++ 6 files changed, 1590 insertions(+), 2 deletions(-) create mode 100644 src/sire/restraints/_boresch_search.py create mode 100644 tests/restraints/test_boresch_search.py diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 2170d6565..d323b624e 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -31,6 +31,11 @@ organisation on `GitHub `__. 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. + `2026.1.0 `__ - June 2026 ----------------------------------------------------------------------------------------- diff --git a/doc/source/tutorial/part06/03_restraints.rst b/doc/source/tutorial/part06/03_restraints.rst index b7a69da47..45eb47657 100644 --- a/doc/source/tutorial/part06/03_restraints.rst +++ b/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,7 +461,7 @@ 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:: @@ -469,6 +469,9 @@ BoreschRestraint( [1574, 1554, 1576] => [4, 3, 5], 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 -------------------------------------------- diff --git a/src/sire/restraints/CMakeLists.txt b/src/sire/restraints/CMakeLists.txt index c6fa94c80..abe7c6865 100644 --- a/src/sire/restraints/CMakeLists.txt +++ b/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 ) diff --git a/src/sire/restraints/__init__.py b/src/sire/restraints/__init__.py index 834f33080..07371e06a 100644 --- a/src/sire/restraints/__init__.py +++ b/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 diff --git a/src/sire/restraints/_boresch_search.py b/src/sire/restraints/_boresch_search.py new file mode 100644 index 000000000..12a25dc1d --- /dev/null +++ b/src/sire/restraints/_boresch_search.py @@ -0,0 +1,1381 @@ +""" +Automatic Boresch restraint generation for ABFE simulations. +""" + +__all__ = ["boresch_search"] + +from collections import deque as _deque + +import numpy as _np + +import sire as _sr + +# --------------------------------------------------------------------------- +# Shared helpers, used by both the "rxrx" and "aldeghi" search protocols. +# --------------------------------------------------------------------------- + + +def _nonH_bonded(connectivity, at_idx, mol, is_lig, ghost_elem, h_elem): + """Non-H atoms bonded to at_idx (lambda=0 non-ghost for ligand atoms).""" + bonded = connectivity.connections_to(at_idx) + result = [] + for b_idx in bonded: + if is_lig: + elem = mol.atom(b_idx).property("element0") + if elem == ghost_elem or elem == h_elem: + continue + else: + if mol.atom(b_idx).property("element") == h_elem: + continue + result.append(b_idx) + return result + + +def _build_triplets(connectivity, anchor_idx, mol, is_lig, ghost_elem, h_elem): + """ + All bonded chains anchor-a2-a3 of non-H (non-ghost) atoms. + + A given anchor can have several bonded heavy-atom neighbours (e.g. ring + branch points), and picking an arbitrary one can point a2/a3 back into + the parent molecule, making the corresponding Boresch angle collinear + with the restraint distance vector. All chains are returned so that the + caller's angle filter/scoring can reject bad geometries in favour of a + better one, rather than committing to a single (possibly degenerate) + choice up front. + """ + triplets = [] + for a2 in _nonH_bonded(connectivity, anchor_idx, mol, is_lig, ghost_elem, h_elem): + for a3 in _nonH_bonded(connectivity, a2, mol, is_lig, ghost_elem, h_elem): + if a3 != anchor_idx: + triplets.append((anchor_idx, a2, a3)) + return triplets + + +def _circular_mean_std(values): + """ + Circular mean/std (radians) for periodic dihedral samples. + + Dihedral angles wrap at +-pi, so a plain arithmetic mean/std is wrong + whenever the true distribution straddles the wrap point (e.g. a mean + dihedral near +-180 degrees): samples flip between +179 and -179 + degrees, giving a mean near 0 and a hugely inflated std. This uses the + standard circular statistics (mean of unit vectors, then wrapped + deviation) instead. + """ + mean = float(_np.arctan2(_np.sin(values).mean(), _np.cos(values).mean())) + dtheta = _np.abs(values - mean) + dtheta = _np.minimum(dtheta, 2.0 * _np.pi - dtheta) + std = float(_np.sqrt(_np.mean(dtheta**2))) + return mean, std + + +def _sample_dofs(system, pert_mol_num, sextuplets, n_frames): + """ + Second trajectory pass: sample the six Boresch DOFs for every candidate + sextuplet. Returns arrays of shape (n_sextuplets, n_frames). + """ + n_sext = len(sextuplets) + + dof_r = _np.zeros((n_sext, n_frames)) + dof_tA = _np.zeros((n_sext, n_frames)) # thetaA: r2-r1-l1 + dof_tB = _np.zeros((n_sext, n_frames)) # thetaB: r1-l1-l2 + dof_pA = _np.zeros((n_sext, n_frames)) # phiA: r3-r2-r1-l1 + dof_pB = _np.zeros((n_sext, n_frames)) # phiB: r2-r1-l1-l2 + dof_pC = _np.zeros((n_sext, n_frames)) # phiC: r1-l1-l2-l3 + + for f_idx, frame in enumerate(system.trajectory()): + space = frame.space() + lig_frame = frame.molecule(pert_mol_num) + rec_frame_cache = {} + + for s_idx, s in enumerate(sextuplets): + r1_mol_num = s["r1_mol_num"] + if r1_mol_num not in rec_frame_cache: + rec_frame_cache[r1_mol_num] = frame.molecule(r1_mol_num) + rec_frame = rec_frame_cache[r1_mol_num] + + l1 = lig_frame.atom(s["l1_idx"]) + l2 = lig_frame.atom(s["l2_idx"]) + l3 = lig_frame.atom(s["l3_idx"]) + r1 = rec_frame.atom(s["r1_idx"]) + r2 = rec_frame.atom(s["r2_idx"]) + r3 = rec_frame.atom(s["r3_idx"]) + + dof_r[s_idx, f_idx] = float( + space.calc_dist(l1.coordinates(), r1.coordinates()) + ) + dof_tA[s_idx, f_idx] = float(_sr.measure(r2, r1, l1).value()) + dof_tB[s_idx, f_idx] = float(_sr.measure(r1, l1, l2).value()) + dof_pA[s_idx, f_idx] = float(_sr.measure(r3, r2, r1, l1).value()) + dof_pB[s_idx, f_idx] = float(_sr.measure(r2, r1, l1, l2).value()) + dof_pC[s_idx, f_idx] = float(_sr.measure(r1, l1, l2, l3).value()) + + return dof_r, dof_tA, dof_tB, dof_pA, dof_pB, dof_pC + + +def _assemble_restraints( + system, + pert_mol_num, + s, + r0, + tA0, + tB0, + pA0, + pB0, + pC0, + kr_str, + kt_str, + kp_str, + temperature, + angle_potential, +): + """Build the sire.mm.BoreschRestraints object and standard state correction.""" + from ._restraints import boresch as _boresch + from ._standard_state_correction import get_standard_state_correction as _get_ssc + + # atomidx is molecule-local, so scope with molnum to disambiguate. Comma- + # separated atomidx values preserve listed order in the returned + # SelectorAtom, giving the required [r1, r2, r3] / [l1, l2, l3] sequence. + receptor_atoms = system[ + f"(molnum {s['r1_mol_num'].value()}) and " + f"(atomidx {s['r1_idx'].value()}, {s['r2_idx'].value()}, {s['r3_idx'].value()})" + ].atoms() + ligand_atoms = system[ + f"(molnum {pert_mol_num.value()}) and " + f"(atomidx {s['l1_idx'].value()}, {s['l2_idx'].value()}, {s['l3_idx'].value()})" + ].atoms() + + restraints = _boresch( + system, + receptor=receptor_atoms, + ligand=ligand_atoms, + kr=kr_str, + ktheta=kt_str, + kphi=kp_str, + r0=f"{r0:.6f} A", + theta0=[f"{tA0:.6f} rad", f"{tB0:.6f} rad"], + phi0=[f"{pA0:.6f} rad", f"{pB0:.6f} rad", f"{pC0:.6f} rad"], + temperature=temperature, + angle_potential=angle_potential, + ) + + correction = _get_ssc(restraints[0], temperature) + + return restraints, correction + + +def _parse_unit(value, name, unit, unit_desc): + """ + Parse 'value' (a unit str or a GeneralUnit) into a GeneralUnit with the + dimensions of 'unit', raising a clear error otherwise. Every str-or- + GeneralUnit parameter in this module (temperature, distance cutoffs, + angle cutoffs, force constants, ...) is validated the same way, so that + they can all be supplied consistently as unit strings, e.g. "298 K", + "3.5 A", "120 degrees", "1 kcal mol-1 A-2". + """ + if isinstance(value, str): + try: + value = _sr.u(value) + except Exception as e: + raise ValueError(f"Could not parse '{name}' as a unit: {e}") from e + else: + try: + float(value.to(unit)) + except Exception as e: + raise TypeError( + f"'{name}' must be a str or a GeneralUnit with {unit_desc} " + f"dimensions, got {type(value)}: {e}" + ) from e + return value + + +def _validate_angle_potential(angle_potential): + """ + Validate 'angle_potential', shared by both search protocols. Must be + either "harmonic" or "restricted_bending" (see + sire.restraints.boresch's angle_potential parameter). Unlike + sire.restraints.boresch's own default of "harmonic" (kept for backwards + compatibility with existing callers), boresch_search defaults to + "restricted_bending" for the restraints it generates itself, since the + collinearity singularity this avoids is exactly the failure mode the + restraint search protocols are designed to steer away from in the first + place. + """ + if angle_potential not in ("harmonic", "restricted_bending"): + raise ValueError( + "'angle_potential' must be either 'harmonic' or 'restricted_bending', " + f"got {angle_potential!r}" + ) + + +# --------------------------------------------------------------------------- +# RXRX protocol (default): H-bond-driven restraint search algorithm named +# and described in: +# "Optimizing Absolute Binding Free Energy Calculations for Production Usage" +# https://doi.org/10.1021/acs.jctc.5c00861 +# +# Atoms from protein-ligand hydrogen bonds rather than from bulk distance +# variance, and scores candidates with a formula that explicitly penalises +# theta angles near 0/180 degrees, in order to avoid the Boresch angle +# singularities that plague automated restraint search. +# --------------------------------------------------------------------------- + + +def _boresch_search_rxrx( + system, + temperature="298 K", + protein_selection="not water", + hbond_distance_cutoff="3.5 A", + hbond_angle_cutoff="120 degrees", + occupancy_cutoff=0.5, + restraint_idx=0, + force_constant_r=None, + force_constant_angle=None, + min_frames=50, + angle_potential="restricted_bending", +): + """ + Generate a Boresch restraint using the RXRX restraint search algorithm. + + The system must contain trajectory frames (i.e. be the result of + ``dynamics.commit()`` after a run with ``frame_frequency > 0``). Ligand + anchor atoms are restricted to non-terminal heavy atoms, and candidate + protein-ligand atom sextuplets are seeded from hydrogen bonds with + occupancy above ``occupancy_cutoff``, with the protein anchor mapped to + the Cα of the hydrogen-bonded residue. Among candidates that survive an + angle-range filter, the ligand anchor closest to the ligand's centre of + mass is chosen, and the lowest-scoring sextuplet using that anchor + (Equation 1 of the RXRX paper) is returned. + + Parameters + ---------- + + system : sire.system.System + A Sire system with embedded trajectory frames. Must contain exactly + one perturbable molecule. + + temperature : str or GeneralUnit, optional + Simulation temperature. Defaults to 298 K. + + protein_selection : str + Sire selection string used to identify protein atoms considered as + hydrogen-bond partners. Defaults to all non-water atoms. + + hbond_distance_cutoff : str or GeneralUnit + Maximum donor-acceptor heavy atom distance for a hydrogen bond. + + hbond_angle_cutoff : str or GeneralUnit + Minimum donor-H...acceptor angle for a hydrogen bond. + + occupancy_cutoff : float + Minimum fraction of frames (0-1) in which a hydrogen bond must be + present for it to be considered as a restraint anchor. + + restraint_idx : int + Index into the candidate list (sharing the chosen ligand anchor + atom) sorted by ascending RXRX score (0 = best). + + force_constant_r : str or GeneralUnit, optional + Distance restraint force constant. Defaults to the RXRX protocol + value of 1 kcal mol-1 A-2. + + force_constant_angle : str or GeneralUnit, optional + Angle and dihedral restraint force constant. Defaults to the RXRX + protocol value of 80 kcal mol-1 rad-2. + + min_frames : int + Minimum number of trajectory frames required. + + angle_potential : str + The functional form used for the two Boresch angle restraint terms, + either "harmonic" or "restricted_bending" (see + sire.restraints.boresch). Defaults to "restricted_bending", which + diverges as the angle approaches 0 or 180 degrees, preventing it + from ever reaching the Boresch collinearity singularity - the exact + failure mode this restraint search is designed to avoid. + + Returns + ------- + + restraints : sire.mm.BoreschRestraints + The generated Boresch restraints object, ready to pass to + ``Config.restraints``. + + correction : sire.units.GeneralUnit + Standard state correction in kcal mol-1. + """ + from ..legacy import Mol as _SireMol + from ..system import System as _System + + if not isinstance(system, _System): + raise TypeError( + f"'system' must be of type 'sire.system.System', got {type(system)}" + ) + + temperature = _parse_unit( + temperature, "temperature", _sr.units.kelvin, "temperature" + ) + + if not isinstance(protein_selection, str): + raise TypeError( + f"'protein_selection' must be a str, got {type(protein_selection)}" + ) + + hbond_distance_cutoff = _parse_unit( + hbond_distance_cutoff, "hbond_distance_cutoff", _sr.units.angstrom, "length" + ) + hbond_cutoff_A = float(hbond_distance_cutoff.to(_sr.units.angstrom)) + + hbond_angle_cutoff = _parse_unit( + hbond_angle_cutoff, "hbond_angle_cutoff", _sr.units.degrees, "angle" + ) + hbond_angle_cutoff_deg = float(hbond_angle_cutoff.to(_sr.units.degrees)) + + if not (0.0 < occupancy_cutoff <= 1.0): + raise ValueError( + f"'occupancy_cutoff' must be in (0, 1], got {occupancy_cutoff!r}" + ) + + if not isinstance(restraint_idx, int) or restraint_idx < 0: + raise ValueError( + f"'restraint_idx' must be a non-negative int, got {restraint_idx!r}" + ) + + r_fc_unit = _sr.units.kcal_per_mol / _sr.units.angstrom**2 + angle_fc_unit = _sr.units.kcal_per_mol / _sr.units.radians**2 + + if force_constant_r is not None: + force_constant_r = _parse_unit( + force_constant_r, "force_constant_r", r_fc_unit, "energy length-2" + ) + if force_constant_angle is not None: + force_constant_angle = _parse_unit( + force_constant_angle, + "force_constant_angle", + angle_fc_unit, + "energy angle-2", + ) + + if not isinstance(min_frames, int) or min_frames < 1: + raise ValueError(f"'min_frames' must be a positive int, got {min_frames!r}") + + _validate_angle_potential(angle_potential) + + n_frames = system.num_frames() + if n_frames < min_frames: + raise ValueError( + f"Trajectory has only {n_frames} frame(s); at least {min_frames} are required " + "for reliable Boresch restraint generation." + ) + + ghost_elem = _SireMol.Element(0) + h_elem = _SireMol.Element("H") + n_elem = _SireMol.Element("N") + o_elem = _SireMol.Element("O") + + # ------------------------------------------------------------------------- + # 1. Locate the perturbable molecule (ligand) and its candidate atoms. + # + # The full heavy-atom set is used as the hydrogen-bond donor/acceptor + # pool; the candidate anchor set is restricted to non-terminal atoms + # (bonded to >= 2 non-H heavy atoms), matching the RXRX paper's Figure 2B. + # ------------------------------------------------------------------------- + pert_mols = system.molecules("property is_perturbable") + if pert_mols.num_molecules() != 1: + raise ValueError( + "System must contain exactly one perturbable molecule for Boresch " + f"restraint generation; found {pert_mols.num_molecules()}." + ) + pert_mol = pert_mols.molecule(0) + pert_mol_num = pert_mol.number() + lig_connectivity = pert_mol.connectivity() + + lig_heavy_idxs = [] + for atom in pert_mol.atoms(): + elem0 = atom.property("element0") + if elem0 != ghost_elem and elem0 != h_elem: + lig_heavy_idxs.append(atom.index()) + + if len(lig_heavy_idxs) < 3: + raise ValueError( + f"Ligand has only {len(lig_heavy_idxs)} non-ghost, non-hydrogen atom(s); " + "need at least 3 for Boresch restraints." + ) + + # Lambda=0 atomic masses, used for the mass-weighted centre of mass below + # (topology doesn't change between frames, so this is computed once). + lig_masses = _np.array( + [float(pert_mol.atom(idx).property("mass0").value()) for idx in lig_heavy_idxs] + ) + + candidate_set = set() + for idx in lig_heavy_idxs: + n_heavy = len( + _nonH_bonded(lig_connectivity, idx, pert_mol, True, ghost_elem, h_elem) + ) + if n_heavy >= 2: + candidate_set.add(idx) + + if not candidate_set: + raise ValueError( + "No non-terminal ligand heavy atoms (bonded to >= 2 heavy atoms) " + "were found for RXRX restraint search." + ) + + def _attached_H(connectivity, at_idx, mol, is_lig): + bonded = connectivity.connections_to(at_idx) + result = [] + for b_idx in bonded: + elem = mol.atom(b_idx).property("element0" if is_lig else "element") + if elem == h_elem: + result.append(b_idx) + return result + + lig_hbond_idxs = [ + idx + for idx in lig_heavy_idxs + if pert_mol.atom(idx).property("element0") in (n_elem, o_elem) + ] + if not lig_hbond_idxs: + raise ValueError( + "Ligand has no N/O heavy atoms to use as hydrogen-bond partners." + ) + lig_donor_H = { + idx: _attached_H(lig_connectivity, idx, pert_mol, True) + for idx in lig_hbond_idxs + } + + # ------------------------------------------------------------------------- + # 2. Locate protein hydrogen-bond partner atoms (N/O) and cache each + # residue's C-alpha atom index for the anchor mapping. + # + # 'protein_selection' defaults to the whole protein, which can be + # thousands of atoms; a hydrogen bond to the ligand can only involve + # atoms within a few Angstrom of it, so the N/O search is narrowed with + # a generous spatial "within" clause up front, using Sire's own (native, + # spatially-indexed) selection engine rather than a manual per-atom + # Python distance loop. This is evaluated once, using the system's + # default (first-frame) coordinates. + # ------------------------------------------------------------------------- + _prefilter_margin_A = 15.0 + try: + prot_hbond_sel = system[ + f"(element N,O) and ({protein_selection}) and not (molnum {pert_mol_num.value()}) and " + f"(atoms within {hbond_cutoff_A + _prefilter_margin_A} of (molnum {pert_mol_num.value()}))" + ] + except Exception as e: + raise ValueError( + f"Could not apply protein selection '{protein_selection}': {e}" + ) from e + + prot_mol_cache = {} + prot_hbond_atoms = [] # dicts: mol_num, idx, res_key, attached_H + res_objs = {} # res_key -> Residue view, for the lazy Cα lookup below + for atom in prot_hbond_sel.atoms(): + mn = atom.molecule().number() + if mn not in prot_mol_cache: + prot_mol_cache[mn] = system.molecule(mn) + mol = prot_mol_cache[mn] + residue = atom.residue() + res_key = (mn, residue.number()) + res_objs.setdefault(res_key, residue) + prot_hbond_atoms.append( + { + "mol_num": mn, + "idx": atom.index(), + "res_key": res_key, + "attached_H": _attached_H(mol.connectivity(), atom.index(), mol, False), + } + ) + + if not prot_hbond_atoms: + raise ValueError( + f"No protein N/O atoms found within {hbond_cutoff_A + _prefilter_margin_A:.1f} A " + f"of the ligand for protein selection '{protein_selection}'. Try adjusting " + "'protein_selection' or the hydrogen-bond cutoffs." + ) + + # Cα lookups are only needed for the (typically small) set of residues + # that turn out to be hydrogen-bonded below, so resolve them lazily by + # scanning just that residue's own (small) atom list, rather than every + # atom of every protein molecule up front. + res_ca_idx = {} + + def _get_ca_idx(res_key): + if res_key not in res_ca_idx: + res_ca_idx[res_key] = None + for a in res_objs[res_key].atoms(): + if a.name().value() == "CA": + res_ca_idx[res_key] = a.index() + break + return res_ca_idx[res_key] + + # ------------------------------------------------------------------------- + # 3. First trajectory pass: hydrogen-bond occupancy per (ligand atom, + # protein residue) pair, plus per-frame ligand centre of mass and + # candidate anchor positions (for the CoM tie-break in step 5). + # ------------------------------------------------------------------------- + def _vec3(v): + return _np.array( + [ + float(v.x().to(_sr.units.angstrom)), + float(v.y().to(_sr.units.angstrom)), + float(v.z().to(_sr.units.angstrom)), + ] + ) + + def _min_image_vec3(space, from_coord, to_coord): + """ + Minimum-image displacement vector (Angstrom, numpy) from from_coord + to to_coord. Sire wraps at the atom level (not per-molecule), so even + directly-bonded atoms cannot be assumed to share the same periodic + image; raw coordinate subtraction is not safe. Mirrors the approach + used in BioSimSpace's SireWrapper._getCenterOfMass. + """ + from ..legacy.Maths import Vector as _SireVector + + return _vec3(_SireVector(space.calc_dist_vector(from_coord, to_coord))) + + def _angle_deg(space, d_coord, h_coord, a_coord): + v1 = _min_image_vec3(space, h_coord, d_coord) + v2 = _min_image_vec3(space, h_coord, a_coord) + cosang = _np.dot(v1, v2) / (_np.linalg.norm(v1) * _np.linalg.norm(v2)) + return _np.degrees(_np.arccos(_np.clip(cosang, -1.0, 1.0))) + + occ_counts = {} # (lig_idx, res_key) -> frame count with a satisfied H-bond + candidate_com_dists = {idx: [] for idx in candidate_set} + + for frame in system.trajectory(): + space = frame.space() + lig_frame = frame.molecule(pert_mol_num) + + lig_coords = {idx: lig_frame.atom(idx).coordinates() for idx in lig_heavy_idxs} + ref_coord = lig_coords[lig_heavy_idxs[0]] + ref_vec3 = _vec3(ref_coord) + # Unwrap every ligand atom relative to a single reference atom before + # averaging, so the mass-weighted centre of mass is computed from a + # mutually consistent (non-split-across-the-box) set of positions. + lig_unwrapped = { + idx: ref_vec3 + _min_image_vec3(space, ref_coord, lig_coords[idx]) + for idx in lig_heavy_idxs + } + lig_positions = _np.array([lig_unwrapped[idx] for idx in lig_heavy_idxs]) + com = _np.average(lig_positions, axis=0, weights=lig_masses) + for idx in candidate_set: + candidate_com_dists[idx].append( + float(_np.linalg.norm(lig_unwrapped[idx] - com)) + ) + + prot_frame_cache = {} + + def _prot_frame(mn): + if mn not in prot_frame_cache: + prot_frame_cache[mn] = frame.molecule(mn) + return prot_frame_cache[mn] + + prot_coords = [ + _prot_frame(p["mol_num"]).atom(p["idx"]).coordinates() + for p in prot_hbond_atoms + ] + + pairs_this_frame = set() + + for li in lig_hbond_idxs: + l_coord = lig_coords[li] + l_H_coords = [lig_frame.atom(h).coordinates() for h in lig_donor_H[li]] + + for p_idx, p in enumerate(prot_hbond_atoms): + p_coord = prot_coords[p_idx] + + dist = float(space.calc_dist(l_coord, p_coord)) + if dist > hbond_cutoff_A: + continue + + bonded = False + # Ligand donor -> protein acceptor. + for h_coord in l_H_coords: + if ( + _angle_deg(space, l_coord, h_coord, p_coord) + >= hbond_angle_cutoff_deg + ): + bonded = True + break + if not bonded: + # Protein donor -> ligand acceptor. + p_frame = _prot_frame(p["mol_num"]) + for h_idx in p["attached_H"]: + h_coord = p_frame.atom(h_idx).coordinates() + if ( + _angle_deg(space, p_coord, h_coord, l_coord) + >= hbond_angle_cutoff_deg + ): + bonded = True + break + + if bonded: + pairs_this_frame.add((li, p["res_key"])) + + for pair in pairs_this_frame: + occ_counts[pair] = occ_counts.get(pair, 0) + 1 + + occupied_pairs = [ + pair + for pair, count in occ_counts.items() + if count / n_frames > occupancy_cutoff + ] + + if not occupied_pairs: + raise ValueError( + "No protein-ligand hydrogen bonds exceeded the occupancy cutoff " + f"({occupancy_cutoff}). Try lowering 'occupancy_cutoff' or widening " + "'hbond_distance_cutoff'/'hbond_angle_cutoff'." + ) + + # ------------------------------------------------------------------------- + # 4. Map each occupied (ligand atom, protein residue) pair to a candidate + # six-atom combination: ligand atom -> nearest candidate anchor (BFS + # over connectivity), protein residue -> C-alpha, then extend both + # into bonded triplets. + # ------------------------------------------------------------------------- + def _nearest_candidate(atom_idx, connectivity): + if atom_idx in candidate_set: + return atom_idx + visited = {atom_idx} + queue = _deque([atom_idx]) + while queue: + cur = queue.popleft() + for nb in connectivity.connections_to(cur): + if nb in visited: + continue + if nb in candidate_set: + return nb + visited.add(nb) + queue.append(nb) + return None + + seen_sextuplets = set() + valid_sextuplets = [] + + for li, res_key in occupied_pairs: + l1_idx = _nearest_candidate(li, lig_connectivity) + if l1_idx is None: + continue + + r1_mol_num = res_key[0] + r1_idx = _get_ca_idx(res_key) + if r1_idx is None: + continue + r1_mol = prot_mol_cache[r1_mol_num] + r1_conn = r1_mol.connectivity() + + ligand_triplets = _build_triplets( + lig_connectivity, l1_idx, pert_mol, True, ghost_elem, h_elem + ) + receptor_triplets = _build_triplets( + r1_conn, r1_idx, r1_mol, False, ghost_elem, h_elem + ) + + for l1_idx, l2_idx, l3_idx in ligand_triplets: + for r1_idx, r2_idx, r3_idx in receptor_triplets: + key = (l1_idx, l2_idx, l3_idx, r1_mol_num, r1_idx, r2_idx, r3_idx) + if key in seen_sextuplets: + continue + seen_sextuplets.add(key) + + valid_sextuplets.append( + { + "l1_idx": l1_idx, + "l2_idx": l2_idx, + "l3_idx": l3_idx, + "r1_mol_num": r1_mol_num, + "r1_idx": r1_idx, + "r2_idx": r2_idx, + "r3_idx": r3_idx, + } + ) + + if not valid_sextuplets: + raise ValueError( + "Could not construct a valid Boresch sextuplet from any hydrogen-bonded " + "candidate. Try adjusting 'protein_selection' or the hydrogen-bond cutoffs." + ) + + # ------------------------------------------------------------------------- + # 5. Second trajectory pass: sample all 6 DOFs for every valid sextuplet. + # ------------------------------------------------------------------------- + dof_r, dof_tA, dof_tB, dof_pA, dof_pB, dof_pC = _sample_dofs( + system, pert_mol_num, valid_sextuplets, n_frames + ) + + # ------------------------------------------------------------------------- + # 6. Filter to the RXRX angle window (45-135 degrees), then reduce to the + # single ligand anchor atom closest to the ligand's centre of mass, + # and score the survivors sharing that anchor with Equation 1. + # ------------------------------------------------------------------------- + _min_angle = _np.deg2rad(45.0) + _max_angle = _np.deg2rad(135.0) + + angle_ok_idxs = [ + s_idx + for s_idx in range(len(valid_sextuplets)) + if _min_angle < dof_tA[s_idx].mean() < _max_angle + and _min_angle < dof_tB[s_idx].mean() < _max_angle + ] + + if not angle_ok_idxs: + raise ValueError( + "All hydrogen-bonded Boresch sextuplets were rejected because one or " + "both anchor angles have a mean outside the 45-135 degree RXRX window." + ) + + candidate_l1_idxs = {valid_sextuplets[s_idx]["l1_idx"] for s_idx in angle_ok_idxs} + mean_com_dist = { + idx: float(_np.mean(candidate_com_dists[idx])) for idx in candidate_l1_idxs + } + best_l1 = min(mean_com_dist, key=mean_com_dist.get) + + scored = [] + for s_idx in angle_ok_idxs: + if valid_sextuplets[s_idx]["l1_idx"] != best_l1: + continue + + mean_tA = dof_tA[s_idx].mean() + mean_tB = dof_tB[s_idx].mean() + _, std_pA = _circular_mean_std(dof_pA[s_idx]) + _, std_pB = _circular_mean_std(dof_pB[s_idx]) + _, std_pC = _circular_mean_std(dof_pC[s_idx]) + + score = ( + dof_r[s_idx].std() + * dof_tA[s_idx].std() + * dof_tB[s_idx].std() + * std_pA + * std_pB + * std_pC + * dof_r[s_idx].mean() ** 2 + * (1.0 - abs(_np.sin(mean_tA))) + * (1.0 - abs(_np.sin(mean_tB))) + ) + scored.append((score, s_idx)) + + scored.sort() + + if restraint_idx >= len(scored): + raise ValueError( + f"restraint_idx={restraint_idx} exceeds the number of valid " + f"candidates ({len(scored)})." + ) + + _, best = scored[restraint_idx] + s = valid_sextuplets[best] + + # ------------------------------------------------------------------------- + # 7. Compute equilibrium values (trajectory means, circular for dihedrals) + # and the RXRX protocol force constants. + # ------------------------------------------------------------------------- + r0 = float(dof_r[best].mean()) + tA0 = float(dof_tA[best].mean()) + tB0 = float(dof_tB[best].mean()) + pA0, _ = _circular_mean_std(dof_pA[best]) + pB0, _ = _circular_mean_std(dof_pB[best]) + pC0, _ = _circular_mean_std(dof_pC[best]) + + if force_constant_r is None: + kr_val = 1.0 + else: + kr_val = float(force_constant_r.to(r_fc_unit)) + + if force_constant_angle is None: + ka_val = 80.0 + else: + ka_val = float(force_constant_angle.to(angle_fc_unit)) + + kr_str = f"{kr_val:.6f} kcal mol-1 A-2" + kt_str = [f"{ka_val:.6f} kcal mol-1 rad-2"] * 2 + kp_str = [f"{ka_val:.6f} kcal mol-1 rad-2"] * 3 + + return _assemble_restraints( + system, + pert_mol_num, + s, + r0, + tA0, + tB0, + pA0, + pB0, + pC0, + kr_str, + kt_str, + kp_str, + temperature, + angle_potential, + ) + + +# --------------------------------------------------------------------------- +# Aldeghi protocol (reference implementation): candidate sextuplets are +# generated from bulk ligand-receptor distance variance rather than +# hydrogen-bond occupancy, matching the MDRestraintsGenerator/BioSimSpace +# approach. Kept for comparison; RXRX is the default protocol. +# --------------------------------------------------------------------------- + + +def _boresch_search_aldeghi( + system, + temperature="298 K", + receptor_selection="(not water) and (atomidx > 1) and (atomname CA, C, N)", + cutoff="10 A", + restraint_idx=0, + force_constant=None, + max_candidates=100, + min_frames=50, + angle_potential="restricted_bending", +): + """ + Generate a Boresch restraint for an ABFE simulation by analysing a + trajectory of the protein-ligand complex, using the Aldeghi-style + (MDRestraintsGenerator/BioSimSpace) restraint search protocol. + + The system must contain trajectory frames (i.e. be the result of + ``dynamics.commit()`` after a run with ``frame_frequency > 0``). The + perturbable molecule is used as the ligand; the receptor anchor atoms are + selected via ``receptor_selection``. + + The six Boresch degrees of freedom are sampled over all frames. Candidate + sextuplets are scored by configurational volume (lower = tighter restraint) + and the winner is used to construct a ``sire.mm.BoreschRestraints`` object. + Force constants are derived from per-DOF trajectory variance via the + equipartition theorem unless ``force_constant`` is given. + + Parameters + ---------- + + system : sire.system.System + A Sire system with embedded trajectory frames. Must contain exactly + one perturbable molecule. + + temperature : str or GeneralUnit, optional + Simulation temperature. Defaults to 298 K. + + receptor_selection : str + Sire selection string for receptor anchor atom candidates. + Defaults to backbone heavy atoms (CA, C, N) in non-water molecules, + which is appropriate for AMBER-format input. + + cutoff : str or GeneralUnit + Maximum mean ligand-receptor anchor distance. Pairs whose mean + distance exceeds this value are excluded. + + restraint_idx : int + Index into the candidate list sorted by ascending configurational + volume (0 = tightest). + + force_constant : str or GeneralUnit, optional + Override for all force constants. If None (default), force constants + are fitted from trajectory variance. + + max_candidates : int + Maximum number of (l1, r1) pairs to evaluate for full DOF sampling + in the second trajectory pass. The pairs with the lowest distance + variance are evaluated first. + + min_frames : int + Minimum number of trajectory frames required. Raise ``ValueError`` + if the trajectory has fewer frames. Default 50. + + angle_potential : str + The functional form used for the two Boresch angle restraint terms, + either "harmonic" or "restricted_bending" (see + sire.restraints.boresch). Defaults to "restricted_bending". + + Returns + ------- + + restraints : sire.mm.BoreschRestraints + The generated Boresch restraints object, ready to pass to + ``Config.restraints``. + + correction : sire.units.GeneralUnit + Standard state correction in kcal mol-1. + """ + + from ..legacy import Mol as _SireMol + from ..system import System as _System + from ..legacy.Units import k_boltz as _k_boltz + + # ------------------------------------------------------------------------- + # Parameter validation. + # ------------------------------------------------------------------------- + if not isinstance(system, _System): + raise TypeError( + f"'system' must be of type 'sire.system.System', got {type(system)}" + ) + + temperature = _parse_unit( + temperature, "temperature", _sr.units.kelvin, "temperature" + ) + + if not isinstance(receptor_selection, str): + raise TypeError( + f"'receptor_selection' must be a str, got {type(receptor_selection)}" + ) + + cutoff = _parse_unit(cutoff, "cutoff", _sr.units.angstrom, "length") + + if not isinstance(restraint_idx, int) or restraint_idx < 0: + raise ValueError( + f"'restraint_idx' must be a non-negative int, got {restraint_idx!r}" + ) + + if force_constant is not None: + force_constant = _parse_unit( + force_constant, + "force_constant", + _sr.units.kcal_per_mol / _sr.units.angstrom**2, + "energy length-2", + ) + + if not isinstance(max_candidates, int) or max_candidates < 1: + raise ValueError( + f"'max_candidates' must be a positive int, got {max_candidates!r}" + ) + + if not isinstance(min_frames, int) or min_frames < 1: + raise ValueError(f"'min_frames' must be a positive int, got {min_frames!r}") + + _validate_angle_potential(angle_potential) + + n_frames = system.num_frames() + if n_frames < min_frames: + raise ValueError( + f"Trajectory has only {n_frames} frame(s); at least {min_frames} are required " + "for reliable Boresch restraint generation." + ) + + # kBT as a plain float in kcal mol-1 (for equipartition). + kBT = float((_k_boltz * temperature).to(_sr.units.kcal_per_mol)) + + # Cutoff as a plain float in Angstroms. + cutoff_A = float(cutoff.to(_sr.units.angstrom)) + + # ------------------------------------------------------------------------- + # 1. Locate the perturbable molecule (ligand). + # ------------------------------------------------------------------------- + pert_mols = system.molecules("property is_perturbable") + if pert_mols.num_molecules() != 1: + raise ValueError( + "System must contain exactly one perturbable molecule for Boresch " + f"restraint generation; found {pert_mols.num_molecules()}." + ) + pert_mol = pert_mols.molecule(0) + pert_mol_num = pert_mol.number() + + ghost_elem = _SireMol.Element(0) + h_elem = _SireMol.Element("H") + + # Collect non-ghost, non-H ligand AtomIdx values (lambda=0 state). + lig_atom_idxs = [] + for atom in pert_mol.atoms(): + elem0 = atom.property("element0") + if elem0 != ghost_elem and elem0 != h_elem: + lig_atom_idxs.append(atom.index()) + + if len(lig_atom_idxs) < 3: + raise ValueError( + f"Ligand has only {len(lig_atom_idxs)} non-ghost, non-hydrogen atom(s); " + "need at least 3 for Boresch restraints." + ) + + lig_connectivity = pert_mol.connectivity() + + # ------------------------------------------------------------------------- + # 2. Locate receptor atoms. + # ------------------------------------------------------------------------- + try: + rec_sel = system[receptor_selection] + except Exception as e: + raise ValueError( + f"Could not apply receptor selection '{receptor_selection}': {e}" + ) from e + + # (mol_num, AtomIdx-within-mol) for each receptor atom. + rec_atom_info = [] + for atom in rec_sel.atoms(): + rec_atom_info.append((atom.molecule().number(), atom.index())) + + if len(rec_atom_info) < 3: + raise ValueError( + f"Receptor selection '{receptor_selection}' matched only " + f"{len(rec_atom_info)} atom(s); need at least 3." + ) + + n_lig = len(lig_atom_idxs) + n_rec = len(rec_atom_info) + + # ------------------------------------------------------------------------- + # 3. First trajectory pass: accumulate all ligand-receptor distances. + # ------------------------------------------------------------------------- + frame_dists = [] + + for frame in system.trajectory(): + space = frame.space() + lig_frame = frame.molecule(pert_mol_num) + frame_rec_mol_cache = {} + d = _np.empty((n_lig, n_rec)) + for i, l_idx in enumerate(lig_atom_idxs): + l_coord = lig_frame.atom(l_idx).coordinates() + for j, (r_mol_num, r_idx) in enumerate(rec_atom_info): + if r_mol_num not in frame_rec_mol_cache: + frame_rec_mol_cache[r_mol_num] = frame.molecule(r_mol_num) + r_coord = frame_rec_mol_cache[r_mol_num].atom(r_idx).coordinates() + d[i, j] = float(space.calc_dist(l_coord, r_coord)) + frame_dists.append(d) + + # Shape: (n_frames, n_lig, n_rec) + dists = _np.array(frame_dists) + mean_dists = dists.mean(axis=0) + var_dists = dists.var(axis=0) + + # Candidate (l1, r1) pairs within cutoff, sorted by ascending variance. + candidate_pairs = sorted( + (var_dists[i, j], i, j) + for i in range(n_lig) + for j in range(n_rec) + if mean_dists[i, j] <= cutoff_A + ) + + if not candidate_pairs: + raise ValueError( + f"No ligand-receptor atom pairs found within cutoff {cutoff}. " + "Try increasing the cutoff or adjusting the receptor selection." + ) + + candidate_pairs = candidate_pairs[:max_candidates] + + # ------------------------------------------------------------------------- + # 4. Build valid sextuplets (l1,l2,l3, r1,r2,r3). + # + # Convention (matching sire.restraints.boresch ordering): + # r = dist(l1, r1) + # θA = angle(r2, r1, l1) + # θB = angle(r1, l1, l2) + # φA = dihedral(r3, r2, r1, l1) -- r3 bonded to r2 + # φB = dihedral(r2, r1, l1, l2) + # φC = dihedral(r1, l1, l2, l3) -- l3 bonded to l2 + # ------------------------------------------------------------------------- + rec_mol_cache = {} + + def _rec_mol(mol_num): + if mol_num not in rec_mol_cache: + rec_mol_cache[mol_num] = system.molecule(mol_num) + return rec_mol_cache[mol_num] + + valid_sextuplets = [] + + for _, l1_pos, r1_pos in candidate_pairs: + l1_idx = lig_atom_idxs[l1_pos] + r1_mol_num, r1_idx = rec_atom_info[r1_pos] + r1_mol = _rec_mol(r1_mol_num) + r1_conn = r1_mol.connectivity() + + ligand_triplets = _build_triplets( + lig_connectivity, l1_idx, pert_mol, True, ghost_elem, h_elem + ) + receptor_triplets = _build_triplets( + r1_conn, r1_idx, r1_mol, False, ghost_elem, h_elem + ) + + for l1_idx, l2_idx, l3_idx in ligand_triplets: + for r1_idx, r2_idx, r3_idx in receptor_triplets: + valid_sextuplets.append( + { + "l1_idx": l1_idx, + "l2_idx": l2_idx, + "l3_idx": l3_idx, + "r1_mol_num": r1_mol_num, + "r1_idx": r1_idx, + "r2_idx": r2_idx, + "r3_idx": r3_idx, + } + ) + + if not valid_sextuplets: + raise ValueError( + "Could not construct a valid Boresch sextuplet from the candidate " + "pairs. Try increasing the cutoff or adjusting the receptor selection." + ) + + # ------------------------------------------------------------------------- + # 5. Second trajectory pass: sample all 6 DOFs for every valid sextuplet. + # ------------------------------------------------------------------------- + dof_r, dof_tA, dof_tB, dof_pA, dof_pB, dof_pC = _sample_dofs( + system, pert_mol_num, valid_sextuplets, n_frames + ) + + # ------------------------------------------------------------------------- + # 6. Score sextuplets by configurational volume; apply stability filter. + # + # score ∝ r₀² · |sin θA₀| · |sin θB₀| · σr · σθA · σθB · σφA · σφB · σφC + # (from Boresch et al. 2003, J. Phys. Chem. B, Eqn 32, lower = tighter) + # + # Dihedral (φ) statistics use circular mean/std since these DOFs wrap at + # +-180 degrees; a plain arithmetic std would be wildly inflated whenever + # the true distribution straddles the wrap point. + # ------------------------------------------------------------------------- + _min_angle = _np.deg2rad(10.0) + _max_angle = _np.deg2rad(170.0) + + scored = [] + for s_idx in range(len(valid_sextuplets)): + mean_tA = dof_tA[s_idx].mean() + mean_tB = dof_tB[s_idx].mean() + + if not ( + _min_angle < mean_tA < _max_angle and _min_angle < mean_tB < _max_angle + ): + continue + + _, std_pA = _circular_mean_std(dof_pA[s_idx]) + _, std_pB = _circular_mean_std(dof_pB[s_idx]) + _, std_pC = _circular_mean_std(dof_pC[s_idx]) + + score = ( + dof_r[s_idx].mean() ** 2 + * abs(_np.sin(mean_tA)) + * abs(_np.sin(mean_tB)) + * dof_r[s_idx].std() + * dof_tA[s_idx].std() + * dof_tB[s_idx].std() + * std_pA + * std_pB + * std_pC + ) + scored.append((score, s_idx)) + + if not scored: + raise ValueError( + "All candidate Boresch sextuplets were rejected because one or both " + "anchor angles have a mean near 0° or 180° (collinearity instability). " + "Try increasing the cutoff or adjusting the receptor selection." + ) + + scored.sort() + + if restraint_idx >= len(scored): + raise ValueError( + f"restraint_idx={restraint_idx} exceeds the number of valid " + f"candidates ({len(scored)})." + ) + + _, best = scored[restraint_idx] + s = valid_sextuplets[best] + + # ------------------------------------------------------------------------- + # 7. Compute equilibrium values (trajectory means) and force constants. + # ------------------------------------------------------------------------- + r0 = float(dof_r[best].mean()) + tA0 = float(dof_tA[best].mean()) + tB0 = float(dof_tB[best].mean()) + pA0, std_pA_best = _circular_mean_std(dof_pA[best]) + pB0, std_pB_best = _circular_mean_std(dof_pB[best]) + pC0, std_pC_best = _circular_mean_std(dof_pC[best]) + + if force_constant is not None: + # Use the numeric magnitude for all DOFs with appropriate units per DOF + # type (A-2 for distance, rad-2 for angles). This matches BSS behaviour. + fc_val = float( + force_constant.to(_sr.units.kcal_per_mol / _sr.units.angstrom**2) + ) + kr_str = f"{fc_val:.6f} kcal mol-1 A-2" + kt_str = [f"{fc_val:.6f} kcal mol-1 rad-2"] * 2 + kp_str = [f"{fc_val:.6f} kcal mol-1 rad-2"] * 3 + else: + # Equipartition: k = kBT / (2σ²). + # Sire's boresch() uses E = k·x² (half-spring-constant convention). + def _k(sigma): + return kBT / (2.0 * float(sigma) ** 2) + + kr_str = f"{_k(dof_r[best].std()):.6f} kcal mol-1 A-2" + kt_str = [ + f"{_k(dof_tA[best].std()):.6f} kcal mol-1 rad-2", + f"{_k(dof_tB[best].std()):.6f} kcal mol-1 rad-2", + ] + kp_str = [ + f"{_k(std_pA_best):.6f} kcal mol-1 rad-2", + f"{_k(std_pB_best):.6f} kcal mol-1 rad-2", + f"{_k(std_pC_best):.6f} kcal mol-1 rad-2", + ] + + # ------------------------------------------------------------------------- + # 8. Build the BoreschRestraints object. + # ------------------------------------------------------------------------- + return _assemble_restraints( + system, + pert_mol_num, + s, + r0, + tA0, + tB0, + pA0, + pB0, + pC0, + kr_str, + kt_str, + kp_str, + temperature, + angle_potential, + ) + + +# --------------------------------------------------------------------------- +# Public entry point. +# --------------------------------------------------------------------------- + + +def boresch_search( + system, + protocol="rxrx", + temperature="298 K", + restraint_idx=0, + min_frames=50, + protein_selection="not water", + hbond_distance_cutoff="3.5 A", + hbond_angle_cutoff="120 degrees", + occupancy_cutoff=0.5, + force_constant_r=None, + force_constant_angle=None, + receptor_selection="(not water) and (atomidx > 1) and (atomname CA, C, N)", + cutoff="10 A", + force_constant=None, + max_candidates=100, + angle_potential="restricted_bending", +): + """ + Generate a Boresch restraint for an ABFE simulation by analysing a + trajectory of the protein-ligand complex. + + Two restraint search protocols are available, selected via ``protocol``; + the parameters below are grouped by which protocol(s) use them. + + - ``"rxrx"`` (default): the RXRX restraint search algorithm, which seeds + candidate anchor atoms from protein-ligand hydrogen bonds and scores + them with a formula that explicitly penalises angles near the Boresch + 0/180 degree singularities. + + - ``"aldeghi"``: a reference implementation of the MDRestraintsGenerator/ + BioSimSpace restraint search, kept for comparison. + + Parameters + ---------- + + system : sire.system.System + A Sire system with embedded trajectory frames. Must contain exactly + one perturbable molecule. Used by both protocols. + + protocol : str + The restraint search protocol to use: ``"rxrx"`` (default) or + ``"aldeghi"``. + + temperature : str or GeneralUnit, optional + Simulation temperature. Defaults to 298 K. Used by both protocols. + + restraint_idx : int + Index into the candidate list sorted by ascending score (0 = best). + Used by both protocols. + + min_frames : int + Minimum number of trajectory frames required. Used by both protocols. + + angle_potential : str + The functional form used for the two Boresch angle restraint terms, + either "harmonic" or "restricted_bending" (see + sire.restraints.boresch). Defaults to "restricted_bending", which + diverges as the angle approaches 0 or 180 degrees, preventing it + from ever reaching the Boresch collinearity singularity. Used by + both protocols. + + protein_selection : str + [rxrx only] Sire selection string used to identify protein atoms + considered as hydrogen-bond partners. Defaults to all non-water + atoms. + + hbond_distance_cutoff : str or GeneralUnit + [rxrx only] Maximum donor-acceptor heavy atom distance for a + hydrogen bond. + + hbond_angle_cutoff : str or GeneralUnit + [rxrx only] Minimum donor-H...acceptor angle for a hydrogen bond. + + occupancy_cutoff : float + [rxrx only] Minimum fraction of frames (0-1) in which a hydrogen + bond must be present for it to be considered as a restraint anchor. + + force_constant_r : str or GeneralUnit, optional + [rxrx only] Distance restraint force constant. Defaults to the RXRX + protocol value of 1 kcal mol-1 A-2. + + force_constant_angle : str or GeneralUnit, optional + [rxrx only] Angle and dihedral restraint force constant. Defaults to + the RXRX protocol value of 80 kcal mol-1 rad-2. + + receptor_selection : str + [aldeghi only] Sire selection string for receptor anchor atom + candidates. Defaults to backbone heavy atoms (CA, C, N) in non-water + molecules, which is appropriate for AMBER-format input. + + cutoff : str or GeneralUnit + [aldeghi only] Maximum mean ligand-receptor anchor distance. Pairs + whose mean distance exceeds this value are excluded. + + force_constant : str or GeneralUnit, optional + [aldeghi only] Override for all force constants. If None (default), + force constants are fitted from trajectory variance. + + max_candidates : int + [aldeghi only] Maximum number of (l1, r1) pairs to evaluate for full + DOF sampling in the second trajectory pass. The pairs with the + lowest distance variance are evaluated first. + + Returns + ------- + + restraints : sire.mm.BoreschRestraints + The generated Boresch restraints object, ready to pass to + ``Config.restraints``. + + correction : sire.units.GeneralUnit + Standard state correction in kcal mol-1. + """ + if not isinstance(protocol, str): + raise TypeError(f"'protocol' must be a str, got {type(protocol)}") + + if protocol == "rxrx": + return _boresch_search_rxrx( + system, + temperature=temperature, + protein_selection=protein_selection, + hbond_distance_cutoff=hbond_distance_cutoff, + hbond_angle_cutoff=hbond_angle_cutoff, + occupancy_cutoff=occupancy_cutoff, + restraint_idx=restraint_idx, + force_constant_r=force_constant_r, + force_constant_angle=force_constant_angle, + min_frames=min_frames, + angle_potential=angle_potential, + ) + elif protocol == "aldeghi": + return _boresch_search_aldeghi( + system, + temperature=temperature, + receptor_selection=receptor_selection, + cutoff=cutoff, + restraint_idx=restraint_idx, + force_constant=force_constant, + max_candidates=max_candidates, + min_frames=min_frames, + angle_potential=angle_potential, + ) + else: + raise ValueError( + f"Unknown 'protocol'={protocol!r}; must be 'rxrx' or 'aldeghi'." + ) diff --git a/tests/restraints/test_boresch_search.py b/tests/restraints/test_boresch_search.py new file mode 100644 index 000000000..4ee06d9e1 --- /dev/null +++ b/tests/restraints/test_boresch_search.py @@ -0,0 +1,160 @@ +""" +Unit tests for boresch_search(). + +Test data (boresch_restraints.prm7 + boresch_restraints.dcd) is a lambda=0 +trajectory of the 1jr5 protein with decoupled ligand01. +""" + +import numpy as np +import pytest +import sire as sr +import sire.legacy.MM as _SireMM + +PROTOCOLS = ["rxrx", "aldeghi"] + + +@pytest.fixture(scope="module") +def abfe_system(): + """ + Load the protein-ligand test system with embedded trajectory frames. + The topology is loaded from prm7; the trajectory from DCD. Molecule 0 + is the protein, molecule 1 is the ligand, the rest is water and ions. + The ligand is then decoupled so that is_perturbable is set, matching + what SOMD2's runner provides to boresch_search(). + """ + mols = sr.load_test_files("boresch_restraints.prm7", "boresch_restraints.dcd") + lig_decoupled = sr.morph.decouple(mols.molecule(1), as_new_molecule=False) + mols.update(lig_decoupled) + return mols + + +@pytest.fixture(scope="module", params=PROTOCOLS) +def boresch_result(request, abfe_system): + from sire.restraints import boresch_search + + return boresch_search(abfe_system, protocol=request.param, temperature="298 K") + + +class TestGenerateBoreschRestraint: + def test_returns_tuple(self, boresch_result): + assert isinstance(boresch_result, tuple) and len(boresch_result) == 2 + + def test_restraints_type(self, boresch_result): + restraints, _ = boresch_result + assert isinstance(restraints, _SireMM.BoreschRestraints) + + def test_single_restraint(self, boresch_result): + restraints, _ = boresch_result + assert restraints.at(0) is not None + + def test_angle_potential_defaults_to_restricted_bending(self, boresch_result): + """ + boresch_search() defaults to "restricted_bending" (unlike + sire.restraints.boresch's own default of "harmonic"), since it + directly avoids the collinearity singularity the restraint search + protocols are designed to steer away from in the first place. + """ + restraints, _ = boresch_result + assert restraints.angle_potential() == "restricted_bending" + + def test_correction_is_negative(self, boresch_result): + """Standard state correction is always negative (costs free energy to restrain).""" + _, correction = boresch_result + assert float(correction.to(sr.units.kcal_per_mol)) < 0 + + def test_distance_positive(self, boresch_result): + restraints, _ = boresch_result + r = restraints.at(0) + assert float(r.r0().value()) > 0 + + def test_angles_not_collinear(self, boresch_result): + """Both anchor angles must be away from 0° and 180° (our stability filter).""" + restraints, _ = boresch_result + r = restraints.at(0) + for theta in r.theta0(): + deg = float(theta.to(sr.units.degrees)) + assert 10.0 < deg < 170.0 + + def test_force_constants_positive(self, boresch_result): + restraints, _ = boresch_result + r = restraints.at(0) + assert float(r.kr().value()) > 0 + for k in r.ktheta(): + assert float(k.value()) > 0 + for k in r.kphi(): + assert float(k.value()) > 0 + + def test_receptor_atoms_count(self, boresch_result): + """Receptor selection must yield exactly 3 anchor atoms.""" + restraints, _ = boresch_result + assert len(list(restraints.at(0).receptor_atoms())) == 3 + + def test_ligand_atoms_count(self, boresch_result): + """Ligand selection must yield exactly 3 anchor atoms.""" + restraints, _ = boresch_result + assert len(list(restraints.at(0).ligand_atoms())) == 3 + + @pytest.mark.parametrize("protocol", PROTOCOLS) + def test_restraint_idx_selects_different_candidate(self, abfe_system, protocol): + from sire.restraints import boresch_search + + r0, _ = boresch_search(abfe_system, protocol=protocol, restraint_idx=0) + r1, _ = boresch_search(abfe_system, protocol=protocol, restraint_idx=1) + # Different candidates must differ in the receptor and/or ligand anchor + # atoms: two candidates can share the same receptor anchor triplet + # while using a different ligand triplet branching off the same l1, + # or vice versa. + r0_sextuplet = ( + list(r0.at(0).receptor_atoms()), + list(r0.at(0).ligand_atoms()), + ) + r1_sextuplet = ( + list(r1.at(0).receptor_atoms()), + list(r1.at(0).ligand_atoms()), + ) + assert r0_sextuplet != r1_sextuplet + + @pytest.mark.parametrize("protocol", PROTOCOLS) + def test_too_few_frames_raises(self, abfe_system, protocol): + from sire.restraints import boresch_search + + with pytest.raises(ValueError, match="frame"): + boresch_search(abfe_system, protocol=protocol, min_frames=10_000) + + @pytest.mark.parametrize("protocol", PROTOCOLS) + def test_angle_potential_harmonic_override(self, abfe_system, protocol): + from sire.restraints import boresch_search + + restraints, _ = boresch_search( + abfe_system, protocol=protocol, angle_potential="harmonic" + ) + assert restraints.angle_potential() == "harmonic" + + @pytest.mark.parametrize("protocol", PROTOCOLS) + def test_angle_potential_invalid_raises(self, abfe_system, protocol): + from sire.restraints import boresch_search + + with pytest.raises(ValueError, match="angle_potential"): + boresch_search(abfe_system, protocol=protocol, angle_potential="nonsense") + + def test_force_constant_override(self, abfe_system): + """'force_constant' is an Aldeghi-only kwarg (RXRX uses two separate + force constants, matching the fixed protocol values from the paper).""" + from sire.restraints import boresch_search + + kval = 10.0 + restraints, _ = boresch_search( + abfe_system, + protocol="aldeghi", + force_constant=f"{kval} kcal mol-1 A-2", + ) + r = restraints.at(0) + assert np.isclose(float(r.kr().value()), kval, atol=1e-3) + + def test_tight_cutoff_raises(self, abfe_system): + """'cutoff' is an Aldeghi-only kwarg (RXRX has no equivalent distance + cutoff; candidates are seeded from hydrogen-bond occupancy instead).""" + from sire.restraints import boresch_search + + with pytest.raises(ValueError, match="cutoff"): + boresch_search(abfe_system, protocol="aldeghi", cutoff="0.1 A") From 7fce3e5ffd31def6453047f2032c18114c47b724 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 2 Jul 2026 10:09:08 +0100 Subject: [PATCH 3/4] Fix print formatting error. --- corelib/src/libs/SireMM/boreschrestraints.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/corelib/src/libs/SireMM/boreschrestraints.cpp b/corelib/src/libs/SireMM/boreschrestraints.cpp index aa8e95599..f172db227 100644 --- a/corelib/src/libs/SireMM/boreschrestraints.cpp +++ b/corelib/src/libs/SireMM/boreschrestraints.cpp @@ -283,7 +283,7 @@ QString BoreschRestraint::toString() const .arg(k.join(", ")) .arg(_r0.toString()) .arg(t.join(", ")) - .arg(p.join(', ')); + .arg(p.join(", ")); } } From d6cb58972eff7291778f51604935d59cbef02511 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 2 Jul 2026 13:25:54 +0100 Subject: [PATCH 4/4] Add restraint_lever option to BoreschRestraints for staged RXRX turn-on. --- corelib/src/libs/SireMM/boreschrestraints.cpp | 41 +++- corelib/src/libs/SireMM/boreschrestraints.h | 12 + doc/source/changelog.rst | 5 + src/sire/restraints/_boresch_search.py | 53 +++++ src/sire/restraints/_restraints.py | 39 +++- tests/restraints/test_boresch.py | 44 ++++ tests/restraints/test_boresch_search.py | 31 +++ .../SireOpenMM/sire_to_openmm_system.cpp | 216 +++++++++++------- wrapper/MM/BoreschRestraints.pypp.cpp | 16 ++ 9 files changed, 367 insertions(+), 90 deletions(-) diff --git a/corelib/src/libs/SireMM/boreschrestraints.cpp b/corelib/src/libs/SireMM/boreschrestraints.cpp index f172db227..4c500d337 100644 --- a/corelib/src/libs/SireMM/boreschrestraints.cpp +++ b/corelib/src/libs/SireMM/boreschrestraints.cpp @@ -347,7 +347,7 @@ QDataStream &operator<<(QDataStream &ds, const BoreschRestraints &borrests) SharedDataStream sds(ds); - sds << borrests.r << borrests.use_pbc << borrests.angle_potential + sds << borrests.r << borrests.use_pbc << borrests.angle_potential << borrests.restraint_lever << static_cast(borrests); return ds; @@ -365,6 +365,7 @@ QDataStream &operator>>(QDataStream &ds, BoreschRestraints &borrests) borrests.use_pbc = false; borrests.angle_potential = "harmonic"; + borrests.restraint_lever = "combined"; } else if (v == 2) { @@ -373,12 +374,14 @@ QDataStream &operator>>(QDataStream &ds, BoreschRestraints &borrests) sds >> borrests.r >> borrests.use_pbc >> static_cast(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 >> static_cast(borrests); + sds >> borrests.r >> borrests.use_pbc >> borrests.angle_potential >> borrests.restraint_lever >> + static_cast(borrests); } else throw version_error(v, "1,2,3", r_borrests, CODELOC); @@ -435,7 +438,7 @@ BoreschRestraints::BoreschRestraints(const QString &name, BoreschRestraints::BoreschRestraints(const BoreschRestraints &other) : ConcreteProperty(other), r(other.r), use_pbc(other.use_pbc), - angle_potential(other.angle_potential) + angle_potential(other.angle_potential), restraint_lever(other.restraint_lever) { } @@ -448,6 +451,7 @@ 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; } @@ -455,7 +459,8 @@ BoreschRestraints &BoreschRestraints::operator=(const BoreschRestraints &other) bool BoreschRestraints::operator==(const BoreschRestraints &other) const { return r == other.r and Restraints::operator==(other) and - use_pbc == other.use_pbc and angle_potential == other.angle_potential; + use_pbc == other.use_pbc and angle_potential == other.angle_potential and + restraint_lever == other.restraint_lever; } bool BoreschRestraints::operator!=(const BoreschRestraints &other) const @@ -509,11 +514,13 @@ QString BoreschRestraints::toString() const } } - return QObject::tr("BoreschRestraints( name=%1, size=%2, use_pbc=%3, angle_potential=%4\n%5\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")); } @@ -644,3 +651,27 @@ 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; +} diff --git a/corelib/src/libs/SireMM/boreschrestraints.h b/corelib/src/libs/SireMM/boreschrestraints.h index 4b54d8938..14be7abba 100644 --- a/corelib/src/libs/SireMM/boreschrestraints.h +++ b/corelib/src/libs/SireMM/boreschrestraints.h @@ -200,6 +200,9 @@ namespace SireMM 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 r; @@ -214,6 +217,15 @@ namespace SireMM * 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"; }; } diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index d323b624e..7ec96c694 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -36,6 +36,11 @@ organisation on `GitHub `__. 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 `__ - June 2026 ----------------------------------------------------------------------------------------- diff --git a/src/sire/restraints/_boresch_search.py b/src/sire/restraints/_boresch_search.py index 12a25dc1d..9f83f9c08 100644 --- a/src/sire/restraints/_boresch_search.py +++ b/src/sire/restraints/_boresch_search.py @@ -128,6 +128,7 @@ def _assemble_restraints( kp_str, temperature, angle_potential, + restraint_lever, ): """Build the sire.mm.BoreschRestraints object and standard state correction.""" from ._restraints import boresch as _boresch @@ -157,6 +158,7 @@ def _assemble_restraints( phi0=[f"{pA0:.6f} rad", f"{pB0:.6f} rad", f"{pC0:.6f} rad"], temperature=temperature, angle_potential=angle_potential, + restraint_lever=restraint_lever, ) correction = _get_ssc(restraints[0], temperature) @@ -208,6 +210,19 @@ def _validate_angle_potential(angle_potential): ) +def _validate_restraint_lever(restraint_lever): + """ + Validate 'restraint_lever', shared by both search protocols. Must be + either "combined" or "split" (see sire.restraints.boresch's + restraint_lever parameter). + """ + if restraint_lever not in ("combined", "split"): + raise ValueError( + "'restraint_lever' must be either 'combined' or 'split', " + f"got {restraint_lever!r}" + ) + + # --------------------------------------------------------------------------- # RXRX protocol (default): H-bond-driven restraint search algorithm named # and described in: @@ -233,6 +248,7 @@ def _boresch_search_rxrx( force_constant_angle=None, min_frames=50, angle_potential="restricted_bending", + restraint_lever="split", ): """ Generate a Boresch restraint using the RXRX restraint search algorithm. @@ -294,6 +310,13 @@ def _boresch_search_rxrx( from ever reaching the Boresch collinearity singularity - the exact failure mode this restraint search is designed to avoid. + restraint_lever : str + How the restraint's six degrees of freedom are grouped into + lambda-addressable OpenMM Forces (see sire.restraints.boresch). + Defaults to "split", matching the RXRX protocol's staged restraint + turn-on, where the dihedral terms and the distance/angle terms are + turned on according to different lambda schedule equations. + Returns ------- @@ -360,6 +383,7 @@ def _boresch_search_rxrx( raise ValueError(f"'min_frames' must be a positive int, got {min_frames!r}") _validate_angle_potential(angle_potential) + _validate_restraint_lever(restraint_lever) n_frames = system.num_frames() if n_frames < min_frames: @@ -809,6 +833,7 @@ def _nearest_candidate(atom_idx, connectivity): kp_str, temperature, angle_potential, + restraint_lever, ) @@ -830,6 +855,7 @@ def _boresch_search_aldeghi( max_candidates=100, min_frames=50, angle_potential="restricted_bending", + restraint_lever="combined", ): """ Generate a Boresch restraint for an ABFE simulation by analysing a @@ -888,6 +914,12 @@ def _boresch_search_aldeghi( either "harmonic" or "restricted_bending" (see sire.restraints.boresch). Defaults to "restricted_bending". + restraint_lever : str + How the restraint's six degrees of freedom are grouped into + lambda-addressable OpenMM Forces (see sire.restraints.boresch). + Defaults to "combined", matching the Aldeghi protocol, where all + six restraint terms are turned on together. + Returns ------- @@ -944,6 +976,7 @@ def _boresch_search_aldeghi( raise ValueError(f"'min_frames' must be a positive int, got {min_frames!r}") _validate_angle_potential(angle_potential) + _validate_restraint_lever(restraint_lever) n_frames = system.num_frames() if n_frames < min_frames: @@ -1224,6 +1257,7 @@ def _k(sigma): kp_str, temperature, angle_potential, + restraint_lever, ) @@ -1249,6 +1283,7 @@ def boresch_search( force_constant=None, max_candidates=100, angle_potential="restricted_bending", + restraint_lever=None, ): """ Generate a Boresch restraint for an ABFE simulation by analysing a @@ -1294,6 +1329,15 @@ def boresch_search( from ever reaching the Boresch collinearity singularity. Used by both protocols. + restraint_lever : str, optional + How the restraint's six degrees of freedom are grouped into + lambda-addressable OpenMM Forces (see sire.restraints.boresch). + Defaults to None, which is matched to ``protocol``: "split" for + "rxrx" (matching the RXRX protocol's staged restraint turn-on) and + "combined" for "aldeghi" (matching the Aldeghi protocol, where all + six restraint terms are turned on together). Used by both + protocols. + protein_selection : str [rxrx only] Sire selection string used to identify protein atoms considered as hydrogen-bond partners. Defaults to all non-water @@ -1349,6 +1393,13 @@ def boresch_search( if not isinstance(protocol, str): raise TypeError(f"'protocol' must be a str, got {type(protocol)}") + if restraint_lever is None: + # Matches the paper: the RXRX protocol turns the dihedral and + # distance/angle restraint terms on according to different lambda + # schedule equations, whereas the Aldeghi protocol turns all six + # terms on together. + restraint_lever = "split" if protocol == "rxrx" else "combined" + if protocol == "rxrx": return _boresch_search_rxrx( system, @@ -1362,6 +1413,7 @@ def boresch_search( force_constant_angle=force_constant_angle, min_frames=min_frames, angle_potential=angle_potential, + restraint_lever=restraint_lever, ) elif protocol == "aldeghi": return _boresch_search_aldeghi( @@ -1374,6 +1426,7 @@ def boresch_search( max_candidates=max_candidates, min_frames=min_frames, angle_potential=angle_potential, + restraint_lever=restraint_lever, ) else: raise ValueError( diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 3c3db78eb..246eadece 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -144,6 +144,7 @@ def boresch( map=None, temperature=u("298 K"), angle_potential=None, + restraint_lever=None, ): """ Create a set of Boresch restraints that will restrain the 6 @@ -243,6 +244,20 @@ def boresch( harmonic potential away from theta0. Default is None, which is equivalent to "harmonic". + restraint_lever : str, optional + How the restraint's six degrees of freedom are grouped into + lambda-addressable OpenMM Forces, either "combined" or "split". + With "combined" (the default), all six terms (distance, two + angles, three dihedrals) share a single scale factor and are + therefore always turned on/off together according to a single + lambda schedule equation. With "split", the distance and two + angle terms share one scale factor, and the three dihedral terms + share a second, independent scale factor, allowing the two + groups to be turned on according to different lambda schedule + equations (e.g. to reproduce the RXRX protocol's staged + restraint turn-on). Default is None, which is equivalent to + "combined". + Returns ------- BoreschRestraints : SireMM::BoreschRestraints @@ -293,11 +308,18 @@ def boresch( if angle_potential is not None else map_dict.get("angle_potential", None) ) + restraint_lever = ( + restraint_lever + if restraint_lever is not None + else map_dict.get("restraint_lever", None) + ) # Values retrieved from the map are wrapped as PropertyName, not a plain - # str, which the strict BoreschRestraints.set_angle_potential(QString) - # signature doesn't accept directly. + # str, which the strict BoreschRestraints.set_angle_potential(QString)/ + # set_restraint_lever(QString) signatures don't accept directly. if angle_potential is not None: angle_potential = str(angle_potential) + if restraint_lever is not None: + restraint_lever = str(restraint_lever) receptor = _to_atoms(mols, receptor) ligand = _to_atoms(mols, ligand) @@ -326,6 +348,15 @@ def boresch( else: angle_potential = "harmonic" + if restraint_lever is not None: + if restraint_lever not in ("combined", "split"): + raise ValueError( + "'restraint_lever' must be either 'combined' or " + f"'split', got {restraint_lever!r}" + ) + else: + restraint_lever = "combined" + from .. import measure default_distance_k = u("5 kcal mol-1 A-2") @@ -482,6 +513,10 @@ def boresch( # Set the functional form used for the two angle restraint terms. b.set_angle_potential(angle_potential) + # Set how the restraint's degrees of freedom are grouped into + # lambda-addressable OpenMM Forces. + b.set_restraint_lever(restraint_lever) + return b diff --git a/tests/restraints/test_boresch.py b/tests/restraints/test_boresch.py index 8365469ba..f676d3f6a 100644 --- a/tests/restraints/test_boresch.py +++ b/tests/restraints/test_boresch.py @@ -361,3 +361,47 @@ def test_boresch_angle_potential_via_map(thrombin_complex): thrombin_complex, map={"angle_potential": "restricted_bending"} ) assert boresch_restraints.angle_potential() == "restricted_bending" + + +def test_boresch_restraint_lever_defaults_to_combined(thrombin_complex): + """ + Check that restraint_lever defaults to "combined" when not specified, + matching the pre-existing behaviour (no change for existing callers). + """ + boresch_restraints = _make_default_boresch(thrombin_complex) + assert boresch_restraints.restraint_lever() == "combined" + + +def test_boresch_restraint_lever_split(thrombin_complex): + """ + Check that restraint_lever="split" is set correctly, and that it + doesn't change any of the other restraint parameters. + """ + boresch_restraints = _make_default_boresch( + thrombin_complex, restraint_lever="split" + ) + assert boresch_restraints.restraint_lever() == "split" + + boresch_restraint = boresch_restraints[0] + assert boresch_restraint.kr().value() == 6.2012 + assert boresch_restraint.theta0()[0].value() == 1.3031 + assert boresch_restraint.theta0()[1].value() == 1.4777 + + +def test_boresch_restraint_lever_invalid_raises(thrombin_complex): + """ + Check that an invalid restraint_lever value raises a ValueError. + """ + with pytest.raises(ValueError, match="restraint_lever"): + _make_default_boresch(thrombin_complex, restraint_lever="not_a_real_lever") + + +def test_boresch_restraint_lever_via_map(thrombin_complex): + """ + Check that restraint_lever can also be set via the map, matching every + other boresch() option. + """ + boresch_restraints = _make_default_boresch( + thrombin_complex, map={"restraint_lever": "split"} + ) + assert boresch_restraints.restraint_lever() == "split" diff --git a/tests/restraints/test_boresch_search.py b/tests/restraints/test_boresch_search.py index 4ee06d9e1..69debde78 100644 --- a/tests/restraints/test_boresch_search.py +++ b/tests/restraints/test_boresch_search.py @@ -57,6 +57,18 @@ def test_angle_potential_defaults_to_restricted_bending(self, boresch_result): restraints, _ = boresch_result assert restraints.angle_potential() == "restricted_bending" + def test_restraint_lever_matches_protocol(self, boresch_result, request): + """ + boresch_search() defaults restraint_lever to match the protocol: + "split" for rxrx (matching the RXRX protocol's staged restraint + turn-on), "combined" for aldeghi (matching the Aldeghi protocol, + where all six restraint terms are turned on together). + """ + restraints, _ = boresch_result + protocol = request.node.callspec.params["boresch_result"] + expected = "split" if protocol == "rxrx" else "combined" + assert restraints.restraint_lever() == expected + def test_correction_is_negative(self, boresch_result): """Standard state correction is always negative (costs free energy to restrain).""" _, correction = boresch_result @@ -137,6 +149,25 @@ def test_angle_potential_invalid_raises(self, abfe_system, protocol): with pytest.raises(ValueError, match="angle_potential"): boresch_search(abfe_system, protocol=protocol, angle_potential="nonsense") + @pytest.mark.parametrize("protocol", PROTOCOLS) + def test_restraint_lever_override(self, abfe_system, protocol): + from sire.restraints import boresch_search + + # Explicitly pick whichever value isn't the protocol's own default, + # to confirm the override actually takes effect. + override = "combined" if protocol == "rxrx" else "split" + restraints, _ = boresch_search( + abfe_system, protocol=protocol, restraint_lever=override + ) + assert restraints.restraint_lever() == override + + @pytest.mark.parametrize("protocol", PROTOCOLS) + def test_restraint_lever_invalid_raises(self, abfe_system, protocol): + from sire.restraints import boresch_search + + with pytest.raises(ValueError, match="restraint_lever"): + boresch_search(abfe_system, protocol=protocol, restraint_lever="nonsense") + def test_force_constant_override(self, abfe_system): """'force_constant' is an Aldeghi-only kwarg (RXRX uses two separate force constants, matching the fixed protocol values from the paper).""" diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 387bfe5de..1fdb73fb5 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -113,78 +113,119 @@ void _add_boresch_restraints(const SireMM::BoreschRestraints &restraints, // the harmonic form with the same ktheta, so no change is needed to // the analytic standard state correction formula. // - QString energy_expression; + // restraints.restraintLever() selects how these terms are grouped into + // lambda-addressable OpenMM Forces: + // + // "combined" (default): all six terms (distance, two angles, three + // dihedrals) share a single scale factor 'rho', in a single Force, + // registered under 'restraints.name()'. + // + // "split": the distance and two angle terms share one scale factor, + // in a Force registered under 'restraints.name() + "_distance_angle"'; + // the three dihedral terms share a second, independent scale factor, + // in a Force registered under 'restraints.name() + "_dihedral"'. This + // allows the two groups to be turned on according to different lambda + // schedule equations (e.g. to reproduce the RXRX protocol's staged + // restraint turn-on). + // + const QString angle_terms = + restraints.anglePotential() == "restricted_bending" + ? QString("e_angle_A=ktheta_A*(cos(theta_A)-cos(theta0_A))^2/sin(theta_A)^2;" + "e_angle_B=ktheta_B*(cos(theta_B)-cos(theta0_B))^2/sin(theta_B)^2;") + : QString("e_angle_B=ktheta_B*(theta_B-theta0_B)^2;" + "e_angle_A=ktheta_A*(theta_A-theta0_A)^2;"); + + const QString distance_angle_expression = + QString("rho*(e_bond + e_angle_A + e_angle_B);" + "e_bond=kr*(r-r0)^2;%1" + "theta_B=angle(p1, p4, p5);" + "theta_A=angle(p2, p1, p4);" + "r=distance(p1, p4);") + .arg(angle_terms); + + const QString dihedral_expression = + QString("rho*(e_torsion_A + e_torsion_B + e_torsion_C);" + "e_torsion_C=kphi_C*(min(dphi_C, two_pi-dphi_C))^2;" + "e_torsion_B=kphi_B*(min(dphi_B, two_pi-dphi_B))^2;" + "e_torsion_A=kphi_A*(min(dphi_A, two_pi-dphi_A))^2;" + "dphi_C=abs(phi_C-phi0_C);" + "dphi_B=abs(phi_B-phi0_B);" + "dphi_A=abs(phi_A-phi0_A);" + "two_pi=6.283185307179586;" + "phi_C=dihedral(p1, p4, p5, p6);" + "phi_B=dihedral(p2, p1, p4, p5);" + "phi_A=dihedral(p3, p2, p1, p4);"); + + const QString combined_expression = + QString("rho*(e_bond + e_angle_A + e_angle_B + e_torsion_A + e_torsion_B + e_torsion_C);" + "e_bond=kr*(r-r0)^2;%1" + "e_torsion_C=kphi_C*(min(dphi_C, two_pi-dphi_C))^2;" + "e_torsion_B=kphi_B*(min(dphi_B, two_pi-dphi_B))^2;" + "e_torsion_A=kphi_A*(min(dphi_A, two_pi-dphi_A))^2;" + "dphi_C=abs(phi_C-phi0_C);" + "dphi_B=abs(phi_B-phi0_B);" + "dphi_A=abs(phi_A-phi0_A);" + "two_pi=6.283185307179586;" + "phi_C=dihedral(p1, p4, p5, p6);" + "phi_B=dihedral(p2, p1, p4, p5);" + "phi_A=dihedral(p3, p2, p1, p4);" + "theta_B=angle(p1, p4, p5);" + "theta_A=angle(p2, p1, p4);" + "r=distance(p1, p4);") + .arg(angle_terms); + + const bool split_lever = restraints.restraintLever() == "split"; - if (restraints.anglePotential() == "restricted_bending") - { - energy_expression = QString( - "rho*(e_bond + e_angle_A + e_angle_B + e_torsion_A + e_torsion_B + e_torsion_C);" - "e_bond=kr*(r-r0)^2;" - "e_angle_A=ktheta_A*(cos(theta_A)-cos(theta0_A))^2/sin(theta_A)^2;" - "e_angle_B=ktheta_B*(cos(theta_B)-cos(theta0_B))^2/sin(theta_B)^2;" - "e_torsion_C=kphi_C*(min(dphi_C, two_pi-dphi_C))^2;" - "e_torsion_B=kphi_B*(min(dphi_B, two_pi-dphi_B))^2;" - "e_torsion_A=kphi_A*(min(dphi_A, two_pi-dphi_A))^2;" - "dphi_C=abs(phi_C-phi0_C);" - "dphi_B=abs(phi_B-phi0_B);" - "dphi_A=abs(phi_A-phi0_A);" - "two_pi=6.283185307179586;" - "phi_C=dihedral(p1, p4, p5, p6);" - "phi_B=dihedral(p2, p1, p4, p5);" - "phi_A=dihedral(p3, p2, p1, p4);" - "theta_B=angle(p1, p4, p5);" - "theta_A=angle(p2, p1, p4);" - "r=distance(p1, p4);"); - } - else + 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); + const double internal_to_ktheta = (1 * SireUnits::kcal_per_mol / (SireUnits::radian2)).to(SireUnits::kJ_per_mol / SireUnits::radian2); + + // Create and register a CustomCompoundBondForce for 'restraints', with + // the passed name, energy expression and per-bond parameter names. + auto make_force = [&](const QString &name, const QString &expression, + const QStringList ¶m_names) { - energy_expression = QString( - "rho*(e_bond + e_angle_A + e_angle_B + e_torsion_A + e_torsion_B + e_torsion_C);" - "e_bond=kr*(r-r0)^2;" - "e_angle_B=ktheta_B*(theta_B-theta0_B)^2;" - "e_angle_A=ktheta_A*(theta_A-theta0_A)^2;" - "e_torsion_C=kphi_C*(min(dphi_C, two_pi-dphi_C))^2;" - "e_torsion_B=kphi_B*(min(dphi_B, two_pi-dphi_B))^2;" - "e_torsion_A=kphi_A*(min(dphi_A, two_pi-dphi_A))^2;" - "dphi_C=abs(phi_C-phi0_C);" - "dphi_B=abs(phi_B-phi0_B);" - "dphi_A=abs(phi_A-phi0_A);" - "two_pi=6.283185307179586;" - "phi_C=dihedral(p1, p4, p5, p6);" - "phi_B=dihedral(p2, p1, p4, p5);" - "phi_A=dihedral(p3, p2, p1, p4);" - "theta_B=angle(p1, p4, p5);" - "theta_A=angle(p2, p1, p4);" - "r=distance(p1, p4);"); - } + auto *ff = new OpenMM::CustomCompoundBondForce(6, expression.toStdString()); + ff->setName("BoreschRestraintForce"); - auto *restraintff = new OpenMM::CustomCompoundBondForce(6, energy_expression.toStdString()); - restraintff->setName("BoreschRestraintForce"); + for (const auto ¶m_name : param_names) + { + ff->addPerBondParameter(param_name.toStdString()); + } - restraintff->addPerBondParameter("rho"); - restraintff->addPerBondParameter("kr"); - restraintff->addPerBondParameter("r0"); - restraintff->addPerBondParameter("ktheta_A"); - restraintff->addPerBondParameter("theta0_A"); - restraintff->addPerBondParameter("ktheta_B"); - restraintff->addPerBondParameter("theta0_B"); - restraintff->addPerBondParameter("kphi_A"); - restraintff->addPerBondParameter("phi0_A"); - restraintff->addPerBondParameter("kphi_B"); - restraintff->addPerBondParameter("phi0_B"); - restraintff->addPerBondParameter("kphi_C"); - restraintff->addPerBondParameter("phi0_C"); + ff->setUsesPeriodicBoundaryConditions(restraints.usesPbc()); - restraintff->setUsesPeriodicBoundaryConditions(restraints.usesPbc()); + ff->setForceGroup(force_group_counter); + lambda_lever.addRestraintIndex(name, system.addForce(ff)); + lambda_lever.setRestraintForceGroup(name, force_group_counter++); - restraintff->setForceGroup(force_group_counter); - lambda_lever.addRestraintIndex(restraints.name(), - system.addForce(restraintff)); - lambda_lever.setRestraintForceGroup(restraints.name(), force_group_counter++); + return ff; + }; - 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); - const double internal_to_ktheta = (1 * SireUnits::kcal_per_mol / (SireUnits::radian2)).to(SireUnits::kJ_per_mol / SireUnits::radian2); + OpenMM::CustomCompoundBondForce *restraintff = nullptr; + OpenMM::CustomCompoundBondForce *distance_angle_ff = nullptr; + OpenMM::CustomCompoundBondForce *dihedral_ff = nullptr; + + if (split_lever) + { + distance_angle_ff = make_force(restraints.name() + "_distance_angle", + distance_angle_expression, + {"rho", "kr", "r0", "ktheta_A", "theta0_A", + "ktheta_B", "theta0_B"}); + + dihedral_ff = make_force(restraints.name() + "_dihedral", + dihedral_expression, + {"rho", "kphi_A", "phi0_A", "kphi_B", "phi0_B", + "kphi_C", "phi0_C"}); + } + else + { + restraintff = make_force( + restraints.name(), + combined_expression, + {"rho", "kr", "r0", "ktheta_A", "theta0_A", "ktheta_B", "theta0_B", + "kphi_A", "phi0_A", "kphi_B", "phi0_B", "kphi_C", "phi0_C"}); + } for (const auto &restraint : restraints.restraints()) { @@ -194,9 +235,6 @@ void _add_boresch_restraints(const SireMM::BoreschRestraints &restraints, std::vector particles; particles.resize(6); - std::vector parameters; - parameters.resize(13); - for (int i = 0; i < 3; ++i) { particles[i] = real_atoms[restraint.receptorAtoms()[i]]; @@ -215,21 +253,33 @@ void _add_boresch_restraints(const SireMM::BoreschRestraints &restraints, } } - parameters[0] = 1.0; // rho - parameters[1] = restraint.kr().value() * internal_to_k; // kr - parameters[2] = restraint.r0().value() * internal_to_nm; // r0 - parameters[3] = restraint.ktheta()[0].value() * internal_to_ktheta; // ktheta_A - parameters[4] = restraint.theta0()[0].value(); // theta0_A (already in radians) - parameters[5] = restraint.ktheta()[1].value() * internal_to_ktheta; // ktheta_B - parameters[6] = restraint.theta0()[1].value(); // theta0_B - parameters[7] = restraint.kphi()[0].value() * internal_to_ktheta; // kphi_A - parameters[8] = restraint.phi0()[0].value(); // phi0_A - parameters[9] = restraint.kphi()[1].value() * internal_to_ktheta; // kphi_B - parameters[10] = restraint.phi0()[1].value(); // phi0_B - parameters[11] = restraint.kphi()[2].value() * internal_to_ktheta; // kphi_C - parameters[12] = restraint.phi0()[2].value(); // phi0_C - - restraintff->addBond(particles, parameters); + const double kr = restraint.kr().value() * internal_to_k; + const double r0 = restraint.r0().value() * internal_to_nm; + const double ktheta_A = restraint.ktheta()[0].value() * internal_to_ktheta; + const double theta0_A = restraint.theta0()[0].value(); // already in radians + const double ktheta_B = restraint.ktheta()[1].value() * internal_to_ktheta; + const double theta0_B = restraint.theta0()[1].value(); + const double kphi_A = restraint.kphi()[0].value() * internal_to_ktheta; + const double phi0_A = restraint.phi0()[0].value(); + const double kphi_B = restraint.kphi()[1].value() * internal_to_ktheta; + const double phi0_B = restraint.phi0()[1].value(); + const double kphi_C = restraint.kphi()[2].value() * internal_to_ktheta; + const double phi0_C = restraint.phi0()[2].value(); + + if (split_lever) + { + distance_angle_ff->addBond( + particles, {1.0, kr, r0, ktheta_A, theta0_A, ktheta_B, theta0_B}); + + dihedral_ff->addBond( + particles, {1.0, kphi_A, phi0_A, kphi_B, phi0_B, kphi_C, phi0_C}); + } + else + { + restraintff->addBond(particles, {1.0, kr, r0, ktheta_A, theta0_A, ktheta_B, + theta0_B, kphi_A, phi0_A, kphi_B, phi0_B, + kphi_C, phi0_C}); + } } } diff --git a/wrapper/MM/BoreschRestraints.pypp.cpp b/wrapper/MM/BoreschRestraints.pypp.cpp index 9f57733b8..585c164fc 100644 --- a/wrapper/MM/BoreschRestraints.pypp.cpp +++ b/wrapper/MM/BoreschRestraints.pypp.cpp @@ -97,6 +97,14 @@ void register_BoreschRestraints_class() BoreschRestraints_exposer.def( "anglePotential", anglePotential_function_value, bp::release_gil_policy(), "Return the functional form used for the two Boresch angle restraint terms,\neither \"harmonic\" or \"restricted_bending\"."); } + { //::SireMM::BoreschRestraints::restraintLever + + typedef ::QString (::SireMM::BoreschRestraints::*restraintLever_function_type)() const; + restraintLever_function_type restraintLever_function_value(&::SireMM::BoreschRestraints::restraintLever); + + BoreschRestraints_exposer.def( + "restraintLever", restraintLever_function_value, bp::release_gil_policy(), "Return how the restraint's six degrees of freedom are grouped into\nlambda-addressable OpenMM Forces, either \"combined\" or \"split\"."); + } { //::SireMM::BoreschRestraints::isEmpty typedef bool (::SireMM::BoreschRestraints::*isEmpty_function_type)() const; @@ -165,6 +173,14 @@ void register_BoreschRestraints_class() BoreschRestraints_exposer.def( "setAnglePotential", setAnglePotential_function_value, (bp::arg("angle_potential")), bp::release_gil_policy(), "Set the functional form used for the two Boresch angle restraint terms.\nMust be either \"harmonic\" or \"restricted_bending\"."); } + { //::SireMM::BoreschRestraints::setRestraintLever + + typedef void (::SireMM::BoreschRestraints::*setRestraintLever_function_type)(::QString const &); + setRestraintLever_function_type setRestraintLever_function_value(&::SireMM::BoreschRestraints::setRestraintLever); + + BoreschRestraints_exposer.def( + "setRestraintLever", setRestraintLever_function_value, (bp::arg("restraint_lever")), bp::release_gil_policy(), "Set how the restraint's six degrees of freedom are grouped into\nlambda-addressable OpenMM Forces. Must be either \"combined\" (the\ndefault) or \"split\"."); + } { //::SireMM::BoreschRestraints::size typedef int (::SireMM::BoreschRestraints::*size_function_type)() const;