diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index a78abcd71..77dde7a7b 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -9,6 +9,25 @@ company supporting open-source development of fostering academic/industrial coll within the biomolecular simulation community. Our software is hosted via the `OpenBioSim` `GitHub `__ organisation. +`2023.4.1 `_ - Dec 14 2023 +------------------------------------------------------------------------------------------------- + +* Make sure ``match_water`` keyword argument is passed to specialised solvation functions (`#190 `__). +* Check perturbable molecules for velocities when combining molecules (`#192 `__). +* Make sure velocities are double counted when searching for velocity properties when combining molecules (`#197 `__). +* Remove redundant ``BioSimSpace.Types.Type.__ne__`` operator (`#201 `__). +* Minor internal updates due to Sire API fixes (`#203 `__). +* Fixed bug in the BSS Boresch restraint search code (`@fjclark `_) (`#204 `__). +* Fixed ``renumber`` option in :meth:`extract ` method (`#210 `__). +* Add workaround for fixing reconstruction of intrascale matrix in :func:`readPerturbableSystem ` function (`#210 `__). +* Remove incorrect ``try_import`` statement in metadynamics driver script and make sure that global parameters in OpenMM script are unique (`#217 `__). +* Ensure the existing trajectory backend is used when getting the number of trajectory frames from a running process (`#219 `__). +* Fixed setting of ``igb`` config parameter for PMEMD simulations (`@annamherz `_) (`#220 `__). +* Make sure AMBER restraint mask matches all hydrogen atoms (`#222 `__). +* Ensure all searches for disulphide bonds are convert to a ``SelectorBond`` object (`#224 `__). +* Fix injection of custom commands into ``LEaP`` script (`#226 `__). + + `2023.4.0 `_ - Oct 13 2023 ------------------------------------------------------------------------------------------------- diff --git a/python/BioSimSpace/IO/_io.py b/python/BioSimSpace/IO/_io.py index fce3193e0..705e8f9d6 100644 --- a/python/BioSimSpace/IO/_io.py +++ b/python/BioSimSpace/IO/_io.py @@ -1100,9 +1100,13 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}): # Flag that the molecule is perturbable. mol.setProperty("is_perturbable", _SireBase.wrap(True)) + # Get the two molecules. + mol0 = system0[idx]._sire_object + mol1 = system1[idx]._sire_object + # Add the molecule0 and molecule1 properties. - mol.setProperty("molecule0", system0[idx]._sire_object) - mol.setProperty("molecule1", system1[idx]._sire_object) + mol.setProperty("molecule0", mol0) + mol.setProperty("molecule1", mol1) # Get the connectivity property name. conn_prop = property_map.get("connectivity", "connectivity") @@ -1121,6 +1125,23 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}): mol = mol.removeProperty(conn_prop + "0").molecule() mol = mol.removeProperty(conn_prop + "1").molecule() + # Reconstruct the intrascale matrices using the GroTop parser. + intra0 = ( + _SireIO.GroTop(_Molecule(mol0).toSystem()._sire_object) + .toSystem()[0] + .property("intrascale") + ) + intra1 = ( + _SireIO.GroTop(_Molecule(mol1).toSystem()._sire_object) + .toSystem()[0] + .property("intrascale") + ) + + # Set the "intrascale" properties. + intrascale_prop = property_map.get("intrascale", "intrascale") + mol.setProperty(intrascale_prop + "0", intra0) + mol.setProperty(intrascale_prop + "1", intra0) + # Commit the changes. mol = _Molecule(mol.commit()) diff --git a/python/BioSimSpace/Metadynamics/_aux/metadynamics.py b/python/BioSimSpace/Metadynamics/_aux/metadynamics.py index 83fe58c91..c6c37be1a 100644 --- a/python/BioSimSpace/Metadynamics/_aux/metadynamics.py +++ b/python/BioSimSpace/Metadynamics/_aux/metadynamics.py @@ -28,10 +28,8 @@ USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from ..._Utils import _try_import - -mm = _try_import("openmm") -unit = _try_import("openmm.unit", "conda install -c conda-forge openmm") +import openmm as mm +from openmm import unit from collections import namedtuple from functools import reduce diff --git a/python/BioSimSpace/Parameters/_Protocol/_amber.py b/python/BioSimSpace/Parameters/_Protocol/_amber.py index d83d016aa..5b1601686 100644 --- a/python/BioSimSpace/Parameters/_Protocol/_amber.py +++ b/python/BioSimSpace/Parameters/_Protocol/_amber.py @@ -122,8 +122,8 @@ def __init__( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, property_map={}, @@ -155,14 +155,17 @@ def __init__( Run 'BioSimSpace.Solvent.waterModels()' to see the supported water models. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -218,34 +221,25 @@ def __init__( else: self._water_model = water_model - # Validate the custom parameter file list. - if custom_parameters is not None: - if not isinstance(custom_parameters, (list, tuple)): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") + # Validate the additional leap commands. + if pre_mol_commands is not None: + if not isinstance(pre_mol_commands, (list, tuple)): + raise TypeError("'pre_mol_commands must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in custom_parameters): + if not all(isinstance(x, str) for x in pre_mol_commands): raise TypeError( - "'custom_parameters' must be a 'list' of 'str' types." + "'pre_mol_commands' must be a 'list' of 'str' types." ) - for x in custom_parameters: - if not os.path.isfile(x): - raise ValueError(f"Custom parameter file does not exist: '{x}'") - - # Convert to absolute paths. - self._custom_parameters = [] - for x in enumerate(custom_parameters): - self._custom_parameters.append(_os.path.abspath(x)) - else: - self._custom_parameters = None - - # Validate the additional leap commands. - if leap_commands is not None: - if not isinstance(leap_commands, (list, tuple)): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + self._pre_mol_commands = pre_mol_commands + if post_mol_commands is not None: + if not isinstance(post_mol_commands, (list, tuple)): + raise TypeError("'post_mol_commands must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in leap_commands): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") - self._leap_commands = leap_commands + if not all(isinstance(x, str) for x in post_mol_commands): + raise TypeError( + "'post_mol_commands' must be a 'list' of 'str' types." + ) + self._post_mol_commands = post_mol_commands # Validate the bond records. if bonds is not None: @@ -273,7 +267,7 @@ def __init__( # Set the compatibility flags. self._tleap = True - if self._custom_parameters is not None or self._leap_commands is not None: + if self._pre_mol_commands is not None or self._post_mol_commands is not None: self._pdb2gmx = False def run(self, molecule, work_dir=None, queue=None): @@ -512,10 +506,10 @@ def _run_tleap(self, molecule, work_dir): file.write("source leaprc.water.tip4pew\n") else: file.write("source leaprc.water.%s\n" % self._water_model) - # Write custom parameters. - if self._custom_parameters is not None: - for param in self._custom_parameters: - file.write("%s\n" % param) + # Write pre-mol user leap commands. + if self._pre_mol_commands is not None: + for command in self._pre_mol_commands: + file.write("%s\n" % command) file.write("mol = loadPdb leap.pdb\n") # Add any disulphide bond records. for bond in disulphide_bonds: @@ -523,9 +517,9 @@ def _run_tleap(self, molecule, work_dir): # Add any additional bond records. for bond in pruned_bond_records: file.write("%s\n" % bond) - # Write user leap commands. - if self._leap_commands is not None: - for command in self._leap_commands: + # Write post-mol user leap commands. + if self._post_mol_commands is not None: + for command in self._post_mol_commands: file.write("%s\n" % command) file.write("saveAmberParm mol leap.top leap.crd\n") file.write("quit") @@ -725,7 +719,7 @@ def _get_disulphide_bonds( # Try searching for disulphide bonds. try: - disulphides = query(mol, property_map) + disulphides = query(mol, property_map).bonds() except: disulphides = [] diff --git a/python/BioSimSpace/Parameters/_parameters.py b/python/BioSimSpace/Parameters/_parameters.py index f62d5630e..945319d99 100644 --- a/python/BioSimSpace/Parameters/_parameters.py +++ b/python/BioSimSpace/Parameters/_parameters.py @@ -135,8 +135,8 @@ def _parameterise_amber_protein( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -172,14 +172,17 @@ def _parameterise_amber_protein( or lone oxygen atoms corresponding to structural (crystal) water molecules. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -241,8 +244,8 @@ def _parameterise_amber_protein( max_distance=max_distance, water_model=water_model, check_ions=True, - custom_parameters=custom_parameters, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, property_map=property_map, @@ -255,8 +258,8 @@ def _parameterise_amber_protein( tolerance=tolerance, water_model=water_model, max_distance=max_distance, - custom_parameters=custom_parameters, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, property_map=property_map, @@ -802,8 +805,8 @@ def _validate( max_distance=None, water_model=None, check_ions=False, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -839,14 +842,17 @@ def _validate( Whether to check for the presence of structural ions. This is only required when parameterising with protein force fields. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -905,22 +911,19 @@ def _validate( "Please choose a 'water_model' for the ion parameters." ) - if custom_parameters is not None: - if not isinstance(custom_parameters, (list, tuple)): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") + if pre_mol_commands is not None: + if not isinstance(pre_mol_commands, (list, tuple)): + raise TypeError("'pre_mol_commands' must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in custom_parameters): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") - for x in custom_parameters: - if not os.path.isfile(x): - raise ValueError(f"Custom parameter file does not exist: '{x}'") - - if leap_commands is not None: - if not isinstance(leap_commands, (list, tuple)): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + if not all(isinstance(x, str) for x in pre_mol_commands): + raise TypeError("'pre_mol_commands' must be a 'list' of 'str' types.") + + if post_mol_commands is not None: + if not isinstance(post_mol_commands, (list, tuple)): + raise TypeError("'post_mol_commands' must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in leap_commands): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + if not all(isinstance(x, str) for x in post_mol_commands): + raise TypeError("'post_mol_commands' must be a 'list' of 'str' types.") if bonds is not None: if molecule is None: @@ -981,7 +984,8 @@ def _function( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -1013,11 +1017,17 @@ def _function( or lone oxygen atoms corresponding to structural (crystal) water molecules. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands and can be used to, e.g., - load additional parameter files. When this option is set, we can no - longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. + + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -1052,7 +1062,8 @@ def _function( tolerance=tolerance, max_distance=max_distance, water_model=water_model, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, work_dir=work_dir, diff --git a/python/BioSimSpace/Process/_openmm.py b/python/BioSimSpace/Process/_openmm.py index 8a913b667..92a40be04 100644 --- a/python/BioSimSpace/Process/_openmm.py +++ b/python/BioSimSpace/Process/_openmm.py @@ -877,13 +877,13 @@ def _generate_config(self): self.addToConfig("\n# Upper wall.") self.addToConfig( - "upper_wall_rest = CustomCentroidBondForce(3, '(k/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" + "upper_wall_rest = CustomCentroidBondForce(3, '(k1/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" ) self.addToConfig("upper_wall_rest.addGroup(lig)") self.addToConfig("upper_wall_rest.addGroup(p1)") self.addToConfig("upper_wall_rest.addGroup(p2)") self.addToConfig("upper_wall_rest.addBond([0,1,2])") - self.addToConfig("upper_wall_rest.addGlobalParameter('k', k1)") + self.addToConfig("upper_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "upper_wall_rest.addGlobalParameter('upper_wall', upper_wall)" ) @@ -903,13 +903,13 @@ def _generate_config(self): ) self.addToConfig( - "dist_restraint = CustomCentroidBondForce(3, '(k/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" + "dist_restraint = CustomCentroidBondForce(3, '(k2/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" ) self.addToConfig("dist_restraint.addGroup(lig)") self.addToConfig("dist_restraint.addGroup(p1)") self.addToConfig("dist_restraint.addGroup(p2)") self.addToConfig("dist_restraint.addBond([0,1,2])") - self.addToConfig("dist_restraint.addGlobalParameter('k', k2)") + self.addToConfig("dist_restraint.addGlobalParameter('k2', k2)") self.addToConfig("dist_restraint.addGlobalParameter('a', wall_width)") self.addToConfig("dist_restraint.addGlobalParameter('b', beta_cent)") self.addToConfig("dist_restraint.addGlobalParameter('c', s_cent)") @@ -919,13 +919,13 @@ def _generate_config(self): self.addToConfig("\n# Lower wall.") self.addToConfig( - "lower_wall_rest = CustomCentroidBondForce(3, '(k/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" + "lower_wall_rest = CustomCentroidBondForce(3, '(k1/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" ) self.addToConfig("lower_wall_rest.addGroup(lig)") self.addToConfig("lower_wall_rest.addGroup(p1)") self.addToConfig("lower_wall_rest.addGroup(p2)") self.addToConfig("lower_wall_rest.addBond([0,1,2])") - self.addToConfig("lower_wall_rest.addGlobalParameter('k', k1)") + self.addToConfig("lower_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "lower_wall_rest.addGlobalParameter('lower_wall', lower_wall)" ) diff --git a/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py b/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py index fce3193e0..705e8f9d6 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py +++ b/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py @@ -1100,9 +1100,13 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}): # Flag that the molecule is perturbable. mol.setProperty("is_perturbable", _SireBase.wrap(True)) + # Get the two molecules. + mol0 = system0[idx]._sire_object + mol1 = system1[idx]._sire_object + # Add the molecule0 and molecule1 properties. - mol.setProperty("molecule0", system0[idx]._sire_object) - mol.setProperty("molecule1", system1[idx]._sire_object) + mol.setProperty("molecule0", mol0) + mol.setProperty("molecule1", mol1) # Get the connectivity property name. conn_prop = property_map.get("connectivity", "connectivity") @@ -1121,6 +1125,23 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}): mol = mol.removeProperty(conn_prop + "0").molecule() mol = mol.removeProperty(conn_prop + "1").molecule() + # Reconstruct the intrascale matrices using the GroTop parser. + intra0 = ( + _SireIO.GroTop(_Molecule(mol0).toSystem()._sire_object) + .toSystem()[0] + .property("intrascale") + ) + intra1 = ( + _SireIO.GroTop(_Molecule(mol1).toSystem()._sire_object) + .toSystem()[0] + .property("intrascale") + ) + + # Set the "intrascale" properties. + intrascale_prop = property_map.get("intrascale", "intrascale") + mol.setProperty(intrascale_prop + "0", intra0) + mol.setProperty(intrascale_prop + "1", intra0) + # Commit the changes. mol = _Molecule(mol.commit()) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Metadynamics/_aux/metadynamics.py b/python/BioSimSpace/Sandpit/Exscientia/Metadynamics/_aux/metadynamics.py index 83fe58c91..c6c37be1a 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Metadynamics/_aux/metadynamics.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Metadynamics/_aux/metadynamics.py @@ -28,10 +28,8 @@ USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from ..._Utils import _try_import - -mm = _try_import("openmm") -unit = _try_import("openmm.unit", "conda install -c conda-forge openmm") +import openmm as mm +from openmm import unit from collections import namedtuple from functools import reduce diff --git a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py index d83d016aa..5b1601686 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py @@ -122,8 +122,8 @@ def __init__( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, property_map={}, @@ -155,14 +155,17 @@ def __init__( Run 'BioSimSpace.Solvent.waterModels()' to see the supported water models. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -218,34 +221,25 @@ def __init__( else: self._water_model = water_model - # Validate the custom parameter file list. - if custom_parameters is not None: - if not isinstance(custom_parameters, (list, tuple)): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") + # Validate the additional leap commands. + if pre_mol_commands is not None: + if not isinstance(pre_mol_commands, (list, tuple)): + raise TypeError("'pre_mol_commands must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in custom_parameters): + if not all(isinstance(x, str) for x in pre_mol_commands): raise TypeError( - "'custom_parameters' must be a 'list' of 'str' types." + "'pre_mol_commands' must be a 'list' of 'str' types." ) - for x in custom_parameters: - if not os.path.isfile(x): - raise ValueError(f"Custom parameter file does not exist: '{x}'") - - # Convert to absolute paths. - self._custom_parameters = [] - for x in enumerate(custom_parameters): - self._custom_parameters.append(_os.path.abspath(x)) - else: - self._custom_parameters = None - - # Validate the additional leap commands. - if leap_commands is not None: - if not isinstance(leap_commands, (list, tuple)): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + self._pre_mol_commands = pre_mol_commands + if post_mol_commands is not None: + if not isinstance(post_mol_commands, (list, tuple)): + raise TypeError("'post_mol_commands must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in leap_commands): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") - self._leap_commands = leap_commands + if not all(isinstance(x, str) for x in post_mol_commands): + raise TypeError( + "'post_mol_commands' must be a 'list' of 'str' types." + ) + self._post_mol_commands = post_mol_commands # Validate the bond records. if bonds is not None: @@ -273,7 +267,7 @@ def __init__( # Set the compatibility flags. self._tleap = True - if self._custom_parameters is not None or self._leap_commands is not None: + if self._pre_mol_commands is not None or self._post_mol_commands is not None: self._pdb2gmx = False def run(self, molecule, work_dir=None, queue=None): @@ -512,10 +506,10 @@ def _run_tleap(self, molecule, work_dir): file.write("source leaprc.water.tip4pew\n") else: file.write("source leaprc.water.%s\n" % self._water_model) - # Write custom parameters. - if self._custom_parameters is not None: - for param in self._custom_parameters: - file.write("%s\n" % param) + # Write pre-mol user leap commands. + if self._pre_mol_commands is not None: + for command in self._pre_mol_commands: + file.write("%s\n" % command) file.write("mol = loadPdb leap.pdb\n") # Add any disulphide bond records. for bond in disulphide_bonds: @@ -523,9 +517,9 @@ def _run_tleap(self, molecule, work_dir): # Add any additional bond records. for bond in pruned_bond_records: file.write("%s\n" % bond) - # Write user leap commands. - if self._leap_commands is not None: - for command in self._leap_commands: + # Write post-mol user leap commands. + if self._post_mol_commands is not None: + for command in self._post_mol_commands: file.write("%s\n" % command) file.write("saveAmberParm mol leap.top leap.crd\n") file.write("quit") @@ -725,7 +719,7 @@ def _get_disulphide_bonds( # Try searching for disulphide bonds. try: - disulphides = query(mol, property_map) + disulphides = query(mol, property_map).bonds() except: disulphides = [] diff --git a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py index f62d5630e..945319d99 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py @@ -135,8 +135,8 @@ def _parameterise_amber_protein( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -172,14 +172,17 @@ def _parameterise_amber_protein( or lone oxygen atoms corresponding to structural (crystal) water molecules. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -241,8 +244,8 @@ def _parameterise_amber_protein( max_distance=max_distance, water_model=water_model, check_ions=True, - custom_parameters=custom_parameters, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, property_map=property_map, @@ -255,8 +258,8 @@ def _parameterise_amber_protein( tolerance=tolerance, water_model=water_model, max_distance=max_distance, - custom_parameters=custom_parameters, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, property_map=property_map, @@ -802,8 +805,8 @@ def _validate( max_distance=None, water_model=None, check_ions=False, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -839,14 +842,17 @@ def _validate( Whether to check for the presence of structural ions. This is only required when parameterising with protein force fields. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -905,22 +911,19 @@ def _validate( "Please choose a 'water_model' for the ion parameters." ) - if custom_parameters is not None: - if not isinstance(custom_parameters, (list, tuple)): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") + if pre_mol_commands is not None: + if not isinstance(pre_mol_commands, (list, tuple)): + raise TypeError("'pre_mol_commands' must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in custom_parameters): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") - for x in custom_parameters: - if not os.path.isfile(x): - raise ValueError(f"Custom parameter file does not exist: '{x}'") - - if leap_commands is not None: - if not isinstance(leap_commands, (list, tuple)): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + if not all(isinstance(x, str) for x in pre_mol_commands): + raise TypeError("'pre_mol_commands' must be a 'list' of 'str' types.") + + if post_mol_commands is not None: + if not isinstance(post_mol_commands, (list, tuple)): + raise TypeError("'post_mol_commands' must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in leap_commands): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + if not all(isinstance(x, str) for x in post_mol_commands): + raise TypeError("'post_mol_commands' must be a 'list' of 'str' types.") if bonds is not None: if molecule is None: @@ -981,7 +984,8 @@ def _function( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -1013,11 +1017,17 @@ def _function( or lone oxygen atoms corresponding to structural (crystal) water molecules. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands and can be used to, e.g., - load additional parameter files. When this option is set, we can no - longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. + + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -1052,7 +1062,8 @@ def _function( tolerance=tolerance, max_distance=max_distance, water_model=water_model, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, work_dir=work_dir, diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py index 374e6700e..07dc8c7a2 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py @@ -878,13 +878,13 @@ def _generate_config(self): self.addToConfig("\n# Upper wall.") self.addToConfig( - "upper_wall_rest = CustomCentroidBondForce(3, '(k/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" + "upper_wall_rest = CustomCentroidBondForce(3, '(k1/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" ) self.addToConfig("upper_wall_rest.addGroup(lig)") self.addToConfig("upper_wall_rest.addGroup(p1)") self.addToConfig("upper_wall_rest.addGroup(p2)") self.addToConfig("upper_wall_rest.addBond([0,1,2])") - self.addToConfig("upper_wall_rest.addGlobalParameter('k', k1)") + self.addToConfig("upper_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "upper_wall_rest.addGlobalParameter('upper_wall', upper_wall)" ) @@ -904,13 +904,13 @@ def _generate_config(self): ) self.addToConfig( - "dist_restraint = CustomCentroidBondForce(3, '(k/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" + "dist_restraint = CustomCentroidBondForce(3, '(k2/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" ) self.addToConfig("dist_restraint.addGroup(lig)") self.addToConfig("dist_restraint.addGroup(p1)") self.addToConfig("dist_restraint.addGroup(p2)") self.addToConfig("dist_restraint.addBond([0,1,2])") - self.addToConfig("dist_restraint.addGlobalParameter('k', k2)") + self.addToConfig("dist_restraint.addGlobalParameter('k2', k2)") self.addToConfig("dist_restraint.addGlobalParameter('a', wall_width)") self.addToConfig("dist_restraint.addGlobalParameter('b', beta_cent)") self.addToConfig("dist_restraint.addGlobalParameter('c', s_cent)") @@ -920,13 +920,13 @@ def _generate_config(self): self.addToConfig("\n# Lower wall.") self.addToConfig( - "lower_wall_rest = CustomCentroidBondForce(3, '(k/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" + "lower_wall_rest = CustomCentroidBondForce(3, '(k1/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" ) self.addToConfig("lower_wall_rest.addGroup(lig)") self.addToConfig("lower_wall_rest.addGroup(p1)") self.addToConfig("lower_wall_rest.addGroup(p2)") self.addToConfig("lower_wall_rest.addBond([0,1,2])") - self.addToConfig("lower_wall_rest.addGlobalParameter('k', k1)") + self.addToConfig("lower_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "lower_wall_rest.addGlobalParameter('lower_wall', lower_wall)" ) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py index e71e1369d..20ac4e3fe 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py @@ -278,7 +278,7 @@ def generateAmberConfig(self, extra_options=None, extra_lines=None): if restraint == "backbone": restraint_mask = "@CA,C,O,N" elif restraint == "heavy": - restraint_mask = "!:WAT & !@H" + restraint_mask = "!:WAT & !@H=" elif restraint == "all": restraint_mask = "!:WAT" diff --git a/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py b/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py index b49a05d7d..1dc7e7d24 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py @@ -921,9 +921,9 @@ def nFrames(self): The number of trajectory frames. """ - # First get the current MDTraj object. + # First get the current trajectory object using the existing backend. if self._process is not None and self._process.isRunning(): - self._trajectory = self.getTrajectory() + self._trajectory = self.getTrajectory(format=self._backend) # There is no trajectory. if self._trajectory is None: diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py index 05ca21681..60f43224d 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py @@ -340,15 +340,10 @@ def extract(self, indices, renumber=False, property_map={}): # Append to the updated atom index. indices_.append(_SireMol.AtomIdx(x)) - if renumber: - mol = self.copy() - else: - mol = self - # Extract a partial molecule. try: # Create an empty atom selection for this molecule. - selection = mol._sire_object.selection() + selection = self._sire_object.selection() selection.selectNone() # Add the atom indices to the selection. @@ -356,7 +351,7 @@ def extract(self, indices, renumber=False, property_map={}): selection.select(idx) partial_mol = ( - _SireMol.PartialMolecule(mol._sire_object, selection) + _SireMol.PartialMolecule(self._sire_object, selection) .extract() .molecule() ) @@ -371,7 +366,7 @@ def extract(self, indices, renumber=False, property_map={}): intrascale = property_map.get("intrascale", "intrascale") # Flag whether the molecule has an intrascale property. - has_intrascale = mol._sire_object.hasProperty(intrascale) + has_intrascale = self._sire_object.hasProperty(intrascale) # Remove the "intrascale" property, since this doesn't correspond to the # extracted molecule. @@ -409,6 +404,15 @@ def extract(self, indices, renumber=False, property_map={}): else: mol = Molecule(partial_mol) + # Keep the same MolNum. + if not renumber: + mol._sire_object = ( + mol._sire_object.edit() + .renumber(self._sire_object.number()) + .commit() + .molecule() + ) + return mol def molecule0(self): diff --git a/python/BioSimSpace/Trajectory/_trajectory.py b/python/BioSimSpace/Trajectory/_trajectory.py index b49a05d7d..1dc7e7d24 100644 --- a/python/BioSimSpace/Trajectory/_trajectory.py +++ b/python/BioSimSpace/Trajectory/_trajectory.py @@ -921,9 +921,9 @@ def nFrames(self): The number of trajectory frames. """ - # First get the current MDTraj object. + # First get the current trajectory object using the existing backend. if self._process is not None and self._process.isRunning(): - self._trajectory = self.getTrajectory() + self._trajectory = self.getTrajectory(format=self._backend) # There is no trajectory. if self._trajectory is None: diff --git a/python/BioSimSpace/_Config/_amber.py b/python/BioSimSpace/_Config/_amber.py index 744de459f..5941bc51a 100644 --- a/python/BioSimSpace/_Config/_amber.py +++ b/python/BioSimSpace/_Config/_amber.py @@ -178,7 +178,7 @@ def createConfig( protocol_dict["cut"] = "999." if is_pmemd: # Use vacuum generalised Born model. - self.addToConfig(" igb=6,") + protocol_dict["igb"] = "6" else: # Non-bonded cut-off. protocol_dict["cut"] = "8.0" @@ -212,7 +212,7 @@ def createConfig( if restraint == "backbone": restraint_mask = "@CA,C,O,N" elif restraint == "heavy": - restraint_mask = "!:WAT & !@H" + restraint_mask = "!:WAT & !@H=" elif restraint == "all": restraint_mask = "!:WAT" diff --git a/python/BioSimSpace/_SireWrappers/_molecule.py b/python/BioSimSpace/_SireWrappers/_molecule.py index d619f1d48..424d6196d 100644 --- a/python/BioSimSpace/_SireWrappers/_molecule.py +++ b/python/BioSimSpace/_SireWrappers/_molecule.py @@ -340,15 +340,10 @@ def extract(self, indices, renumber=False, property_map={}): # Append to the updated atom index. indices_.append(_SireMol.AtomIdx(x)) - if renumber: - mol = self.copy() - else: - mol = self - # Extract a partial molecule. try: # Create an empty atom selection for this molecule. - selection = mol._sire_object.selection() + selection = self._sire_object.selection() selection.selectNone() # Add the atom indices to the selection. @@ -356,7 +351,7 @@ def extract(self, indices, renumber=False, property_map={}): selection.select(idx) partial_mol = ( - _SireMol.PartialMolecule(mol._sire_object, selection) + _SireMol.PartialMolecule(self._sire_object, selection) .extract() .molecule() ) @@ -371,7 +366,7 @@ def extract(self, indices, renumber=False, property_map={}): intrascale = property_map.get("intrascale", "intrascale") # Flag whether the molecule has an intrascale property. - has_intrascale = mol._sire_object.hasProperty(intrascale) + has_intrascale = self._sire_object.hasProperty(intrascale) # Remove the "intrascale" property, since this doesn't correspond to the # extracted molecule. @@ -409,6 +404,15 @@ def extract(self, indices, renumber=False, property_map={}): else: mol = Molecule(partial_mol) + # Keep the same MolNum. + if not renumber: + mol._sire_object = ( + mol._sire_object.edit() + .renumber(self._sire_object.number()) + .commit() + .molecule() + ) + return mol def molecule0(self): diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index 0a5bf5157..cc41b3020 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -662,7 +662,7 @@ def addMolecules(self, molecules): _warnings.warn( "Failed to remove 'velocity' property from all molecules!" ) - if has_perturnable: + if has_perturbable: try: self._sire_object = _SireIO.removeProperty( self._sire_object, "velocity0" diff --git a/requirements.txt b/requirements.txt index 39fffe478..b501e7840 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ # BioSimSpace runtime requirements. # main -sire~=2023.4.0 +sire~=2023.4.2 # devel #sire==2023.5.0.dev diff --git a/tests/Parameters/test_parameters.py b/tests/Parameters/test_parameters.py index 55ab15afb..8acd41ff7 100644 --- a/tests/Parameters/test_parameters.py +++ b/tests/Parameters/test_parameters.py @@ -1,4 +1,5 @@ import pytest +import tempfile import BioSimSpace as BSS @@ -77,3 +78,62 @@ def test_molecule_rename(): # Make sure the name is correct. assert mol._sire_object.name().value() == "smiles:[C@@H](C(F)(F)F)(OC(F)F)Cl" + + +@pytest.mark.skipif( + has_antechamber is False or has_tleap is False, + reason="Requires AmberTools/antechamber and tLEaP to be installed.", +) +def test_leap_commands(molecule0): + """ + Test that custom commands are correctly inserted into the LEaP input + script. + """ + + # Create lists of pre- and post-commands. + pre_mol_commands = ["command1", "command2"] + post_mol_commands = ["command3", "command4"] + + with tempfile.TemporaryDirectory() as tmp: + # This will fail, but we only want to check the LEaP script. + try: + mol = BSS.Parameters.ff14SB( + molecule0, + work_dir=tmp, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, + ).getMolecule() + except: + pass + + # Load the script and check that the custom parameters are present and + # in the correct order. + with open(f"{tmp}/leap.txt", "r") as f: + script = f.readlines() + + # Create lists to store the indices of the custom commands. + line_pre = [-1 for _ in range(len(pre_mol_commands))] + line_post = [-1 for _ in range(len(post_mol_commands))] + + # Loop over the lines in the script and store the line numbers + # where the custom commands are found. + for x, line in enumerate(script): + for y, command in enumerate(pre_mol_commands): + if command in line: + line_pre[y] = x + for y, command in enumerate(post_mol_commands): + if command in line: + line_post[y] = x + + # Make sure the lines are found. + for line in line_pre: + assert line != -1 + for line in line_post: + assert line != -1 + + # Make sure the lines are in the correct order. + for x in range(len(line_pre) - 1): + assert line_pre[x] < line_pre[x + 1] + assert line_pre[-1] < line_post[0] + for x in range(len(line_post) - 1): + assert line_post[x] < line_post[x + 1] diff --git a/tests/Sandpit/Exscientia/Parameters/test_parameters.py b/tests/Sandpit/Exscientia/Parameters/test_parameters.py index db89a8025..f42f2d539 100644 --- a/tests/Sandpit/Exscientia/Parameters/test_parameters.py +++ b/tests/Sandpit/Exscientia/Parameters/test_parameters.py @@ -1,4 +1,5 @@ import pytest +import tempfile import BioSimSpace.Sandpit.Exscientia as BSS @@ -82,3 +83,62 @@ def test_molecule_rename(): # Make sure the name is correct. assert mol._sire_object.name().value() == "smiles:[C@@H](C(F)(F)F)(OC(F)F)Cl" + + +@pytest.mark.skipif( + has_antechamber is False or has_tleap is False, + reason="Requires AmberTools/antechamber and tLEaP to be installed.", +) +def test_leap_commands(molecule0): + """ + Test that custom commands are correctly inserted into the LEaP input + script. + """ + + # Create lists of pre- and post-commands. + pre_mol_commands = ["command1", "command2"] + post_mol_commands = ["command3", "command4"] + + with tempfile.TemporaryDirectory() as tmp: + # This will fail, but we only want to check the LEaP script. + try: + mol = BSS.Parameters.ff14SB( + molecule0, + work_dir=tmp, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, + ).getMolecule() + except: + pass + + # Load the script and check that the custom parameters are present and + # in the correct order. + with open(f"{tmp}/leap.txt", "r") as f: + script = f.readlines() + + # Create lists to store the indices of the custom commands. + line_pre = [-1 for _ in range(len(pre_mol_commands))] + line_post = [-1 for _ in range(len(post_mol_commands))] + + # Loop over the lines in the script and store the line numbers + # where the custom commands are found. + for x, line in enumerate(script): + for y, command in enumerate(pre_mol_commands): + if command in line: + line_pre[y] = x + for y, command in enumerate(post_mol_commands): + if command in line: + line_post[y] = x + + # Make sure the lines are found. + for line in line_pre: + assert line != -1 + for line in line_post: + assert line != -1 + + # Make sure the lines are in the correct order. + for x in range(len(line_pre) - 1): + assert line_pre[x] < line_pre[x + 1] + assert line_pre[-1] < line_post[0] + for x in range(len(line_post) - 1): + assert line_post[x] < line_post[x + 1] diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py b/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py index d19b385f8..7f6bc0aa3 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py +++ b/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py @@ -102,3 +102,27 @@ def test_hydrogen_mass_repartitioning(system, ignore_waters): # Assert the the masses are approximately the same. assert final_mass == pytest.approx(initial_mass) + + +def test_extract(system): + """Test the extract method.""" + + # A list of atom indices to extract. + idxs = [0, 1, 2, 3] + + # Extract the first molecule from the system. + mol = system[0] + + # Extract the atoms. + partial_mol = mol.extract(idxs) + + assert partial_mol.nAtoms() == 4 + + # Make sure the numbers match. + assert partial_mol.number() == mol.number() + + # Extract and renumber. + partial_mol = mol.extract(idxs, renumber=True) + + # Make sure the numbers are different. + assert partial_mol.number() != mol.number() diff --git a/tests/_SireWrappers/test_molecule.py b/tests/_SireWrappers/test_molecule.py index 867bac2c6..28b6c28df 100644 --- a/tests/_SireWrappers/test_molecule.py +++ b/tests/_SireWrappers/test_molecule.py @@ -96,3 +96,27 @@ def test_hydrogen_mass_repartitioning(system, ignore_waters): # Assert the the masses are approximately the same. assert final_mass == pytest.approx(initial_mass) + + +def test_extract(system): + """Test the extract method.""" + + # A list of atom indices to extract. + idxs = [0, 1, 2, 3] + + # Extract the first molecule from the system. + mol = system[0] + + # Extract the atoms. + partial_mol = mol.extract(idxs) + + assert partial_mol.nAtoms() == 4 + + # Make sure the numbers match. + assert partial_mol.number() == mol.number() + + # Extract and renumber. + partial_mol = mol.extract(idxs, renumber=True) + + # Make sure the numbers are different. + assert partial_mol.number() != mol.number()