Switch to simplified patchIntrascale wrapper function.

Author: lohedges

src/BioSimSpace/Align/_align.py

@@ -2035,6 +2035,7 @@ def merge(
     roi=None,
     property_map0={},
     property_map1={},
+    **kwargs,
 ):
     """
     Create a merged molecule from 'molecule0' and 'molecule1' based on the
@@ -2181,6 +2182,7 @@ def merge(
         roi=roi,
         property_map0=property_map0,
         property_map1=property_map1,
+        **kwargs,
     )
 
 

src/BioSimSpace/Align/_merge.py

@@ -38,6 +38,7 @@ def merge(
     roi=None,
     property_map0={},
     property_map1={},
+    **kwargs,
 ):
     """
     Merge this molecule with 'other'.
@@ -1313,34 +1314,34 @@ def merge(
     else:
         edit_mol.set_property("connectivity0", conn0)
         edit_mol.set_property("connectivity1", conn1)
-    # Merge the intrascale properties of the two molecules.
+    # Build the intrascale matrices from the per-state connectivity.
     ff = molecule0.property(ff0)
     sf14 = _SireMM.CLJScaleFactor(
         ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor()
     )
-    if roi is not None:
-        # For ROI protein merges, build the intrascale matrices directly from
-        # the per-state connectivity, bypassing mergeIntrascale. Protein
-        # mutations involve no ring-breaking and no GLYCAM-style per-pair
-        # overrides, so the overrideIntrascale step in mergeIntrascale is a
-        # no-op. Bypassing it avoids a Windows-specific performance issue
-        # where the O(n²) override loop is very slow for large proteins.
-        # TODO: investigate the root cause and remove this workaround.
-        merged_intrascale = [
-            _SireMM.CLJNBPairs(conn0, sf14),
-            _SireMM.CLJNBPairs(conn1, sf14),
-        ]
-    else:
-        merged_intrascale = _SireIO.mergeIntrascale(
+    intra0 = _SireMM.CLJNBPairs(conn0, sf14)
+    intra1 = _SireMM.CLJNBPairs(conn1, sf14)
+
+    # For non-ROI merges, patch with any non-default per-pair scale factors
+    # from the individual molecule intrascales (e.g. GLYCAM funct=2 overrides).
+    # This step can be skipped via apply_scale_factors=False when the caller
+    # knows no non-default scale factors are present.
+    if "apply_scale_factors" in kwargs:
+        if not isinstance(kwargs["apply_scale_factors"], bool):
+            raise TypeError("'apply_scale_factors' must be of type 'bool'")
+
+    if kwargs.get("apply_scale_factors", roi is None):
+        intra0, intra1 = _SireIO.patchIntrascale(
             molecule0.property("intrascale"),
             molecule1.property("intrascale"),
-            conn0,
-            conn1,
-            sf14,
+            intra0,
+            intra1,
             mol0_merged_mapping,
             mol1_merged_mapping,
         )
 
+    merged_intrascale = [intra0, intra1]
+
     # Store the two molecular components.
     edit_mol.set_property("molecule0", molecule0)
     edit_mol.set_property("molecule1", molecule1)

src/BioSimSpace/Sandpit/Exscientia/Align/_align.py

@@ -1370,6 +1370,7 @@ def merge(
     roi=None,
     property_map0={},
     property_map1={},
+    **kwargs,
 ):
     """
     Create a merged molecule from 'molecule0' and 'molecule1' based on the
@@ -1504,6 +1505,7 @@ def merge(
         roi=roi,
         property_map0=property_map0,
         property_map1=property_map1,
+        **kwargs,
     )
 
 

src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py

@@ -37,6 +37,7 @@ def merge(
     roi=None,
     property_map0={},
     property_map1={},
+    **kwargs,
 ):
     """
     Merge this molecule with 'other'.
@@ -1306,30 +1307,30 @@ def merge(
         ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor()
     )
 
-    # Merge the intrascale properties of the two molecules.
-    if roi is not None:
-        # For ROI protein merges, build the intrascale matrices directly from
-        # the per-state connectivity, bypassing mergeIntrascale. Protein
-        # mutations involve no ring-breaking and no GLYCAM-style per-pair
-        # overrides, so the overrideIntrascale step in mergeIntrascale is a
-        # no-op. Bypassing it avoids a Windows-specific performance issue
-        # where the O(n²) override loop is very slow for large proteins.
-        # TODO: investigate the root cause and remove this workaround.
-        merged_intrascale = [
-            _SireMM.CLJNBPairs(conn0, sf14),
-            _SireMM.CLJNBPairs(conn1, sf14),
-        ]
-    else:
-        merged_intrascale = _SireIO.mergeIntrascale(
+    # Build the intrascale matrices from the per-state connectivity.
+    intra0 = _SireMM.CLJNBPairs(conn0, sf14)
+    intra1 = _SireMM.CLJNBPairs(conn1, sf14)
+
+    # For non-ROI merges, patch with any non-default per-pair scale factors
+    # from the individual molecule intrascales (e.g. GLYCAM funct=2 overrides).
+    # This step can be skipped via apply_scale_factors=False when the caller
+    # knows no non-default scale factors are present.
+    if "apply_scale_factors" in kwargs:
+        if not isinstance(kwargs["apply_scale_factors"], bool):
+            raise TypeError("'apply_scale_factors' must be of type 'bool'")
+
+    if kwargs.get("apply_scale_factors", roi is None):
+        intra0, intra1 = _SireIO.patchIntrascale(
             molecule0.property("intrascale"),
             molecule1.property("intrascale"),
-            conn0,
-            conn1,
-            sf14,
+            intra0,
+            intra1,
             mol0_merged_mapping,
             mol1_merged_mapping,
         )
 
+    merged_intrascale = [intra0, intra1]
+
     # Store the two molecular components.
     edit_mol.set_property("molecule0", molecule0)
     edit_mol.set_property("molecule1", molecule1)

tests/Align/test_align.py

@@ -857,13 +857,14 @@ def test_ring_opening_and_size_change(ligands, mapping):
 )
 def test_ring_breaking_intrascale():
     """
-    Regression test for mergeIntrascale with ring-breaking perturbations.
-
-    Before the fix, CLJNBPairs pairs that were excluded/1-4 in one end state
-    but fully interacting (1,1) in the other were incorrectly left at the
-    non-(1,1) value in the merged intrascale. For cyclopentane->cyclohexane
-    this produced only 6 changed OpenMM exceptions instead of the correct 18,
-    causing simulation instabilities due to missing non-bonded interactions.
+    Test that ring-breaking merges produce correct intrascale matrices for a
+    standard force field (GAFF2) with no non-default per-pair scale factors.
+
+    The intrascale matrices are built from per-state connectivity, which
+    correctly captures bonded distances across the ring-closure bond in the
+    merged atom space. Since GAFF2 has no non-default per-pair scale factors,
+    patchIntrascale is a no-op — verified by checking that apply_scale_factors=False
+    gives the same changed bond and exception counts as the default path.
     """
     # Parameterise both molecules with GAFF2.
     cyclopentane = BSS.Parameters.gaff2("C1CCCC1").getMolecule()
@@ -928,6 +929,22 @@ def test_ring_breaking_intrascale():
     assert len(omm_rev.changed_bonds()) == len(ref_bonds)
     assert len(omm_rev.changed_exceptions()) == len(ref_exceptions)
 
+    # Verify patchIntrascale is a no-op for GAFF2: skipping it with
+    # apply_scale_factors=False must give identical results.
+    merged_nopatch = BSS.Align.merge(
+        cyclopentane_aligned,
+        cyclohexane,
+        mapping,
+        allow_ring_size_change=True,
+        allow_ring_breaking=True,
+        apply_scale_factors=False,
+    )
+    omm_nopatch = merged_nopatch._sire_object.perturbation().to_openmm(
+        map={"coordinates": "coordinates0"}
+    )
+    assert len(omm_nopatch.changed_bonds()) == len(ref_bonds)
+    assert len(omm_nopatch.changed_exceptions()) == len(ref_exceptions)
+
 
 @pytest.mark.skipif(
     not has_antechamber or not has_tleap,
@@ -937,24 +954,20 @@ def test_ring_breaking_intrascale():
     not has_openff,
     reason="Requires OpenFF to be installed.",
 )
-def test_ring_breaking_intrascale_connectivity():
+def test_ring_breaking_intrascale_m338():
     """
-    Regression test for mergeIntrascale using a real-world ring-breaking
-    perturbation (int1 -> m338).
-
-    Validates that the merged intrascale matrices produced by mergeIntrascale
-    (which builds CLJNBPairs from per-state connectivity and overrides with
-    per-pair values) exactly match those produced directly from
-    CLJNBPairs(conn0/conn1, sf14). For ring-breaking perturbations the two
-    approaches must agree: discrepancies indicate that bonded distances in the
-    merged topology are not being correctly captured, which causes missing or
-    spurious OpenMM exceptions and simulation instabilities.
+    Test that ring-breaking merges produce correct intrascale matrices for a
+    real-world perturbation (int1 -> m338) with a standard force field (OpenFF).
+
+    Since OpenFF has no non-default per-pair scale factors, patchIntrascale is
+    a no-op and the merged intrascale matrices must exactly match those built
+    directly from CLJNBPairs(conn0/conn1, sf14).
     """
     from sire.legacy import CAS as _SireCAS
     from sire.legacy import MM as _SireMM
     from sire.legacy import Mol as _SireMol
 
-    # Atom mapping: {int1_idx: m338_idx}, read from m338_int1_MCS.txt.
+    # Atom mapping: {int1_idx: m338_idx}
     mapping = {
         21: 0,
         0: 1,
@@ -1013,11 +1026,11 @@ def test_ring_breaking_intrascale_connectivity():
 
     sire_mol = merged._sire_object
 
-    # Extract the intrascale matrices produced by mergeIntrascale.
     intra0 = sire_mol.property("intrascale0")
     intra1 = sire_mol.property("intrascale1")
 
-    # Build the reference intrascale matrices from per-state connectivity,
+    # Build reference matrices directly from per-state connectivity. For OpenFF
+    # these must be identical to the merge output since patchIntrascale is a no-op.
     ff = mol0._sire_object.property("forcefield")
     sf14 = _SireMM.CLJScaleFactor(
         ff.electrostatic14_scale_factor(), ff.vdw14_scale_factor()