From 77fc745264fa4577e8d07314ef03dfa0be5c2a17 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 21 Jan 2024 13:52:57 +0000 Subject: [PATCH 01/13] Improve flexibility of boresch function This allows the user to pass in equilibrium values as well as force constants. In addition, the likely stability of the restraints is assessed. --- src/sire/restraints/_restraints.py | 363 ++++++++++++++++++++++------- 1 file changed, 283 insertions(+), 80 deletions(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 8689e8b7f..5a2dafd7f 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -23,34 +23,115 @@ def boresch( kr=None, ktheta=None, kphi=None, + r0=None, + theta0=None, + phi0=None, name: str = None, map=None, + temperature=sr.u("298 K"), ): """ - Create a set of Boresch restraints that will hold the passed - ligand in a its relative binding mode relative to the - passed receptor. All of the atoms in both 'ligand' and - 'receptor' must be contained in 'mols'. + Create a set of Boresch restraints that will restrain the 6 + external degrees of freedom of the ligand relative to the receptor. + All of the atoms in both 'ligand' and 'receptor' must be contained in + 'mols'. The BoreschRestraint will be a set of six restraints between three identified ligand atoms, and three identified receptor - atoms. + atoms: - 1. A single distance restraint, with specified kr and r0 parameters - 2. Two angle restraints, with the specified two ktheta and theta0 - parameters - 3. Three torsion restraints, with the specified three kphi and phi0 - parameters + 1. A single distance restraint, with specified force constant (kr) + and equilibrium distance (r0) parameters. + 2. Two angle restraints, with specified force constants (ktheta) + and equilibrium angles (theta0) parameters. + 3. Three torsion restraints, with specified force constants (kphi) + and equilibrium angles (phi0) parameters. This will create a single BoreschRestraint, which will be passed back in a BoreschRestraints object. - If the force constants (kr, ktheta and kphi) are None, then they - will have default values of 150 kcal mol-1 A-2 and - 150 kcal mol-1 rad-2 - - The equilibium distances and angles are based on the current coordinates - of that atoms + Parameters + ---------- + mols : sire.system._system.System + The system containing the receptor and ligand. + + receptor : SireMol::Selector + The receptor atoms to restrain. + + ligand : SireMol::Selector + The ligand atoms to restrain. + + kr : str or SireUnits::Dimension::GeneralUnit, optional + The force constant for the distance restraint. If None, this will + default to 10 kcal mol-1 A-2. Default is None. + + ktheta : str or SireUnits::Dimension::GeneralUnit or list of str or SireUnits::Dimension::GeneralUnit, optional + The force constants for the angle restraints, in the order kthetaA, + kthetaB If None, this will default to 100 kcal mol-1 rad-2 for + both angle restraints. If a list, then this should be a list of + length 2 containing the force constants for the two angle + restraints. If a single value, then this will be used for both + angle restraints. Default is None. + + kphi : str or SireUnits::Dimension::GeneralUnit or list of str or SireUnits::Dimension::GeneralUnit, optional + The force constants for the torsion restraints, in the order kthetaA, + kthetaB, kthetaC. If None, this will default to 100 kcal mol-1 rad-2 + for all three torsion restraints. If a list, then this should be a + list of length 3 containing the force constants for the three + torsion restraints. If a single value, then this will be used for + all three torsion restraints. Default is None. + + r0 : str or SireUnits::Dimension::GeneralUnit, optional + The equilibrium distance for the distance restraint. If None, this + will be measured from the current coordinates of the atoms. Default + is None. + + theta0 : list of str or SireUnits::Dimension::GeneralUnit, optional + The equilibrium angles for the angle restraints. If None, these + will be measured from the current coordinates of the atoms. If a + list, then this should be a list of length 2 containing the + equilibrium angles for the two angle restraints. Default is None. + + phi0 : list of str or SireUnits::Dimension::GeneralUnit, optional + The equilibrium angles for the torsion restraints. If None, these + will be measured from the current coordinates of the atoms. If a + list, then this should be a list of length 3 containing the + equilibrium angles for the three torsion restraints. Default is None. + + name : str, optional + The name of the restraint. If None, then a default name will be + used. Default is None. + + map : dict, optional + A property map. Default is None. + + temperature : str or SireUnits::Dimension::GeneralUnit, optional + The temperature to use when checking for unstable restraints. If + None, then this will default to 298 K. Default is None. + + Returns + ------- + BoreschRestraints : SireMM::BoreschRestraints + The Boresch restraints. + + Examples + -------- + Create a set of Boresch restraints for the ligand in the system + 'system', specifying all of the force constants and equilibrium + values: + ``` + my_boresch_restraint = boresch( + system, + receptor=system["protein"]["atomidx 1538 or atomidx 1518 or atomidx 1540"], + ligand=system["resname LIG"]["atomidx 4 or atomidx 3 or atomidx 5"], + kr="6.2012 kcal mol-1 A-2", + ktheta=["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], + kphi=["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2", "55.1775 kcal mol-1 rad-2"], + r0="7.687 A", + theta0=["1.3031 rad", "1.4777 rad"], + phi0=["2.5569 rad", "2.9359 rad", "1.4147 rad"], + ) + ``` """ from .. import u from ..base import create_map @@ -64,73 +145,148 @@ def boresch( if len(receptor) != 3 or len(ligand) != 3: # Eventually will choose the best atoms from the receptor # and ligand... - raise ValueError( - "You need to provide 3 receptor atoms and 3 ligand atoms" - ) + raise ValueError("You need to provide 3 receptor atoms and 3 ligand atoms") - if kr is None: - kr = u("150 kcal mol-1 A-2") - else: - kr = u(kr) + from .. import measure - _default_k = u("150 kcal mol-1 rad-2") + default_distance_k = u("10 kcal mol-1 A-2") + default_angle_k = u("100 kcal mol-1 rad-2") - if ktheta is None: - ktheta = [_default_k, _default_k] - elif type(ktheta) is not list: - ktheta = 2 * [u(ktheta)] - else: - if len(ktheta) == 0: - ktheta = [_default_k, _default_k] - if len(ktheta) < 2: - ktheta = 2 * [u(ktheta[0])] - else: - ktheta = [u(x) for x in ktheta[0:2]] + # Use the user-specified equilibrium values if they are provided. + distance = [[ligand[0], receptor[0]]] + angles = [ + [receptor[1], receptor[0], ligand[0]], + [receptor[0], ligand[0], ligand[1]], + ] + dihedrals = [ + [receptor[2], receptor[1], receptor[0], ligand[0]], + [receptor[1], receptor[0], ligand[0], ligand[1]], + [receptor[0], ligand[0], ligand[1], ligand[2]], + ] - if kphi is None: - kphi = [_default_k, _default_k, _default_k] - elif type(kphi) is not list: - kphi = 3 * [u(kphi)] - else: - if len(kphi) == 0: - kphi = [_default_k, _default_k, _default_k] - elif len(kphi) == 1: - kphi = 3 * [u(kphi[0])] - elif len(kphi) == 2: - kphi = [u(kphi[0]), u(kphi[1]), u(kphi[1])] + restraint_components = { + "distance": { + "input_k": kr, + "validated_k": [], + "input_equil": r0, + "measure": distance, + "validated_equil": [], + }, + "angle": { + "input_k": ktheta, + "validated_k": [], + "input_equil": theta0, + "measure": angles, + "validated_equil": [], + }, + "dihedral": { + "input_k": kphi, + "validated_k": [], + "input_equil": phi0, + "measure": dihedrals, + "validated_equil": [], + }, + } + + for restraint_component in restraint_components: + n_measures = len(restraint_components[restraint_component]["measure"]) + + # Get the force constants. + if restraint_components[restraint_component]["input_k"] is None: + restraint_components[restraint_component]["validated_k"] = n_measures * [ + default_distance_k + if restraint_component == "distance" + else default_angle_k + ] + elif type(restraint_components[restraint_component]["input_k"]) is not list: + # Populate the list with the single specified value. + restraint_components[restraint_component]["validated_k"] = n_measures * [ + u(restraint_components[restraint_component]["input_k"]) + ] else: - kphi = [u(x) for x in kphi[0:3]] - - # r is | Ligand1 - Receptor1 | = distance(P1, P4) - # thetaA = angle(R2, R1, L1) = angle(P2, P1, P4) - # thetaB = angle(R1, L1, L2) = angle(P1, P4, P5) - # phiA = dihedral(R3, R2, R1, L1) = dihedral(P3, P2, P1, P4) - # phiB = dihedral(R2, R1, L1, L2) = dihedral(P2, P1, P4, P5) - # phiC = dihedral(R1, L1, L2, L3) = dihedral(P1, P4, P5, P6) - from .. import measure + if len(restraint_components[restraint_component]["input_k"]) == 0: + # Empty list - populate with default values. + restraint_components[restraint_component][ + "validated_k" + ] = n_measures * [ + default_distance_k + if restraint_component == "distance" + else default_angle_k + ] + elif len(restraint_components[restraint_component]["input_k"]) == 1: + # List of length 1 - populate with that value. + restraint_components[restraint_component][ + "validated_k" + ] = n_measures * [ + u(restraint_components[restraint_component]["input_k"][0]) + ] + elif ( + len(restraint_components[restraint_component]["input_k"]) == n_measures + ): + # List of the correct length for this restraint component. + restraint_components[restraint_component]["validated_k"] = [ + u(x) for x in restraint_components[restraint_component]["input_k"] + ] + else: + # List of the wrong length. + raise ValueError( + f"Input force constants for {restraint_component} must be a single value or a list of length {n_measures}" + ) - r0 = measure(ligand[0], receptor[0]) - theta0 = [ - measure(receptor[1], receptor[0], ligand[0]), - measure(receptor[0], ligand[0], ligand[1]), - ] - phi0 = [ - measure(receptor[2], receptor[1], receptor[0], ligand[0]), - measure(receptor[1], receptor[0], ligand[0], ligand[1]), - measure(receptor[0], ligand[0], ligand[1], ligand[2]), - ] + # Get the equilibrium values. + if restraint_components[restraint_component]["input_equil"] is None: + # Set all components to None - these will be measured from the structure later. + restraint_components[restraint_component]["input_equil"] = [ + None for i in range(n_measures) + ] + + if type(restraint_components[restraint_component]["input_equil"]) is not list: + # Only allow this if we are dealing with the distance component, as this is the only one that can be a single value. + if restraint_component == "distance": + restraint_components[restraint_component][ + "input_equil" + ] = n_measures * [ + u(restraint_components[restraint_component]["input_equil"]) + ] + else: + raise ValueError( + f"Input equilibrium values for {restraint_component} must be a list of length {n_measures} of values or Nones" + ) + + elif ( + len(restraint_components[restraint_component]["input_equil"]) != n_measures + ): + raise ValueError( + f"If specified, equilibrium values for {restraint_component} must be a list of length {n_measures} of values or Nones" + ) + + # Now validate the input equilibrium values, replacing Nones with measured values. + for i, equil in enumerate( + restraint_components[restraint_component]["input_equil"] + ): + if equil is not None: + restraint_components[restraint_component]["validated_equil"].append( + u(equil) + ) + else: + restraint_components[restraint_component]["validated_equil"].append( + measure(*restraint_components[restraint_component]["measure"][i]) + ) + + # Warn the user if the restraint is likely to be unstable. + _check_stability_boresch_restraint(restraint_components, temperature) mols = mols.atoms() b = BoreschRestraint( receptor=mols.find(receptor), ligand=mols.find(ligand), - r0=r0, - theta0=theta0, - phi0=phi0, - kr=kr, - ktheta=ktheta, - kphi=kphi, + r0=restraint_components["distance"]["validated_equil"][0], + theta0=restraint_components["angle"]["validated_equil"], + phi0=restraint_components["dihedral"]["validated_equil"], + kr=restraint_components["distance"]["validated_k"][0], + ktheta=restraint_components["angle"]["validated_k"], + kphi=restraint_components["dihedral"]["validated_k"], ) if name is None: @@ -139,9 +295,60 @@ def boresch( return BoreschRestraint(name, b) -def distance( - mols, atoms0, atoms1, r0=None, k=None, name: str = None, map=None -): +def _check_stability_boresch_restraint(restraint_components, temperature=_u("298 K")): + """ + Internal function to check for likely unstable Boresch restraints. + """ + import warnings as _warnings + + from .. import u + + # Check for unstable combinations of force constants. + if restraint_components["distance"]["validated_k"][0].value() == 0: + raise ValueError('"kr" cannot be zero') + + if restraint_components["angle"]["validated_k"][0].value() == 0: + if ( + restraint_components["dihedral"]["validated_k"][0].value() != 0 + or restraint_components["dihedral"]["validated_k"][1].value() != 0 + ): + raise ValueError( + "Restraining phiA or phiB without restraining thetaA " + "will produce unstable Boresch restraints." + ) + + if restraint_components["angle"]["validated_k"][1].value() == 0: + if ( + restraint_components["dihedral"]["validated_k"][1].value() != 0 + or restraint_components["dihedral"]["validated_k"][2].value() != 0 + ): + raise ValueError( + "Restraining phiB or phiC without restraining thetaB " + "will produce unstable Boresch restraints." + ) + + # Ensure that restrained angles are >= 10 kT from collinear. + for equil_angle, k_angle in zip( + restraint_components["angle"]["validated_equil"], + restraint_components["angle"]["validated_k"], + ): + if k_angle.value() != 0: + # Find the minimum stable angle "distance". We use the squared values as sire units don't support + # taking the square root. + min_stable_dist_sq = (2 * 10 * sr.units.k_boltz * temperature) / k_angle + min_dist_sq = min( + [abs(equil_angle - 0), abs(equil_angle - 180 * sr.units.degrees)] + ).pow(2) + if min_dist_sq < min_stable_dist_sq: + _warnings.warn( + f"The equilibrium value angle value of {equil_angle} is within 10 kT of" + "collinearity, which may result in unstable Boresch restraints." + " Consider increasing the force constants or selecting equilibrium" + " values further from 0 or pi radians." + ) + + +def distance(mols, atoms0, atoms1, r0=None, k=None, name: str = None, map=None): """ Create a set of distance restraints from all of the atoms in 'atoms0' to all of the atoms in 'atoms1' where all atoms are @@ -251,9 +458,7 @@ def bond(*args, **kwargs): return distance(*args, **kwargs) -def positional( - mols, atoms, k=None, r0=None, position=None, name: str = None, map=None -): +def positional(mols, atoms, k=None, r0=None, position=None, name: str = None, map=None): """ Create a set of position restraints for the atoms specified in 'atoms' that are contained in the container 'mols', using the @@ -335,9 +540,7 @@ def positional( if position is None: restraints.add( - PositionalRestraint( - idxs[0], atom.property(coords_prop), ik, ir0 - ) + PositionalRestraint(idxs[0], atom.property(coords_prop), ik, ir0) ) else: restraints.add(PositionalRestraint(idxs[0], position[i], ik, ir0)) From ccb4107aa21807d5c87eb2a3acd6b001f44480d7 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 21 Jan 2024 13:59:26 +0000 Subject: [PATCH 02/13] Clear up import and type errors in _restraints.py --- src/sire/restraints/_restraints.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 5a2dafd7f..d32cc301d 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -5,6 +5,8 @@ "positional", ] +from .. import u + def _to_atoms(mols, atoms): """ @@ -26,9 +28,9 @@ def boresch( r0=None, theta0=None, phi0=None, - name: str = None, + name=None, map=None, - temperature=sr.u("298 K"), + temperature=u("298 K"), ): """ Create a set of Boresch restraints that will restrain the 6 @@ -133,7 +135,6 @@ def boresch( ) ``` """ - from .. import u from ..base import create_map from ..mm import BoreschRestraint, BoreschRestraints @@ -295,13 +296,13 @@ def boresch( return BoreschRestraint(name, b) -def _check_stability_boresch_restraint(restraint_components, temperature=_u("298 K")): +def _check_stability_boresch_restraint(restraint_components, temperature=u("298 K")): """ Internal function to check for likely unstable Boresch restraints. """ import warnings as _warnings - from .. import u + from .. import units # Check for unstable combinations of force constants. if restraint_components["distance"]["validated_k"][0].value() == 0: @@ -335,9 +336,9 @@ def _check_stability_boresch_restraint(restraint_components, temperature=_u("298 if k_angle.value() != 0: # Find the minimum stable angle "distance". We use the squared values as sire units don't support # taking the square root. - min_stable_dist_sq = (2 * 10 * sr.units.k_boltz * temperature) / k_angle + min_stable_dist_sq = (2 * 10 * units.k_boltz * temperature) / k_angle min_dist_sq = min( - [abs(equil_angle - 0), abs(equil_angle - 180 * sr.units.degrees)] + [abs(equil_angle - 0), abs(equil_angle - 180 * units.degrees)] ).pow(2) if min_dist_sq < min_stable_dist_sq: _warnings.warn( @@ -348,7 +349,7 @@ def _check_stability_boresch_restraint(restraint_components, temperature=_u("298 ) -def distance(mols, atoms0, atoms1, r0=None, k=None, name: str = None, map=None): +def distance(mols, atoms0, atoms1, r0=None, k=None, name=None, map=None): """ Create a set of distance restraints from all of the atoms in 'atoms0' to all of the atoms in 'atoms1' where all atoms are @@ -458,7 +459,7 @@ def bond(*args, **kwargs): return distance(*args, **kwargs) -def positional(mols, atoms, k=None, r0=None, position=None, name: str = None, map=None): +def positional(mols, atoms, k=None, r0=None, position=None, name=None, map=None): """ Create a set of position restraints for the atoms specified in 'atoms' that are contained in the container 'mols', using the From cea5f0a31bc7d2748056ef9652ef8af8e69e0ca8 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 12:29:58 +0000 Subject: [PATCH 03/13] Fix bug in boresch return value If a name was supplied, BoreschRestraint was called instead of BoreschRestraints. --- src/sire/restraints/_restraints.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index d32cc301d..d493d93e3 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -293,7 +293,7 @@ def boresch( if name is None: return BoreschRestraints(b) else: - return BoreschRestraint(name, b) + return BoreschRestraints(name, b) def _check_stability_boresch_restraint(restraint_components, temperature=u("298 K")): From 83d9ad6fc4463fdda5aba49b09504f73af99a65e Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 12:33:50 +0000 Subject: [PATCH 04/13] Add tests for boresch plus minor tweaks to boresch --- src/sire/restraints/_restraints.py | 53 +++--- tests/conftest.py | 8 +- tests/restraints/test_boresch.py | 251 +++++++++++++++++++++++++++++ 3 files changed, 289 insertions(+), 23 deletions(-) create mode 100644 tests/restraints/test_boresch.py diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index d493d93e3..90fd8c59e 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -114,7 +114,9 @@ def boresch( Returns ------- BoreschRestraints : SireMM::BoreschRestraints - The Boresch restraints. + A container of Boresch restraints, where the first restraint is + the BoreschRestraint created. The Boresch restraint created can be + extracted with BoreschRestraints[0]. Examples -------- @@ -146,7 +148,11 @@ def boresch( if len(receptor) != 3 or len(ligand) != 3: # Eventually will choose the best atoms from the receptor # and ligand... - raise ValueError("You need to provide 3 receptor atoms and 3 ligand atoms") + raise ValueError( + "You need to provide 3 receptor atoms and 3 ligand atoms" + f"but only {len(receptor)} receptor atoms and {len(ligand)} " + f"ligand atoms were provided." + ) from .. import measure @@ -195,9 +201,11 @@ def boresch( # Get the force constants. if restraint_components[restraint_component]["input_k"] is None: restraint_components[restraint_component]["validated_k"] = n_measures * [ - default_distance_k - if restraint_component == "distance" - else default_angle_k + ( + default_distance_k + if restraint_component == "distance" + else default_angle_k + ) ] elif type(restraint_components[restraint_component]["input_k"]) is not list: # Populate the list with the single specified value. @@ -207,20 +215,22 @@ def boresch( else: if len(restraint_components[restraint_component]["input_k"]) == 0: # Empty list - populate with default values. - restraint_components[restraint_component][ - "validated_k" - ] = n_measures * [ - default_distance_k - if restraint_component == "distance" - else default_angle_k - ] + restraint_components[restraint_component]["validated_k"] = ( + n_measures + * [ + ( + default_distance_k + if restraint_component == "distance" + else default_angle_k + ) + ] + ) elif len(restraint_components[restraint_component]["input_k"]) == 1: # List of length 1 - populate with that value. - restraint_components[restraint_component][ - "validated_k" - ] = n_measures * [ - u(restraint_components[restraint_component]["input_k"][0]) - ] + restraint_components[restraint_component]["validated_k"] = ( + n_measures + * [u(restraint_components[restraint_component]["input_k"][0])] + ) elif ( len(restraint_components[restraint_component]["input_k"]) == n_measures ): @@ -244,11 +254,10 @@ def boresch( if type(restraint_components[restraint_component]["input_equil"]) is not list: # Only allow this if we are dealing with the distance component, as this is the only one that can be a single value. if restraint_component == "distance": - restraint_components[restraint_component][ - "input_equil" - ] = n_measures * [ - u(restraint_components[restraint_component]["input_equil"]) - ] + restraint_components[restraint_component]["input_equil"] = ( + n_measures + * [u(restraint_components[restraint_component]["input_equil"])] + ) else: raise ValueError( f"Input equilibrium values for {restraint_component} must be a list of length {n_measures} of values or Nones" diff --git a/tests/conftest.py b/tests/conftest.py index 3232687fc..321d3edd6 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,6 +1,7 @@ -import sire as sr import pytest +import sire as sr + def pytest_addoption(parser): parser.addoption( @@ -202,3 +203,8 @@ def openmm_platform(): pass return "Reference" + + +@pytest.fixture(scope="session") +def thrombin_complex(): + return sr.load_test_files("thrombin.top", "thrombin.rst7") diff --git a/tests/restraints/test_boresch.py b/tests/restraints/test_boresch.py new file mode 100644 index 000000000..6d2b7b565 --- /dev/null +++ b/tests/restraints/test_boresch.py @@ -0,0 +1,251 @@ +import pytest + +import sire as sr +from sire.restraints import boresch + +# Valid Boresch restraint parameters. +BORESCH_PARAMS_DEFAULT = { + "receptor_selection": "atomidx 1538 or atomidx 1518 or atomidx 1540", + "ligand_selection": "atomidx 4 or atomidx 3 or atomidx 5", + "kr": "6.2012 kcal mol-1 A-2", + "ktheta": ["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], + "kphi": [ + "59.8626 kcal mol-1 rad-2", + "0.7923 kcal mol-1 rad-2", + "55.1775 kcal mol-1 rad-2", + ], + "r0": "7.687 A", + "theta0": ["1.3031 rad", "1.4777 rad"], + "phi0": ["2.5569 rad", "2.9359 rad", "1.4147 rad"], + "name": "boresch_restraint_1", +} +BORESCH_PARAMS_UNNAMED = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_UNNAMED["name"] = None +BORESCH_PARAMS_NON_EXPLICIT = BORESCH_PARAMS_DEFAULT.copy() +for param in ("kr", "ktheta", "kphi", "r0", "theta0", "phi0"): + BORESCH_PARAMS_NON_EXPLICIT[param] = None +BORESCH_PARAMS_SINGLE_FORCE_CONSTS = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_SINGLE_FORCE_CONSTS["ktheta"] = ["28.7685 kcal mol-1 rad-2"] +BORESCH_PARAMS_SINGLE_FORCE_CONSTS["kphi"] = ["59.8626 kcal mol-1 rad-2"] + +# Invalid Boresch restraint parameters. +BORESCH_PARAMS_2_RECEPTOR_ATOMS = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_2_RECEPTOR_ATOMS["receptor_selection"] = "atomidx 1538 or atomidx 1518" +BORESCH_PARAMS_2_LIGAND_ATOMS = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_2_LIGAND_ATOMS["ligand_selection"] = "atomidx 4 or atomidx 3" +BORESCH_PARAMS_WRONG_NUM_KPHI = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_WRONG_NUM_KPHI["kphi"] = [ + "28.7685 kcal mol-1 rad-2", + "24.8204 kcal mol-1 rad-2", +] +BORESCH_PARAMS_KR_0 = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_KR_0["kr"] = "0.0 kcal mol-1 A-2" +BORESCH_PARAMS_KTHETA_A_0 = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_KTHETA_A_0["ktheta"] = [ + "0.0 kcal mol-1 rad-2", + "24.8204 kcal mol-1 rad-2", +] +BORESCH_PARAMS_KTHETA_B_0 = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_KTHETA_B_0["ktheta"] = [ + "28.7685 kcal mol-1 rad-2", + "0.0 kcal mol-1 rad-2", +] + +# Restraints likely to be unstable due to thetaA0 or thetaB0 being too close to 0 or pi +# or the associated force constants being too low. +BORESCH_PARAMS_UNSTABLE_FORCE_CONSTS_LOW = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_UNSTABLE_FORCE_CONSTS_LOW["ktheta"] = [ + "0.1 kcal mol-1 rad-2", + "24.8204 kcal mol-1 rad-2", +] +BORESCH_PARAMS_UNSTABLE_THETA_A_0 = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_UNSTABLE_THETA_A_0["theta0"] = ["0.1 rad", "1.4777 rad"] +BORESCH_PARAMS_UNSTABLE_THETA_B_0 = BORESCH_PARAMS_DEFAULT.copy() +BORESCH_PARAMS_UNSTABLE_THETA_B_0["theta0"] = ["1.3031 rad", "3.1 rad"] + + +# Parameterise the test with many combinations of parameters which should pass. +@pytest.mark.parametrize( + ( + "receptor_selection", + "ligand_selection", + "kr", + "ktheta", + "kphi", + "r0", + "theta0", + "phi0", + "name", + ), + [ + [val for val in BORESCH_PARAMS_DEFAULT.values()], + [val for val in BORESCH_PARAMS_UNNAMED.values()], + [val for val in BORESCH_PARAMS_NON_EXPLICIT.values()], + [val for val in BORESCH_PARAMS_SINGLE_FORCE_CONSTS.values()], + ], +) +def test_create_boresch_restraint( + thrombin_complex, + receptor_selection, + ligand_selection, + kr, + ktheta, + kphi, + r0, + theta0, + phi0, + name, +): + """ + Check that we can create a Boresch restraint with valid parameters. + """ + boresch_restraints = boresch( + thrombin_complex, + receptor=thrombin_complex["protein"][receptor_selection], + ligand=thrombin_complex["resname LIG"][ligand_selection], + kr=kr, + ktheta=ktheta, + kphi=kphi, + r0=r0, + theta0=theta0, + phi0=phi0, + name=name, + ) + + +# Parameterise the test with many combinations of parameters which should fail. +@pytest.mark.parametrize( + ( + "receptor_selection", + "ligand_selection", + "kr", + "ktheta", + "kphi", + "r0", + "theta0", + "phi0", + "name", + ), + [ + [val for val in BORESCH_PARAMS_2_RECEPTOR_ATOMS.values()], + [val for val in BORESCH_PARAMS_2_LIGAND_ATOMS.values()], + [val for val in BORESCH_PARAMS_WRONG_NUM_KPHI.values()], + [val for val in BORESCH_PARAMS_KR_0.values()], + [val for val in BORESCH_PARAMS_KTHETA_A_0.values()], + [val for val in BORESCH_PARAMS_KTHETA_B_0.values()], + ], +) +def test_create_boresch_restraint_fails( + thrombin_complex, + receptor_selection, + ligand_selection, + kr, + ktheta, + kphi, + r0, + theta0, + phi0, + name, +): + """ + Check that we can create a Boresch restraint with valid parameters. + """ + with pytest.raises(ValueError): + boresch_restraints = boresch( + thrombin_complex, + receptor=thrombin_complex["protein"][receptor_selection], + ligand=thrombin_complex["resname LIG"][ligand_selection], + kr=kr, + ktheta=ktheta, + kphi=kphi, + r0=r0, + theta0=theta0, + phi0=phi0, + name=name, + ) + + +def test_boresch_restraint_params(thrombin_complex): + """ + Check that the parameters of the created Boresch restraint are as expected. + """ + boresch_restraints = 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"], + ) + boresch_restraint = boresch_restraints[0] + + assert boresch_restraint.receptor_atoms() == [1574, 1554, 1576] + assert boresch_restraint.ligand_atoms() == [4, 3, 5] + assert boresch_restraint.kr().value() == 6.2012 + assert boresch_restraint.ktheta()[0].value() == 28.7685 + assert boresch_restraint.ktheta()[1].value() == 24.8204 + assert boresch_restraint.kphi()[0].value() == 59.8626 + assert boresch_restraint.kphi()[1].value() == 0.7923 + assert boresch_restraint.kphi()[2].value() == 55.1775 + assert boresch_restraint.r0().value() == 7.687 + assert boresch_restraint.theta0()[0].value() == 1.3031 + assert boresch_restraint.theta0()[1].value() == 1.4777 + assert boresch_restraint.phi0()[0].value() == 2.5569 + assert boresch_restraint.phi0()[1].value() == 2.9359 + assert boresch_restraint.phi0()[2].value() == 1.4147 + + +@pytest.mark.parametrize( + ( + "receptor_selection", + "ligand_selection", + "kr", + "ktheta", + "kphi", + "r0", + "theta0", + "phi0", + "name", + ), + [ + [val for val in BORESCH_PARAMS_UNSTABLE_FORCE_CONSTS_LOW.values()], + [val for val in BORESCH_PARAMS_UNSTABLE_THETA_A_0.values()], + [val for val in BORESCH_PARAMS_UNSTABLE_THETA_B_0.values()], + ], +) +def test_boresch_restraint_params_unstable( + thrombin_complex, + receptor_selection, + ligand_selection, + kr, + ktheta, + kphi, + r0, + theta0, + phi0, + name, +): + """ + Check that a warning is raised when creating a Boresch restraint with unstable parameters. + """ + with pytest.warns(UserWarning): + boresch_restraints = boresch( + thrombin_complex, + receptor=thrombin_complex["protein"][receptor_selection], + ligand=thrombin_complex["resname LIG"][ligand_selection], + kr=kr, + ktheta=ktheta, + kphi=kphi, + r0=r0, + theta0=theta0, + phi0=phi0, + name=name, + ) From d82b2a28980baf0da21993a7630f38e0daee6707 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 12:35:13 +0000 Subject: [PATCH 05/13] Add Boresch analytical correction and test --- src/sire/restraints/CMakeLists.txt | 1 + src/sire/restraints/__init__.py | 5 +- .../restraints/_standard_state_correction.py | 165 ++++++++++++++++++ .../test_standard_state_correction.py | 51 ++++++ 4 files changed, 220 insertions(+), 2 deletions(-) create mode 100644 src/sire/restraints/_standard_state_correction.py create mode 100644 tests/restraints/test_standard_state_correction.py diff --git a/src/sire/restraints/CMakeLists.txt b/src/sire/restraints/CMakeLists.txt index 6411a8319..c6fa94c80 100644 --- a/src/sire/restraints/CMakeLists.txt +++ b/src/sire/restraints/CMakeLists.txt @@ -8,6 +8,7 @@ set ( SCRIPTS __init__.py _restraints.py + _standard_state_correction.py ) # installation diff --git a/src/sire/restraints/__init__.py b/src/sire/restraints/__init__.py index 4b59de758..5f5af1465 100644 --- a/src/sire/restraints/__init__.py +++ b/src/sire/restraints/__init__.py @@ -1,3 +1,4 @@ -__all__ = ["positional", "bond", "distance", "boresch"] +__all__ = ["positional", "bond", "distance", "boresch", "get_standard_state_correction"] -from ._restraints import positional, bond, distance, boresch +from ._restraints import bond, boresch, distance, positional +from ._standard_state_correction import get_standard_state_correction diff --git a/src/sire/restraints/_standard_state_correction.py b/src/sire/restraints/_standard_state_correction.py new file mode 100644 index 000000000..17da780c2 --- /dev/null +++ b/src/sire/restraints/_standard_state_correction.py @@ -0,0 +1,165 @@ +__all__ = ["get_standard_state_correction"] + +import numpy as _np + +import sire as sr + +from .. import units as _units + + +def get_standard_state_correction(restraint, temperature=300.0 * _units.kelvin): + """ + Get the entropic correction for releasing a given restraint to + the standard state volume at a given temperature. + + Parameters + ---------- + restraint : sire.legacy.MM._MM.BoreschRestraint + The restraint for which to calculate the standard state correction. + + Returns + ------- + correction : sire.legacy.Units._Units.GeneralUnit + The standard state correction for the given restraint. + + Examples + -------- + To create a Boresch restraint and calculate the standard state correction: + + ``` + # Create the Boresch restraints object. + my_boresch_restraints = boresch( + system, + receptor=system["protein"]["atomidx 1538 or atomidx 1518 or atomidx 1540"], + ligand=system["resname LIG"]["atomidx 4 or atomidx 3 or atomidx 5"], + kr="6.2012 kcal mol-1 A-2", + ktheta=["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], + kphi=["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2", "55.1775 kcal mol-1 rad-2"], + r0="7.687 A", + theta0=["1.3031 rad", "1.4777 rad"], + phi0=["2.5569 rad", "2.9359 rad", "1.4147 rad"], + ) + + # Extract the single Boresch restraint from the Boresch restraints. + my_boresch_restraint = my_boresch_restraints[0] + + # Calculate the standard state correction for the Boresch restraint. + standard_state_correction = get_standard_state_correction(my_boresch_restraint) + print(standard_state_correction) + ``` + """ + if isinstance(restraint, sr.legacy.MM._MM.BoreschRestraint): + return _get_boresch_standard_state_correction(restraint, temperature) + else: + raise NotImplementedError( + f"Standard state correction for restraint type {type(restraint)} is not implemented. " + "This function currently only supports Boresch restraints." + ) + + +def _get_boresch_standard_state_correction(restraint, temperature): + """ + Get the entropic correction for releasing a given Boresch restraint to + the standard state volume at a given temperature. + + Parameters + ---------- + restraint : sire.legacy.MM._MM.BoreschRestraint + The Boresch restraint for which to calculate the standard state correction. + temperature : float + The temperature at which to calculate the standard state correction. + + Returns + ------- + correction : sire.legacy.Units._Units.GeneralUnit + The standard state correction for the given Boresch restraint. + + Examples + -------- + To create a Boresch restraint and calculate the standard state correction: + ``` + # Create the Boresch restraints object. + my_boresch_restraints = boresch( + system, + receptor=system["protein"]["atomidx 1538 or atomidx 1518 or atomidx 1540"], + ligand=system["resname LIG"]["atomidx 4 or atomidx 3 or atomidx 5"], + kr="6.2012 kcal mol-1 A-2", + ktheta=["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], + kphi=["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2", "55.1775 kcal mol-1 rad-2"], + r0="7.687 A", + theta0=["1.3031 rad", "1.4777 rad"], + phi0=["2.5569 rad", "2.9359 rad", "1.4147 rad"], + ) + + # Extract the single Boresch restraint from the Boresch restraints. + my_boresch_restraint = my_boresch_restraints[0] + + # Calculate the standard state correction for the Boresch restraint. + standard_state_correction = get_standard_state_correction(my_boresch_restraint) + print(standard_state_correction) + ``` + """ + # Remove units so that we can raise to non-integer powers and take sine. + + # Params. + T = (temperature / _units.kelvin).value() # K + r0 = (restraint.r0() / _units.angstrom).value() # A + thetaA0 = ( + restraint.theta0()[0] / _units.radians + ).value() # Remove units so we can take sin. + thetaB0 = ( + restraint.theta0()[1] / _units.radians + ).value() # Remove units so we can take sin. + + # Constants. + v0 = ( + (_units.meter3 / 1000) / _units.mole.value() + ).value() # A^3, the standard state volume. + R = _units.gasr.value() # kcal mol-1, the molar gas constant. + prefactor = 8 * (_np.pi**2) * v0 # Divide this to account for force constants of 0. + force_constants = [] + + # Correct for force constants of zero which break the analytical correction. + # kr + kr = ( + restraint.kr() / (_units.kcal * _units.mole.pow(-1) * _units.angstrom.pow(-2)) + ).value() + force_constants.append(kr) + + # kthetas + for i in range(2): + if restraint.ktheta()[i] == 0: + prefactor /= 2 / _np.sin((restraint.theta0()[i] / _units.radians).value()) + else: + force_const = ( + restraint.ktheta()[i] + / (_units.kcal * _units.mole.pow(-1) * _units.radians.pow(-2)) + ).value() + force_constants.append(force_const) + + # kphis + for i in range(3): + if restraint.kphi()[i] == 0: + prefactor /= 2 * _np.pi + else: + force_const = ( + restraint.kphi()[i] + / (_units.kcal * _units.mole.pow(-1) * _units.radians.pow(-2)) + ).value() + force_constants.append(force_const) + + n_nonzero_k = len(force_constants) + + # Calculation + numerator = prefactor * _np.sqrt(_np.prod(force_constants)) + denominator = ( + r0**2 + * _np.sin(thetaA0) + * _np.sin(thetaB0) + * (2 * _np.pi * R * T) ** (n_nonzero_k / 2) + ) + breakpoint() + + # Use values with units to return a sire.legacy.Units._Units.GeneralUnit object. + dg = -_units.gasr * temperature * _np.log(numerator / denominator) + return dg diff --git a/tests/restraints/test_standard_state_correction.py b/tests/restraints/test_standard_state_correction.py new file mode 100644 index 000000000..977e064dc --- /dev/null +++ b/tests/restraints/test_standard_state_correction.py @@ -0,0 +1,51 @@ +import pytest + +# Boresch parameters from old test_boresch_corr SOMD test so we can compare +# to previous results. +BORESCH_PARAMS_DEFAULT = { + "receptor_selection": "atomidx 1538 or atomidx 1518 or atomidx 1540", + "ligand_selection": "atomidx 4 or atomidx 3 or atomidx 5", + "kr": "50.98 kcal mol-1 A-2", + "ktheta": ["133.48 kcal mol-1 rad-2", "76.78 kcal mol-1 rad-2"], + "kphi": [ + "430.72 kcal mol-1 rad-2", + "98.46 kcal mol-1 rad-2", + "99.58 kcal mol-1 rad-2", + ], + "r0": "5.92 A", + "theta0": ["1.85 rad", "1.59 rad"], + "phi0": ["-0.30 rad", "-1.55 rad", "2.90 rad"], + "name": "boresch_restraint_1", +} + + +@pytest.mark.slow +def test_standard_state_correction_boresch(thrombin_complex): + """ + Check that the parameters of the created Boresch restraint are as expected. + """ + import sire as sr + from sire.restraints import boresch, get_standard_state_correction + + boresch_restraints = 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"], + ) + boresch_restraint = boresch_restraints[0] + + std_state_correction = get_standard_state_correction( + boresch_restraint, temperature=298.15 * sr.units.kelvin + ) + assert std_state_correction.value() == pytest.approx(-10.98, abs=1e-2) From 21e52baf7a2fbe6880d6331009b679774f3e965e Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 12:46:36 +0000 Subject: [PATCH 06/13] Improve formatting of examples in boresch and get_standard_state_correction doc strings --- src/sire/restraints/_restraints.py | 25 ++++++------ .../restraints/_standard_state_correction.py | 40 +++++++++---------- 2 files changed, 31 insertions(+), 34 deletions(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 90fd8c59e..ab0846961 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -123,19 +123,18 @@ def boresch( Create a set of Boresch restraints for the ligand in the system 'system', specifying all of the force constants and equilibrium values: - ``` - my_boresch_restraint = boresch( - system, - receptor=system["protein"]["atomidx 1538 or atomidx 1518 or atomidx 1540"], - ligand=system["resname LIG"]["atomidx 4 or atomidx 3 or atomidx 5"], - kr="6.2012 kcal mol-1 A-2", - ktheta=["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], - kphi=["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2", "55.1775 kcal mol-1 rad-2"], - r0="7.687 A", - theta0=["1.3031 rad", "1.4777 rad"], - phi0=["2.5569 rad", "2.9359 rad", "1.4147 rad"], - ) - ``` + >>> my_boresch_restraints = boresch( + >>> system, + >>> receptor=system["protein"]["atomidx 1538 or atomidx 1518 or atomidx 1540"], + >>> ligand=system["resname LIG"]["atomidx 4 or atomidx 3 or atomidx 5"], + >>> kr="6.2012 kcal mol-1 A-2", + >>> ktheta=["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], + >>> kphi=["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2", "55.1775 kcal mol-1 rad-2"], + >>> r0="7.687 A", + >>> theta0=["1.3031 rad", "1.4777 rad"], + >>> phi0=["2.5569 rad", "2.9359 rad", "1.4147 rad"], + >>> ) + >>> my_boresch_restraint = my_boresch_restraits[0] """ from ..base import create_map from ..mm import BoreschRestraint, BoreschRestraints diff --git a/src/sire/restraints/_standard_state_correction.py b/src/sire/restraints/_standard_state_correction.py index 17da780c2..ee1812ba4 100644 --- a/src/sire/restraints/_standard_state_correction.py +++ b/src/sire/restraints/_standard_state_correction.py @@ -26,27 +26,25 @@ def get_standard_state_correction(restraint, temperature=300.0 * _units.kelvin): -------- To create a Boresch restraint and calculate the standard state correction: - ``` - # Create the Boresch restraints object. - my_boresch_restraints = boresch( - system, - receptor=system["protein"]["atomidx 1538 or atomidx 1518 or atomidx 1540"], - ligand=system["resname LIG"]["atomidx 4 or atomidx 3 or atomidx 5"], - kr="6.2012 kcal mol-1 A-2", - ktheta=["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], - kphi=["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2", "55.1775 kcal mol-1 rad-2"], - r0="7.687 A", - theta0=["1.3031 rad", "1.4777 rad"], - phi0=["2.5569 rad", "2.9359 rad", "1.4147 rad"], - ) - - # Extract the single Boresch restraint from the Boresch restraints. - my_boresch_restraint = my_boresch_restraints[0] - - # Calculate the standard state correction for the Boresch restraint. - standard_state_correction = get_standard_state_correction(my_boresch_restraint) - print(standard_state_correction) - ``` + >>> # Create the Boresch restraints object. + >>> my_boresch_restraints = boresch( + >>> system, + >>> receptor=system["protein"]["atomidx 1538 or atomidx 1518 or atomidx 1540"], + >>> ligand=system["resname LIG"]["atomidx 4 or atomidx 3 or atomidx 5"], + >>> kr="6.2012 kcal mol-1 A-2", + >>> ktheta=["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], + >>> kphi=["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2", "55.1775 kcal mol-1 rad-2"], + >>> r0="7.687 A", + >>> theta0=["1.3031 rad", "1.4777 rad"], + >>> phi0=["2.5569 rad", "2.9359 rad", "1.4147 rad"], + >>> ) + >>> + >>> # Extract the single Boresch restraint from the Boresch restraints. + >>> my_boresch_restraint = my_boresch_restraints[0] + >>> + >>> # Calculate the standard state correction for the Boresch restraint. + >>> standard_state_correction = get_standard_state_correction(my_boresch_restraint) + >>> print(standard_state_correction) """ if isinstance(restraint, sr.legacy.MM._MM.BoreschRestraint): return _get_boresch_standard_state_correction(restraint, temperature) From d0a0268e271b0d300223115bbf21a50757565876 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 12:49:19 +0000 Subject: [PATCH 07/13] Fix typo in name of sphinx-issues in development guide --- doc/source/contributing/development.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/contributing/development.rst b/doc/source/contributing/development.rst index 275f5c4f7..50b9ea699 100644 --- a/doc/source/contributing/development.rst +++ b/doc/source/contributing/development.rst @@ -444,7 +444,7 @@ additional packages as described in the .. code-block:: bash - conda install sphinx sphinxcontrib-programoutput sphinx_issues furo + conda install sphinx sphinxcontrib-programoutput sphinx-issues furo Then move to the ``doc`` directory and run: From cdc53e08579e1110e27426d09ad51eaefab0ede9 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 13:10:35 +0000 Subject: [PATCH 08/13] Add missing space in boresch example --- src/sire/restraints/_restraints.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index ab0846961..237b8c094 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -123,6 +123,7 @@ def boresch( Create a set of Boresch restraints for the ligand in the system 'system', specifying all of the force constants and equilibrium values: + >>> my_boresch_restraints = boresch( >>> system, >>> receptor=system["protein"]["atomidx 1538 or atomidx 1518 or atomidx 1540"], @@ -134,7 +135,7 @@ def boresch( >>> theta0=["1.3031 rad", "1.4777 rad"], >>> phi0=["2.5569 rad", "2.9359 rad", "1.4147 rad"], >>> ) - >>> my_boresch_restraint = my_boresch_restraits[0] + >>> my_boresch_restraint = my_boresch_restraints[0] """ from ..base import create_map from ..mm import BoreschRestraint, BoreschRestraints From 2a538fe5c0c545f6a8789e329626632286726f3d Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 14:48:48 +0000 Subject: [PATCH 09/13] Allow boresch arguments to be passed using a map --- src/sire/restraints/_restraints.py | 19 ++++++++++++++++++- tests/restraints/test_boresch.py | 18 ++++++++++++++++++ 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/src/sire/restraints/_restraints.py b/src/sire/restraints/_restraints.py index 237b8c094..d4c64774e 100644 --- a/src/sire/restraints/_restraints.py +++ b/src/sire/restraints/_restraints.py @@ -105,7 +105,10 @@ def boresch( used. Default is None. map : dict, optional - A property map. Default is None. + A dictionary of additional options. Note that any options + set in this dictionary that are also specified via one of + the arguments above will be overridden by the argument + value temperature : str or SireUnits::Dimension::GeneralUnit, optional The temperature to use when checking for unstable restraints. If @@ -140,7 +143,21 @@ def boresch( from ..base import create_map from ..mm import BoreschRestraint, BoreschRestraints + # If an argument is not specified in the function + # arguments, but is specified in the map, update + # the argument from the map. map = create_map(map) + map_dict = map.to_dict() + kr = kr if kr is not None else map_dict.get("kr", None) + ktheta = ktheta if ktheta is not None else map_dict.get("ktheta", None) + kphi = kphi if kphi is not None else map_dict.get("kphi", None) + r0 = r0 if r0 is not None else map_dict.get("r0", None) + theta0 = theta0 if theta0 is not None else map_dict.get("theta0", None) + phi0 = phi0 if phi0 is not None else map_dict.get("phi0", None) + name = name if name is not None else map_dict.get("name", None) + temperature = ( + temperature if temperature is not None else map_dict.get("temperature", None) + ) receptor = _to_atoms(mols, receptor) ligand = _to_atoms(mols, ligand) diff --git a/tests/restraints/test_boresch.py b/tests/restraints/test_boresch.py index 6d2b7b565..b134ee975 100644 --- a/tests/restraints/test_boresch.py +++ b/tests/restraints/test_boresch.py @@ -249,3 +249,21 @@ def test_boresch_restraint_params_unstable( phi0=phi0, name=name, ) + + +def test_boresch_creation_with_map(thrombin_complex): + """ + Test that we can provide options through the map. + """ + boresch_restraints = boresch( + thrombin_complex, + receptor=thrombin_complex["protein"][ + BORESCH_PARAMS_DEFAULT["receptor_selection"] + ], + ligand=thrombin_complex["resname LIG"][ + BORESCH_PARAMS_DEFAULT["ligand_selection"] + ], + map={"kr": "44 kcal mol-1 A-2"}, + ) + boresch_restraint = boresch_restraints[0] + assert boresch_restraint.kr().value() == 44.0 From 155d41993b80b31a7fb107df11a323d2feeea580 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 15:51:08 +0000 Subject: [PATCH 10/13] Tiny tweak to how default temperature is specified --- src/sire/restraints/_standard_state_correction.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/sire/restraints/_standard_state_correction.py b/src/sire/restraints/_standard_state_correction.py index ee1812ba4..37b45711c 100644 --- a/src/sire/restraints/_standard_state_correction.py +++ b/src/sire/restraints/_standard_state_correction.py @@ -4,10 +4,11 @@ import sire as sr +from .. import u as _u from .. import units as _units -def get_standard_state_correction(restraint, temperature=300.0 * _units.kelvin): +def get_standard_state_correction(restraint, temperature=_u("300 K")): """ Get the entropic correction for releasing a given restraint to the standard state volume at a given temperature. @@ -17,6 +18,9 @@ def get_standard_state_correction(restraint, temperature=300.0 * _units.kelvin): restraint : sire.legacy.MM._MM.BoreschRestraint The restraint for which to calculate the standard state correction. + temperature : sire.legacy.Units._Units.GeneralUnit, optional + The temperature at which to calculate the standard state correction. + Returns ------- correction : sire.legacy.Units._Units.GeneralUnit From 8c773062001b6b31a307e08227d631ff11099ab9 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 16:03:28 +0000 Subject: [PATCH 11/13] Add Boresch restraints to the tutorial --- doc/source/tutorial/part06/03_restraints.rst | 79 ++++++++++++++++++++ 1 file changed, 79 insertions(+) diff --git a/doc/source/tutorial/part06/03_restraints.rst b/doc/source/tutorial/part06/03_restraints.rst index e715aca9e..0b35d49fd 100644 --- a/doc/source/tutorial/part06/03_restraints.rst +++ b/doc/source/tutorial/part06/03_restraints.rst @@ -229,6 +229,85 @@ BondRestraints( size=630 the first atom in the first molecule and the oxygen atoms in all of the water molecules. +Boresch Restraints +--------------------------- + +The :func:`sire.restraints.boresch` function is used to create Boresch restraints, +which are commonly used for restraining the relative positions and orientations +of a receptor and ligand during alchemical absolute binding free energy calculations. +They restrain the six relative external degrees of freedom of the receptor and ligand +by restraining one distance, two angles, and three dihedrals angles which are +defined according to three anchor points in the receptor and three anchor points +in the ligand. For more detail, please see J. Phys. Chem. B 2003, 107, 35, 9535–9551. + +To create a Boresch restraint, you need to specify the receptor and ligand anchor +atoms (note that the order of the atoms is important). Like the distance restraints, +the atoms can be specified using a search string, passing lists of atom indexes, or +molecule views holding the atoms. You can also specify the force constants and equilibrium +values for the restraints. If not supplied, default force constants of 10 kcal mol-1 Å-2 +and 100 kcal mol-1 rad-2 are used for the distance and angle restraints, respectively, +and the equilibrium values are set to the current values of the distances and angles in +the system supplied. For example, + +>>> boresch_restraints = sr.restraints.boresch(mols, receptor=[1574, 1554, 1576], ligand=[4,3,5]) +>>> boresch_restraint = boresch_restraints[0] +>>> print(boresch_restraint) +BoreschRestraint( [1574, 1554, 1576] => [4, 3, 5], + k=[10 kcal mol-1 Å-2, 0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2, + 0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2, 0.0304617 kcal mol-1 °-2] + r0=15.1197 Å, theta0=[80.5212°, 59.818°], + 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 force constants have been set +and the equilibrium values have been set to the current values of the distances and angles in the +system supplied. + +.. note:: + + Boresch restraints can be subject to instabilities if any three contiguous anchor points + approach collinearity (J. Chem. Theory Comput. 2023, 19, 12, 3686–3704). It is important to + prevent this by ensuring the associated angles are sufficiently far from 0 or 180 degrees, + and that the `ktheta` force constants are high enough. Sire will raise a warning if the + `theta0` values are too close to 0 or 180 degrees for the given temperature and force constants. + +Alternatively, we could have explicitly set the force constants and equilibrium values, e.g. + +>>> boresch_restraints = sr.restraints.boresch(mols, receptor=[1574, 1554, 1576], ligand=[4,3,5], + kr = "6.2012 kcal mol-1 A-2", + ktheta = ["28.7685 kcal mol-1 rad-2", "24.8204 kcal mol-1 rad-2"], + kphi = ["59.8626 kcal mol-1 rad-2", "0.7923 kcal mol-1 rad-2", + "55.1775 kcal mol-1 rad-2"], + r0 = "16 A", + theta0 = ["1.2 rad", "1.3 rad"], + phi0 = ["2.2 rad", "2.5 rad", "1.5 rad"], + name = "boresch_restraint_1") +>>> boresch_restraint = boresch_restraints[0] +>>> print(boresch_restraint) +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°] ) + +.. note:: + + :func:`sire.restraints.boresch` returns a list of Boresch restraints. If you are only + interested in a single Boresch restraint, you can extract it with the index, e.g. + ``boresch_restraint = boresch_restraints[0]``. + +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 +selected (see 10.26434/chemrxiv-2023-8s9dz-v2). This can be calculated with +:func:`sire.restraints.standard_state_correction`: + +>>> boresch_restraint = boresch_restraints[0] +>>> from sire import u +>>> correction = sr.restraints.standard_state_correction(boresch_restraint, temperature=u("298 K")) +>>> print(correction) +-6.2399 kcal mol-1 + + Using restraints in minimisation or dynamics -------------------------------------------- From deaff85ad2b8232e7382a80abe4e844b19b83d89 Mon Sep 17 00:00:00 2001 From: finlayclark Date: Sun, 11 Feb 2024 16:10:54 +0000 Subject: [PATCH 12/13] Update changelog --- doc/source/changelog.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 5b3473017..221b99d50 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -17,6 +17,13 @@ organisation on `GitHub `__. * Please add an item to this changelog when you create your PR +* Added more support for Boresch restraints. Specifically, :func:`sire.restraints.boresch` + now supports the specification of equilibrium values, uses different default force + constants, and warns the user if the restraints are likely to be unstable. + :func:`sire.restraints.get_standard_state_correction` was implemented for Boresch + restraints. Tests were added for restraint creation and for the standard state + correction. Boresch restraints were added to :doc:`tutorial `. + `2023.5.1 `__ - January 2024 -------------------------------------------------------------------------------------------- From f4e440fac6766126f9c30db13efe0151477ce32f Mon Sep 17 00:00:00 2001 From: finlayclark Date: Tue, 20 Feb 2024 20:18:35 +0000 Subject: [PATCH 13/13] Autoformatting changes [ci skip] Also remove forgotten breakpoint. --- .../restraints/_standard_state_correction.py | 1 - tests/somd/test_boresch_corr.py | 3 +- wrapper/Tools/BoreschAnalyticalCorrection.py | 53 +++++++++++-------- 3 files changed, 33 insertions(+), 24 deletions(-) diff --git a/src/sire/restraints/_standard_state_correction.py b/src/sire/restraints/_standard_state_correction.py index 37b45711c..3420b5eb2 100644 --- a/src/sire/restraints/_standard_state_correction.py +++ b/src/sire/restraints/_standard_state_correction.py @@ -160,7 +160,6 @@ def _get_boresch_standard_state_correction(restraint, temperature): * _np.sin(thetaB0) * (2 * _np.pi * R * T) ** (n_nonzero_k / 2) ) - breakpoint() # Use values with units to return a sire.legacy.Units._Units.GeneralUnit object. dg = -_units.gasr * temperature * _np.log(numerator / denominator) diff --git a/tests/somd/test_boresch_corr.py b/tests/somd/test_boresch_corr.py index 2aaf7e604..fdf488146 100644 --- a/tests/somd/test_boresch_corr.py +++ b/tests/somd/test_boresch_corr.py @@ -1,9 +1,10 @@ """Test the Boresch analytical and numerical correction methods.""" +import os import shlex import subprocess -import os import sys + import pytest try: diff --git a/wrapper/Tools/BoreschAnalyticalCorrection.py b/wrapper/Tools/BoreschAnalyticalCorrection.py index 5472fb638..86e91c868 100644 --- a/wrapper/Tools/BoreschAnalyticalCorrection.py +++ b/wrapper/Tools/BoreschAnalyticalCorrection.py @@ -7,22 +7,25 @@ # @author: Finlay Clark # Based very closely on StandardState.py by Stefano Bosisio and Julien Michel -import numpy as np import os -from math import pi, sin, log -from Sire import Units +from math import log, pi, sin +import numpy as np +from Sire import Units from Sire.Tools import Parameter, resolveParameters from Sire.Tools.OpenMMMD import * # Constants -v0 = ((Units.meter3/1000)/Units.mole.value()).value() # A^3, the standard state volume -R = Units.gasr.value() # kcal mol-1, the molar gas constant +v0 = ( + (Units.meter3 / 1000) / Units.mole.value() +).value() # A^3, the standard state volume +R = Units.gasr.value() # kcal mol-1, the molar gas constant + @resolveParameters def run(): try: - host = os.environ['HOSTNAME'] + host = os.environ["HOSTNAME"] except KeyError: host = "unknown" print("### Running Standard state correction calculation on %s ###" % host) @@ -31,23 +34,23 @@ def run(): print("###====================Utilised Parameters=====================###") print(temperature) print(boresch_restraints_dict) - print ("###===========================================================###\n") + print("###===========================================================###\n") # Get Boresch restraint dict in dict form boresch_dict = dict(boresch_restraints_dict.val) # Params - T = temperature.val.value() # K - r0 = boresch_dict['equilibrium_values']['r0'] # A - thetaA0 = boresch_dict["equilibrium_values"]["thetaA0"] # rad - thetaB0 = boresch_dict["equilibrium_values"]["thetaB0"] # rad + T = temperature.val.value() # K + r0 = boresch_dict["equilibrium_values"]["r0"] # A + thetaA0 = boresch_dict["equilibrium_values"]["thetaA0"] # rad + thetaB0 = boresch_dict["equilibrium_values"]["thetaB0"] # rad # Force constants defined as E = k*x**2, so need to multiply all force constants # by 2 to correct for original definition (E= 0.5*k*x**2) for k in boresch_dict["force_constants"]: boresch_dict["force_constants"][k] *= 2 - prefactor = 8*(pi**2)*v0 # Divide this to account for force constants of 0 + prefactor = 8 * (pi**2) * v0 # Divide this to account for force constants of 0 force_constants = [] # Loop through and correct for force constants of zero, @@ -58,23 +61,29 @@ def run(): print("Error: kr must not be zero") sys.exit(-1) if k == "kthetaA": - prefactor /= 2/sin(thetaA0) + prefactor /= 2 / sin(thetaA0) if k == "kthetaB": - prefactor /= 2/sin(thetaB0) + prefactor /= 2 / sin(thetaB0) if k[:4] == "kphi": - prefactor /= 2*pi + prefactor /= 2 * pi else: force_constants.append(val) n_nonzero_k = len(force_constants) # Calculation - numerator = prefactor*np.sqrt(np.prod(force_constants)) - denominator = (r0**2)*sin(thetaA0)*sin(thetaB0)*(2*pi*R*T)**(n_nonzero_k/2) + numerator = prefactor * np.sqrt(np.prod(force_constants)) + denominator = ( + (r0**2) * sin(thetaA0) * sin(thetaB0) * (2 * pi * R * T) ** (n_nonzero_k / 2) + ) - dg = -R*T*log(numerator/denominator) + dg = -R * T * log(numerator / denominator) - print(f"Analytical correction for releasing Boresch restraints = {dg:.2f} kcal mol-1") - print("WARNING !!! The analytical correction is only reliable when restraints are " - "sufficiently strong, r is sufficiently far from zero, and r2, r1, l1, and " - "r1, l1, l2 are sufficiently far from collinear") + print( + f"Analytical correction for releasing Boresch restraints = {dg:.2f} kcal mol-1" + ) + print( + "WARNING !!! The analytical correction is only reliable when restraints are " + "sufficiently strong, r is sufficiently far from zero, and r2, r1, l1, and " + "r1, l1, l2 are sufficiently far from collinear" + )