diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py index 8a6c0e9..d4482be 100644 --- a/src/somd2/config/_config.py +++ b/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 diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py index 739d45a..0948780 100644 --- a/src/somd2/runner/_base.py +++ b/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