Add option to control epsilon scaling for Beutler soft-core

src/somd2/config/_config.py

@@ -158,6 +158,7 @@ def __init__(
         softcore_form="zacharias",
         taylor_power=1,
         beutler_alpha=0.5,
+        beutler_fix_epsilon=True,
         output_directory="output",
         restart=False,
         use_backup=False,
@@ -477,6 +478,15 @@ def __init__(
             form. Must be >= 0. The default is 0.5. Only used when softcore_form is
             "beutler".
 
+        beutler_fix_epsilon: bool
+            Whether to hold LJ epsilon fixed at its real-atom value for
+            ghost-decoupling molecules when softcore_form is "beutler", so that the
+            Beutler (1-alpha) prefactor provides the sole LJ decay pathway. The
+            default is True. This is automatically disabled (regardless of this
+            setting) for any alchemical ion added to maintain a constant charge,
+            since an ion's persisting atom is a real (non-ghost) mutation and needs
+            its LJ epsilon to interpolate normally rather than being held fixed.
+
         output_directory: str
             Path to a directory to store output files.
 
@@ -623,6 +633,7 @@ def __init__(
         self.softcore_form = softcore_form
         self.taylor_power = taylor_power
         self.beutler_alpha = beutler_alpha
+        self.beutler_fix_epsilon = beutler_fix_epsilon
         self.somd1_compatibility = somd1_compatibility
         self.pert_file = pert_file
         self.auto_fix_minimise = auto_fix_minimise
@@ -2141,6 +2152,16 @@ def beutler_alpha(self, beutler_alpha):
             raise ValueError("'beutler_alpha' must be >= 0")
         self._beutler_alpha = beutler_alpha
 
+    @property
+    def beutler_fix_epsilon(self):
+        return self._beutler_fix_epsilon
+
+    @beutler_fix_epsilon.setter
+    def beutler_fix_epsilon(self, beutler_fix_epsilon):
+        if not isinstance(beutler_fix_epsilon, bool):
+            raise TypeError("'beutler_fix_epsilon' must be of type 'bool'")
+        self._beutler_fix_epsilon = beutler_fix_epsilon
+
     @property
     def use_backup(self):
         return self._use_backup

src/somd2/runner/_base.py

@@ -253,31 +253,6 @@ def __init__(self, system, config):
             self._config._extra_args["use_gcmc_lrc"] = True
             self._config._extra_args["num_gcmc_waters"] = self._config.gcmc_num_waters
 
-        # Set the soft-core form.
-        if self._config.softcore_form == "taylor":
-            self._config._extra_args["use_taylor_softening"] = True
-            self._config._extra_args["taylor_power"] = self._config.taylor_power
-        elif self._config.softcore_form == "beutler":
-            schedule_name = self._config._lambda_schedule_name
-            if schedule_name not in (None, "annihilate", "decouple"):
-                raise ValueError(
-                    "The Beutler soft-core form is only supported with the 'annihilate' "
-                    "or 'decouple' lambda schedules, or a custom schedule."
-                )
-            self._config._extra_args["use_beutler_softening"] = True
-            self._config._extra_args["beutler_alpha"] = self._config.beutler_alpha
-
-        # Build deferred schedules now that the softcore form is known.
-        fix_epsilon = self._config.softcore_form == "beutler"
-        if self._config._lambda_schedule_name == "annihilate":
-            from .._utils._schedules import annihilate as _annihilate
-
-            self._config._lambda_schedule = _annihilate(fix_epsilon=fix_epsilon)
-        elif self._config._lambda_schedule_name == "decouple":
-            from .._utils._schedules import decouple as _decouple
-
-            self._config._lambda_schedule = _decouple(fix_epsilon=fix_epsilon)
-
         # We're running in SOMD1 compatibility mode.
         if self._config.somd1_compatibility:
             from .._utils._somd1 import make_compatible
@@ -407,6 +382,47 @@ def __init__(self, system, config):
                     coalchemical_restraints
                 )
 
+        # Set the soft-core form.
+        if self._config.softcore_form == "taylor":
+            self._config._extra_args["use_taylor_softening"] = True
+            self._config._extra_args["taylor_power"] = self._config.taylor_power
+        elif self._config.softcore_form == "beutler":
+            schedule_name = self._config._lambda_schedule_name
+            if schedule_name not in (None, "annihilate", "decouple"):
+                raise ValueError(
+                    "The Beutler soft-core form is only supported with the 'annihilate' "
+                    "or 'decouple' lambda schedules, or a custom schedule."
+                )
+            self._config._extra_args["use_beutler_softening"] = True
+            self._config._extra_args["beutler_alpha"] = self._config.beutler_alpha
+
+        # Build deferred schedules now that the softcore form is known. Epsilon is
+        # only held fixed (with LJ decay handled entirely by the Beutler soft-core
+        # prefactor) for molecules undergoing a ghost-atom decoupling/annihilation.
+        # An alchemical ion is a real (non-ghost) atom mutating identity (e.g. a
+        # water oxygen turning into Na+), so its LJ epsilon needs to interpolate
+        # normally; fixing it would leave the ion's persisting atom stuck at its
+        # initial LJ parameters for the whole stage. Disable fix_epsilon whenever
+        # an alchemical ion has been added, regardless of the configured value.
+        fix_epsilon = (
+            self._config.softcore_form == "beutler" and self._config.beutler_fix_epsilon
+        )
+        if fix_epsilon and charge_diff != 0:
+            _logger.info(
+                "Disabling Beutler 'fix_epsilon' since an alchemical ion has been "
+                "added: the ion's persisting atom is a real (non-ghost) mutation "
+                "and needs its LJ epsilon to interpolate normally."
+            )
+            fix_epsilon = False
+        if self._config._lambda_schedule_name == "annihilate":
+            from .._utils._schedules import annihilate as _annihilate
+
+            self._config._lambda_schedule = _annihilate(fix_epsilon=fix_epsilon)
+        elif self._config._lambda_schedule_name == "decouple":
+            from .._utils._schedules import decouple as _decouple
+
+            self._config._lambda_schedule = _decouple(fix_epsilon=fix_epsilon)
+
         # Set the lambda values.
         if self._config.lambda_values:
             self._lambda_values = self._config.lambda_values