Backport fix from PR #514

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,31 +864,45 @@ 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
     )
 
     # -------------------------------------------------------------------------
     # 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."
         )

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)