Repository: openforcefield/openff-interchange

Issue:
GROMACS writer assumes unique molecules are sorted
**Description**

Working with a data source which is not sorted by chemically uniques molecule, i.e.

```

ATOM    194 Cl-  Cl- A   2      37.846  20.981  12.009  1.00  0.00          Cl  
ATOM    195 Na+  Na+ A   3      12.107   6.695   0.364  1.00  0.00          Na  
ATOM    196 Na+  Na+ A   4      26.880  32.136   0.363  1.00  0.00          Na  
ATOM    197 Na+  Na+ A   5      10.502  16.107  26.742  1.00  0.00          Na  
ATOM    198 Na+  Na+ A   6      14.887  38.236  20.253  1.00  0.00          Na  
ATOM    199 Na+  Na+ A   7      17.994  22.653  20.731  1.00  0.00          Na  
ATOM    200 Na+  Na+ A   8      26.623  28.181  29.279  1.00  0.00          Na  
ATOM    201 Cl-  Cl- A   9      30.356  36.332  33.324  1.00  0.00          Cl  
ATOM    202 Cl-  Cl- A  10      31.594  29.048  28.622  1.00  0.00          Cl  
ATOM    203 Cl-  Cl- A  11      15.932  22.505  11.696  1.00  0.00          Cl  
ATOM    204 Cl-  Cl- A  12      21.690   5.088   3.933  1.00  0.00          Cl  
ATOM    205 Cl-  Cl- A  13      22.319  37.825  21.847  1.00  0.00          Cl  
ATOM    206 Cl-  Cl- A  14      21.755  40.314  18.121  1.00  0.00          Cl  
```

causes issues in GROMACS files. The exporter finds each instance of a chemically uniques molecule and counts the number of them that are found, but assumes they are sorted. Molecules from this snippet are effectively re-ordered in the written GROMACS files, i.e.

```
    3Cl-    Cl-  194   3.785   2.098   1.201
    4Cl-    Cl-  195   1.211   0.670   0.036
    5Cl-    Cl-  196   2.688   3.214   0.036
    6Cl-    Cl-  197   1.050   1.611   2.674
    7Cl-    Cl-  198   1.489   3.824   2.025
    8Cl-    Cl-  199   1.799   2.265   2.073
    9Cl-    Cl-  200   2.662   2.818   2.928
    4Na+    Na+  201   3.036   3.633   3.332
    5Na+    Na+  202   3.159   2.905   2.862
    6Na+    Na+  203   1.593   2.250   1.170
    7Na+    Na+  204   2.169   0.509   0.393
    8Na+    Na+  205   2.232   3.783   2.185
    9Na+    Na+  206   2.176   4.031   1.812
```

```

[ molecules ]
;name	number
; ...
MOL2	7   ; Cl- molecule
MOL3	6   ; Na+ molecule
; ...
```

This is self-consistent but bad - atoms should not be re-ordered at export time.

**Reproduction**

This can be more easily reproduced with a simple water-ammonia-water topology. The molecules should be ordered identically in the GROMACS files and input topology, but that is not what happens.


```python
from openff.interchange._tests import MoleculeWithConformer
from openff.toolkit import ForceField, Topology


ForceField("openff-2.2.1.offxml").create_interchange(
    Topology.from_molecules(
        [
            MoleculeWithConformer.from_smiles("O"),
            MoleculeWithConformer.from_smiles("N"),
            MoleculeWithConformer.from_smiles("O"),
        ]
    )
).to_top("ordering.top")

!tail ordering.top
```
**Output**

```

[ system ]
;name
FOO

[ molecules ]
;name	number
MOL0	2
MOL1	1
```

**Software versions**

- Which operating system and version are you using?
- How did you install Interchange?
- What is the output of running `conda list`?


Discussion:

Patch:
diff --git a/docs/releasehistory.md b/docs/releasehistory.md
index 3e62319f..d51da967 100644
--- a/docs/releasehistory.md
+++ b/docs/releasehistory.md
@@ -16,6 +16,11 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 ### Behavior changes
 
 * #1137 Internally use `openmm.LangevinMiddleIntegrator`, including in the OpenMM driver.
+* #1186 `GROMACSSystem.molecules` is changed from a dictionary to a list of tuples.
+
+### Bug fixes
+
+* #1186 Fixes an issue in which GROMACS topology files were mangled when molecules were not grouped within a topology.
 
 ## 0.4.2 - 2025-02-26
 
diff --git a/openff/interchange/interop/gromacs/_import/_import.py b/openff/interchange/interop/gromacs/_import/_import.py
index a49a7c4b..efb8cc1e 100644
--- a/openff/interchange/interop/gromacs/_import/_import.py
+++ b/openff/interchange/interop/gromacs/_import/_import.py
@@ -115,7 +115,7 @@ def from_files(top_file, gro_file, cls=GROMACSSystem) -> GROMACSSystem:
         elif current_directive == "molecules":
             molecule_name, number_of_copies = _process_molecule(line)
 
-            system.molecules[molecule_name] = number_of_copies
+            system.molecules.append((molecule_name, number_of_copies))
 
         else:
             raise ValueError(f"Invalid directive {current_directive}")
diff --git a/openff/interchange/interop/gromacs/_import/_topology.py b/openff/interchange/interop/gromacs/_import/_topology.py
index 6c60aeb0..711811d5 100644
--- a/openff/interchange/interop/gromacs/_import/_topology.py
+++ b/openff/interchange/interop/gromacs/_import/_topology.py
@@ -9,7 +9,9 @@ def _create_topology_from_system(system: GROMACSSystem) -> Topology:
 
     topology = Topology()
 
-    for molecule_name, gromacs_molecule in system.molecule_types.items():
+    for molecule_name, n_copies in system.molecules:
+        gromacs_molecule = system.molecule_types[molecule_name]
+
         molecule = _SimpleMolecule()
         molecule.name = molecule_name
 
@@ -38,7 +40,7 @@ def _create_topology_from_system(system: GROMACSSystem) -> Topology:
             molecule.add_bond(molecule.atom(0), molecule.atom(1))
             molecule.add_bond(molecule.atom(0), molecule.atom(2))
 
-        for _ in range(system.molecules[molecule_name]):
+        for _ in range(n_copies):
             topology.add_molecule(molecule)
 
         # TODO: Add residue metadata
diff --git a/openff/interchange/interop/gromacs/_interchange.py b/openff/interchange/interop/gromacs/_interchange.py
index 0ea8f3ec..3931630c 100644
--- a/openff/interchange/interop/gromacs/_interchange.py
+++ b/openff/interchange/interop/gromacs/_interchange.py
@@ -81,8 +81,10 @@ def to_interchange(
 
     molecule_start_index = 0
 
-    for molecule_name, molecule_type in system.molecule_types.items():
-        for _ in range(system.molecules[molecule_name]):
+    for molecule_name, n_copies in system.molecules:
+        molecule_type = system.molecule_types[molecule_name]
+
+        for _ in range(n_copies):
             for atom in molecule_type.atoms:
                 topology_atom_index = molecule_start_index + atom.index - 1
                 topology_key = TopologyKey(
@@ -324,10 +326,10 @@ def _convert_topology(
 
     topology = Topology()
 
-    for molecule_name, molecule_type in system.molecule_types.items():
-        graph = networkx.Graph()
+    for molecule_name, n_molecules in system.molecules:
+        molecule_type = system.molecule_types[molecule_name]
 
-        n_molecules = system.molecules[molecule_name]
+        graph = networkx.Graph()
 
         for atom in molecule_type.atoms:
             graph.add_node(
diff --git a/openff/interchange/interop/gromacs/export/_export.py b/openff/interchange/interop/gromacs/export/_export.py
index 130b4a51..d35b5e5a 100644
--- a/openff/interchange/interop/gromacs/export/_export.py
+++ b/openff/interchange/interop/gromacs/export/_export.py
@@ -432,7 +432,7 @@ class GROMACSWriter(_BaseModel):
         top.write("[ molecules ]\n")
         top.write(";name\tnumber\n")
 
-        for molecule_name, n_molecules in self.system.molecules.items():
+        for molecule_name, n_molecules in self.system.molecules:
             top.write(f"{molecule_name.replace(' ', '_')}\t{n_molecules}\n")
 
         top.write("\n")
@@ -449,8 +449,8 @@ class GROMACSWriter(_BaseModel):
             )
 
         n_particles = sum(
-            len(molecule_type.atoms) * self.system.molecules[molecule_name]
-            for molecule_name, molecule_type in self.system.molecule_types.items()
+            len(self.system.molecule_types[molecule_name].atoms) * n_copies
+            for molecule_name, n_copies in self.system.molecules
         )
 
         assert n_particles == self.system.positions.shape[0], (
@@ -465,8 +465,8 @@ class GROMACSWriter(_BaseModel):
         gro.write(f"{n_particles}\n")
 
         count = 0
-        for molecule_name, molecule in self.system.molecule_types.items():
-            n_copies = self.system.molecules[molecule_name]
+        for molecule_name, n_copies in self.system.molecules:
+            molecule = self.system.molecule_types[molecule_name]
 
             for copy_index in range(n_copies):
                 for atom in molecule.atoms:
diff --git a/openff/interchange/interop/gromacs/models/models.py b/openff/interchange/interop/gromacs/models/models.py
index 8b9337dd..c087ae5f 100644
--- a/openff/interchange/interop/gromacs/models/models.py
+++ b/openff/interchange/interop/gromacs/models/models.py
@@ -341,9 +341,9 @@ class GROMACSSystem(_BaseModel):
         dict(),
         description="Molecule types, keyed by name.",
     )
-    molecules: dict[str, int] = Field(
-        dict(),
-        description="The number of each molecule type in the system, keyed by the name of each molecule.",
+    molecules: list[tuple[str, int]] = Field(
+        list(),
+        description="The number of each molecule type in the system, ordered topologically.",
     )
 
     @classmethod
@@ -395,18 +395,17 @@ class GROMACSSystem(_BaseModel):
                 f"The molecule type {molecule_name} is not present in this system.",
             )
 
-        if n_copies > self.molecules[molecule_name]:
-            raise ValueError(
-                f"Cannot remove {n_copies} copies of {molecule_name} from this system "
-                f"because only {self.molecules[molecule_name]} are present.",
-            )
-
-        if n_copies != 1 or self.molecules[molecule_name] != 1:
+        if n_copies != 1 or sum(_n_copies for name, _n_copies in self.molecules if name == molecule_name) != 1:
             raise NotImplementedError()
 
-        molecule_names = [*self.molecules.keys()]
-        molecules_before = molecule_names[: molecule_names.index(molecule_name)]
-        n_atoms_before = sum(len(self.molecule_types[name].atoms) * self.molecules[name] for name in molecules_before)
+        n_atoms_before = sum(
+            len(self.molecule_types[molecule_name].atoms) * n_copies
+            for molecule_name, n_copies in self.molecules[
+                : self.molecules.index(
+                    (molecule_name, 1),
+                )
+            ]
+        )
 
         if self.positions is not None:
             row_indices_to_delete = [
@@ -423,10 +422,7 @@ class GROMACSSystem(_BaseModel):
             )
 
         self.molecule_types.pop(molecule_name)
-        self.molecules[molecule_name] -= n_copies
-
-        if self.molecules[molecule_name] == 0:
-            self.molecules.pop(molecule_name)
+        self.molecules = [(name, n_copies) for name, n_copies in self.molecules if name != molecule_name]
 
     def add_molecule_type(self, molecule: GROMACSMolecule, n_copies: int):
         """Add a molecule type to the system."""
@@ -449,4 +445,4 @@ class GROMACSSystem(_BaseModel):
             self.atom_types[atom_type_name] = atom_type
 
         self.molecule_types[molecule.name] = molecule
-        self.molecules[molecule.name] = n_copies
+        self.molecules.append((molecule.name, n_copies))
diff --git a/openff/interchange/smirnoff/_gromacs.py b/openff/interchange/smirnoff/_gromacs.py
index 0eab0548..9e68cbc4 100644
--- a/openff/interchange/smirnoff/_gromacs.py
+++ b/openff/interchange/smirnoff/_gromacs.py
@@ -189,6 +189,8 @@ def _convert(
         else:
             raise RuntimeError()
 
+    _ordered_molecules = [None] * interchange.topology.n_molecules
+
     for unique_molecule_index in unique_molecule_map:
         unique_molecule = interchange.topology.molecule(unique_molecule_index)
 
@@ -273,9 +275,8 @@ def _convert(
 
         system.molecule_types[unique_molecule.name] = molecule
 
-        system.molecules[unique_molecule.name] = len(
-            interchange.topology.identical_molecule_groups[unique_molecule_index],
-        )
+        for duplicate_molecule_index, _ in unique_molecule_map[unique_molecule_index]:
+            _ordered_molecules[duplicate_molecule_index] = unique_molecule.name
 
     if "VirtualSites" in interchange.collections:
         # TODO: Some say to skip this if the user only wants a topology file?
@@ -291,6 +292,8 @@ def _convert(
     else:
         system.positions = interchange.positions
 
+    system.molecules = [(key, len(list(group))) for key, group in itertools.groupby(_ordered_molecules)]
+
     system.box = interchange.box
 
     return system


Test patch:
diff --git a/openff/interchange/_tests/unit_tests/interop/gromacs/test_interchange.py b/openff/interchange/_tests/unit_tests/interop/gromacs/test_interchange.py
index a5bc361d..56528753 100644
--- a/openff/interchange/_tests/unit_tests/interop/gromacs/test_interchange.py
+++ b/openff/interchange/_tests/unit_tests/interop/gromacs/test_interchange.py
@@ -3,6 +3,7 @@ import pytest
 from openff.toolkit import Molecule, Topology, unit
 
 from openff.interchange import Interchange
+from openff.interchange._tests import MoleculeWithConformer
 from openff.interchange.interop.gromacs._interchange import (
     _convert_topology,
     to_interchange,
@@ -95,3 +96,26 @@ class TestConvertTopology:
         assert converted.n_bonds == 4
 
         assert [atom.atomic_number for atom in converted.atoms] == [8, 1, 1, 8, 1, 1]
+
+    def test_nongrouped_molecules(self, sage):
+        """Test that a system containing un-grouped identical molecules is handled correctly."""
+
+        pytest.importorskip("openmm")
+
+        import openmm.app
+
+        sage.create_interchange(
+            Topology.from_molecules(
+                [
+                    MoleculeWithConformer.from_smiles("O"),
+                    MoleculeWithConformer.from_smiles("N"),
+                    MoleculeWithConformer.from_smiles("O"),
+                ],
+            ),
+        ).to_top("ungrouped.top")
+
+        assert [
+            "MOL0",
+            "MOL1",
+            "MOL0",
+        ] == [residue.name for residue in openmm.app.GromacsTopFile("ungrouped.top").topology.residues()]