Merge pull request #138 from OpenBioSim/feature_terminal_flip

Add support for Terminal Flip Monte Carlo

Author: lohedges

README.md

@@ -219,6 +219,29 @@ require a different `nvcc` to that provided by conda, you can set the
 Depending on your setup, you may also need to install the `cuda-nvvm` package
 from `conda-forge`.
 
+## Terminal ring flip Monte Carlo
+
+SOMD2 supports terminal ring flip Monte Carlo (MC) moves to improve sampling
+of terminal aromatic rings in perturbable ligands, as described in
+[this paper](https://chemrxiv.org/doi/full/10.26434/chemrxiv-2025-2zkx5).
+Each move attempts a discrete rotation of a terminal ring around the bond
+connecting it to the rest of the molecule, accepted or rejected via the
+Metropolis criterion. Terminal ring groups are detected automatically from
+the molecular connectivity of perturbable molecules.
+
+To enable terminal flip MC, set the frequency at which moves are attempted:
+
+```
+somd2 perturbable_system.bss --terminal-flip-frequency "1 ps"
+```
+
+The flip angle for each group is determined automatically from the ring
+geometry. To override this for all groups:
+
+```
+somd2 perturbable_system.bss --terminal-flip-frequency "1 ps" --terminal-flip-angle "180 degrees"
+```
+
 ## Analysis
 
 Simulation output will be written to the directory specified using the

src/somd2/_utils/__init__.py

@@ -19,14 +19,19 @@
 # along with SOMD2. If not, see <http://www.gnu.org/licenses/>.
 #####################################################################
 
-import platform as _platform
+import sys as _sys
 
-if _platform.system() == "Windows":
-    _lam_sym = "lambda"
-else:
+try:
+    "λΔ°".encode(_sys.stdout.encoding or "utf-8")
     _lam_sym = "λ"
+    _delta_sym = "ΔE"
+    _degree_sym = "°"
+except (UnicodeEncodeError, LookupError):
+    _lam_sym = "lambda"
+    _delta_sym = "delta"
+    _degree_sym = "deg"
 
-del _platform
+del _sys
 
 
 def _has_ghost(mol, idxs, is_lambda1=False):

src/somd2/config/_config.py

@@ -139,6 +139,9 @@ def __init__(
         replica_exchange=False,
         randomise_velocities=False,
         perturbed_system=None,
+        terminal_flip_frequency=None,
+        terminal_flip_angle=None,
+        terminal_flip_max_mobile_atoms=None,
         gcmc=False,
         gcmc_frequency=None,
         gcmc_selection=None,
@@ -370,13 +373,30 @@ def __init__(
             GPU resources are available.
 
         randomise_velocities: bool
-            Whether to randomise velocities at the start of each replica exchange cycle.
+            Whether to randomise velocities at the start of each replica exchange cycle
+            or following a terminal flip Monte Carlo move.
 
         perturbed_system: str
             The path to a stream file containing a Sire system for the equilibrated perturbed
             end state (lambda = 1). This will be used as the starting conformation all lambda
             windows > 0.5 when performing a replica exchange simulation.
 
+        terminal_flip_frequency: str
+            Frequency at which to attempt terminal ring flip Monte Carlo moves. If None
+            (the default), no terminal flip moves will be performed. When set, terminal
+            ring groups in perturbable molecules are detected automatically using Sire's
+            native connectivity. This must be a multiple of 'energy_frequency'.
+
+        terminal_flip_angle: str
+            Override the flip angle used for all terminal ring groups, e.g.
+            ``"180 degrees"``. If None (the default), the angle is determined
+            automatically for each group from its geometry.
+
+        terminal_flip_max_mobile_atoms: int or None
+            Maximum number of mobile atoms allowed in a terminal ring group.
+            Groups with more mobile atoms than this threshold are skipped during
+            detection. Defaults to None (no limit).
+
         gcmc: bool
             Whether to perform Grand Canonical Monte Carlo (GCMC) water insertions/deletions.
 
@@ -559,6 +579,9 @@ def __init__(
         self.replica_exchange = replica_exchange
         self.randomise_velocities = randomise_velocities
         self.perturbed_system = perturbed_system
+        self.terminal_flip_frequency = terminal_flip_frequency
+        self.terminal_flip_angle = terminal_flip_angle
+        self.terminal_flip_max_mobile_atoms = terminal_flip_max_mobile_atoms
         self.gcmc = gcmc
         self.gcmc_frequency = gcmc_frequency
         self.gcmc_selection = gcmc_selection
@@ -1994,6 +2017,80 @@ def perturbed_system(self, perturbed_system):
             self._perturbed_system = None
             self._perturbed_system_file = None
 
+    @property
+    def terminal_flip_frequency(self):
+        return self._terminal_flip_frequency
+
+    @terminal_flip_frequency.setter
+    def terminal_flip_frequency(self, terminal_flip_frequency):
+        if terminal_flip_frequency is not None:
+            if not isinstance(terminal_flip_frequency, str):
+                raise TypeError("'terminal_flip_frequency' must be of type 'str'")
+
+            from sire.units import picosecond
+
+            try:
+                t = _sr.u(terminal_flip_frequency)
+            except Exception:
+                raise ValueError(
+                    f"Unable to parse 'terminal_flip_frequency' as a Sire GeneralUnit: "
+                    f"{terminal_flip_frequency}"
+                )
+
+            if t.value() != 0 and not t.has_same_units(picosecond):
+                raise ValueError("'terminal_flip_frequency' units are invalid.")
+
+            self._terminal_flip_frequency = t
+        else:
+            self._terminal_flip_frequency = None
+
+    @property
+    def terminal_flip_angle(self):
+        return self._terminal_flip_angle
+
+    @terminal_flip_angle.setter
+    def terminal_flip_angle(self, terminal_flip_angle):
+        if terminal_flip_angle is not None:
+            if not isinstance(terminal_flip_angle, str):
+                raise TypeError("'terminal_flip_angle' must be of type 'str'")
+
+            from sire.units import degrees
+
+            try:
+                a = _sr.u(terminal_flip_angle)
+            except Exception:
+                raise ValueError(
+                    f"Unable to parse 'terminal_flip_angle' as a Sire GeneralUnit: "
+                    f"{terminal_flip_angle}"
+                )
+
+            if not a.has_same_units(degrees):
+                raise ValueError("'terminal_flip_angle' units are invalid.")
+
+            self._terminal_flip_angle = a
+        else:
+            self._terminal_flip_angle = None
+
+    @property
+    def terminal_flip_max_mobile_atoms(self):
+        return self._terminal_flip_max_mobile_atoms
+
+    @terminal_flip_max_mobile_atoms.setter
+    def terminal_flip_max_mobile_atoms(self, terminal_flip_max_mobile_atoms):
+        if terminal_flip_max_mobile_atoms is not None:
+            if not isinstance(terminal_flip_max_mobile_atoms, int):
+                try:
+                    terminal_flip_max_mobile_atoms = int(terminal_flip_max_mobile_atoms)
+                except:
+                    raise ValueError(
+                        "'terminal_flip_max_mobile_atoms' must be of type 'int'"
+                    )
+            if terminal_flip_max_mobile_atoms < 1:
+                raise ValueError(
+                    "'terminal_flip_max_mobile_atoms' must be greater than 0"
+                )
+        self._terminal_flip_max_mobile_atoms = terminal_flip_max_mobile_atoms
+
     @property
     def gcmc(self):
         return self._gcmc

src/somd2/runner/_base.py

@@ -664,6 +664,59 @@ def __init__(self, system, config):
             # Store the excess chemcical potential value.
             self._mu_ex = self._config.gcmc_excess_chemical_potential.value()
 
+        # Terminal flip specific validation and setup.
+        if self._config.terminal_flip_frequency is not None:
+            from math import isclose
+
+            # Make sure the terminal flip frequency is a multiple of the
+            # energy frequency.
+            ratio = (
+                self._config.terminal_flip_frequency / self._config.energy_frequency
+            ).value()
+
+            if not isclose(ratio, round(ratio), abs_tol=1e-4):
+                msg = "'terminal_flip_frequency' must be a multiple of 'energy_frequency'."
+                _logger.error(msg)
+                raise ValueError(msg)
+
+            # Auto-detect terminal ring groups using Sire connectivity.
+            from ._samplers import detect_terminal_groups
+
+            if isinstance(self._system, list):
+                mols = self._system[0]
+            else:
+                mols = self._system
+
+            flip_angle = (
+                self._config.terminal_flip_angle.to("degrees").value()
+                if self._config.terminal_flip_angle is not None
+                else None
+            )
+            self._terminal_groups = detect_terminal_groups(
+                mols,
+                flip_angle=flip_angle,
+                max_mobile_atoms=self._config.terminal_flip_max_mobile_atoms,
+            )
+
+            if not self._terminal_groups:
+                _logger.warning(
+                    "No terminal ring groups detected. Terminal flip moves will not "
+                    "be performed."
+                )
+            else:
+                _logger.info(
+                    f"Detected {len(self._terminal_groups)} terminal ring group(s) "
+                    f"for terminal flip MC."
+                )
+                for i, (angle, indices) in enumerate(self._terminal_groups):
+                    _logger.info(
+                        f"  Group {i}: flip angle = {angle}°, "
+                        f"anchor = {indices[0]}, pivot = {indices[1]}, "
+                        f"{len(indices) - 2} mobile atom(s)"
+                    )
+        else:
+            self._terminal_groups = []
+
         # Store the initial system time.
         if isinstance(self._system, list):
             self._initial_time = []
@@ -1148,6 +1201,7 @@ def increment_filename(base_filename, suffix):
             output_directory / f"energy_components_{lam}.txt"
         )
         filenames["gcmc_ghosts"] = str(output_directory / f"gcmc_ghosts_{lam}.txt")
+        filenames["sampler_stats"] = str(output_directory / f"sampler_stats_{lam}.pkl")
         if restart:
             filenames["config"] = str(
                 output_directory / increment_filename("config", "yaml")