From 99b90c6291116694237e54667c09d3006eca4838 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 7 Sep 2026 09:06:23 +0100 Subject: [PATCH] Backport fix from PR #476. [ci skip] --- corelib/src/libs/SireMM/amberparams.cpp | 25 +++++------ tests/convert/test_openmm_zero_k_torsions.py | 47 ++++++++++++++++++++ 2 files changed, 59 insertions(+), 13 deletions(-) create mode 100644 tests/convert/test_openmm_zero_k_torsions.py diff --git a/corelib/src/libs/SireMM/amberparams.cpp b/corelib/src/libs/SireMM/amberparams.cpp index fdfedb950..402886b89 100644 --- a/corelib/src/libs/SireMM/amberparams.cpp +++ b/corelib/src/libs/SireMM/amberparams.cpp @@ -66,6 +66,8 @@ #include +#include + using namespace SireMol; using namespace SireCAS; using namespace SireMM; @@ -1453,8 +1455,10 @@ QStringList AmberParams::validateAndFix() if (not exc_atoms.isEmpty()) { + // the connectivity implied by the bonds - only built if a worker below + // actually needs it, as many molecules (e.g. water) have no 1-4 pairs Connectivity conn; - bool has_connectivity = false; + std::once_flag conn_flag; auto new_dihedrals = amber_dihedrals; auto new_nb14s = amber_nb14s; @@ -1504,17 +1508,9 @@ QStringList AmberParams::validateAndFix() const auto atm3 = molinfo.atomIdx(CGAtomIdx(CGIdx(jcg), Index(j))); - if (not has_connectivity) - { - // have to use the connectivity that is implied by the - // bonds - QMutexLocker lkr(&mutex); - if (not has_connectivity) - { - conn = this->connectivity(); - has_connectivity = true; - } - } + std::call_once(conn_flag, + [&]() + { conn = this->connectivity(); }); // find the shortest bonded paths between these two atoms const auto paths = conn.findPaths(atm0, atm3, 4); @@ -1577,6 +1573,10 @@ QStringList AmberParams::validateAndFix() auto dih = this->convert( DihedralID(path[0], path[1], path[2], path[3])); + // the check and the updates below must be a single + // atomic operation, as workers share these containers + QMutexLocker lkr(&mutex); + // skip if we already have this dihedral if (new_dihedrals.contains(dih)) continue; @@ -1603,7 +1603,6 @@ QStringList AmberParams::validateAndFix() // create a null dihedral parameter and add this to the // set - QMutexLocker lkr(&mutex); new_dihedrals.insert( dih, qMakePair(AmberDihedral(Expression(0), Symbol("phi")), diff --git a/tests/convert/test_openmm_zero_k_torsions.py b/tests/convert/test_openmm_zero_k_torsions.py new file mode 100644 index 000000000..df6e98835 --- /dev/null +++ b/tests/convert/test_openmm_zero_k_torsions.py @@ -0,0 +1,47 @@ +import sire as sr +import pytest + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_openmm_zero_k_torsions(openmm_platform): + """ + GROMACS topologies can contain torsions with a zero force constant, which + the GroTop reader drops. In AMBER-derived topologies those torsions are what + carry the 1-4 scaling, so AmberParams::validateAndFix() has to rebuild a null + dihedral for every 1-4 pair that no longer has one. Check that every 1-4 pair + in the file survives that repair, and that repeating the conversion is stable + (the repair loop is run in parallel and used to race). + """ + import openmm + + mols = sr.load_test_files("zero_k_torsions.gro", "zero_k_torsions.top") + + # the number of entries in the [ pairs ] section of the topology + num_pairs = 8575 + + for _ in range(5): + omm = sr.convert.to(mols, "openmm", map={"platform": openmm_platform}) + + nonbonded = None + + for force in omm.getSystem().getForces(): + if isinstance(force, openmm.NonbondedForce): + nonbonded = force + break + + assert nonbonded is not None + + # count the exceptions that are real 1-4 pairs, i.e. those that are + # scaled rather than fully excluded + num_14 = 0 + + for i in range(nonbonded.getNumExceptions()): + _, _, chg, _, lj = nonbonded.getExceptionParameters(i) + + if chg._value != 0.0 or lj._value != 0.0: + num_14 += 1 + + assert num_14 == num_pairs