Fix issue #243

CHANGELOG.md

@@ -21,6 +21,7 @@ Changelog
 * Fixed alchemical ion creation when an existing counter-ion is perturbed to water. The water was recorded as the alchemical ion, so it got the wrong lambda schedule, and the same counter-ion was reused for every unit of charge. The furthest counter-ions from the perturbable molecule are now used, once each, whenever the system has no free ion of the required type [#239](https://github.com/OpenBioSim/somd2/pull/239).
 * Raise an error when the `annihilate` or `decouple` lambda schedule is used with a perturbable molecule created by `sire.morph.annihilate()`, which also removes the bonded terms. ABFE ligands should be created with `sire.morph.decouple()`, leaving the choice between decoupling and annihilation to the lambda schedule [#239](https://github.com/OpenBioSim/somd2/pull/239).
 * Only fall back to the Aldeghi Boresch restraint search protocol when the system can't be used with the RXRX protocol. Systems that can't be used with either protocol, or restraint search settings that would save too few frames, now raise an error before the restraint search trajectory is run, rather than after a misleading fallback warning [#241](https://github.com/OpenBioSim/somd2/pull/241).
+* Save the GCMC end-state topologies, which include the additional waters, at the first checkpoint of each run as well as the last, and keep them on restart, so that they match the trajectory and checkpoints while a simulation is running, if it crashes, or if a restart never reaches a checkpoint [#244](https://github.com/OpenBioSim/somd2/pull/244).
 
 [2026.2.0](https://github.com/openbiosim/somd2/compare/2026.1.0...2026.2.0) - Sep 2026
 --------------------------------------------------------------------------------------

src/somd2/runner/_base.py

@@ -927,15 +927,21 @@ def __init__(self, system, config):
         # Create the lock file name.
         self._lock_file = str(self._config.output_directory / "somd2.lock")
 
-        # Write the end-state topologies to the output directory.
+        # Write the end-state topologies. A GCMC restart keeps those already
+        # written, since they include the additional waters.
         if isinstance(self._system, list):
             mols = self._system[0]
         else:
             mols = self._system
-        mols0 = _sr.morph.link_to_reference(mols)
-        mols1 = _sr.morph.link_to_perturbed(mols)
-        _sr.save(mols0, self._filenames["topology0"])
-        _sr.save(mols1, self._filenames["topology1"])
+        topologies = (self._filenames["topology0"], self._filenames["topology1"])
+        if not (
+            self._config.restart
+            and self._config.gcmc
+            and all(_Path(t).exists() for t in topologies)
+        ):
+            _sr.save(_sr.morph.link_to_reference(mols), topologies[0])
+            _sr.save(_sr.morph.link_to_perturbed(mols), topologies[1])
+        self._gcmc_topology_saved = False
 
         # Saved with checkpoints for the viewer, so a failure is only a warning.
         try:
@@ -1091,7 +1097,9 @@ def _set_tuned_pme_parameters(self):
             if pme_file.exists():
                 with open(pme_file) as f:
                     params = _yaml.safe_load(f)
-                _logger.info(f"Using PME parameters from {pme_file}: {params}")
+                _logger.info(
+                    f"Using PME parameters from {pme_file}: {self._describe_pme(params)}"
+                )
                 self._config._extra_args["pme_alpha"] = params["pme_alpha"]
                 self._config._extra_args["pme_grid"] = params["pme_grid"]
             elif self._config.tune_pme and not has_pme_options:
@@ -1142,15 +1150,39 @@ def _set_tuned_pme_parameters(self):
             _logger.info(f"The default PME parameters are already the fastest, {error}")
             return
 
-        _logger.info(
-            f"Using tuned PME parameters: alpha={params['pme_alpha']:.4f} nm^-1, "
-            f"grid={params['pme_grid']}, {error}"
-        )
+        _logger.info(f"Using tuned PME parameters: {self._describe_pme(params)}")
 
         _dict_to_yaml(params, str(pme_file))
         self._config._extra_args["pme_alpha"] = params["pme_alpha"]
         self._config._extra_args["pme_grid"] = params["pme_grid"]
 
+    @staticmethod
+    def _describe_pme(params):
+        """
+        A description of tuned PME parameters for the log.
+
+        Parameters
+        ----------
+
+        params : dict
+            The tuned PME parameters, and their relative force errors if known.
+
+        Returns
+        -------
+
+        description : str
+            The parameters, rounded for reading.
+        """
+        description = (
+            f"alpha={params['pme_alpha']:.4f} nm^-1, grid={params['pme_grid']}"
+        )
+        if "pme_error" in params and "pme_target_error" in params:
+            description += (
+                f", relative force error {params['pme_error']:.2e} "
+                f"(target {params['pme_target_error']:.2e})"
+            )
+        return description
+
     @property
     def _is_abfe_bound(self):
         """
@@ -2538,28 +2570,13 @@ def _checkpoint(
                         f"{self._standard_state_correction:.6f}"
                     )
 
-            if is_final_block:
-                # Save the end-state GCMC topologies for trajectory analysis and visualisation.
-                # This topology contains additional water molecules that are used for GCMC
-                # insertion moves.
-                if self._config.gcmc:
-                    mols0 = _sr.morph.link_to_reference(system)
-                    mols1 = _sr.morph.link_to_perturbed(system)
-
-                    # Save to AMBER format.
-                    _sr.save(mols0, self._filenames["topology0"])
-                    _sr.save(mols1, self._filenames["topology1"])
-
-                    # Save to PDB format.
-                    _sr.save(
-                        mols0,
-                        self._filenames["topology0"].replace(".prm7", ".pdb"),
-                    )
-                    _sr.save(
-                        mols1,
-                        self._filenames["topology1"].replace(".prm7", ".pdb"),
-                    )
+            # Saved at the first checkpoint of each run, and the final one for the
+            # final coordinates in the PDB files.
+            if self._config.gcmc and (is_final_block or not self._gcmc_topology_saved):
+                self._save_gcmc_topology(system)
+                self._gcmc_topology_saved = True
 
+            if is_final_block:
                 # Assemble and save the final trajectory.
                 if self._config.save_trajectories:
                     # Save the final trajectory chunk to file.
@@ -2948,6 +2965,27 @@ def _flush_energy_components(self, index):
             )
             _pq_local.write_table(table, filepath)
 
+    def _save_gcmc_topology(self, system):
+        """
+        Save the end-state topologies, including the additional waters used for
+        GCMC insertion moves, in AMBER and PDB format.
+
+        Parameters
+        ----------
+
+        system : :class: `System <sr.system.System>`
+            The system being simulated.
+        """
+        mols0 = _sr.morph.link_to_reference(system)
+        mols1 = _sr.morph.link_to_perturbed(system)
+
+        for mols, topology in (
+            (mols0, self._filenames["topology0"]),
+            (mols1, self._filenames["topology1"]),
+        ):
+            _sr.save(mols, topology)
+            _sr.save(mols, topology.replace(".prm7", ".pdb"))
+
     def _save_ghost_residues(self, index, gcmc_sampler):
         """
         Record the current GCMC ghost residue indices for a window. This must