Skip to content

Commit

Permalink
Adds aromaticity and Gasteiger charges guessers (#2926)
Browse files Browse the repository at this point in the history
Towards #2468 

## Work done in this PR
* Add aromaticity and Gasteiger charges guessers which work via the RDKIT converter.
  • Loading branch information
cbouy authored Apr 23, 2021
1 parent b82635a commit 32e3254
Show file tree
Hide file tree
Showing 2 changed files with 96 additions and 0 deletions.
47 changes: 47 additions & 0 deletions package/MDAnalysis/topology/guessers.py
Original file line number Diff line number Diff line change
Expand Up @@ -475,3 +475,50 @@ def guess_atom_charge(atomname):
"""
# TODO: do something slightly smarter, at least use name/element
return 0.0


def guess_aromaticities(atomgroup):
"""Guess aromaticity of atoms using RDKit
Parameters
----------
atomgroup : mda.core.groups.AtomGroup
Atoms for which the aromaticity will be guessed
Returns
-------
aromaticities : numpy.ndarray
Array of boolean values for the aromaticity of each atom
.. versionadded:: 2.0.0
"""
mol = atomgroup.convert_to("RDKIT")
atoms = sorted(mol.GetAtoms(),
key=lambda a: a.GetIntProp("_MDAnalysis_index"))
return np.array([atom.GetIsAromatic() for atom in atoms])


def guess_gasteiger_charges(atomgroup):
"""Guess Gasteiger partial charges using RDKit
Parameters
----------
atomgroup : mda.core.groups.AtomGroup
Atoms for which the charges will be guessed
Returns
-------
charges : numpy.ndarray
Array of float values representing the charge of each atom
.. versionadded:: 2.0.0
"""
mol = atomgroup.convert_to("RDKIT")
from rdkit.Chem.rdPartialCharges import ComputeGasteigerCharges
ComputeGasteigerCharges(mol, throwOnParamFailure=True)
atoms = sorted(mol.GetAtoms(),
key=lambda a: a.GetIntProp("_MDAnalysis_index"))
return np.array([atom.GetDoubleProp("_GasteigerCharge") for atom in atoms],
dtype=np.float32)
49 changes: 49 additions & 0 deletions testsuite/MDAnalysisTests/topology/test_guessers.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,19 @@
from MDAnalysisTests.core.test_fragments import make_starshape
import MDAnalysis.tests.datafiles as datafiles

from MDAnalysisTests.util import import_not_available


try:
from rdkit import Chem
from rdkit.Chem.rdPartialCharges import ComputeGasteigerCharges
except ImportError:
pass

requires_rdkit = pytest.mark.skipif(import_not_available("rdkit"),
reason="requires RDKit")


class TestGuessMasses(object):
def test_guess_masses(self):
out = guessers.guess_masses(['C', 'C', 'H'])
Expand Down Expand Up @@ -139,3 +152,39 @@ def test_guess_bonds_peptide():
bonds = guessers.guess_bonds(u.atoms, u.atoms.positions)
assert_equal(np.sort(u.bonds.indices, axis=0),
np.sort(bonds, axis=0))


@pytest.mark.parametrize("smi", [
"c1ccccc1",
"C1=CC=CC=C1",
"CCO",
"c1ccccc1Cc1ccccc1",
"CN1C=NC2=C1C(=O)N(C(=O)N2C)C",
])
@requires_rdkit
def test_guess_aromaticities(smi):
mol = Chem.MolFromSmiles(smi)
mol = Chem.AddHs(mol)
expected = np.array([atom.GetIsAromatic() for atom in mol.GetAtoms()])
u = mda.Universe(mol)
values = guessers.guess_aromaticities(u.atoms)
assert_equal(values, expected)


@pytest.mark.parametrize("smi", [
"c1ccccc1",
"C1=CC=CC=C1",
"CCO",
"c1ccccc1Cc1ccccc1",
"CN1C=NC2=C1C(=O)N(C(=O)N2C)C",
])
@requires_rdkit
def test_guess_gasteiger_charges(smi):
mol = Chem.MolFromSmiles(smi)
mol = Chem.AddHs(mol)
ComputeGasteigerCharges(mol, throwOnParamFailure=True)
expected = np.array([atom.GetDoubleProp("_GasteigerCharge")
for atom in mol.GetAtoms()], dtype=np.float32)
u = mda.Universe(mol)
values = guessers.guess_gasteiger_charges(u.atoms)
assert_equal(values, expected)

0 comments on commit 32e3254

Please sign in to comment.