Added simple rounding error solution

openff/interchange/components/interchange.py

@@ -332,6 +332,7 @@ def to_gromacs(
         hydrogen_mass: PositiveFloat = 1.007947,
         monolithic: bool = True,
         _merge_atom_types: bool = False,
+        _normalize_charges: bool = False,
     ):
         """
         Export this Interchange object to GROMACS files.
@@ -353,6 +354,14 @@ def to_gromacs(
         _merge_atom_types: bool, default = False
             The flag to define behaviour of GROMACSWriter. If True, then similar atom types will be merged.
             If False, each atom will have its own atom type.
+        _normalize_charges: bool, default = False
+            If True charges written to the topology are normalized per molecule to ensure that each molecule has a
+            charge which is exactly integer. Charges that are exactly 0 are not touched.
+            If False, charges are untouched.
+        _normalize_charges: bool, default = False
+            If True charges written to the topology are normalized per molecule to ensure that each molecule has a
+            charge which is exactly integer. Charges that are exactly 0 are not touched.
+            If False, charges are untouched.
 
         Notes
         -----
@@ -371,7 +380,11 @@ def to_gromacs(
             gro_file=prefix + ".gro",
         )
 
-        writer.to_top(monolithic=monolithic, _merge_atom_types=_merge_atom_types)
+        writer.to_top(
+            monolithic=monolithic,
+            _merge_atom_types=_merge_atom_types,
+            _normalize_charges=_normalize_charges,
+        )
         writer.to_gro(decimal=decimal)
 
         self.to_mdp(prefix + "_pointenergy.mdp")
@@ -407,6 +420,7 @@ def to_top(
         hydrogen_mass: PositiveFloat = 1.007947,
         monolithic: bool = True,
         _merge_atom_types: bool = False,
+        _normalize_charges: bool = False,
     ):
         """
         Export this Interchange to a GROMACS topology file.
@@ -425,12 +439,19 @@ def to_top(
         _merge_atom_types: book, default=False
             The flag to define behaviour of GROMACSWriter. If True, then similar atom types will be merged.
             If False, each atom will have its own atom type.
+        _normalize_charges: bool, default = False
+            If True charges written to the topology are normalized per molecule to ensure that each molecule has a
+            charge which is exactly integer. Charges that are exactly 0 are not touched.
+            If False, charges are untouched.
 
         Notes
         -----
         Molecule names in written files are not guaranteed to match the `Moleclue.name` attribute of the
         molecules in the topology, especially if they are empty strings or not unique.
 
+        Charges written to the topology are normalized per molecule to ensure that each molecule has a charge which is
+        exactly integer. Charges that are exactly 0 are not touched.
+
         """
         from openff.interchange.interop.gromacs.export._export import GROMACSWriter
         from openff.interchange.smirnoff._gromacs import _convert
@@ -442,6 +463,7 @@ def to_top(
         ).to_top(
             monolithic=monolithic,
             _merge_atom_types=_merge_atom_types,
+            _normalize_charges=_normalize_charges,
         )
 
     def to_gro(self, file_path: Path | str, decimal: int = 3):

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

@@ -27,7 +27,12 @@
     top_file: pathlib.Path | str
     gro_file: pathlib.Path | str
 
-    def to_top(self, monolithic: bool = True, _merge_atom_types: bool = False):
+    def to_top(
+        self,
+        monolithic: bool = True,
+        _merge_atom_types: bool = False,
+        _normalize_charges: bool = False,
+    ):
         """Write a GROMACS topology file."""
         with open(self.top_file, "w") as top:
             self._write_defaults(top)
@@ -52,6 +57,7 @@
                 monolithic=monolithic,
                 mapping_to_reduced_atom_types=mapping_to_reduced_atom_types,
                 merge_atom_types=_merge_atom_types,
+                normalize_charges=_normalize_charges,
             )
 
             self._write_system(top)
@@ -173,6 +179,7 @@
         monolithic: bool,
         mapping_to_reduced_atom_types,
         merge_atom_types: bool,
+        normalize_charges: bool,
     ):
         for molecule_name, molecule_type in self.system.molecule_types.items():
             # this string needs to be something that plays nicely in file paths
@@ -198,6 +205,7 @@
                 molecule_type,
                 mapping_to_reduced_atom_types,
                 merge_atom_types,
+                normalize_charges,
             )
             self._write_pairs(molecule_file, molecule_type)
             self._write_bonds(molecule_file, molecule_type)
@@ -218,11 +226,51 @@
         molecule_type,
         mapping_to_reduced_atom_types,
         merge_atom_types: bool,
+        normalize_charges: bool,
     ):
         top.write("[ atoms ]\n")
         top.write(";index, atom type, resnum, resname, name, cgnr, charge, mass\n")
 
-        for atom in molecule_type.atoms:
+        def _adjust_charges(charges: numpy.array, tolerance=8) -> numpy.array:
+            """
+            Adjust charges so that written charge for a molecule type is integer.
+
+            We first get the initial charges and round them.
+            Note, that charges that are exactly equal to zero are never touched.
+            """
+            # placeholder for output charges
+            rounded_charges = numpy.round(charges, tolerance)
+            # integer total charge
+            total_charge = numpy.round(numpy.sum(charges), 0)
+            # non-zero charge indices
+            indices = numpy.where(rounded_charges != 0)[0]
+
+            def _rounding_error(_arr, _sum, _tolerance):
+                return numpy.round(numpy.sum(_arr) - _sum, _tolerance)
+
+            # Initial error due to rounding
+            rounding_error = _rounding_error(rounded_charges, total_charge, tolerance)
+
+            # We correct rounded_charges to achieve tolerance
+            if rounding_error != 0:
+                rounded_charges[indices] += numpy.round(
+                    -rounding_error / len(indices),
+                    tolerance,
+                )
+                diff = _rounding_error(rounded_charges, total_charge, tolerance)
+
+                # If there's a remainder of a charge after spreading everything out evenly,
+                # dump the remainder on the first nonzero atom.
+                rounded_charges[indices[0]] = numpy.round(
+                    rounded_charges[indices[0]] - diff,
+                    tolerance,
+                )
+            return rounded_charges
+
+        charges_to_write = numpy.array([atom.charge.m for atom in molecule_type.atoms])
+        rounded_charges = _adjust_charges(charges_to_write) if normalize_charges else charges_to_write
+
+        for atom, charge in zip(molecule_type.atoms, rounded_charges):
             if merge_atom_types:
                 top.write(
                     f"{atom.index:6d} "
@@ -231,7 +279,7 @@
                     f"{atom.residue_name:8s} "
                     f"{atom.name:6s}"
                     f"{atom.charge_group_number:6d}"
-                    f"{atom.charge.m:20.12f}"
+                    f"{charge:20.12f}"
                     f"{atom.mass.m:20.12f}\n",
                 )
             else:
@@ -242,7 +290,7 @@
                     f"{atom.residue_name:8s} "
                     f"{atom.name:6s}"
                     f"{atom.charge_group_number:6d}"
-                    f"{atom.charge.m:20.12f}"
+                    f"{charge:20.12f}"
                     f"{atom.mass.m:20.12f}\n",
                 )