From d32e415197163cc3320a8c002e61a7807f23beaf Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 27 May 2025 09:52:10 +0100 Subject: [PATCH 1/7] Allow pert file without modifications to dummy bonded terms. --- python/BioSimSpace/Process/_somd.py | 131 ++++++++++++++++------------ 1 file changed, 75 insertions(+), 56 deletions(-) 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: From 41534e618a612987f07241af9c8d6cf72918ebbf Mon Sep 17 00:00:00 2001 From: Matthew Date: Fri, 30 May 2025 13:52:42 +0100 Subject: [PATCH 2/7] Fix for is_lambda1 in Somd1 processes --- python/BioSimSpace/Process/_process.py | 2 +- python/BioSimSpace/Sandpit/Exscientia/Process/_process.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/BioSimSpace/Process/_process.py b/python/BioSimSpace/Process/_process.py index 020558e0..28437edd 100644 --- a/python/BioSimSpace/Process/_process.py +++ b/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 diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py index fdcca27c..47a09a82 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py +++ b/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 From db03a53130b748c715fa8c3ceb6321e1a2a0ef86 Mon Sep 17 00:00:00 2001 From: Matthew Date: Mon, 2 Jun 2025 12:04:20 +0100 Subject: [PATCH 3/7] Import proper errors in somd config --- python/BioSimSpace/_Config/_somd.py | 1 + 1 file changed, 1 insertion(+) 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 From 5d748ce9a25bcd51623363b3cdd0705ef36b84d0 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 2 Jun 2025 17:08:32 +0100 Subject: [PATCH 4/7] Typo. [ci skip] --- python/BioSimSpace/FreeEnergy/_relative.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/BioSimSpace/FreeEnergy/_relative.py b/python/BioSimSpace/FreeEnergy/_relative.py index 1ba69e79..b852a1a5 100644 --- a/python/BioSimSpace/FreeEnergy/_relative.py +++ b/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. From 44eeb7c643c582e66dc74d3bb0481174e2d29bca Mon Sep 17 00:00:00 2001 From: Matthew Date: Tue, 3 Jun 2025 11:29:52 +0100 Subject: [PATCH 5/7] Amber config typo fix [ci skip] --- python/BioSimSpace/_Config/_amber.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 580c91b3be4d7cfa6cf77a28bb4af0b1fbd3317b Mon Sep 17 00:00:00 2001 From: Matthew Date: Tue, 3 Jun 2025 15:10:36 +0100 Subject: [PATCH 6/7] Fix for AMBER runs on perturbable systems with restraints --- python/BioSimSpace/Process/_amber.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/python/BioSimSpace/Process/_amber.py b/python/BioSimSpace/Process/_amber.py index c585c670..6f5a35c8 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,6 +354,17 @@ def _setup(self, **kwargs): # Reference file for position restraints. try: + 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) file = _os.path.splitext(self._ref_file)[0] _IO.saveMolecules( file, reference_system, "rst7", property_map=self._property_map From 8fe807ad16b77356b7764fe10fe31f593835e50e Mon Sep 17 00:00:00 2001 From: Matthew Date: Tue, 3 Jun 2025 15:56:29 +0100 Subject: [PATCH 7/7] Fix for previous commit - removes old file save [ci skip] --- python/BioSimSpace/Process/_amber.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/python/BioSimSpace/Process/_amber.py b/python/BioSimSpace/Process/_amber.py index 6f5a35c8..44f20929 100644 --- a/python/BioSimSpace/Process/_amber.py +++ b/python/BioSimSpace/Process/_amber.py @@ -365,10 +365,6 @@ def _setup(self, **kwargs): ) else: _shutil.copy(self._rst_file, self._ref_file) - file = _os.path.splitext(self._ref_file)[0] - _IO.saveMolecules( - file, reference_system, "rst7", property_map=self._property_map - ) except Exception as e: msg = "Failed to write reference system to 'RST7' format." if _isVerbose():