diff --git a/src/somd2/_utils/_somd1.py b/src/somd2/_utils/_somd1.py index d95046c..a4fa418 100644 --- a/src/somd2/_utils/_somd1.py +++ b/src/somd2/_utils/_somd1.py @@ -84,7 +84,7 @@ def apply_pert(system, pert_file): return system -def make_compatible(system): +def make_compatible(system, fix_perturbable_zero_sigams=False): """ Makes a perturbation SOMD1 compatible. @@ -94,6 +94,9 @@ def make_compatible(system): system : sire.system.System, sire.legacy.System.System The system containing the molecules to be perturbed. + fix_perturbable_zero_sigams : bool + Whether to prevent LJ sigma values being perturbed to zero. + Returns ------- @@ -107,6 +110,9 @@ def make_compatible(system): "'system' must of type 'sire.system.System' or 'sire.legacy.System.System'" ) + if not isinstance(fix_perturbable_zero_sigams, bool): + raise TypeError("'fix_perturbable_zero_sigams' must be of type 'bool'.") + # Extract the legacy system. if isinstance(system, _LegacySystem): system = _System(system) @@ -132,36 +138,36 @@ def make_compatible(system): ################################## # First fix zero LJ sigmas values. ################################## - - # Tolerance for zero sigma values. - null_lj_sigma = 1e-9 - - for atom in mol.atoms(): - # Get the end state LJ sigma values. - lj0 = atom.property("LJ0") - lj1 = atom.property("LJ1") - - # Lambda = 0 state has a zero sigma value. - if abs(lj0.sigma().value()) <= null_lj_sigma: - # Use the sigma value from the lambda = 1 state. - edit_mol = ( - edit_mol.atom(atom.index()) - .set_property( - "LJ0", _SireMM.LJParameter(lj1.sigma(), lj0.epsilon()) + if fix_perturbable_zero_sigams: + # Tolerance for zero sigma values. + null_lj_sigma = 1e-9 + + for atom in mol.atoms(): + # Get the end state LJ sigma values. + lj0 = atom.property("LJ0") + lj1 = atom.property("LJ1") + + # Lambda = 0 state has a zero sigma value. + if abs(lj0.sigma().value()) <= null_lj_sigma: + # Use the sigma value from the lambda = 1 state. + edit_mol = ( + edit_mol.atom(atom.index()) + .set_property( + "LJ0", _SireMM.LJParameter(lj1.sigma(), lj0.epsilon()) + ) + .molecule() ) - .molecule() - ) - - # Lambda = 1 state has a zero sigma value. - if abs(lj1.sigma().value()) <= null_lj_sigma: - # Use the sigma value from the lambda = 0 state. - edit_mol = ( - edit_mol.atom(atom.index()) - .set_property( - "LJ1", _SireMM.LJParameter(lj0.sigma(), lj1.epsilon()) + + # Lambda = 1 state has a zero sigma value. + if abs(lj1.sigma().value()) <= null_lj_sigma: + # Use the sigma value from the lambda = 0 state. + edit_mol = ( + edit_mol.atom(atom.index()) + .set_property( + "LJ1", _SireMM.LJParameter(lj0.sigma(), lj1.epsilon()) + ) + .molecule() ) - .molecule() - ) ######################## # Now process the bonds. diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py index 2ed5307..e75f1ef 100644 --- a/src/somd2/config/_config.py +++ b/src/somd2/config/_config.py @@ -109,11 +109,13 @@ def __init__( include_constrained_energies=False, dynamic_constraints=True, ghost_modifications=True, + fix_perturbable_zero_sigmas=False, charge_difference=None, coalchemical_restraint_dist=None, com_reset_frequency=10, minimise=True, minimisation_constraints=False, + minimisation_errors=False, equilibration_time="0 ps", equilibration_timestep="2 fs", equilibration_constraints=True, @@ -258,6 +260,9 @@ def __init__( sampling of non-physical conformations. We implement the recommended modifcations from https://pubs.acs.org/doi/10.1021/acs.jctc.0c01328 + fix_perturbable_zero_sigmas: bool + Whether to prevent LJ sigma values being perturbed to zero. + charge_difference: int The charge difference between the two end states. (Perturbed minus reference.) If None, then alchemical ions will automatically be @@ -285,6 +290,9 @@ def __init__( constraints will be used. If True, then the use of constraints will be determined based on the value of 'equilibration_constraints'. + minimisation_errors: bool + Whether to raise an exception if a minimisation fails to converge. + equilibration_time: str Time interval for equilibration. Only simulations starting from scratch will be equilibrated. @@ -498,11 +506,13 @@ def __init__( self.include_constrained_energies = include_constrained_energies self.dynamic_constraints = dynamic_constraints self.ghost_modifications = ghost_modifications + self.fix_perturbable_zero_sigmas = fix_perturbable_zero_sigmas self.charge_difference = charge_difference self.coalchemical_restraint_dist = coalchemical_restraint_dist self.com_reset_frequency = com_reset_frequency self.minimise = minimise self.minimisation_constraints = minimisation_constraints + self.minimisation_errors = minimisation_errors self.equilibration_time = equilibration_time self.equilibration_timestep = equilibration_timestep self.equilibration_constraints = equilibration_constraints @@ -1145,6 +1155,16 @@ def ghost_modifications(self, ghost_modifications): raise ValueError("'ghost_modifications' must be of type 'bool'") self._ghost_modifications = ghost_modifications + @property + def fix_perturbable_zero_sigmas(self): + return self._fix_perturbable_zero_sigmas + + @fix_perturbable_zero_sigmas.setter + def fix_perturbable_zero_sigmas(self, fix_perturbable_zero_sigmas): + if not isinstance(fix_perturbable_zero_sigmas, bool): + raise ValueError("'fix_perturbable_zero_sigmas' must be of type 'bool'") + self._fix_perturbable_zero_sigmas = fix_perturbable_zero_sigmas + @property def charge_difference(self): return self._charge_difference @@ -1222,6 +1242,16 @@ def minimisation_constraints(self, minimisation_constraints): raise ValueError("'minimisation_constraints' must be of type 'bool'") self._minimisation_constraints = minimisation_constraints + @property + def minimisation_errors(self): + return self._minimisation_errors + + @minimisation_errors.setter + def minimisation_errors(self, minimisation_errors): + if not isinstance(minimisation_errors, bool): + raise ValueError("'minimisation_errors' must be of type 'bool'") + self._minimisation_errors = minimisation_errors + @property def equilibration_time(self): return self._equilibration_time diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py index 44ca71a..2b8bbc3 100644 --- a/src/somd2/runner/_base.py +++ b/src/somd2/runner/_base.py @@ -190,7 +190,9 @@ def __init__(self, system, config): self._config._extra_args["check_for_h_by_ambertype"] = False # Make sure that perturbable LJ sigmas aren't scaled to zero. - self._config._extra_args["fix_perturbable_zero_sigmas"] = True + self._config._extra_args["fix_perturbable_zero_sigmas"] = ( + config.fix_perturbable_zero_sigmas + ) # We're running in SOMD1 compatibility mode. if self._config.somd1_compatibility: @@ -199,7 +201,10 @@ def __init__(self, system, config): # First, try to make the perturbation SOMD1 compatible. _logger.info("Applying SOMD1 perturbation compatibility.") - self._system = make_compatible(self._system) + self._system = make_compatible( + self._system, + fix_perturbable_zero_sigmas=self.config.fix_perturbable_zero_sigmas, + ) self._system = _sr.morph.link_to_reference(self._system) # Next, swap the water topology so that it is in AMBER format. diff --git a/src/somd2/runner/_repex.py b/src/somd2/runner/_repex.py index 00718a0..01694ca 100644 --- a/src/somd2/runner/_repex.py +++ b/src/somd2/runner/_repex.py @@ -824,10 +824,12 @@ def run(self): replica_list[i * num_workers : (i + 1) * num_workers], ): if not success: - _logger.error( - f"Minimisation failed for {_lam_sym} = {self._lambda_values[index]:.5f}: {e}" - ) - raise e + msg = f"Minimisation failed for {_lam_sym} = {self._lambda_values[index]:.5f}: {e}" + if self.config.minimisation_errors: + _logger.error(msg) + raise e + else: + _logger.warning(msg) except KeyboardInterrupt: _logger.error("Minimisation cancelled. Exiting.") _sys.exit(1) @@ -1387,6 +1389,8 @@ def _equilibrate(self, index): # Reset the timer. if self._initial_time[index].value() != 0: system.set_time(self._initial_time[index]) + else: + system.set_time(_sr.u("0ps")) # Delete the dynamics object. self._dynamics_cache.delete(index) diff --git a/src/somd2/runner/_runner.py b/src/somd2/runner/_runner.py index 0e42473..88bfbcd 100644 --- a/src/somd2/runner/_runner.py +++ b/src/somd2/runner/_runner.py @@ -460,7 +460,11 @@ def generate_lam_vals(lambda_base, increment=0.001): perturbable_constraint=perturbable_constraint, ) except Exception as e: - raise RuntimeError(f"Minimisation failed: {e}") + msg = f"Minimisation failed for {_lam_sym} = {lambda_value:.5f}: {e}" + if self.confing.minimisation_errors: + raise RuntimeError(msg) + else: + _logger.warning(msg) # Equilibration. is_equilibrated = False @@ -523,6 +527,8 @@ def generate_lam_vals(lambda_base, increment=0.001): # Reset the timer. if self._initial_time[index].value() != 0: system.set_time(self._initial_time[index]) + else: + system.set_time(_sr.u("0ps")) except Exception as e: try: