Backport fix from PR #476

corelib/src/libs/SireMM/amberparams.cpp

@@ -66,6 +66,8 @@
 
 #include <QDebug>
 
+#include <mutex>
+
 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")),

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