Merge pull request #94 from OpenBioSim/fix_93

Fix issue #93

Author: lohedges

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.

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

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.

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)

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: