Skip to content

Commit 41fb514

Browse files
lohedgeslohedges
committed
Presever SMILES molecule properties when parameterising. [closes #331]
1 parent e72f913 commit 41fb514

File tree

6 files changed

+248
-368
lines changed

6 files changed

+248
-368
lines changed

‎python/BioSimSpace/Parameters/_Protocol/_amber.py

Lines changed: 51 additions & 110 deletions
Original file line numberDiff line numberDiff line change
@@ -310,15 +310,24 @@ def run(self, molecule, work_dir=None, queue=None):
310310
if work_dir is None:
311311
work_dir = _os.getcwd()
312312

313-
# Flag whether the molecule is a SMILES string.
313+
# Try to create a molecule from the SMILES string.
314314
if isinstance(molecule, str):
315-
is_smiles = True
316-
else:
317-
is_smiles = False
315+
try:
316+
smiles = molecule
317+
molecule = _smiles(molecule)
318+
edit_mol = molecule._sire_object.edit()
319+
edit_mol = edit_mol.rename(f"smiles:{smiles}").molecule()
320+
molecule._sire_object = edit_mol.commit()
321+
except Exception as e:
322+
msg = "Unable to convert SMILES to Molecule using RDKit."
323+
if _isVerbose():
324+
msg += ": " + getattr(e, "message", repr(e))
325+
raise _ThirdPartyError(msg) from e
326+
else:
327+
raise _ThirdPartyError(msg) from None
318328

319-
if not is_smiles:
320-
# Create a copy of the molecule.
321-
new_mol = molecule.copy()
329+
# Create a copy of the molecule.
330+
new_mol = molecule.copy()
322331

323332
# Choose the program to run with depending on the force field compatibility.
324333
# If tLEaP and pdb2gmx are supported, default to tLEaP, then use pdb2gmx if
@@ -328,8 +337,7 @@ def run(self, molecule, work_dir=None, queue=None):
328337
if self._tleap:
329338
if _tleap_exe is not None:
330339
output = self._run_tleap(molecule, str(work_dir))
331-
if not is_smiles:
332-
new_mol._ion_water_model = self._water_model
340+
new_mol._ion_water_model = self._water_model
333341
# Otherwise, try using pdb2gmx.
334342
elif self._pdb2gmx:
335343
if _gmx_exe is not None:
@@ -371,41 +379,23 @@ def run(self, molecule, work_dir=None, queue=None):
371379
# Make the molecule 'mol' compatible with 'par_mol'. This will create
372380
# a mapping between atom indices in the two molecules and add all of
373381
# the new properties from 'par_mol' to 'mol'.
374-
if is_smiles:
375-
new_mol = par_mol
376-
377-
# We'll now add MolName and ResName info to the molecule, since
378-
# this will be missing.
379-
380-
# Rename the molecule with the original SMILES string.
381-
edit_mol = new_mol._sire_object.edit()
382-
edit_mol = edit_mol.rename(molecule).molecule()
383-
384-
# Rename the residue LIG.
385-
resname = _SireMol.ResName("LIG")
386-
edit_mol = edit_mol.residue(_SireMol.ResIdx(0)).rename(resname).molecule()
387-
388-
# Commit the changes.
389-
new_mol._sire_object = edit_mol.commit()
390-
382+
if self._ensure_compatible:
383+
new_mol.makeCompatibleWith(
384+
par_mol,
385+
property_map=self._property_map,
386+
overwrite=True,
387+
verbose=False,
388+
)
391389
else:
392-
if self._ensure_compatible:
390+
try:
393391
new_mol.makeCompatibleWith(
394392
par_mol,
395393
property_map=self._property_map,
396394
overwrite=True,
397395
verbose=False,
398396
)
399-
else:
400-
try:
401-
new_mol.makeCompatibleWith(
402-
par_mol,
403-
property_map=self._property_map,
404-
overwrite=True,
405-
verbose=False,
406-
)
407-
except:
408-
new_mol = par_mol
397+
except:
398+
new_mol = par_mol
409399

410400
# Record the forcefield used to parameterise the molecule.
411401
new_mol._forcefield = self._forcefield
@@ -421,33 +411,18 @@ def _run_tleap(self, molecule, work_dir):
421411
Parameters
422412
----------
423413
424-
molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`, str
425-
The molecule to parameterise, either as a Molecule object or SMILES
426-
string.
414+
molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
415+
The molecule to parameterise.
427416
428417
work_dir : str
429418
The working directory.
430419
"""
431420

432-
# Convert SMILES to a molecule.
433-
if isinstance(molecule, str):
434-
try:
435-
_molecule = _smiles(molecule)
436-
except Exception as e:
437-
msg = "Unable to convert SMILES to Molecule using RDKit."
438-
if _isVerbose():
439-
msg += ": " + getattr(e, "message", repr(e))
440-
raise _ThirdPartyError(msg) from e
441-
else:
442-
raise _ThirdPartyError(msg) from None
443-
else:
444-
_molecule = molecule
445-
446421
# Write the system to a PDB file.
447422
try:
448423
# LEaP expects residue numbering to be ascending and continuous.
449424
renumbered_molecule = _SireIO.renumberConstituents(
450-
_molecule.toSystem()._sire_object
425+
molecule.toSystem()._sire_object
451426
)[0]
452427
renumbered_molecule = _Molecule(renumbered_molecule)
453428
_IO.saveMolecules(
@@ -570,9 +545,8 @@ def _run_pdb2gmx(self, molecule, work_dir):
570545
Parameters
571546
----------
572547
573-
molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`, str
574-
The molecule to parameterise, either as a Molecule object or SMILES
575-
string.
548+
molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
549+
The molecule to parameterise.
576550
577551
work_dir : str
578552
The working directory.
@@ -587,25 +561,11 @@ def _run_pdb2gmx(self, molecule, work_dir):
587561
"'pdb2gmx' does not support the '%s' force field." % self._forcefield
588562
)
589563

590-
# Convert SMILES to a molecule.
591-
if isinstance(molecule, str):
592-
try:
593-
_molecule = _smiles(molecule)
594-
except Exception as e:
595-
msg = "Unable to convert SMILES to Molecule using RDKit."
596-
if _isVerbose():
597-
msg += ": " + getattr(e, "message", repr(e))
598-
raise _ThirdPartyError(msg) from e
599-
else:
600-
raise _ThirdPartyError(msg) from None
601-
else:
602-
_molecule = molecule
603-
604564
# Write the system to a PDB file.
605565
try:
606566
_IO.saveMolecules(
607567
_os.path.join(str(work_dir), "input"),
608-
_molecule,
568+
molecule,
609569
"pdb",
610570
property_map=self._property_map,
611571
)
@@ -1007,21 +967,24 @@ def run(self, molecule, work_dir=None, queue=None):
1007967
if work_dir is None:
1008968
work_dir = _os.getcwd()
1009969

1010-
# Convert SMILES to a molecule.
970+
# Try to create a molecule from the SMILES string.
1011971
if isinstance(molecule, str):
1012-
is_smiles = True
1013972
try:
1014-
new_mol = _smiles(molecule)
973+
smiles = molecule
974+
molecule = _smiles(molecule)
975+
edit_mol = molecule._sire_object.edit()
976+
edit_mol = edit_mol.rename(f"smiles:{smiles}").molecule()
977+
molecule._sire_object = edit_mol.commit()
1015978
except Exception as e:
1016979
msg = "Unable to convert SMILES to Molecule using RDKit."
1017980
if _isVerbose():
1018981
msg += ": " + getattr(e, "message", repr(e))
1019982
raise _ThirdPartyError(msg) from e
1020983
else:
1021984
raise _ThirdPartyError(msg) from None
1022-
else:
1023-
is_smiles = False
1024-
new_mol = molecule.copy()
985+
986+
# Create a copy of the molecule.
987+
new_mol = molecule.copy()
1025988

1026989
# Use the net molecular charge passed as an option.
1027990
if self._net_charge is not None:
@@ -1251,45 +1214,23 @@ def run(self, molecule, work_dir=None, queue=None):
12511214
# Make the molecule 'mol' compatible with 'par_mol'. This will create
12521215
# a mapping between atom indices in the two molecules and add all of
12531216
# the new properties from 'par_mol' to 'mol'.
1254-
if is_smiles:
1255-
new_mol = par_mol
1256-
1257-
# We'll now add MolName and ResName info to the molecule, since
1258-
# this will be missing.
1259-
1260-
# Rename the molecule with the original SMILES string.
1261-
edit_mol = new_mol._sire_object.edit()
1262-
edit_mol = edit_mol.rename(molecule).molecule()
1263-
1264-
# Rename the residue LIG.
1265-
resname = _SireMol.ResName("LIG")
1266-
edit_mol = (
1267-
edit_mol.residue(_SireMol.ResIdx(0))
1268-
.rename(resname)
1269-
.molecule()
1217+
if self._ensure_compatible:
1218+
new_mol.makeCompatibleWith(
1219+
par_mol,
1220+
property_map=self._property_map,
1221+
overwrite=True,
1222+
verbose=False,
12701223
)
1271-
1272-
# Commit the changes.
1273-
new_mol._sire_object = edit_mol.commit()
1274-
12751224
else:
1276-
if self._ensure_compatible:
1225+
try:
12771226
new_mol.makeCompatibleWith(
12781227
par_mol,
12791228
property_map=self._property_map,
12801229
overwrite=True,
12811230
verbose=False,
12821231
)
1283-
else:
1284-
try:
1285-
new_mol.makeCompatibleWith(
1286-
par_mol,
1287-
property_map=self._property_map,
1288-
overwrite=True,
1289-
verbose=False,
1290-
)
1291-
except:
1292-
new_mol = par_mol
1232+
except:
1233+
new_mol = par_mol
12931234

12941235
# Record the forcefield used to parameterise the molecule.
12951236
new_mol._forcefield = ff

0 commit comments

Comments
 (0)