Merge pull request #50 from OpenBioSim/fix_45-48_main

Backport fixes from PR #49

Author: lohedges

doc/source/tutorials/hydration_freenrg.rst

@@ -99,12 +99,25 @@ dynamics engine called SOMD. (This is a wrapper around the excellent
 `OpenMM <https://openmm.org>`__ package and is the default engine if no other
 packages are present.)
 
+>>> protocol = BSS.Protocol.FreeEnergyProduction()
+>>> free_somd = BSS.FreeEnergy.Relative(solvated, protocol, engine="somd", work_dir="freenrg_somd/free")
+>>> vac_somd = BSS.FreeEnergy.Relative(merged.toSystem(), protocol, engine="somd", work_dir="freenrg_somd/vacuum")
+
 .. note ::
 
-   It's possible to use a different protocol or molecular dynamics engine for each leg.
+   For GROMACS, you can use :data:`BioSimSpace.Protocol.FreeEnergyMinimisation` and
+   :data:`BioSimSpace.Protocol.FreeEnergyEquilibration` to minimise and equilibrate at each
+   :math:`{\lambda}` value prior to setting up to the production run above. For SOMD this
+   isn't possible, so you can either use GROMACS to prepare production input for each window,
+   or use the regular :data:`BioSimSpace.Protocol.Minimisation` and
+   :data:`BioSimSpace.Protocol.Equilibration` protocols to minimise and equilibrate
+   the :math:`{\lambda=0}` and :math:`{\lambda=1}` states only using any supported
+   engine from :data:`BioSimSpace.Process`.
 
->>> free_somd = BSS.FreeEnergy.Relative(solvated, protocol, engine="somd", work_dir="freenrg_somd/free")
->>> vac_somd = BSS.FreeEnergy.Relative(merged.toSystem(), protocol, engine="somd", work_dir="freenrg_somd/vacuum")
+.. note ::
+
+   It's possible to use a different protocol or molecular dynamics engine for each leg,
+   e.g. if you want to use a different :math:`{\lambda}` schedule.
 
 When complete, BioSimSpace will have set up a folder hierarchy containing
 everything that is needed to run the hydration free energy calculation

python/BioSimSpace/FreeEnergy/_relative.py

@@ -95,6 +95,8 @@ def __init__(
         setup_only=False,
         ignore_warnings=False,
         show_errors=True,
+        extra_options={},
+        extra_lines=[],
         property_map={},
     ):
         """
@@ -139,6 +141,13 @@ def __init__(
             run file. This option is specific to GROMACS and will be ignored
             when a different molecular dynamics engine is chosen.
 
+        extra_options : dict
+            A dictionary containing extra options. Overrides the defaults generated
+            by the protocol.
+
+        extra_lines : [str]
+            A list of extra lines to put at the end of the configuration file.
+
         property_map : dict
             A dictionary that maps system "properties" to their user defined
             values. This allows the user to refer to properties with their
@@ -249,6 +258,23 @@ def __init__(
             raise ValueError("'show_errors' must be of type 'bool.")
         self._show_errors = show_errors
 
+        # Check the extra options.
+        if not isinstance(extra_options, dict):
+            raise TypeError("'extra_options' must be of type 'dict'.")
+        else:
+            keys = extra_options.keys()
+            if not all(isinstance(k, str) for k in keys):
+                raise TypeError("Keys of 'extra_options' must be of type 'str'.")
+        self._extra_options = extra_options
+
+        # Check the extra lines.
+        if not isinstance(extra_lines, list):
+            raise TypeError("'extra_lines' must be of type 'list'.")
+        else:
+            if not all(isinstance(line, str) for line in extra_lines):
+                raise TypeError("Lines in 'extra_lines' must be of type 'str'.")
+        self._extra_lines = extra_lines
+
         # Check that the map is valid.
         if not isinstance(property_map, dict):
             raise TypeError("'property_map' must be of type 'dict'")
@@ -856,6 +882,8 @@ def _initialise_runner(self, system):
                 self._protocol,
                 platform=platform,
                 work_dir=first_dir,
+                extra_options=self._extra_options,
+                extra_lines=self._extra_lines,
                 property_map=self._property_map,
             )
             if self._setup_only:
@@ -871,6 +899,9 @@ def _initialise_runner(self, system):
                 work_dir=first_dir,
                 ignore_warnings=self._ignore_warnings,
                 show_errors=self._show_errors,
+                extra_options=self._extra_options,
+                extra_lines=self._extra_lines,
+                property_map=self._property_map,
             )
             if self._setup_only:
                 del first_process

python/BioSimSpace/_Config/_somd.py

@@ -232,8 +232,11 @@ def createConfig(self, extra_options={}, extra_lines=[]):
                     % self._protocol.getStartTemperature().kelvin().value()
                 )
 
-            # Friction coefficient (1 / ps).
-            protocol_dict["inverse friction"] = "{:.5f}".format(
+            # Friction coefficient (1 / ps). Note that the unit is written
+            # to the configuration file as "picosecond". This is because SOMD
+            # applies the value to OpenMM without inverting it, so has things the
+            # wrong way round, i.e. it uses the inverse friction as the friction.
+            protocol_dict["inverse friction"] = "{:.5f} picosecond".format(
                 1 / self._protocol.getThermostatTimeConstant().picoseconds().value()
             )
 

test/Process/test_somd.py

@@ -172,10 +172,8 @@ def run_process(system, protocol):
                 )
                 res = False
 
-    return not res
-
     # Make sure the process didn't error.
-    assert not process.isError()
+    assert not res
 
     # Make sure that we get a molecular system back.
     assert process.getSystem() is not None

test/Sandpit/Exscientia/Process/test_somd.py

@@ -155,10 +155,8 @@ def run_process(system, protocol):
                 )
                 res = False
 
-    return not res
-
     # Make sure the process didn't error.
-    assert not process.isError()
+    assert not res
 
     # Make sure that we get a molecular system back.
     assert process.getSystem() is not None