diff --git a/python/BioSimSpace/FreeEnergy/_relative.py b/python/BioSimSpace/FreeEnergy/_relative.py
index 1ba69e79..d6e1ce8f 100644
--- a/python/BioSimSpace/FreeEnergy/_relative.py
+++ b/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
diff --git a/python/BioSimSpace/Node/_node.py b/python/BioSimSpace/Node/_node.py
index 2e047df2..3d63cfda 100644
--- a/python/BioSimSpace/Node/_node.py
+++ b/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
diff --git a/python/BioSimSpace/Parameters/_parameters.py b/python/BioSimSpace/Parameters/_parameters.py
index dade5ac5..66834bbc 100644
--- a/python/BioSimSpace/Parameters/_parameters.py
+++ b/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
diff --git a/python/BioSimSpace/Process/_amber.py b/python/BioSimSpace/Process/_amber.py
index 59d14b36..5095b0fc 100644
--- a/python/BioSimSpace/Process/_amber.py
+++ b/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():
diff --git a/python/BioSimSpace/Process/_process.py b/python/BioSimSpace/Process/_process.py
index 020558e0..ca57b898 100644
--- a/python/BioSimSpace/Process/_process.py
+++ b/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
diff --git a/python/BioSimSpace/Process/_somd.py b/python/BioSimSpace/Process/_somd.py
index cababe0f..b606042b 100644
--- a/python/BioSimSpace/Process/_somd.py
+++ b/python/BioSimSpace/Process/_somd.py
@@ -238,6 +238,9 @@ def __init__(
         self._zero_dummy_impropers = kwargs.get("zero_dummy_impropers", False)
         if not isinstance(self._zero_dummy_impropers, bool):
             self._zero_dummy_impropers = False
+        self._no_dummy_modifications = kwargs.get("no_dummy_modifications", False)
+        if not isinstance(self._no_dummy_modifications, bool):
+            self._no_dummy_modifications = False
 
         # The names of the input files.
         self._rst_file = _os.path.join(str(self._work_dir), f"{name}.rst7")
@@ -340,6 +343,7 @@ def _setup(self):
                     zero_dummy_impropers=self._zero_dummy_impropers,
                     property_map=self._property_map,
                     perturbation_type=self._protocol.getPerturbationType(),
+                    no_dummy_modifications=self._no_dummy_modifications,
                 )
 
                 self._input_files.append(self._pert_file)
@@ -922,6 +926,7 @@ def _to_pert_file(
     print_all_atoms=False,
     property_map={},
     perturbation_type="full",
+    no_dummy_modifications=False,
 ):
     """
     Write a perturbation file for a perturbable molecule.
@@ -961,6 +966,10 @@ def _to_pert_file(
         "grow_soft" : Perturb all growing soft atom LJ terms (i.e. 0.0->value).
         "charge_soft" : Perturb all charging soft atom LJ terms (i.e. 0.0->value).
 
+    no_dummy_modifications : bool
+        Whether to skip modifications to dummy atoms. This is useful when
+        modifications to dummy atoms have already been applied.
+
     Returns
     -------
 
@@ -997,6 +1006,9 @@ def _to_pert_file(
     if not isinstance(perturbation_type, str):
         raise TypeError("'perturbation_type' must be of type 'str'")
 
+    if not isinstance(no_dummy_modifications, bool):
+        raise TypeError("'no_dummy_modifications' must be of type 'bool'")
+
     # Convert to lower case and strip whitespace.
     perturbation_type = perturbation_type.lower().replace(" ", "")
 
@@ -1973,19 +1985,22 @@ def sort_angles(angles, idx):
                 initial_dummy = _has_dummy(mol, [idx0, idx1, idx2])
                 final_dummy = _has_dummy(mol, [idx0, idx1, idx2], True)
 
-                # Set the angle parameters of the dummy state to those of the non-dummy end state.
-                if initial_dummy and final_dummy:
-                    has_dummy = True
-                    amber_angle0 = _SireMM.AmberAngle()
-                    amber_angle1 = _SireMM.AmberAngle()
-                elif initial_dummy or final_dummy:
-                    has_dummy = True
-                    if initial_dummy:
-                        amber_angle0 = amber_angle1
+                # Modifications to dummy states.
+                if not no_dummy_modifications:
+                    # Set the angle parameters of the dummy state to those of
+                    # the non-dummy end state.
+                    if initial_dummy and final_dummy:
+                        has_dummy = True
+                        amber_angle0 = _SireMM.AmberAngle()
+                        amber_angle1 = _SireMM.AmberAngle()
+                    elif initial_dummy or final_dummy:
+                        has_dummy = True
+                        if initial_dummy:
+                            amber_angle0 = amber_angle1
+                        else:
+                            amber_angle1 = amber_angle0
                     else:
-                        amber_angle1 = amber_angle0
-                else:
-                    has_dummy = False
+                        has_dummy = False
 
                 # Only write record if the angle parameters change.
                 if has_dummy or amber_angle0 != amber_angle1:
@@ -2313,28 +2328,30 @@ def sort_dihedrals(dihedrals, idx):
                 all_dummy_initial = all(_is_dummy(mol, [idx0, idx1, idx2, idx3]))
                 all_dummy_final = all(_is_dummy(mol, [idx0, idx1, idx2, idx3], True))
 
-                # Dummies are present in both end states, use null potentials.
-                if has_dummy_initial and has_dummy_final:
-                    amber_dihedral0 = _SireMM.AmberDihedral()
-                    amber_dihedral1 = _SireMM.AmberDihedral()
-
-                # Dummies in the initial state.
-                elif has_dummy_initial:
-                    if all_dummy_initial and not zero_dummy_dihedrals:
-                        # Use the potential at lambda = 1 and write to the pert file.
-                        amber_dihedral0 = amber_dihedral1
-                        force_write = True
-                    else:
-                        zero_k = True
-
-                # Dummies in the final state.
-                elif has_dummy_final:
-                    if all_dummy_final and not zero_dummy_dihedrals:
-                        # Use the potential at lambda = 0 and write to the pert file.
-                        amber_dihedral1 = amber_dihedral0
-                        force_write = True
-                    else:
-                        zero_k = True
+                # Perform dummy modifications if required.
+                if not no_dummy_modifications:
+                    # Dummies are present in both end states, use null potentials.
+                    if has_dummy_initial and has_dummy_final:
+                        amber_dihedral0 = _SireMM.AmberDihedral()
+                        amber_dihedral1 = _SireMM.AmberDihedral()
+
+                    # Dummies in the initial state.
+                    elif has_dummy_initial:
+                        if all_dummy_initial and not zero_dummy_dihedrals:
+                            # Use the potential at lambda = 1 and write to the pert file.
+                            amber_dihedral0 = amber_dihedral1
+                            force_write = True
+                        else:
+                            zero_k = True
+
+                    # Dummies in the final state.
+                    elif has_dummy_final:
+                        if all_dummy_final and not zero_dummy_dihedrals:
+                            # Use the potential at lambda = 0 and write to the pert file.
+                            amber_dihedral1 = amber_dihedral0
+                            force_write = True
+                        else:
+                            zero_k = True
 
                 # Only write record if the dihedral parameters change.
                 if zero_k or force_write or amber_dihedral0 != amber_dihedral1:
@@ -2711,28 +2728,30 @@ def sort_dihedrals(dihedrals, idx):
                 all_dummy_initial = all(_is_dummy(mol, [idx0, idx1, idx2, idx3]))
                 all_dummy_final = all(_is_dummy(mol, [idx0, idx1, idx2, idx3], True))
 
-                # Dummies are present in both end states, use null potentials.
-                if has_dummy_initial and has_dummy_final:
-                    amber_dihedral0 = _SireMM.AmberDihedral()
-                    amber_dihedral1 = _SireMM.AmberDihedral()
-
-                # Dummies in the initial state.
-                elif has_dummy_initial:
-                    if all_dummy_initial and not zero_dummy_impropers:
-                        # Use the potential at lambda = 1 and write to the pert file.
-                        amber_dihedral0 = amber_dihedral1
-                        force_write = True
-                    else:
-                        zero_k = True
-
-                # Dummies in the final state.
-                elif has_dummy_final:
-                    if all_dummy_final and not zero_dummy_impropers:
-                        # Use the potential at lambda = 0 and write to the pert file.
-                        amber_dihedral1 = amber_dihedral0
-                        force_write = True
-                    else:
-                        zero_k = True
+                # Perform dummy modifications if required.
+                if not no_dummy_modifications:
+                    # Dummies are present in both end states, use null potentials.
+                    if has_dummy_initial and has_dummy_final:
+                        amber_dihedral0 = _SireMM.AmberDihedral()
+                        amber_dihedral1 = _SireMM.AmberDihedral()
+
+                    # Dummies in the initial state.
+                    elif has_dummy_initial:
+                        if all_dummy_initial and not zero_dummy_impropers:
+                            # Use the potential at lambda = 1 and write to the pert file.
+                            amber_dihedral0 = amber_dihedral1
+                            force_write = True
+                        else:
+                            zero_k = True
+
+                    # Dummies in the final state.
+                    elif has_dummy_final:
+                        if all_dummy_final and not zero_dummy_impropers:
+                            # Use the potential at lambda = 0 and write to the pert file.
+                            amber_dihedral1 = amber_dihedral0
+                            force_write = True
+                        else:
+                            zero_k = True
 
                 # Only write record if the improper parameters change.
                 if zero_k or force_write or amber_dihedral0 != amber_dihedral1:
diff --git a/python/BioSimSpace/Sandpit/Exscientia/Node/_node.py b/python/BioSimSpace/Sandpit/Exscientia/Node/_node.py
index 0a244da0..3d63cfda 100644
--- a/python/BioSimSpace/Sandpit/Exscientia/Node/_node.py
+++ b/python/BioSimSpace/Sandpit/Exscientia/Node/_node.py
@@ -22,6 +22,7 @@
 from glob import glob as _glob
 
 from .._Utils import _try_import
+
 from .. import _Utils
 
 import os as _os
@@ -35,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():
@@ -89,7 +90,7 @@ def help(name):
     print(proc.stdout)
 
 
-def run(name, args={}):
+def run(name, args={}, work_dir=None):
     """
     Run a node.
 
@@ -102,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
     -------
 
@@ -111,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
 
@@ -122,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)
+    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)
 
-        # Create the command.
-        command = "%s/python %s --config input.yaml" % (
-            _SireBase.getBinDir(),
-            full_name,
-        )
-
-    # 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):
@@ -179,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
diff --git a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py
index dade5ac5..66834bbc 100644
--- a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py
+++ b/python/BioSimSpace/Sandpit/Exscientia/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
diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py
index fdcca27c..cce3479d 100644
--- a/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py
+++ b/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py
@@ -736,7 +736,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
diff --git a/python/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py b/python/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py
index 7464f352..56f5a931 100644
--- a/python/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py
+++ b/python/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py
@@ -24,12 +24,13 @@
 __author__ = "Lester Hedges"
 __email__ = "lester.hedges@gmail.com"
 
-__all__ = ["solvate", "spc", "spce", "tip3p", "tip4p", "tip5p", "waterModels"]
+__all__ = ["solvate", "spc", "spce", "tip3p", "tip4p", "tip5p", "opc", "waterModels"]
 
 import os as _os
 import re as _re
 import subprocess as _subprocess
 import shlex as _shlex
+import shutil as _shutil
 import sys as _sys
 import warnings as _warnings
 
@@ -658,6 +659,108 @@ def tip5p(
     )
 
 
+def opc(
+    molecule=None,
+    box=None,
+    angles=3 * [_Angle(90, "degrees")],
+    shell=None,
+    ion_conc=0,
+    is_neutral=True,
+    is_aligned=False,
+    match_water=True,
+    work_dir=None,
+    property_map={},
+):
+    """
+    Add OPC solvent.
+
+    Parameters
+    ----------
+
+    molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`, \
+               :class:`Molecule <BioSimSpace._SireWrappers.Molecules>`, \
+               :class:`System <BioSimSpace._SireWrappers.System>`
+        A molecule, or container/system of molecules.
+
+    box : [:class:`Length <BioSimSpace.Types.Length>`]
+        A list containing the box size in each dimension.
+
+    angles : [:class:`Angle <BioSimSpace.Types.Angle>`]
+        A list containing the angles between the box vectors: yz, xz, and xy.
+
+    shell : :class:`Length` <BioSimSpace.Types.Length>`
+        Thickness of the water shell around the solute. Note that the
+        base length of the resulting box must be at least twice as large
+        as the cutoff used by the chosen molecular dynamics engine. As such,
+        the shell option is often unsuitable for small molecules.
+
+    ion_conc : float
+        The ion concentration in (mol per litre).
+
+    is_neutral : bool
+        Whether to neutralise the system.
+
+    is_aligned : bool
+        Whether to align the principal axes of the molecule to those of the
+        solvent box.
+
+    match_water : bool
+        Whether to update the naming of existing water molecules to match the
+        expected convention for GROMACS, which is used as the solvation engine.
+
+    work_dir : str
+        The working directory for the process.
+
+    property_map : dict
+        A dictionary that maps system "properties" to their user defined
+        values. This allows the user to refer to properties with their
+        own naming scheme, e.g. { "charge" : "my-charge" }
+
+    Returns
+    -------
+
+    system : :class:`System <BioSimSpace._SireWrappers.System>`
+        The solvated molecular system.
+    """
+
+    if _gmx_exe is None:
+        raise _MissingSoftwareError(
+            "'BioSimSpace.Solvent.opc' is not supported. "
+            "Please install GROMACS (http://www.gromacs.org)."
+        )
+
+    # Validate arguments.
+    molecule, box, angles, shell, work_dir, property_map = _validate_input(
+        "opc",
+        molecule,
+        box,
+        angles,
+        shell,
+        ion_conc,
+        is_neutral,
+        is_aligned,
+        match_water,
+        work_dir,
+        property_map,
+    )
+
+    # Return the solvated system.
+    return _solvate(
+        molecule,
+        box,
+        angles,
+        shell,
+        "opc",
+        4,
+        ion_conc,
+        is_neutral,
+        is_aligned,
+        match_water,
+        work_dir=work_dir,
+        property_map=property_map,
+    )
+
+
 def _validate_input(
     model,
     molecule,
@@ -1094,12 +1197,22 @@ def _solvate(
             mod = model
         command = "%s solvate -cs %s" % (_gmx_exe, mod)
 
+        # OPC is a special case using a local structure file.
+        if mod == "opc":
+            command += ".gro"
+
         # Add the shell information.
         if molecule is not None and shell is not None:
             command += " -shell %f" % shell.nanometers().value()
 
         command += " -cp box.gro -o output.gro"
 
+        if mod == "opc":
+            template = _SireBase.getShareDir() + "/templates/water/opc"
+            _shutil.copyfile(template + ".itp", "opc.top")
+            _shutil.copyfile(template + "_solvate.gro", "opc.gro")
+            command += " -p opc.top"
+
         with open("README.txt", "a") as file:
             # Write the command to file.
             file.write("\n# gmx solvate was run with the following command:\n")
@@ -1151,14 +1264,26 @@ def _solvate(
                 "the 'box' size or 'shell' thickness."
             )
 
+        # Delete the defaults section and final line from the OPC topology file.
+        if model == "opc":
+            with open("opc.top", "r") as file:
+                lines = file.readlines()
+            with open("opc.top", "w") as file:
+                for line in lines[6:-1]:
+                    file.write(line)
+
         # Create a TOP file for the water model. By default we use the Amber03
         # force field to generate a dummy topology for the water model.
         with open("water_ions.top", "w") as file:
             file.write("#define FLEXIBLE 1\n\n")
             file.write("; Include AmberO3 force field\n")
             file.write('#include "amber03.ff/forcefield.itp"\n\n')
+            # Special handling for OPC, which uses a local topology file.
             file.write("; Include %s water topology\n" % model.upper())
-            file.write('#include "amber03.ff/%s.itp"\n\n' % model)
+            if model == "opc":
+                file.write('#include "opc.top"\n\n')
+            else:
+                file.write('#include "amber03.ff/%s.itp"\n\n' % model)
             file.write("; Include ions\n")
             file.write('#include "amber03.ff/ions.itp"\n\n')
             file.write("[ system ] \n")
@@ -1389,8 +1514,12 @@ def _solvate(
                             file.write("#define FLEXIBLE 1\n\n")
                             file.write("; Include AmberO3 force field\n")
                             file.write('#include "amber03.ff/forcefield.itp"\n\n')
+                            # Special handling for OPC, which uses a local topology file.
                             file.write("; Include %s water topology\n" % model.upper())
-                            file.write('#include "amber03.ff/%s.itp"\n\n' % model)
+                            if model == "opc":
+                                file.write('#include "opc.top"\n\n')
+                            else:
+                                file.write('#include "amber03.ff/%s.itp"\n\n' % model)
                             file.write("; Include ions\n")
                             file.write('#include "amber03.ff/ions.itp"\n\n')
                             file.write("[ system ] \n")
diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py
index e3733017..fc5dc4ca 100644
--- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py
+++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py
@@ -2522,7 +2522,17 @@ def _set_water_topology(self, format, is_crystal=False, property_map={}):
                     # TODO: Assume TIP3P. Not sure how to detect SPC/E at present.
                     water_model = "TIP3P"
                 elif num_point == 4:
-                    water_model = "TIP4P"
+                    # Check for OPC water.
+                    try:
+                        if (
+                            waters[0].search("element Xx").atoms()[0].charge().value()
+                            < -1.1
+                        ):
+                            water_model = "OPC"
+                        else:
+                            water_model = "TIP4P"
+                    except:
+                        water_model = "TIP4P"
                 elif num_point == 5:
                     water_model = "TIP5P"
 
diff --git a/python/BioSimSpace/Solvent/_solvent.py b/python/BioSimSpace/Solvent/_solvent.py
index 7464f352..56f5a931 100644
--- a/python/BioSimSpace/Solvent/_solvent.py
+++ b/python/BioSimSpace/Solvent/_solvent.py
@@ -24,12 +24,13 @@
 __author__ = "Lester Hedges"
 __email__ = "lester.hedges@gmail.com"
 
-__all__ = ["solvate", "spc", "spce", "tip3p", "tip4p", "tip5p", "waterModels"]
+__all__ = ["solvate", "spc", "spce", "tip3p", "tip4p", "tip5p", "opc", "waterModels"]
 
 import os as _os
 import re as _re
 import subprocess as _subprocess
 import shlex as _shlex
+import shutil as _shutil
 import sys as _sys
 import warnings as _warnings
 
@@ -658,6 +659,108 @@ def tip5p(
     )
 
 
+def opc(
+    molecule=None,
+    box=None,
+    angles=3 * [_Angle(90, "degrees")],
+    shell=None,
+    ion_conc=0,
+    is_neutral=True,
+    is_aligned=False,
+    match_water=True,
+    work_dir=None,
+    property_map={},
+):
+    """
+    Add OPC solvent.
+
+    Parameters
+    ----------
+
+    molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`, \
+               :class:`Molecule <BioSimSpace._SireWrappers.Molecules>`, \
+               :class:`System <BioSimSpace._SireWrappers.System>`
+        A molecule, or container/system of molecules.
+
+    box : [:class:`Length <BioSimSpace.Types.Length>`]
+        A list containing the box size in each dimension.
+
+    angles : [:class:`Angle <BioSimSpace.Types.Angle>`]
+        A list containing the angles between the box vectors: yz, xz, and xy.
+
+    shell : :class:`Length` <BioSimSpace.Types.Length>`
+        Thickness of the water shell around the solute. Note that the
+        base length of the resulting box must be at least twice as large
+        as the cutoff used by the chosen molecular dynamics engine. As such,
+        the shell option is often unsuitable for small molecules.
+
+    ion_conc : float
+        The ion concentration in (mol per litre).
+
+    is_neutral : bool
+        Whether to neutralise the system.
+
+    is_aligned : bool
+        Whether to align the principal axes of the molecule to those of the
+        solvent box.
+
+    match_water : bool
+        Whether to update the naming of existing water molecules to match the
+        expected convention for GROMACS, which is used as the solvation engine.
+
+    work_dir : str
+        The working directory for the process.
+
+    property_map : dict
+        A dictionary that maps system "properties" to their user defined
+        values. This allows the user to refer to properties with their
+        own naming scheme, e.g. { "charge" : "my-charge" }
+
+    Returns
+    -------
+
+    system : :class:`System <BioSimSpace._SireWrappers.System>`
+        The solvated molecular system.
+    """
+
+    if _gmx_exe is None:
+        raise _MissingSoftwareError(
+            "'BioSimSpace.Solvent.opc' is not supported. "
+            "Please install GROMACS (http://www.gromacs.org)."
+        )
+
+    # Validate arguments.
+    molecule, box, angles, shell, work_dir, property_map = _validate_input(
+        "opc",
+        molecule,
+        box,
+        angles,
+        shell,
+        ion_conc,
+        is_neutral,
+        is_aligned,
+        match_water,
+        work_dir,
+        property_map,
+    )
+
+    # Return the solvated system.
+    return _solvate(
+        molecule,
+        box,
+        angles,
+        shell,
+        "opc",
+        4,
+        ion_conc,
+        is_neutral,
+        is_aligned,
+        match_water,
+        work_dir=work_dir,
+        property_map=property_map,
+    )
+
+
 def _validate_input(
     model,
     molecule,
@@ -1094,12 +1197,22 @@ def _solvate(
             mod = model
         command = "%s solvate -cs %s" % (_gmx_exe, mod)
 
+        # OPC is a special case using a local structure file.
+        if mod == "opc":
+            command += ".gro"
+
         # Add the shell information.
         if molecule is not None and shell is not None:
             command += " -shell %f" % shell.nanometers().value()
 
         command += " -cp box.gro -o output.gro"
 
+        if mod == "opc":
+            template = _SireBase.getShareDir() + "/templates/water/opc"
+            _shutil.copyfile(template + ".itp", "opc.top")
+            _shutil.copyfile(template + "_solvate.gro", "opc.gro")
+            command += " -p opc.top"
+
         with open("README.txt", "a") as file:
             # Write the command to file.
             file.write("\n# gmx solvate was run with the following command:\n")
@@ -1151,14 +1264,26 @@ def _solvate(
                 "the 'box' size or 'shell' thickness."
             )
 
+        # Delete the defaults section and final line from the OPC topology file.
+        if model == "opc":
+            with open("opc.top", "r") as file:
+                lines = file.readlines()
+            with open("opc.top", "w") as file:
+                for line in lines[6:-1]:
+                    file.write(line)
+
         # Create a TOP file for the water model. By default we use the Amber03
         # force field to generate a dummy topology for the water model.
         with open("water_ions.top", "w") as file:
             file.write("#define FLEXIBLE 1\n\n")
             file.write("; Include AmberO3 force field\n")
             file.write('#include "amber03.ff/forcefield.itp"\n\n')
+            # Special handling for OPC, which uses a local topology file.
             file.write("; Include %s water topology\n" % model.upper())
-            file.write('#include "amber03.ff/%s.itp"\n\n' % model)
+            if model == "opc":
+                file.write('#include "opc.top"\n\n')
+            else:
+                file.write('#include "amber03.ff/%s.itp"\n\n' % model)
             file.write("; Include ions\n")
             file.write('#include "amber03.ff/ions.itp"\n\n')
             file.write("[ system ] \n")
@@ -1389,8 +1514,12 @@ def _solvate(
                             file.write("#define FLEXIBLE 1\n\n")
                             file.write("; Include AmberO3 force field\n")
                             file.write('#include "amber03.ff/forcefield.itp"\n\n')
+                            # Special handling for OPC, which uses a local topology file.
                             file.write("; Include %s water topology\n" % model.upper())
-                            file.write('#include "amber03.ff/%s.itp"\n\n' % model)
+                            if model == "opc":
+                                file.write('#include "opc.top"\n\n')
+                            else:
+                                file.write('#include "amber03.ff/%s.itp"\n\n' % model)
                             file.write("; Include ions\n")
                             file.write('#include "amber03.ff/ions.itp"\n\n')
                             file.write("[ system ] \n")
diff --git a/python/BioSimSpace/_Config/_amber.py b/python/BioSimSpace/_Config/_amber.py
index 9c64d359..1d34d4ba 100644
--- a/python/BioSimSpace/_Config/_amber.py
+++ b/python/BioSimSpace/_Config/_amber.py
@@ -240,7 +240,7 @@ def createConfig(
                 # Convert to a squashed representation, if needed
                 if isinstance(self._protocol, _FreeEnergyMixin):
                     atom_mapping0 = _squashed_atom_mapping(
-                        self.system, is_lambda1=False
+                        self._system, is_lambda1=False
                     )
                     atom_mapping1 = _squashed_atom_mapping(
                         self._system, is_lambda1=True
diff --git a/python/BioSimSpace/_Config/_somd.py b/python/BioSimSpace/_Config/_somd.py
index 308c8d69..417b70f6 100644
--- a/python/BioSimSpace/_Config/_somd.py
+++ b/python/BioSimSpace/_Config/_somd.py
@@ -32,6 +32,7 @@
 from .. import Protocol as _Protocol
 from ..Protocol._free_energy_mixin import _FreeEnergyMixin
 from ..Protocol._position_restraint_mixin import _PositionRestraintMixin
+from .._Exceptions import IncompatibleError as _IncompatibleError
 
 from ._config import Config as _Config
 
diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py
index 3b8d2194..b548e126 100644
--- a/python/BioSimSpace/_SireWrappers/_system.py
+++ b/python/BioSimSpace/_SireWrappers/_system.py
@@ -2438,7 +2438,17 @@ def _set_water_topology(self, format, is_crystal=False, property_map={}):
                     # TODO: Assume TIP3P. Not sure how to detect SPC/E at present.
                     water_model = "TIP3P"
                 elif num_point == 4:
-                    water_model = "TIP4P"
+                    # Check for OPC water.
+                    try:
+                        if (
+                            waters[0].search("element Xx").atoms()[0].charge().value()
+                            < -1.1
+                        ):
+                            water_model = "OPC"
+                        else:
+                            water_model = "TIP4P"
+                    except:
+                        water_model = "TIP4P"
                 elif num_point == 5:
                     water_model = "TIP5P"
 
diff --git a/tests/Parameters/test_parameters.py b/tests/Parameters/test_parameters.py
index 9b2d68ed..8e159519 100644
--- a/tests/Parameters/test_parameters.py
+++ b/tests/Parameters/test_parameters.py
@@ -208,3 +208,18 @@ def test_broken_sdf_formal_charge():
     from math import isclose
 
     assert isclose(charge.value(), 0.0, abs_tol=1e-6)
+
+
+def test_ff19SB():
+    """
+    Test that the ff19SB force field can be used to parameterise a molecule.
+    """
+
+    # Load the molecule.
+    mol = BSS.IO.readMolecules(f"{url}/4LYT_Fixed.pdb.bz2")[0]
+
+    # Parameterise the molecule with ff19SB.
+    mol = BSS.Parameters.ff19SB(mol).getMolecule()
+
+    # Make sure the molecule has CMAP terms.
+    assert mol._sire_object.hasProperty("cmap")
diff --git a/tests/Sandpit/Exscientia/Parameters/test_parameters.py b/tests/Sandpit/Exscientia/Parameters/test_parameters.py
index e23f7f4e..dc95d428 100644
--- a/tests/Sandpit/Exscientia/Parameters/test_parameters.py
+++ b/tests/Sandpit/Exscientia/Parameters/test_parameters.py
@@ -213,3 +213,18 @@ def test_broken_sdf_formal_charge():
     from math import isclose
 
     assert isclose(charge.value(), 0.0, abs_tol=1e-6)
+
+
+def test_ff19SB():
+    """
+    Test that the ff19SB force field can be used to parameterise a molecule.
+    """
+
+    # Load the molecule.
+    mol = BSS.IO.readMolecules(f"{url}/4LYT_Fixed.pdb.bz2")[0]
+
+    # Parameterise the molecule with ff19SB.
+    mol = BSS.Parameters.ff19SB(mol).getMolecule()
+
+    # Make sure the molecule has CMAP terms.
+    assert mol._sire_object.hasProperty("cmap")