Skip to content

Commit 974cdd1

Browse files
lohedgeslohedges
authored
Merge pull request #393 from OpenBioSim/fix_392
Fix issue #392
2 parents de6629c + 5fec8d7 commit 974cdd1

5 files changed

Lines changed: 165 additions & 31 deletions

File tree

‎python/BioSimSpace/Parameters/_utils.py‎

Lines changed: 63 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -28,14 +28,14 @@
2828

2929
import tempfile as _tempfile
3030

31-
from .. import _is_notebook
31+
from .. import _isVerbose
3232
from .. import IO as _IO
3333
from .. import _Utils
3434
from ..Units.Charge import electron_charge as _electron_charge
3535
from .._SireWrappers import Molecule as _Molecule
3636

3737

38-
def formalCharge(molecule):
38+
def formalCharge(molecule, property_map={}):
3939
"""
4040
Compute the formal charge on a molecule. This function requires that
4141
the molecule has explicit hydrogen atoms.
@@ -46,6 +46,11 @@ def formalCharge(molecule):
4646
molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
4747
A molecule object.
4848
49+
property_map : dict
50+
A dictionary that maps system "properties" to their user defined
51+
values. This allows the user to refer to properties with their
52+
own naming scheme, e.g. { "charge" : "my-charge" }
53+
4954
Returns
5055
-------
5156
@@ -58,6 +63,9 @@ def formalCharge(molecule):
5863
"'molecule' must be of type 'BioSimSpace._SireWrappers.Molecule'"
5964
)
6065

66+
if not isinstance(property_map, dict):
67+
raise TypeError("'property_map' must be of type 'dict'")
68+
6169
from rdkit import Chem as _Chem
6270
from rdkit import RDLogger as _RDLogger
6371

@@ -71,15 +79,56 @@ def formalCharge(molecule):
7179
# Zero the total formal charge.
7280
formal_charge = 0
7381

74-
# Run in the working directory.
75-
with _Utils.cd(work_dir):
76-
# Save the molecule to a PDB file.
77-
_IO.saveMolecules("tmp", molecule, "PDB")
78-
79-
# Read the ligand PDB into an RDKit molecule.
80-
mol = _Chem.MolFromPDBFile("tmp.pdb")
81-
82-
# Compute the formal charge.
83-
formal_charge = _Chem.rdmolops.GetFormalCharge(mol)
84-
85-
return formal_charge * _electron_charge
82+
# Get the fileformat property name.
83+
property = property_map.get("fileformat", "fileformat")
84+
85+
# Preferentially use the file format that the molecule was loaded from.
86+
try:
87+
# Get the raw list of formats.
88+
raw_formats = molecule._sire_object.property(property).value().split(",")
89+
90+
# Remove all formats other than PDB and SDF.
91+
formats = [f for f in raw_formats if f in ["PDB", "SDF"]]
92+
93+
if len(formats) == 0:
94+
formats = ["PDB", "SDF"]
95+
elif len(formats) == 1:
96+
if formats[0] == "PDB":
97+
formats.append("SDF")
98+
else:
99+
formats.append("PDB")
100+
except:
101+
formats = ["PDB", "SDF"]
102+
103+
# List of exceptions.
104+
exceptions = []
105+
106+
# Try going via each format in turn.
107+
for format in formats:
108+
try:
109+
with _Utils.cd(work_dir):
110+
# Save the molecule in the given format.
111+
_IO.saveMolecules("tmp", molecule, format)
112+
113+
# Load with RDKit.
114+
if format == "SDF":
115+
rdmol = _Chem.MolFromMolFile("tmp.sdf")
116+
else:
117+
rdmol = _Chem.MolFromPDBFile("tmp.pdb")
118+
119+
# Compute the formal charge.
120+
formal_charge = _Chem.rdmolops.GetFormalCharge(rdmol)
121+
122+
return formal_charge * _electron_charge
123+
124+
except Exception as e:
125+
exceptions.append(e)
126+
127+
# If we got this far, then we failed to compute the formal charge.
128+
msg = "Failed to compute the formal charge on the molecule."
129+
if _isVerbose():
130+
for e in exceptions:
131+
msg += "\n\n" + str(e)
132+
raise RuntimeError(msg)
133+
else:
134+
raise RuntimeError(msg) from None

‎python/BioSimSpace/Sandpit/Exscientia/Parameters/_utils.py‎

Lines changed: 63 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -28,14 +28,14 @@
2828

2929
import tempfile as _tempfile
3030

31-
from .. import _is_notebook
31+
from .. import _isVerbose
3232
from .. import IO as _IO
3333
from .. import _Utils
3434
from ..Units.Charge import electron_charge as _electron_charge
3535
from .._SireWrappers import Molecule as _Molecule
3636

3737

38-
def formalCharge(molecule):
38+
def formalCharge(molecule, property_map={}):
3939
"""
4040
Compute the formal charge on a molecule. This function requires that
4141
the molecule has explicit hydrogen atoms.
@@ -46,6 +46,11 @@ def formalCharge(molecule):
4646
molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
4747
A molecule object.
4848
49+
property_map : dict
50+
A dictionary that maps system "properties" to their user defined
51+
values. This allows the user to refer to properties with their
52+
own naming scheme, e.g. { "charge" : "my-charge" }
53+
4954
Returns
5055
-------
5156
@@ -58,6 +63,9 @@ def formalCharge(molecule):
5863
"'molecule' must be of type 'BioSimSpace._SireWrappers.Molecule'"
5964
)
6065

66+
if not isinstance(property_map, dict):
67+
raise TypeError("'property_map' must be of type 'dict'")
68+
6169
from rdkit import Chem as _Chem
6270
from rdkit import RDLogger as _RDLogger
6371

@@ -71,15 +79,56 @@ def formalCharge(molecule):
7179
# Zero the total formal charge.
7280
formal_charge = 0
7381

74-
# Run in the working directory.
75-
with _Utils.cd(work_dir):
76-
# Save the molecule to a PDB file.
77-
_IO.saveMolecules("tmp", molecule, "PDB")
78-
79-
# Read the ligand PDB into an RDKit molecule.
80-
mol = _Chem.MolFromPDBFile("tmp.pdb")
81-
82-
# Compute the formal charge.
83-
formal_charge = _Chem.rdmolops.GetFormalCharge(mol)
84-
85-
return formal_charge * _electron_charge
82+
# Get the fileformat property name.
83+
property = property_map.get("fileformat", "fileformat")
84+
85+
# Preferentially use the file format that the molecule was loaded from.
86+
try:
87+
# Get the raw list of formats.
88+
raw_formats = molecule._sire_object.property(property).value().split(",")
89+
90+
# Remove all formats other than PDB and SDF.
91+
formats = [f for f in raw_formats if f in ["PDB", "SDF"]]
92+
93+
if len(formats) == 0:
94+
formats = ["PDB", "SDF"]
95+
elif len(formats) == 1:
96+
if formats[0] == "PDB":
97+
formats.append("SDF")
98+
else:
99+
formats.append("PDB")
100+
except:
101+
formats = ["PDB", "SDF"]
102+
103+
# List of exceptions.
104+
exceptions = []
105+
106+
# Try going via each format in turn.
107+
for format in formats:
108+
try:
109+
with _Utils.cd(work_dir):
110+
# Save the molecule in the given format.
111+
_IO.saveMolecules("tmp", molecule, format)
112+
113+
# Load with RDKit.
114+
if format == "SDF":
115+
rdmol = _Chem.MolFromMolFile("tmp.sdf")
116+
else:
117+
rdmol = _Chem.MolFromPDBFile("tmp.pdb")
118+
119+
# Compute the formal charge.
120+
formal_charge = _Chem.rdmolops.GetFormalCharge(rdmol)
121+
122+
return formal_charge * _electron_charge
123+
124+
except Exception as e:
125+
exceptions.append(e)
126+
127+
# If we got this far, then we failed to compute the formal charge.
128+
msg = "Failed to compute the formal charge on the molecule."
129+
if _isVerbose():
130+
for e in exceptions:
131+
msg += "\n\n" + str(e)
132+
raise RuntimeError(msg)
133+
else:
134+
raise RuntimeError(msg) from None

‎tests/Parameters/test_parameters.py‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -188,3 +188,20 @@ def test_acdoctor():
188188

189189
# Make sure parameterisation works when acdoctor is disabled.
190190
mol = BSS.Parameters.gaff(mol, acdoctor=False).getMolecule()
191+
192+
193+
def test_broken_sdf_formal_charge():
194+
"""
195+
Test that the PDB fallback works when using a broken SDF file
196+
as input for calculating the total formal charge.
197+
"""
198+
199+
# Load the molecule.
200+
mol = BSS.IO.readMolecules(f"{url}/broken.sdf")[0]
201+
202+
# Compute the formal charge.
203+
charge = BSS.Parameters._utils.formalCharge(mol)
204+
205+
from math import isclose
206+
207+
assert isclose(charge.value(), 0.0, abs_tol=1e-6)

‎tests/Sandpit/Exscientia/Parameters/test_parameters.py‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -193,3 +193,20 @@ def test_acdoctor():
193193

194194
# Make sure parameterisation works when acdoctor is disabled.
195195
mol = BSS.Parameters.gaff(mol, acdoctor=False).getMolecule()
196+
197+
198+
def test_broken_sdf_formal_charge():
199+
"""
200+
Test that the PDB fallback works when using a broken SDF file
201+
as input for calculating the total formal charge.
202+
"""
203+
204+
# Load the molecule.
205+
mol = BSS.IO.readMolecules(f"{url}/broken.sdf")[0]
206+
207+
# Compute the formal charge.
208+
charge = BSS.Parameters._utils.formalCharge(mol)
209+
210+
from math import isclose
211+
212+
assert isclose(charge.value(), 0.0, abs_tol=1e-6)

‎tests/Sandpit/Exscientia/Protocol/test_config.py‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -480,9 +480,11 @@ def system_and_boresch_restraint():
480480
def test_turn_on_restraint_boresch(self, system_and_boresch_restraint):
481481
"""Test for turning on multiple distance restraints"""
482482
system, restraint = system_and_boresch_restraint
483-
protocol = FreeEnergy(perturbation_type="restraint",
484-
runtime=1*BSS.Units.Time.nanosecond,
485-
timestep=2*BSS.Units.Time.femtosecond)
483+
protocol = FreeEnergy(
484+
perturbation_type="restraint",
485+
runtime=1 * BSS.Units.Time.nanosecond,
486+
timestep=2 * BSS.Units.Time.femtosecond,
487+
)
486488
freenrg = BSS.FreeEnergy.AlchemicalFreeEnergy(
487489
system, protocol, engine="SOMD", restraint=restraint
488490
)

0 commit comments

Comments
 (0)