Improve GROMACS parsing

openff/interchange/_tests/data/README.md

@@ -65,10 +65,19 @@ struct.save('packed-argon.pdb')
 
 - [Source](https://github.com/samplchallenges/SAMPL6/blob/c661d3985af7fa0ba8c64a1774cfb2363cd31bda/host_guest/CB8AndGuests/CB8.mol2)
 
-`complex.top`
+`gromacs/complex.top`
 
 - [Source](https://raw.githubusercontent.com/samplchallenges/SAMPL6/master/host_guest/SAMPLing/CB8-G3-0/GROMACS/complex.top)
 
-`complex.gro`
+`gromacs/complex.gro`
 
 - [Source](https://raw.githubusercontent.com/samplchallenges/SAMPL6/master/host_guest/SAMPLing/CB8-G3-0/GROMACS/complex.gro)
+
+`gromacs/water-dimer.top`
+`gromacs/water-dimer.gro`
+
+```python
+ForceField("openff-2.0.0.offxml").create_interchange(
+    Topology.from_pdb("water-dimer.pdb")
+).to_gromacs(prefix="water-dimer")
+```

openff/interchange/interop/gromacs/_import/_import.py

@@ -367,7 +367,7 @@ def _process_dihedral(
 
     elif dihedral_function == 4:
         return PeriodicImproperDihedral(
-            atom1=atom1,
+            atom1=atom1,  # central atom?
             atom2=atom2,
             atom3=atom3,
             atom4=atom4,

openff/interchange/interop/gromacs/export/_export.py

@@ -184,13 +184,24 @@ def _write_dihedrals(self, top, molecule_type):
         for dihedral in molecule_type.dihedrals:
             function = functions[type(dihedral)]
 
-            top.write(
-                f"{dihedral.atom1 :6d}"
-                f"{dihedral.atom2 :6d}"
-                f"{dihedral.atom3 :6d}"
-                f"{dihedral.atom4 :6d}"
-                f"{functions[type(dihedral)] :6d}",
-            )
+            if function in [4]:
+                # GROMACS seems to prefer the central atom first
+                top.write(
+                    f"{dihedral.atom2 :6d}"
+                    f"{dihedral.atom1 :6d}"
+                    f"{dihedral.atom3 :6d}"
+                    f"{dihedral.atom4 :6d}"
+                    f"{functions[type(dihedral)] :6d}",
+                )
+
+            else:
+                top.write(
+                    f"{dihedral.atom1 :6d}"
+                    f"{dihedral.atom2 :6d}"
+                    f"{dihedral.atom3 :6d}"
+                    f"{dihedral.atom4 :6d}"
+                    f"{functions[type(dihedral)] :6d}",
+                )
 
             if function in [1, 4]:
                 top.write(

openff/interchange/smirnoff/_gromacs.py

@@ -396,51 +396,34 @@ def _convert_dihedrals(
 
     # TODO: Ensure number of torsions written matches what is expected
     if improper_torsion_handler:
-        # Molecule/Topology.impropers lists the central atom **second** ...
-        for improper in unique_molecule.smirnoff_impropers:
-            topology_indices = tuple(
-                interchange.topology.atom_index(a) for a in improper
-            )
-            # ... so the tuple must be modified to list the central atom **first**,
-            # which is how the improper handler's slot map is built up
-            indices_to_match = (
-                topology_indices[1],
-                topology_indices[0],
-                topology_indices[2],
-                topology_indices[3],
+        for top_key in improper_torsion_handler.key_map:
+            atoms = tuple(
+                interchange.topology.atom(topology_atom_index)
+                for topology_atom_index in top_key.atom_indices
             )
 
-            molecule_indices = tuple(unique_molecule.atom_index(a) for a in improper)
+            if not all(atom in unique_molecule.atoms for atom in atoms):
+                continue
 
-            # Now, indices_to_match has the central atom listed **first**,
-            # but it's still listed second in molecule_indices
+            key = improper_torsion_handler.key_map[top_key]
+            params = improper_torsion_handler.potentials[key].parameters
 
-            for top_key in improper_torsion_handler.key_map:
-                if top_key.atom_indices[0] != indices_to_match[0]:
-                    continue
-                if top_key.atom_indices[1] != indices_to_match[1]:
-                    continue
-                if top_key.atom_indices[2] != indices_to_match[2]:
-                    continue
-                if top_key.atom_indices[3] != indices_to_match[3]:
-                    continue
-                if indices_to_match == top_key.atom_indices:
-                    key = improper_torsion_handler.key_map[top_key]
-                    params = improper_torsion_handler.potentials[key].parameters
+            idivf = int(params["idivf"])
 
-                    idivf = int(params["idivf"])
+            molecule_indices = tuple(atom.molecule_atom_index for atom in atoms)
 
-                    molecule.dihedrals.append(
-                        PeriodicImproperDihedral(
-                            atom1=molecule_indices[1] + 1,
-                            atom2=molecule_indices[0] + 1,
-                            atom3=molecule_indices[2] + 1,
-                            atom4=molecule_indices[3] + 1,
-                            phi=params["phase"].to(unit.degree),
-                            k=params["k"].to(unit.kilojoule_per_mole) / idivf,
-                            multiplicity=int(params["periodicity"]),
-                        ),
-                    )
+            dihedral = PeriodicImproperDihedral(
+                atom1=molecule_indices[0] + 1,
+                atom2=molecule_indices[1] + 1,
+                atom3=molecule_indices[2] + 1,
+                atom4=molecule_indices[3] + 1,
+                phi=params["phase"].to(unit.degree),
+                k=params["k"].to(unit.kilojoule_per_mole) / idivf,
+                multiplicity=int(params["periodicity"]),
+            )
+
+            if dihedral not in molecule.dihedrals:
+                molecule.dihedrals.append(dihedral)
 
 
 def _convert_settles(