From 315f6350ffd6ede23787a40e9e6ffd56d0938882 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 28 Oct 2026 15:27:19 +0000 Subject: [PATCH] Backport fix from PR #514. [ci skip] --- src/sire/restraints/_boresch_search.py | 61 ++++++++++++++++++++++--- tests/restraints/test_boresch_search.py | 22 +++++++++ 2 files changed, 76 insertions(+), 7 deletions(-) diff --git a/src/sire/restraints/_boresch_search.py b/src/sire/restraints/_boresch_search.py index 0674e8aff..fae93065f 100644 --- a/src/sire/restraints/_boresch_search.py +++ b/src/sire/restraints/_boresch_search.py @@ -144,10 +144,24 @@ def _circular_mean_std(values): return mean, std +# Bonded sp2/sp3 angles average well below this, while sp centres average +# ~170 degrees or more. +_MAX_ANCHOR_ANGLE = _np.deg2rad(160.0) + + +def _anchors_collinear(ang_r, ang_l, s_idx): + """ + Whether the r1-r2-r3 or l1-l2-l3 angle of a sextuplet is near collinear, + which leaves the phiA or phiC dihedral undefined. + """ + return max(ang_r[s_idx].mean(), ang_l[s_idx].mean()) >= _MAX_ANCHOR_ANGLE + + 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). + sextuplet, plus the anchor-internal angles that phiA and phiC depend on. + Returns arrays of shape (n_sextuplets, n_frames). """ n_sext = len(sextuplets) @@ -157,11 +171,19 @@ def _sample_dofs(system, pert_mol_num, sextuplets, n_frames): 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 + ang_r = _np.zeros((n_sext, n_frames)) # r1-r2-r3 + ang_l = _np.zeros((n_sext, n_frames)) # l1-l2-l3 for f_idx, frame in enumerate(system.trajectory()): space = frame.space() lig_frame = frame.molecule(pert_mol_num) rec_frame_cache = {} + angle_cache = {} + + def _angle(key, a, b, c): + if key not in angle_cache: + angle_cache[key] = float(_sr.measure(a, b, c).value()) + return angle_cache[key] for s_idx, s in enumerate(sextuplets): r1_mol_num = s["r1_mol_num"] @@ -184,8 +206,14 @@ def _sample_dofs(system, pert_mol_num, sextuplets, n_frames): 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()) + ang_r[s_idx, f_idx] = _angle( + (r1_mol_num, s["r1_idx"], s["r2_idx"], s["r3_idx"]), r1, r2, r3 + ) + ang_l[s_idx, f_idx] = _angle( + (pert_mol_num, s["l1_idx"], s["l2_idx"], s["l3_idx"]), l1, l2, l3 + ) - return dof_r, dof_tA, dof_tB, dof_pA, dof_pB, dof_pC + return dof_r, dof_tA, dof_tB, dof_pA, dof_pB, dof_pC, ang_r, ang_l def _best_starting_frame( @@ -836,7 +864,7 @@ def _nearest_candidate(atom_idx, connectivity): # ------------------------------------------------------------------------- # 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( + dof_r, dof_tA, dof_tB, dof_pA, dof_pB, dof_pC, ang_r, ang_l = _sample_dofs( system, pert_mol_num, valid_sextuplets, n_frames ) @@ -844,23 +872,37 @@ def _nearest_candidate(atom_idx, connectivity): # 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. + # Sextuplets with a near-collinear r1-r2-r3 or l1-l2-l3 angle are also + # rejected, since phiA or phiC is then undefined. # ------------------------------------------------------------------------- _min_angle = _np.deg2rad(45.0) _max_angle = _np.deg2rad(135.0) - angle_ok_idxs = [ + theta_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: + if not theta_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." ) + angle_ok_idxs = [ + s_idx for s_idx in theta_ok_idxs if not _anchors_collinear(ang_r, ang_l, s_idx) + ] + + if not angle_ok_idxs: + raise ValueError( + "All hydrogen-bonded Boresch sextuplets were rejected because their " + "r1-r2-r3 or l1-l2-l3 anchor angle is near collinear, e.g. through a " + "nitrile or alkyne. Try adjusting 'protein_selection' or the " + "hydrogen-bond cutoffs." + ) + 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 @@ -1266,12 +1308,14 @@ def _rec_mol(mol_num): # ------------------------------------------------------------------------- # 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( + dof_r, dof_tA, dof_tB, dof_pA, dof_pB, dof_pC, ang_r, ang_l = _sample_dofs( system, pert_mol_num, valid_sextuplets, n_frames ) # ------------------------------------------------------------------------- # 6. Score sextuplets by configurational volume; apply stability filter. + # Sextuplets with a near-collinear r1-r2-r3 or l1-l2-l3 angle are also + # rejected, since phiA or phiC is then undefined. # # 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) @@ -1293,6 +1337,9 @@ def _rec_mol(mol_num): ): continue + if _anchors_collinear(ang_r, ang_l, s_idx): + 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]) @@ -1312,7 +1359,7 @@ def _rec_mol(mol_num): if not scored: raise ValueError( - "All candidate Boresch sextuplets were rejected because one or both " + "All candidate Boresch sextuplets were rejected because one or more " "anchor angles have a mean near 0° or 180° (collinearity instability). " "Try increasing the cutoff or adjusting the receptor selection." ) diff --git a/tests/restraints/test_boresch_search.py b/tests/restraints/test_boresch_search.py index 91923eb71..a8d7b362f 100644 --- a/tests/restraints/test_boresch_search.py +++ b/tests/restraints/test_boresch_search.py @@ -139,6 +139,15 @@ def test_angles_not_collinear(self, boresch_result): deg = float(theta.to(sr.units.degrees)) assert 10.0 < deg < 170.0 + def test_anchor_triplets_not_collinear(self, boresch_result): + """The r1-r2-r3 and l1-l2-l3 angles must be away from 180°.""" + restraints, _, start = boresch_result + r = restraints.at(0) + a = start.atoms() + for idxs in (list(r.receptor_atoms()), list(r.ligand_atoms())): + deg = float(sr.measure(*[a[i] for i in idxs]).to(sr.units.degrees)) + assert deg < 160.0 + def test_force_constants_positive(self, boresch_result): restraints, _, _ = boresch_result r = restraints.at(0) @@ -322,3 +331,16 @@ def test_invalid_protocol_raises(self, abfe_system): with pytest.raises(ValueError, match="protocol"): check_boresch_search(abfe_system, protocol="nonsense") + + +def test_anchors_collinear(): + """Exocyclic ring angles pass; a thermally bent nitrile angle is rejected.""" + from sire.restraints._boresch_search import _anchors_collinear + + rng = np.random.default_rng(42) + ring = np.deg2rad(rng.normal(135.0, 4.0, (1, 200))) + nitrile = np.deg2rad(180.0 - np.abs(rng.normal(0.0, 6.0, (1, 200)))) + + assert not _anchors_collinear(ring, ring, 0) + assert _anchors_collinear(ring, nitrile, 0) + assert _anchors_collinear(nitrile, ring, 0)