From 3da5fd3052b2bd72c8a94e6a7ad3fe20c19fda40 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Mon, 2 Mar 2026 10:05:54 +0000
Subject: [PATCH] Fix AMBER CMAP indexing bug. [closes #400]

---
 corelib/src/libs/SireIO/amberprm.cpp | 18 +++++++----
 doc/source/changelog.rst             |  4 +++
 tests/conftest.py                    |  5 +++
 tests/io/test_ambercmap.py           | 47 ++++++++++++++++++++++++++++
 4 files changed, 68 insertions(+), 6 deletions(-)

diff --git a/corelib/src/libs/SireIO/amberprm.cpp b/corelib/src/libs/SireIO/amberprm.cpp
index 4ca5b5127..031f85e00 100644
--- a/corelib/src/libs/SireIO/amberprm.cpp
+++ b/corelib/src/libs/SireIO/amberprm.cpp
@@ -337,11 +337,14 @@ QVector<QVector<int>> indexCMAPs(const QVector<qint64> &cmaps, const QVector<int
     for (int i = 0; i < cmaps.count(); i += 6)
     {
         // format is atom0-atom1-atom2-atom3-atom4-parameter
-        int mol0 = atom_to_mol_data[cmaps_data[i]];     // DO NOT DIVIDE BY THREE
-        int mol1 = atom_to_mol_data[cmaps_data[i + 1]]; // THIS IS THE RAW ATOM INDEX
-        int mol2 = atom_to_mol_data[cmaps_data[i + 2]];
-        int mol3 = atom_to_mol_data[cmaps_data[i + 3]];
-        int mol4 = atom_to_mol_data[cmaps_data[i + 4]];
+        // CMAP indices are plain 1-based atom numbers (unlike bonds/angles/dihedrals
+        // which are 3×(atom_0based)).  Subtract 1 to convert to the 0-based index
+        // used by atom_to_mol.
+        int mol0 = atom_to_mol_data[cmaps_data[i] - 1];
+        int mol1 = atom_to_mol_data[cmaps_data[i + 1] - 1];
+        int mol2 = atom_to_mol_data[cmaps_data[i + 2] - 1];
+        int mol3 = atom_to_mol_data[cmaps_data[i + 3] - 1];
+        int mol4 = atom_to_mol_data[cmaps_data[i + 4] - 1];
 
         if (mol0 != mol1 or mol0 != mol2 or mol0 != mol3 or mol0 != mol4)
             throw SireIO::parse_error(
@@ -2249,7 +2252,10 @@ std::tuple<QVector<qint64>, QHash<CMAPParameter, qint64>> getCMAPData(const Ambe
     QHash<CMAPParameter, qint64> param_to_idx;
     QVector<qint64> cmap_idxs;
 
-    start_idx /= 3;
+    // NOTE: start_idx is in atom units (0-based count of atoms before this
+    // molecule).  CMAP indices in the prmtop file are plain 1-based atom
+    // numbers (NOT multiplied by 3 as bonds/angles/dihedrals are).
+    // Do NOT divide by 3 here.
 
     const auto info = params.info();
     const auto cmaps = params.cmapFunctions().parameters();
diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst
index 5a231d1a9..6f8d469d7 100644
--- a/doc/source/changelog.rst
+++ b/doc/source/changelog.rst
@@ -17,6 +17,10 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR.
 
+* Fixed a bug in the AMBER prmtop writer where CMAP atom indices were calculated
+  incorrectly for systems containing more than one molecule with CMAP terms (e.g.
+  multi-chain glycoproteins).
+
 `2025.4.0 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 
diff --git a/tests/conftest.py b/tests/conftest.py
index 4a46587c8..49f3b42c8 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -260,6 +260,11 @@ def amber_cmap():
     return sr.load_test_files("amber_cmap.prm7")
 
 
+@pytest.fixture(scope="session")
+def multichain_cmap():
+    return sr.load_test_files("multichain_cmap.prm7", "multichain_cmap.rst7")
+
+
 @pytest.fixture(scope="session")
 def gromacs_cmap():
     return sr.load_test_files("1aki.gro", "1aki.top")
diff --git a/tests/io/test_ambercmap.py b/tests/io/test_ambercmap.py
index e2ac451a6..39a83b72e 100644
--- a/tests/io/test_ambercmap.py
+++ b/tests/io/test_ambercmap.py
@@ -60,6 +60,53 @@ def test_amber_cmap(tmpdir, amber_cmap):
                     assert line1 == line2
 
 
+def test_amber_multichain_cmap(tmpdir, multichain_cmap):
+    """Regression test for multi-chain CMAP write bug.
+
+    When a topology has more than one molecule with CMAP terms, the write path
+    previously applied a spurious '/= 3' to the per-molecule atom offset,
+    producing wrong global atom indices for molecules 2, 3, … .  The re-parse
+    of the generated text then raised:
+      "there is a cmap between more than one different molecule"
+    """
+    mols = multichain_cmap
+
+    # Collect per-molecule CMAP term counts indexed by position in the
+    # molecule list. There must be at least two molecules with CMAP
+    # to trigger the multi-chain bug.
+    cmap_counts = {}
+    for i, mol in enumerate(mols.molecules()):
+        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"
+
+    dir = tmpdir.mkdir("test_amber_multichain_cmap")
+
+    # Write the topology back to file. (This is where the bug used to trigger.)
+    f = sr.save(mols, dir.join("output"), format="prm7")
+
+    # Reead the file back in.
+    mols2 = sr.load(f)
+
+    # Each molecule must still have the same number of CMAP terms after the
+    # 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"
+        count2 = len(mol2.property("cmap").parameters())
+        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")
+
+
 def test_amber_cmap_grotop(tmpdir, amber_cmap):
     """Testing reading and writing from amber to gromacs and back to amber."""
     mols = amber_cmap.clone()