Fix issue #447

corelib/src/libs/SireIO/biosimspace.cpp

@@ -1227,6 +1227,16 @@ namespace SireIO
                 // Reduce the mass.
                 mass -= delta_mass;
 
+                if (mass.value() < 0)
+                {
+                    qWarning().noquote() << QObject::tr(
+                                                "Hydrogen mass repartitioning with factor %1 has resulted "
+                                                "in a negative mass for atom index %2. Consider using a "
+                                                "smaller repartitioning factor.")
+                                                .arg(factor)
+                                                .arg(idx.value());
+                }
+
                 // Set the new mass.
                 edit_mol = edit_mol.atom(idx).setProperty(mass_prop, mass).molecule();
 

doc/source/changelog.rst

@@ -88,6 +88,20 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Add softcore ``CustomBondForce`` for ring-breaking and ring-making pairs.
 
+* Add ``max_h_mass`` map option (default ``3.5`` g/mol) to control the mass threshold
+  used when identifying hydrogen atoms in the OpenMM conversion layer. Previously this
+  was hardcoded to ``2.5`` g/mol, which caused hydrogen constraints to be skipped when
+  an HMR factor greater than ~2.5 was applied.
+
+* Warn if hydrogen mass repartiting produces a negative heavy atom mass.
+
+* Update merge code to handle ``kartograf`` API changes in version 2.0.
+
+* Allocate ``ghost-14`` slot if either end-state exception scale is nonzero.
+
+* Constrain perturbable bonds with unchanged parameters regardless of atom mass, matching
+  SOMD1, instead of requiring a light (hydrogen) atom to be present.
+
 `2025.4.0 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 

src/sire/_match.py

@@ -91,7 +91,13 @@ def match_atoms(
             elif "KartografAtomMapper" in str(match.__class__):
                 # use Kartograf to get the mapping - convert to RDKit then Kartograf
                 from kartograf.atom_aligner import align_mol_shape
-                from kartograf import KartografAtomMapper, SmallMoleculeComponent
+                from kartograf import KartografAtomMapper
+
+                try:
+                    # kartograf >= 2.0 moved SmallMoleculeComponent to gufe
+                    from gufe import SmallMoleculeComponent
+                except ImportError:
+                    from kartograf import SmallMoleculeComponent
 
                 if not isinstance(match, KartografAtomMapper):
                     raise TypeError("match must be a KartografAtomMapper")

tests/conftest.py

@@ -228,6 +228,16 @@ def solvated_neopentane_methane():
     return sr.load_test_files("neo_meth_solv.bss")
 
 
+@pytest.fixture(scope="session")
+def ghost_14_bug():
+    return sr.load_test_files("ghost_14_bug.bss")
+
+
+@pytest.fixture(scope="session")
+def cyclopropyl_constraint_bug():
+    return sr.load_test_files("ejm_31_jmc_28.bss")
+
+
 @pytest.fixture(scope="session")
 def zero_lj_mols():
     return sr.load_test_files("zero_lj.prm7", "zero_lj.rst7")

tests/convert/test_openmm.py

@@ -489,3 +489,85 @@ def test_openmm_membrane_barostat(ala_mols, openmm_platform):
     assert barostat.getXYMode() == MonteCarloMembraneBarostat.XYIsotropic
     assert barostat.getZMode() == MonteCarloMembraneBarostat.ZFree
     assert barostat.getFrequency() == 50
+
+
+@pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+def test_openmm_ghost_14_bug(ghost_14_bug, openmm_platform):
+    """
+    Regression test for a bug where a pair involving a ghost atom whose
+    intrascale exclusion differs between end states (excluded at lambda=0,
+    a real 1-4 scale at lambda=1, because the connectivity path between the
+    two atoms only exists at lambda=1) never had a slot allocated in the
+    Ghost14BondForce. The construction-time check that decides whether a
+    pair needs a ghost-14 slot only looked at the lambda=0 exception scale,
+    so a pair that's excluded at lambda=0 but not at lambda=1 was silently
+    skipped, then permanently excluded from every nonbonded force for the
+    lifetime of the simulation - even at lambda=1, where it should have a
+    real interaction.
+
+    The test molecule (cyclopropane -> propane, a minimal "ring-breaking"
+    perturbation with no protein or water) reproduces this exactly: atom 2
+    (a ring/chain carbon) and atoms 9, 10, 11 (new hydrogens that only
+    exist on atom 0 at lambda=1, once the propane chain needs an extra H
+    that the cyclopropane ring didn't) are 1-3 excluded at lambda=0 and a
+    real 1-4 pair, (0.8333, 0.5), at lambda=1.
+    """
+    mols = sr.morph.link_to_reference(ghost_14_bug)
+    mol = mols[0]
+
+    d = mol.dynamics(
+        schedule=sr.cas.LambdaSchedule.standard_morph(),
+        platform=openmm_platform,
+        cutoff="none",
+    )
+
+    def get_force(d, name):
+        for force in d.context().getSystem().getForces():
+            if force.getName() == name:
+                return force
+        return None
+
+    # the bug-triggering pairs - atom 2 (C3) with the three new hydrogens
+    # on atom 0 (C1) that only exist at lambda=1
+    pairs = {(2, 9), (2, 10), (2, 11)}
+
+    def get_ghost14_params(d):
+        ghost14 = get_force(d, "Ghost14BondForce")
+        assert ghost14 is not None
+
+        found = {}
+        for i in range(ghost14.getNumBonds()):
+            p1, p2, params = ghost14.getBondParameters(i)
+            key = (min(p1, p2), max(p1, p2))
+            if key in pairs:
+                found[key] = list(params)
+        return found
+
+    # at every lambda, all three pairs must have a slot in Ghost14BondForce
+    # at all (this is what the bug broke - they were silently missing)
+    d.set_lambda(0.0)
+    params0 = get_ghost14_params(d)
+    assert set(params0.keys()) == pairs
+
+    d.set_lambda(0.5)
+    params_mid = get_ghost14_params(d)
+    assert set(params_mid.keys()) == pairs
+
+    d.set_lambda(1.0)
+    params1 = get_ghost14_params(d)
+    assert set(params1.keys()) == pairs
+
+    # at lambda=0 the pair is fully excluded, matching intrascale0 = (0,0)
+    # (four_epsilon is parameter index 2)
+    for params in params0.values():
+        assert params[2] == pytest.approx(0.0, abs=1e-6)
+
+    # the interaction should grow smoothly and monotonically as lambda
+    # goes from 0 to 1, ending up clearly nonzero (matching the real
+    # intrascale1 = (0.8333, 0.5) 1-4 scale)
+    for key in pairs:
+        assert params0[key][2] < params_mid[key][2] < params1[key][2]
+        assert params1[key][2] > 0.1

tests/convert/test_openmm_constraints.py

@@ -454,6 +454,72 @@ def test_auto_constraints(ala_mols, openmm_platform):
                 assert bond in constrained
 
 
+@pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+def test_cyclopropyl_anchor_bond_constraint(
+    cyclopropyl_constraint_bug, openmm_platform
+):
+    """
+    Regression test for an unperturbed anchor bond (to two alternate
+    heavy-atom substituents) not being constrained under
+    "h-bonds-not-heavy-perturbed", unlike SOMD1.
+    """
+    mols = sr.morph.link_to_reference(cyclopropyl_constraint_bug)
+    mol = mols[0]
+
+    # restrict light-atom detection to the maximum end-state mass, as SOMD2
+    # does, rather than Sire's broader (and individually buggy) defaults
+    map = {
+        "check_for_h_by_max_mass": True,
+        "check_for_h_by_mass": False,
+        "check_for_h_by_element": False,
+        "check_for_h_by_ambertype": False,
+    }
+
+    d = mol.dynamics(
+        constraint="h-bonds-not-heavy-perturbed",
+        platform=openmm_platform,
+        cutoff="none",
+        map=map,
+    )
+
+    anchor_idx = sr.mol.AtomIdx(17)
+    substituent_idxs = {sr.mol.AtomIdx(19), sr.mol.AtomIdx(32)}
+
+    def get_anchor_constraints(d):
+        found = {}
+
+        for constraint in d.get_constraints():
+            idx0 = constraint[0][0].index()
+            idx1 = constraint[0][1].index()
+
+            if anchor_idx in (idx0, idx1):
+                other = idx1 if idx0 == anchor_idx else idx0
+
+                if other in substituent_idxs:
+                    found[other] = constraint[1].value()
+
+        return found
+
+    found = get_anchor_constraints(d)
+
+    assert set(found.keys()) == substituent_idxs
+
+    # both anchor bonds should have the same, static length, matching
+    # SOMD1's value for this generic aromatic-C/sp3-C bond
+    for length in found.values():
+        assert length == pytest.approx(1.51234, abs=1e-3)
+
+    # the constraint should remain static (not lambda-dependent), since the
+    # bond's own parameters don't change between end states
+    d.set_lambda(1.0)
+    found_at_1 = get_anchor_constraints(d)
+
+    assert found_at_1 == found
+
+
 @pytest.mark.skipif(
     "openmm" not in sr.convert.supported_formats(),
     reason="openmm support is not available",

tests/morph/test_merge.py

@@ -53,7 +53,11 @@ def test_merge(ose_mols, zan_mols, openmm_platform):
     zan = zan_mols[0]
 
     merged = sr.morph.merge(
-        ose, zan, match=KartografAtomMapper(atom_map_hydrogens=True)
+        ose,
+        zan,
+        match=KartografAtomMapper(
+            atom_map_hydrogens=True, map_hydrogens_on_hydrogens_only=False
+        ),
     )
 
     ose_nrg = ose.dynamics(platform=openmm_platform).current_potential_energy()
@@ -74,7 +78,11 @@ def test_merge(ose_mols, zan_mols, openmm_platform):
     assert extracted_zan_nrg.value() == pytest.approx(zan_nrg.value())
 
     merged = sr.morph.merge(
-        zan, ose, match=KartografAtomMapper(atom_map_hydrogens=True)
+        zan,
+        ose,
+        match=KartografAtomMapper(
+            atom_map_hydrogens=True, map_hydrogens_on_hydrogens_only=False
+        ),
     )
 
     extracted_zan = merged.perturbation().extract_reference()
@@ -152,7 +160,11 @@ def test_merge_neopentane_methane(neopentane_methane, openmm_platform):
     nrg_met = methane.dynamics(platform=openmm_platform).current_potential_energy()
 
     merged = sr.morph.merge(
-        neopentane, methane, match=KartografAtomMapper(atom_map_hydrogens=True)
+        neopentane,
+        methane,
+        match=KartografAtomMapper(
+            atom_map_hydrogens=True, map_hydrogens_on_hydrogens_only=False
+        ),
     )
 
     nrg_merged_0 = merged.dynamics(