From 8c344663e45a48a7d0aca6733de2837a709c13be Mon Sep 17 00:00:00 2001 From: Rocco Meli Date: Fri, 6 Dec 2024 17:42:30 +0100 Subject: [PATCH 1/2] cleanup --- testsuite/MDAnalysisTests/converters/test_rdkit.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/testsuite/MDAnalysisTests/converters/test_rdkit.py b/testsuite/MDAnalysisTests/converters/test_rdkit.py index 16793a4484..c3440c72f8 100644 --- a/testsuite/MDAnalysisTests/converters/test_rdkit.py +++ b/testsuite/MDAnalysisTests/converters/test_rdkit.py @@ -44,6 +44,7 @@ _standardize_patterns) from rdkit import Chem from rdkit.Chem import AllChem + from rdkit.Chem import rdDistGeom except ImportError: pass @@ -682,6 +683,7 @@ def test_ions(self, smi): def test_reorder_atoms(self, smi): mol = Chem.MolFromSmiles(smi) mol = Chem.AddHs(mol) + rdDistGeom.EmbedMolecule(mol) # remove bond order and charges info pdb = Chem.MolToPDBBlock(mol) u = mda.Universe(StringIO(pdb), format="PDB") From 53fffa63a125c94c5d082a773993ce2460900af0 Mon Sep 17 00:00:00 2001 From: Rocco Meli Date: Fri, 6 Dec 2024 22:06:56 +0100 Subject: [PATCH 2/2] add safeguard --- package/CHANGELOG | 3 ++- package/MDAnalysis/converters/RDKit.py | 9 ++++++--- .../MDAnalysisTests/converters/test_rdkit.py | 16 +++++++++++++--- 3 files changed, 21 insertions(+), 7 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 5bce69eaec..9a40eb53cf 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,13 +14,14 @@ The rules for this file: ------------------------------------------------------------------------------- -??/??/?? IAlibay +??/??/?? IAlibay, RMeli * 2.9.0 Fixes Enhancements +* Add check and warning for empty (all zero) coordinates in RDKit converter (PR #4824) Changes diff --git a/package/MDAnalysis/converters/RDKit.py b/package/MDAnalysis/converters/RDKit.py index 85f55b7900..aa78a7ea0c 100644 --- a/package/MDAnalysis/converters/RDKit.py +++ b/package/MDAnalysis/converters/RDKit.py @@ -363,9 +363,12 @@ def convert(self, obj, cache=True, NoImplicit=True, max_iter=200, # add a conformer for the current Timestep if hasattr(ag, "positions"): - if np.isnan(ag.positions).any(): - warnings.warn("NaN detected in coordinates, the output " - "molecule will not have 3D coordinates assigned") + if np.isnan(ag.positions).any() or np.allclose( + ag.positions, 0.0, rtol=0.0, atol=1e-12 + ): + warnings.warn("NaN or empty coordinates detected in coordinates, " + "the output molecule will not have 3D coordinates " + "assigned") else: # assign coordinates conf = Chem.Conformer(mol.GetNumAtoms()) diff --git a/testsuite/MDAnalysisTests/converters/test_rdkit.py b/testsuite/MDAnalysisTests/converters/test_rdkit.py index c3440c72f8..1d56c4c5f6 100644 --- a/testsuite/MDAnalysisTests/converters/test_rdkit.py +++ b/testsuite/MDAnalysisTests/converters/test_rdkit.py @@ -44,7 +44,6 @@ _standardize_patterns) from rdkit import Chem from rdkit.Chem import AllChem - from rdkit.Chem import rdDistGeom except ImportError: pass @@ -332,7 +331,7 @@ def test_nan_coords(self): xyz = u.atoms.positions xyz[0][2] = np.nan u.atoms.positions = xyz - with pytest.warns(UserWarning, match="NaN detected"): + with pytest.warns(UserWarning, match="NaN .* detected"): mol = u.atoms.convert_to("RDKIT") with pytest.raises(ValueError, match="Bad Conformer Id"): mol.GetConformer() @@ -683,7 +682,6 @@ def test_ions(self, smi): def test_reorder_atoms(self, smi): mol = Chem.MolFromSmiles(smi) mol = Chem.AddHs(mol) - rdDistGeom.EmbedMolecule(mol) # remove bond order and charges info pdb = Chem.MolToPDBBlock(mol) u = mda.Universe(StringIO(pdb), format="PDB") @@ -694,6 +692,18 @@ def test_reorder_atoms(self, smi): expected = [a.GetSymbol() for a in mol.GetAtoms()] assert values == expected + @pytest.mark.parametrize("smi", [ + "O=S(C)(C)=NC", + ]) + def test_warn_empty_coords(self, smi): + mol = Chem.MolFromSmiles(smi) + mol = Chem.AddHs(mol) + # remove bond order and charges info + pdb = Chem.MolToPDBBlock(mol) + u = mda.Universe(StringIO(pdb), format="PDB") + with pytest.warns(match="NaN or empty coordinates detected"): + u.atoms.convert_to.rdkit() + def test_pdb_names(self): u = mda.Universe(PDB_helix) mol = u.atoms.convert_to.rdkit()