From 5c9bd4b2cdc11da99d8cc26735daa8aa6f0ea4f5 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 7 Sep 2026 09:54:31 +0100 Subject: [PATCH] Fix GROMACS CMAP units. --- corelib/src/libs/SireIO/grotop.cpp | 19 +++++--- doc/source/changelog.rst | 3 ++ tests/io/test_ambercmap.py | 73 ++++++++++++++++++++++++++---- 3 files changed, 80 insertions(+), 15 deletions(-) diff --git a/corelib/src/libs/SireIO/grotop.cpp b/corelib/src/libs/SireIO/grotop.cpp index fd44d6ed9..c333848a8 100644 --- a/corelib/src/libs/SireIO/grotop.cpp +++ b/corelib/src/libs/SireIO/grotop.cpp @@ -457,7 +457,10 @@ static QList cmap_id_to_atomtypes(const QString &cmap_id) return parts.mid(0, 5); } -static QString cmap_to_string(const CMAPParameter &cmap) +/** Serialise the CMAP grid, scaling the values by 'scale'. This is used both for + the in-memory string representation (scale of 1) and for writing to file, + where the grid must be converted from kcal mol-1 to kJ mol-1 */ +static QString cmap_to_string(const CMAPParameter &cmap, double scale = 1.0) { // format is "1 nRows nCols param param param..." QStringList params; @@ -471,7 +474,7 @@ static QString cmap_to_string(const CMAPParameter &cmap) for (int i = 0; i < vals.size(); ++i) { - line.append(QString::number(vals[i], 'f', 8)); + line.append(QString::number(vals[i] * scale, 'f', 8)); if (line.count() == 10) { @@ -3717,10 +3720,11 @@ static QStringList writeCMAPTypes(const QHash &cmap_para const auto &cmap = cmap_params[key]; key = key.replace(";", " "); - // Create the line with the parameters. + // Create the line with the parameters, converting the grid from + // kcal mol-1 to the kJ mol-1 expected by gromacs. lines.append(QString("%1 %2") .arg(key) - .arg(cmap_to_string(cmap))); + .arg(cmap_to_string(cmap, (1 * kcal_per_mol).to(kJ_per_mol)))); } lines.append(""); @@ -7406,10 +7410,13 @@ QStringList GroTop::processDirectives(const QMap &taglocs, const Q continue; } - // we can now read in the cmap values + // we can now read in the cmap values, converting the grid from the + // kJ mol-1 used by gromacs to the kcal mol-1 used by sire QVector cmap_values(nrows * ncols); auto *cmap_values_data = cmap_values.data(); + const double to_kcal_per_mol = (1 * kJ_per_mol).to(kcal_per_mol); + ok = true; for (int i = 0; i < nrows * ncols; ++i) @@ -7426,7 +7433,7 @@ QStringList GroTop::processDirectives(const QMap &taglocs, const Q break; } - cmap_values_data[i] = value; + cmap_values_data[i] = value * to_kcal_per_mol; } if (not ok) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 751cd27c3..e76b10d75 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -99,6 +99,9 @@ organisation on `GitHub `__. rebuilt for every affected 1-4 pair, and the parallel loop that did this read the shared dihedral hash without holding the mutex that guarded the inserts. +* Fixed CMAP grids not being converted between kcal mol-1 and kJ mol-1 when writing + to, and reading from, GROMACS topology files. + `2026.1.0 `__ - June 2026 ----------------------------------------------------------------------------------------- diff --git a/tests/io/test_ambercmap.py b/tests/io/test_ambercmap.py index 39a83b72e..51ae069ce 100644 --- a/tests/io/test_ambercmap.py +++ b/tests/io/test_ambercmap.py @@ -79,9 +79,9 @@ def test_amber_multichain_cmap(tmpdir, multichain_cmap): if mol.has_property("cmap"): cmap_counts[i] = len(mol.property("cmap").parameters()) - assert ( - len(cmap_counts) >= 2 - ), "Expected at least two molecules with CMAP terms in this topology" + assert len(cmap_counts) >= 2, ( + "Expected at least two molecules with CMAP terms in this topology" + ) dir = tmpdir.mkdir("test_amber_multichain_cmap") @@ -95,13 +95,13 @@ def test_amber_multichain_cmap(tmpdir, multichain_cmap): # roundtrip. for i, count in cmap_counts.items(): mol2 = mols2[i] - assert mol2.has_property( - "cmap" - ), f"Molecule at index {i} lost its cmap property after roundtrip" + assert mol2.has_property("cmap"), ( + f"Molecule at index {i} lost its cmap property after roundtrip" + ) count2 = len(mol2.property("cmap").parameters()) - assert ( - count2 == count - ), f"Molecule at index {i}: CMAP count changed from {count} to {count2}" + assert count2 == count, ( + f"Molecule at index {i}: CMAP count changed from {count} to {count2}" + ) # Verify a second write also succeeds without error. sr.save(mols2, dir.join("output2"), format="prm7") @@ -169,3 +169,58 @@ def test_amber_cmap_grotop(tmpdir, amber_cmap): found = True assert found + + +def test_amber_cmap_grotop_units(tmpdir, amber_cmap): + """Testing that CMAP grids are converted to kJ mol-1 when written to gromacs.""" + mols = amber_cmap.clone() + + dir = tmpdir.mkdir("test_amber_cmap_grotop_units") + + # Save to a temporary file in GroTop format. + f = sr.save(mols, dir.join("output"), format="GroTop")[0] + + # Read the values from the [ cmaptypes ] section of the file. + file_values = [] + in_cmaptypes = False + + for line in open(f): + line = line.strip() + + if line.startswith("["): + in_cmaptypes = line.replace(" ", "") == "[cmaptypes]" + continue + + if not in_cmaptypes or not line or line.startswith(";"): + continue + + # strip the line continuation, then drop the leading + # "atm0 atm1 atm2 atm3 atm4 func nrows ncols" of a header line + parts = line.rstrip("\\").split() + + try: + float(parts[0]) + except ValueError: + parts = parts[8:] + + file_values += [float(x) for x in parts] + + # Gather the grid values held in memory, which are in kcal mol-1. Only the + # unique grids are written to the file, so deduplicate to match. + mol_values = [] + seen = set() + + for cmap in mols[0].property("cmap").parameters(): + values = tuple(cmap.parameter().values()) + + if values not in seen: + seen.add(values) + mol_values += list(values) + + assert len(file_values) == len(mol_values) + + # The written values must be the in-memory values converted to kJ mol-1. + kcal_to_kj = sr.u("1 kcal mol-1").to("kJ mol-1") + + for written, expected in zip(sorted(file_values), sorted(mol_values)): + assert written == pytest.approx(expected * kcal_to_kj, rel=1e-5)