Skip to content

Commit

Permalink
Preserve SMILES molecule properties when parameterising. [closes #330]
Browse files Browse the repository at this point in the history
  • Loading branch information
@lohedges
lohedges committed Aug 16, 2024
1 parent c837b81 commit afb6b05
Show file tree
Hide file tree
Showing 6 changed files with 192 additions and 206 deletions.
110 changes: 40 additions & 70 deletions python/BioSimSpace/Parameters/_Protocol/_amber.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -310,15 +310,24 @@ def run(self, molecule, work_dir=None, queue=None):
if work_dir is None:
work_dir = _os.getcwd()

# Flag whether the molecule is a SMILES string.
# Try to create a molecule from the SMILES string.
if isinstance(molecule, str):
is_smiles = True
else:
is_smiles = False
try:
smiles = molecule
molecule = _smiles(molecule)
edit_mol = molecule._sire_object.edit()
edit_mol = edit_mol.rename(f"smiles:{smiles}").molecule()
molecule._sire_object = edit_mol.commit()
except Exception as e:
msg = "Unable to convert SMILES to Molecule using RDKit."
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None

if not is_smiles:
# Create a copy of the molecule.
new_mol = molecule.copy()
# Create a copy of the molecule.
new_mol = molecule.copy()

# Choose the program to run with depending on the force field compatibility.
# If tLEaP and pdb2gmx are supported, default to tLEaP, then use pdb2gmx if
Expand All @@ -328,8 +337,7 @@ def run(self, molecule, work_dir=None, queue=None):
if self._tleap:
if _tleap_exe is not None:
output = self._run_tleap(molecule, str(work_dir))
if not is_smiles:
new_mol._ion_water_model = self._water_model
new_mol._ion_water_model = self._water_model
# Otherwise, try using pdb2gmx.
elif self._pdb2gmx:
if _gmx_exe is not None:
Expand Down Expand Up @@ -371,41 +379,23 @@ def run(self, molecule, work_dir=None, queue=None):
# Make the molecule 'mol' compatible with 'par_mol'. This will create
# a mapping between atom indices in the two molecules and add all of
# the new properties from 'par_mol' to 'mol'.
if is_smiles:
new_mol = par_mol

# We'll now add MolName and ResName info to the molecule, since
# this will be missing.

# Rename the molecule with the original SMILES string.
edit_mol = new_mol._sire_object.edit()
edit_mol = edit_mol.rename(molecule).molecule()

# Rename the residue LIG.
resname = _SireMol.ResName("LIG")
edit_mol = edit_mol.residue(_SireMol.ResIdx(0)).rename(resname).molecule()

# Commit the changes.
new_mol._sire_object = edit_mol.commit()

if self._ensure_compatible:
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)
else:
if self._ensure_compatible:
try:
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)
else:
try:
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)
except:
new_mol = par_mol
except:
new_mol = par_mol

# Record the forcefield used to parameterise the molecule.
new_mol._forcefield = self._forcefield
Expand Down Expand Up @@ -1009,9 +999,12 @@ def run(self, molecule, work_dir=None, queue=None):

# Convert SMILES to a molecule.
if isinstance(molecule, str):
is_smiles = True
try:
smiles = molecule
new_mol = _smiles(molecule)
edit_mol = new_mol._sire_object.edit()
edit_mol = edit_mol.rename(f"smiles:{smiles}").molecule()
new_mol._sire_object = edit_mol.commit()
except Exception as e:
msg = "Unable to convert SMILES to Molecule using RDKit."
if _isVerbose():
Expand All @@ -1020,7 +1013,6 @@ def run(self, molecule, work_dir=None, queue=None):
else:
raise _ThirdPartyError(msg) from None
else:
is_smiles = False
new_mol = molecule.copy()

# Use the net molecular charge passed as an option.
Expand Down Expand Up @@ -1251,45 +1243,23 @@ def run(self, molecule, work_dir=None, queue=None):
# Make the molecule 'mol' compatible with 'par_mol'. This will create
# a mapping between atom indices in the two molecules and add all of
# the new properties from 'par_mol' to 'mol'.
if is_smiles:
new_mol = par_mol

# We'll now add MolName and ResName info to the molecule, since
# this will be missing.

# Rename the molecule with the original SMILES string.
edit_mol = new_mol._sire_object.edit()
edit_mol = edit_mol.rename(molecule).molecule()

# Rename the residue LIG.
resname = _SireMol.ResName("LIG")
edit_mol = (
edit_mol.residue(_SireMol.ResIdx(0))
.rename(resname)
.molecule()
if self._ensure_compatible:
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)

# Commit the changes.
new_mol._sire_object = edit_mol.commit()

else:
if self._ensure_compatible:
try:
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)
else:
try:
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)
except:
new_mol = par_mol
except:
new_mol = par_mol

# Record the forcefield used to parameterise the molecule.
new_mol._forcefield = ff
Expand Down
59 changes: 26 additions & 33 deletions python/BioSimSpace/Parameters/_Protocol/_openforcefield.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,7 @@

from ... import _isVerbose
from ... import IO as _IO
from ...Convert import smiles as _smiles
from ..._Exceptions import IncompatibleError as _IncompatibleError
from ..._Exceptions import ThirdPartyError as _ThirdPartyError
from ..._SireWrappers import Molecule as _Molecule
Expand Down Expand Up @@ -415,48 +416,40 @@ def run(self, molecule, work_dir=None, queue=None):

# Make the parameterised molecule compatible with the original topology.
if is_smiles:
new_mol = par_mol

# We'll now add MolName and ResName info to the molecule, since
# this will be missing.

# Rename the molecule with the original SMILES string.
# Since the name is written to topology file formats, we
# need to ensure that it doesn't start with an [ character,
# which would break GROMACS.
name = molecule
if name.startswith("["):
name = f"smiles:{name}"

edit_mol = new_mol._sire_object.edit()
edit_mol = edit_mol.rename(name).molecule()

# Rename the residue LIG.
resname = _SireMol.ResName("LIG")
edit_mol = edit_mol.residue(_SireMol.ResIdx(0)).rename(resname).molecule()

# Commit the changes.
new_mol._sire_object = edit_mol.commit()
# Try to create BioSimSpace molecule from the SMILES string.
try:
smiles = molecule
molecule = _smiles(molecule)
edit_mol = molecule._sire_object.edit()
edit_mol = edit_mol.rename(f"smiles:{smiles}").molecule()
molecule._sire_object = edit_mol.commit()
except Exception as e:
msg = "Unable to convert SMILES to Molecule using RDKit."
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None

if self._ensure_compatible:
new_mol = molecule.copy()
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)
else:
if self._ensure_compatible:
try:
new_mol = molecule.copy()
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)
else:
try:
new_mol.makeCompatibleWith(
par_mol,
property_map=self._property_map,
overwrite=True,
verbose=False,
)
except:
new_mol = par_mol
except:
new_mol = par_mol

# Record the forcefield used to parameterise the molecule.
new_mol._forcefield = self._forcefield
Expand Down
Loading

0 comments on commit afb6b05

Please sign in to comment.