Allow SOMD1 pert file to be created without modifications to ghost atom bonded terms

python/BioSimSpace/FreeEnergy/_relative.py

@@ -1085,7 +1085,7 @@ def _somd2_extract(parquet_file, T=None, estimator="MBAR"):
         df = table.to_pandas()
 
         if is_mbar:
-            # Extract the columns correspodning to the lambda array.
+            # Extract the columns corresponding to the lambda array.
             df = df[[x for x in lambda_array]]
 
             # Subtract the potential at the simulated lambda.

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,7 @@ def _checkPerturbable(self, system):
                     "in the 'property_map' argument."
                 )
             else:
-                is_lambda1 = self._property_map["is_lambda1"].value()
+                is_lambda1 = self._property_map["is_lambda1"]
                 self._property_map.pop("is_lambda1")
 
             # Loop over all perturbable molecules in the system and replace them

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:

python/BioSimSpace/Sandpit/Exscientia/Process/_process.py

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

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

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