Add support for ff19SB force field

python/BioSimSpace/FreeEnergy/_relative.py

@@ -1064,15 +1064,16 @@ def _somd2_extract(parquet_file, T=None, estimator="MBAR"):
             lam = float(metadata["lambda"])
         except:
             raise ValueError("Parquet metadata does not contain 'lambda'.")
-        try:
-            lambda_array = metadata["lambda_array"]
-        except:
-            raise ValueError("Parquet metadata does not contain 'lambda array'")
         if not is_mbar:
             try:
                 lambda_grad = metadata["lambda_grad"]
             except:
                 raise ValueError("Parquet metadata does not contain 'lambda grad'")
+        else:
+            try:
+                lambda_grad = metadata["lambda_grad"]
+            except:
+                lambda_grad = []
 
         # Make sure that the temperature is correct.
         if not T == temperature:
@@ -1085,8 +1086,8 @@ def _somd2_extract(parquet_file, T=None, estimator="MBAR"):
         df = table.to_pandas()
 
         if is_mbar:
-            # Extract the columns correspodning to the lambda array.
-            df = df[[x for x in lambda_array]]
+            # Extract all columns other than those used for the gradient.
+            df = df[[x for x in df.columns if x not in lambda_grad]]
 
             # Subtract the potential at the simulated lambda.
             df = df.subtract(df[lam], axis=0)
@@ -1103,13 +1104,13 @@ def _somd2_extract(parquet_file, T=None, estimator="MBAR"):
 
                 # Forward difference.
                 if lam_delta > lam:
-                    double_incr = (lam_delta - lam) * 2
-                    grad = (df[lam_delta] - df[lam]) * 2 / double_incr
+                    incr = lam_delta - lam
+                    grad = (df[lam_delta] - df[lam]) / incr
 
                 # Backward difference.
                 else:
-                    double_incr = (lam - lam_delta) * 2
-                    grad = (df[lam] - df[lam_delta]) * 2 / double_incr
+                    incr = lam - lam_delta
+                    grad = (df[lam] - df[lam_delta]) / incr
 
             # Central difference.
             else:
@@ -1311,10 +1312,21 @@ def _analyse_internal(files, temperatures, lambdas, engine, estimator, **kwargs)
 
         # Preprocess the data.
         try:
-            processed_data = Relative._preprocess_data(data, estimator, **kwargs)
-            processed_data = _alchemlyb.concat(processed_data)
-        except:
-            _warnings.warn("Could not preprocess the data!")
+            preprocess = kwargs.pop("preprocess", True)
+        except KeyError:
+            preprocess = True
+
+        if not isinstance(preprocess, bool):
+            raise TypeError("'preprocess' must be of type 'bool'.")
+
+        if preprocess:
+            try:
+                processed_data = Relative._preprocess_data(data, estimator, **kwargs)
+                processed_data = _alchemlyb.concat(processed_data)
+            except:
+                _warnings.warn("Could not preprocess the data!")
+                processed_data = _alchemlyb.concat(data)
+        else:
             processed_data = _alchemlyb.concat(data)
 
         mbar_method = None

python/BioSimSpace/Node/_node.py

@@ -36,7 +36,7 @@
 # Set the default node directory.
 _node_dir = _os.path.dirname(__file__) + "/_nodes"
 
-__all__ = ["list", "help", "run", "setNodeDirectory"]
+__all__ = ["list", "help", "run", "setNodeDirectory", "getNodeDirectory"]
 
 
 def list():
@@ -90,7 +90,7 @@ def help(name):
     print(proc.stdout)
 
 
-def run(name, args={}):
+def run(name, args={}, work_dir=None):
     """
     Run a node.
 
@@ -103,6 +103,10 @@ def run(name, args={}):
     args : dict
         A dictionary of arguments to be passed to the node.
 
+    work_dir : str, optional
+        The working directory in which to run the node. If not specified,
+        the current working directory is used.
+
     Returns
     -------
 
@@ -112,9 +116,18 @@ def run(name, args={}):
 
     # Validate the input.
 
+    if not isinstance(name, str):
+        raise TypeError("'name' must be of type 'str'.")
+
     if not isinstance(args, dict):
         raise TypeError("'args' must be of type 'dict'.")
 
+    if work_dir is not None:
+        if not isinstance(work_dir, str):
+            raise TypeError("'work_dir' must be of type 'str'.")
+    else:
+        work_dir = _os.getcwd()
+
     # Apped the node directory name.
     full_name = _node_dir + "/" + name
 
@@ -123,45 +136,50 @@ def run(name, args={}):
         if not _os.path.isfile(full_name + ".py"):
             raise ValueError(
                 "Cannot find node: '%s'. " % name
+                + "in directory '%s'. " % _node_dir
                 + "Run 'Node.list()' to see available nodes!"
             )
         else:
             full_name += ".py"
 
-    # Write a YAML configuration file for the BioSimSpace node.
-    if len(args) > 0:
-        with open("input.yaml", "w") as file:
-            _yaml.dump(args, file, default_flow_style=False)
-
-        # Create the command.
-        command = "%s/python %s --config input.yaml" % (
-            _SireBase.getBinDir(),
-            full_name,
-        )
+    with _Utils.chdir(work_dir):
+        # Write a YAML configuration file for the BioSimSpace node.
+        if len(args) > 0:
+            with open("input.yaml", "w") as file:
+                _yaml.dump(args, file, default_flow_style=False)
 
-    # No arguments.
-    else:
-        command = "%s/python %s" % (_SireBase.getBinDir(), full_name)
+            # Create the command.
+            command = "%s/python %s --config input.yaml" % (
+                _SireBase.getBinDir(),
+                full_name,
+            )
 
-    # Run the node as a subprocess.
-    proc = _subprocess.run(
-        _Utils.command_split(command), shell=False, text=True, stderr=_subprocess.PIPE
-    )
+        # No arguments.
+        else:
+            command = "%s/python %s" % (_SireBase.getBinDir(), full_name)
+
+        # Run the node as a subprocess.
+        proc = _subprocess.run(
+            _Utils.command_split(command),
+            shell=False,
+            text=True,
+            stderr=_subprocess.PIPE,
+        )
 
-    if proc.returncode == 0:
-        # Read the output YAML file into a dictionary.
-        with open("output.yaml", "r") as file:
-            output = _yaml.safe_load(file)
+        if proc.returncode == 0:
+            # Read the output YAML file into a dictionary.
+            with open("output.yaml", "r") as file:
+                output = _yaml.safe_load(file)
 
-        # Delete the redundant YAML files.
-        _os.remove("input.yaml")
-        _os.remove("output.yaml")
+            # Delete the redundant YAML files.
+            _os.remove("input.yaml")
+            _os.remove("output.yaml")
 
-        return output
+            return output
 
-    else:
-        # Print the standard error, decoded as UTF-8.
-        print(proc.stderr)
+        else:
+            # Print the standard error, decoded as UTF-8.
+            print(proc.stderr)
 
 
 def setNodeDirectory(dir):
@@ -180,3 +198,16 @@ def setNodeDirectory(dir):
 
     global _node_dir
     _node_dir = dir
+
+
+def getNodeDirectory():
+    """
+    Get the directory of the node library.
+
+    Returns
+    -------
+
+    dir : str
+        The path to the node library.
+    """
+    return _node_dir

python/BioSimSpace/Parameters/_parameters.py

@@ -37,16 +37,14 @@
 
 # A dictionary mapping AMBER protein force field names to their pdb2gmx
 # compatibility. Note that the names specified below will be used for the
-# parameterisation functions, so they should be suitably formatted. Once we
-# have CMAP support we should be able to determine the available force fields
-# by scanning the AmberTools installation directory, as we do for those from
-# OpenFF.
+# parameterisation functions, so they should be suitably formatted.
 _amber_protein_forcefields = {
     "ff03": True,
     "ff99": True,
     "ff99SB": False,
     "ff99SBildn": False,
     "ff14SB": False,
+    "ff19SB": False,
 }
 
 from .. import _amber_home, _gmx_exe, _gmx_path, _isVerbose

python/BioSimSpace/Process/_amber.py

@@ -172,6 +172,13 @@ def __init__(
         # Flag to indicate whether the original system has a box.
         self._has_box = _AmberConfig.hasBox(self._system, self._property_map)
 
+        # Take note of whether the original reference system was None
+        # This will be used later to avoid duplication
+        if reference_system is not None:
+            self._is_real_reference = True
+        else:
+            self._is_real_reference = False
+
         # If the path to the executable wasn't specified, then search
         # for it in AMBERHOME and the PATH.
         if exe is None:
@@ -333,7 +340,6 @@ def _setup(self, **kwargs):
         else:
             # Check for perturbable molecules and convert to the chosen end state.
             system = self._checkPerturbable(system)
-            reference_system = self._checkPerturbable(reference_system)
 
         # RST file (coordinates).
         try:
@@ -348,10 +354,17 @@ def _setup(self, **kwargs):
 
         # Reference file for position restraints.
         try:
-            file = _os.path.splitext(self._ref_file)[0]
-            _IO.saveMolecules(
-                file, reference_system, "rst7", property_map=self._property_map
-            )
+            if self._is_real_reference:
+                reference_system, _ = _squash(
+                    system, explicit_dummies=self._explicit_dummies
+                )
+                reference_system = self._checkPerturbable(reference_system)
+                file = _os.path.splitext(self._ref_file)[0]
+                _IO.saveMolecules(
+                    file, reference_system, "rst7", property_map=self._property_map
+                )
+            else:
+                _shutil.copy(self._rst_file, self._ref_file)
         except Exception as e:
             msg = "Failed to write reference system to 'RST7' format."
             if _isVerbose():

python/BioSimSpace/Process/_process.py

@@ -731,7 +731,10 @@ def _checkPerturbable(self, system):
                     "in the 'property_map' argument."
                 )
             else:
-                is_lambda1 = self._property_map["is_lambda1"].value()
+                try:
+                    is_lambda1 = self._property_map["is_lambda1"].value()
+                except:
+                    is_lambda1 = self._property_map["is_lambda1"]
                 self._property_map.pop("is_lambda1")
 
             # Loop over all perturbable molecules in the system and replace them