From a19745cf2725f75fa4b4a536463a9a35588915b9 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 16 Jun 2026 16:00:00 +0100 Subject: [PATCH 1/8] Update max hydrogen mass to 4.0 g/mol and expose map option to set. --- doc/source/changelog.rst | 5 +++++ wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 21 ++++++++++++------- 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 053f4b6a8..3dc9a59f3 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -88,6 +88,11 @@ organisation on `GitHub `__. * Add softcore ``CustomBondForce`` for ring-breaking and ring-making pairs. +* Add ``max_h_mass`` map option (default ``4.0`` 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. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index aa19a231f..2e134cc6e 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -614,6 +614,13 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, check_for_h_by_max_mass = map["check_for_h_by_max_mass"].value().asABoolean(); } + double max_h_mass = 4.0; + + if (map.specified("max_h_mass")) + { + max_h_mass = map["max_h_mass"].value().asADouble(); + } + bool check_for_h_by_element = true; if (map.specified("check_for_h_by_element")) @@ -650,8 +657,8 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, const double mass0 = params_masses.at(cgatomidx).to(SireUnits::g_per_mol); const double mass1 = params1_masses.at(cgatomidx).to(SireUnits::g_per_mol); - const bool mass0_is_light = (mass0 >= 1) and (mass0 < 2.5); - const bool mass1_is_light = (mass1 >= 1) and (mass1 < 2.5); + const bool mass0_is_light = (mass0 >= 1) and (mass0 < max_h_mass); + const bool mass1_is_light = (mass1 >= 1) and (mass1 < max_h_mass); double mass = std::max(mass0, mass1); @@ -665,9 +672,9 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, // this must be a ghost in both end states? light_atoms.insert(i); } - else if (check_for_h_by_max_mass and mass < 2.5) + else if (check_for_h_by_max_mass and mass < max_h_mass) { - // the maximum mass is less than 2.5, so this is a H + // the maximum mass is less than max_h_mass, so this is a H light_atoms.insert(i); } else if (check_for_h_by_mass and (mass0_is_light or mass1_is_light)) @@ -725,12 +732,12 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, // this must be a ghost light_atoms.insert(i); } - else if (check_for_h_by_max_mass and mass < 2.5) + else if (check_for_h_by_max_mass and mass < max_h_mass) { - // the maximum mass is less than 2.5, so this is a H + // the maximum mass is less than max_h_mass, so this is a H light_atoms.insert(i); } - else if (check_for_h_by_mass and (mass >= 1 and mass < 2.5)) + else if (check_for_h_by_mass and (mass >= 1 and mass < max_h_mass)) { // one of the atoms is H or He light_atoms.insert(i); From cb67f717683828e79a7b44a8b8b85a6adb37ef5b Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 16 Jun 2026 16:19:23 +0100 Subject: [PATCH 2/8] Update threshold to handle SOMD2 upper bound. --- doc/source/changelog.rst | 2 +- wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 3dc9a59f3..c06071e4f 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -88,7 +88,7 @@ organisation on `GitHub `__. * Add softcore ``CustomBondForce`` for ring-breaking and ring-making pairs. -* Add ``max_h_mass`` map option (default ``4.0`` g/mol) to control the mass threshold +* Add ``max_h_mass`` map option (default ``4.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. diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 2e134cc6e..d9ec717bb 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -614,7 +614,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, check_for_h_by_max_mass = map["check_for_h_by_max_mass"].value().asABoolean(); } - double max_h_mass = 4.0; + double max_h_mass = 4.5; if (map.specified("max_h_mass")) { From 776b5c4cf4d8c766b8b76fedebb0e079fc4c77aa Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 16 Jun 2026 16:37:17 +0100 Subject: [PATCH 3/8] Reduce threshold to handle minimal realistic value. --- doc/source/changelog.rst | 2 +- wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index c06071e4f..9eefa96a0 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -88,7 +88,7 @@ organisation on `GitHub `__. * Add softcore ``CustomBondForce`` for ring-breaking and ring-making pairs. -* Add ``max_h_mass`` map option (default ``4.5`` g/mol) to control the mass threshold +* 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. diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index d9ec717bb..4c38a72c8 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -614,7 +614,7 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, check_for_h_by_max_mass = map["check_for_h_by_max_mass"].value().asABoolean(); } - double max_h_mass = 4.5; + double max_h_mass = 3.5; if (map.specified("max_h_mass")) { From 80cf65bf4f359930681ae4092b919a6075da16aa Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 16 Jun 2026 16:45:59 +0100 Subject: [PATCH 4/8] Warn if HMR produces a negative heavy atom mass. --- corelib/src/libs/SireIO/biosimspace.cpp | 10 ++++++++++ doc/source/changelog.rst | 2 ++ 2 files changed, 12 insertions(+) diff --git a/corelib/src/libs/SireIO/biosimspace.cpp b/corelib/src/libs/SireIO/biosimspace.cpp index 7dbc040f8..fa293afac 100644 --- a/corelib/src/libs/SireIO/biosimspace.cpp +++ b/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(); diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 9eefa96a0..98e5a893f 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -93,6 +93,8 @@ organisation on `GitHub `__. 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. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- From 0df37ae28616d8b18238f205991ad169726521e2 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 17 Jun 2026 09:00:40 +0100 Subject: [PATCH 5/8] Handle kartograf API change. --- src/sire/_match.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/sire/_match.py b/src/sire/_match.py index de83f9003..7ecffb8d4 100644 --- a/src/sire/_match.py +++ b/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") From 87af3857a7b987a535313e68f1a064b3d8801e70 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 18 Jun 2026 09:52:53 +0100 Subject: [PATCH 6/8] Pin map_hydrogens_on_hydrogens_only to keep mapping kartograf-version-stable --- doc/source/changelog.rst | 2 ++ tests/morph/test_merge.py | 18 +++++++++++++++--- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 98e5a893f..449d8f13f 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -95,6 +95,8 @@ organisation on `GitHub `__. * Warn if hydrogen mass repartiting produces a negative heavy atom mass. +* Update merge code to handle ``kartograf`` API changes in version 2.0. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/tests/morph/test_merge.py b/tests/morph/test_merge.py index 25f9550b6..cc8e130b0 100644 --- a/tests/morph/test_merge.py +++ b/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( From dde13763063d717a3363e1deadb31217f9ce5058 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 18 Jun 2026 13:21:14 +0100 Subject: [PATCH 7/8] Allocate ghost-14 slot if either end-state exception scale is nonzero. --- doc/source/changelog.rst | 2 + tests/conftest.py | 5 ++ tests/convert/test_openmm.py | 82 +++++++++++++++++++ .../SireOpenMM/sire_to_openmm_system.cpp | 20 ++++- 4 files changed, 108 insertions(+), 1 deletion(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 449d8f13f..5cb65c1c4 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -97,6 +97,8 @@ organisation on `GitHub `__. * Update merge code to handle ``kartograf`` API changes in version 2.0. +* Allocate ``ghost-14`` slot if either end-state exception scale is nonzero. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/tests/conftest.py b/tests/conftest.py index 49f3b42c8..2b8af0da0 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -228,6 +228,11 @@ 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 zero_lj_mols(): return sr.load_test_files("zero_lj.prm7", "zero_lj.rst7") diff --git a/tests/convert/test_openmm.py b/tests/convert/test_openmm.py index d2efb7c15..c5335ccea 100644 --- a/tests/convert/test_openmm.py +++ b/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 diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 0b5032b2d..7e1f498b5 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -2641,7 +2641,25 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, auto to_from_ghost = (from_ghost_idxs.contains(atom0) and to_ghost_idxs.contains(atom1)) or (from_ghost_idxs.contains(atom1) and to_ghost_idxs.contains(atom0)); - if (not to_from_ghost and (coul_14_scl != 0 or lj_14_scl != 0)) + // Whether a ghost-14 slot is needed is decided here, once, + // for the lifetime of the OpenMM::CustomBondForce - bonds + // can't be added to it later, only updated. coul_14_scl/ + // lj_14_scl are this pair's state0 scale only; also check + // the aligned state1 entry at the same index, since a pair + // can be fully excluded in one end state's connectivity + // and fully (or partially) unmasked in the other - e.g. + // a 1-3 pair that becomes unmasked once a ring-breaking + // bond removes the path between them. Missing this means + // the pair is permanently excluded from ghost/non-ghost + // and ghost/ghost too (their exclusion lists are built + // from this same state0-only loop below), with no way to + // recover the real interaction at the end where it exists. + const auto &pert_param = mol.perturbed->exception_params[j]; + const auto coul_14_scl1 = boost::get<2>(pert_param); + const auto lj_14_scl1 = boost::get<3>(pert_param); + + if (not to_from_ghost and (coul_14_scl != 0 or lj_14_scl != 0 or + coul_14_scl1 != 0 or lj_14_scl1 != 0)) { // this is a 1-4 interaction that should be added // to the ghost-14 forcefield From feba446ef4f59a6f27780e745ebd0d1f434eec45 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 18 Jun 2026 16:11:37 +0100 Subject: [PATCH 8/8] Constrain unperturbed heavy-atom bonds regardless of mass, matching SOMD1 --- doc/source/changelog.rst | 3 + tests/conftest.py | 5 + tests/convert/test_openmm_constraints.py | 66 ++++++++++++++ wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 91 ++++++++++++------- 4 files changed, 133 insertions(+), 32 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 5cb65c1c4..dd764d9ed 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -99,6 +99,9 @@ organisation on `GitHub `__. * 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 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/tests/conftest.py b/tests/conftest.py index 2b8af0da0..816684d5e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -233,6 +233,11 @@ 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") diff --git a/tests/convert/test_openmm_constraints.py b/tests/convert/test_openmm_constraints.py index ef98bb276..e856864e2 100644 --- a/tests/convert/test_openmm_constraints.py +++ b/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", diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 4c38a72c8..78b53999f 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -895,6 +895,19 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, auto bonds = params.bonds(); + // Whether an atom's amber type changes between end states - used below + // to identify bonds that are part of the alchemical junction. Amber + // type, not charge, since charges can be redistributed across the whole + // ligand without affecting bonded parameters. + const auto params1_ambertypes = params1.amberTypes(); + + auto is_perturbing_atom = [&](int atom_idx) -> bool + { + const auto &cgatomidx = idx_to_cgatomidx_data[atom_idx]; + + return ambertypes.at(cgatomidx) != params1_ambertypes.at(cgatomidx); + }; + if (is_perturbable) { // add in any bonds that exist only in the other state @@ -952,6 +965,34 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, bool bond_is_not_constrained = true; bool should_constrain_bond = false; + // Whether the bond's own k/r0 change between end states. + bool bond_is_perturbing = false; + + // Whether the bond is part of the alchemical junction (touches a + // ghost or transmuting atom). SOMD1 constrains an unperturbing bond + // unconditionally only if it's part of the perturbation; ordinary + // unperturbed bonds elsewhere in the molecule are not affected. + bool bond_is_ghost_junction = false; + + if (is_perturbable) + { + const auto bondparam1 = params1.bonds().value(it.key()).first; + + double k_1 = bondparam1.k() * bond_k_to_openmm; + r0_1 = bondparam1.r0() * bond_r0_to_openmm; + + if (r0_1 == 0) + { + // we cannot shrink the bond to 0 - this must be + // a bond that is disappearing - we should simply + // keep it the same length + r0_1 = r0; + } + + bond_is_perturbing = (std::abs(k_1 - k) > 1e-3 or std::abs(r0_1 - r0) > 1e-3); + bond_is_ghost_junction = is_perturbing_atom(atom0) or is_perturbing_atom(atom1); + } + if (this_constraint_type == CONSTRAIN_AUTO_BONDS) { // constrain the bond if its predicted vibrational frequency is less than @@ -973,45 +1014,31 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, should_constrain_bond = vibrational_period < auto_constraints_factor * timestep_in_fs; } else if ((not has_massless_atom) and ((this_constraint_type & CONSTRAIN_BONDS) or - (has_light_atom and (this_constraint_type & CONSTRAIN_HBONDS)))) + (has_light_atom and (this_constraint_type & CONSTRAIN_HBONDS)) or + (is_perturbable and not bond_is_perturbing and + bond_is_ghost_junction and + (this_constraint_type & CONSTRAIN_HBONDS)))) { should_constrain_bond = true; - if (is_perturbable) + if (is_perturbable and bond_is_perturbing) { - // we need to see if this bond is being perturbed - const auto bondparam1 = params1.bonds().value(it.key()).first; - - double k_1 = bondparam1.k() * bond_k_to_openmm; - r0_1 = bondparam1.r0() * bond_r0_to_openmm; - - if (r0_1 == 0) + // we need to check against the "NOT_PERTURBED"-style constraints + if (this_constraint_type & CONSTRAIN_NOT_PERTURBED) { - // we cannot shrink the bond to 0 - this must be - // a bond that is disappearing - we should simply - // keep it the same length - r0_1 = r0; + // DO NOT CONSTRAIN any perturbing bonds + should_constrain_bond = false; } - - if (std::abs(k_1 - k) > 1e-3 or std::abs(r0_1 - r0) > 1e-3) + else if (this_constraint_type & CONSTRAIN_NOT_HEAVY_PERTURBED) { - // we need to check against the "NOT_PERTURBED"-style constraints - if (this_constraint_type & CONSTRAIN_NOT_PERTURBED) - { - // DO NOT CONSTRAIN any perturbing bonds - should_constrain_bond = false; - } - else if (this_constraint_type & CONSTRAIN_NOT_HEAVY_PERTURBED) - { - // DO NOT CONSTRAIN any perturbing bonds that DON'T contain hydrogen - should_constrain_bond = has_light_atom; - } - else - { - // DO CONSTRAIN any perturbing bonds - we only don't constrain - // perturbing bonds if we have one of the "NOT_PERTURBED" constraints - should_constrain_bond = true; - } + // DO NOT CONSTRAIN any perturbing bonds that DON'T contain hydrogen + should_constrain_bond = has_light_atom; + } + else + { + // DO CONSTRAIN any perturbing bonds - we only don't constrain + // perturbing bonds if we have one of the "NOT_PERTURBED" constraints + should_constrain_bond = true; } } }