Skip to content

Commit

Permalink
Write directly to SDF format where possible to preserve stereochemistry.
Browse files Browse the repository at this point in the history 
[closes #389]
  • Loading branch information
@lohedges @annamherz
lohedges authored and annamherz committed Oct 6, 2022
1 parent 2e4e6ca commit 15b7c01
Show file tree
Hide file tree
Showing 2 changed files with 108 additions and 78 deletions.
93 changes: 54 additions & 39 deletions python/BioSimSpace/Parameters/Protocol/_openforcefield.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,40 +197,59 @@ def run(self, molecule, work_dir=None, queue=None):
else:
raise IOError(msg) from None
else:
# Write the molecule to a PDB file.
try:
_IO.saveMolecules(prefix + "molecule", molecule, "pdb")
except Exception as e:
msg = "Failed to write the molecule to 'PDB' format."
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise IOError(msg) from e
else:
raise IOError(msg) from None

# Create an RDKit molecule from the PDB file.
try:
rdmol = _Chem.MolFromPDBFile(prefix + "molecule.pdb", removeHs=False)
except Exception as e:
msg = "RDKit was unable to read the molecular PDB file!"
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None

# Use RDKit to write back to SDF format.
try:
writer = _Chem.SDWriter(prefix + "molecule.sdf")
writer.write(rdmol)
writer.close()
except Exception as e:
msg = "RDKit was unable to write the molecular SDF file!"
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None
# If the molecule was originally loaded from an SDF format file,
# then write back to the same format.
fileformat_prop = self._property_map.get("fileformat", "fileformat")
if (molecule._sire_object.hasProperty(fileformat_prop) and
"SDF" in molecule._sire_object.property("fileformat").value()):
# Write the molecule to SDF format.
try:
_IO.saveMolecules(prefix + "molecule", molecule, "sdf")
except Exception as e:
msg = "Failed to write the molecule to 'SDF' format."
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise IOError(msg) from e
else:
raise IOError(msg) from None

# Otherwise, go via an intermediate PDB file and use RDKit to try
# to recover stereochemistry.
else:
# Write the molecule to a PDB file.
try:
_IO.saveMolecules(prefix + "molecule", molecule, "pdb")
except Exception as e:
msg = "Failed to write the molecule to 'PDB' format."
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise IOError(msg) from e
else:
raise IOError(msg) from None

# Create an RDKit molecule from the PDB file.
try:
rdmol = _Chem.MolFromPDBFile(prefix + "molecule.pdb", removeHs=False)
except Exception as e:
msg = "RDKit was unable to read the molecular PDB file!"
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None

# Use RDKit to write back to SDF format.
try:
writer = _Chem.SDWriter(prefix + "molecule.sdf")
writer.write(rdmol)
writer.close()
except Exception as e:
msg = "RDKit was unable to write the molecular SDF file!"
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None

# Create the Open Forcefield Molecule from the intermediate SDF file,
# as recommended by @j-wags and @mattwthompson.
Expand Down Expand Up @@ -291,11 +310,7 @@ def run(self, molecule, work_dir=None, queue=None):

# Convert the OpenMM System to a ParmEd structure.
try:
if is_smiles:
parmed_structure = _parmed.openmm.load_topology(omm_topology, omm_system, off_molecule.conformers[0])
else:
pdbfile = _PDBFile(prefix + "molecule.pdb")
parmed_structure = _parmed.openmm.load_topology(omm_topology, omm_system, pdbfile.positions)
parmed_structure = _parmed.openmm.load_topology(omm_topology, omm_system, off_molecule.conformers[0])
except Exception as e:
msg = "Unable to convert OpenMM System to ParmEd structure!"
if _isVerbose():
Expand Down
93 changes: 54 additions & 39 deletions python/BioSimSpace/Sandpit/Exscientia/Parameters/Protocol/_openforcefield.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,40 +197,59 @@ def run(self, molecule, work_dir=None, queue=None):
else:
raise IOError(msg) from None
else:
# Write the molecule to a PDB file.
try:
_IO.saveMolecules(prefix + "molecule", molecule, "pdb")
except Exception as e:
msg = "Failed to write the molecule to 'PDB' format."
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise IOError(msg) from e
else:
raise IOError(msg) from None

# Create an RDKit molecule from the PDB file.
try:
rdmol = _Chem.MolFromPDBFile(prefix + "molecule.pdb", removeHs=False)
except Exception as e:
msg = "RDKit was unable to read the molecular PDB file!"
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None

# Use RDKit to write back to SDF format.
try:
writer = _Chem.SDWriter(prefix + "molecule.sdf")
writer.write(rdmol)
writer.close()
except Exception as e:
msg = "RDKit was unable to write the molecular SDF file!"
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None
# If the molecule was originally loaded from an SDF format file,
# then write back to the same format.
fileformat_prop = self._property_map.get("fileformat", "fileformat")
if (molecule._sire_object.hasProperty(fileformat_prop) and
"SDF" in molecule._sire_object.property("fileformat").value()):
# Write the molecule to SDF format.
try:
_IO.saveMolecules(prefix + "molecule", molecule, "sdf")
except Exception as e:
msg = "Failed to write the molecule to 'SDF' format."
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise IOError(msg) from e
else:
raise IOError(msg) from None

# Otherwise, go via an intermediate PDB file and use RDKit to try
# to recover stereochemistry.
else:
# Write the molecule to a PDB file.
try:
_IO.saveMolecules(prefix + "molecule", molecule, "pdb")
except Exception as e:
msg = "Failed to write the molecule to 'PDB' format."
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise IOError(msg) from e
else:
raise IOError(msg) from None

# Create an RDKit molecule from the PDB file.
try:
rdmol = _Chem.MolFromPDBFile(prefix + "molecule.pdb", removeHs=False)
except Exception as e:
msg = "RDKit was unable to read the molecular PDB file!"
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None

# Use RDKit to write back to SDF format.
try:
writer = _Chem.SDWriter(prefix + "molecule.sdf")
writer.write(rdmol)
writer.close()
except Exception as e:
msg = "RDKit was unable to write the molecular SDF file!"
if _isVerbose():
msg += ": " + getattr(e, "message", repr(e))
raise _ThirdPartyError(msg) from e
else:
raise _ThirdPartyError(msg) from None

# Create the Open Forcefield Molecule from the intermediate SDF file,
# as recommended by @j-wags and @mattwthompson.
Expand Down Expand Up @@ -291,11 +310,7 @@ def run(self, molecule, work_dir=None, queue=None):

# Convert the OpenMM System to a ParmEd structure.
try:
if is_smiles:
parmed_structure = _parmed.openmm.load_topology(omm_topology, omm_system, off_molecule.conformers[0])
else:
pdbfile = _PDBFile(prefix + "molecule.pdb")
parmed_structure = _parmed.openmm.load_topology(omm_topology, omm_system, pdbfile.positions)
parmed_structure = _parmed.openmm.load_topology(omm_topology, omm_system, off_molecule.conformers[0])
except Exception as e:
msg = "Unable to convert OpenMM System to ParmEd structure!"
if _isVerbose():
Expand Down

0 comments on commit 15b7c01

Please sign in to comment.