9265907197

Committed 28 May 2024 08:19AM UTC coverage: 98.317% (+2.0%) from 96.296%

Build # 9265907197

Build Type

Pull #42

travis-ci

Committed by

web-flow

Commit Message

Merge 0b9305451 into a1684afb9

Pull Request Pull Request #42: Use charges for determining valence

Run Details

244 of 253 branches covered (96.44%)

Branch coverage included in aggregate %.

55 of 59 new or added lines in 2 files covered. (93.22%)

1 existing line in 1 file now uncovered.

457 of 460 relevant lines covered (99.35%)

3.93 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

99.43

/pysmiles/smiles_helper.py

# -*- coding: utf-8 -*-
# Copyright 2018 Peter C Kroon

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Contains helper functions for parsing and writing SMILES strings, as well as
some convenience functions for adding hydrogens, and detecting aromaticity.
"""

import logging
import re
import operator

import networkx as nx
from . import PTE

LOGGER = logging.getLogger(__name__)

ISOTOPE_PATTERN = r'(?P<isotope>[\d]+)?'
ELEMENT_PATTERN = r'(?P<element>b|c|n|o|s|p|as|se|\*|[A-Z][a-z]{0,2})'
STEREO_PATTERN = r'(?P<stereo>@|@@|@TH[1-2]|@AL[1-2]|@SP[1-3]|@OH[\d]{1,2}|'\
                  r'@TB[\d]{1,2})?'
HCOUNT_PATTERN = r'(?P<hcount>H[\d]?)?'
CHARGE_PATTERN = r'(?P<charge>(-|\+)(\++|-+|[\d]{1,2})?)?'
CLASS_PATTERN = r'(?::(?P<class>[\d]+))?'
ATOM_PATTERN = re.compile(r'^\[' + ISOTOPE_PATTERN + ELEMENT_PATTERN +
                          STEREO_PATTERN + HCOUNT_PATTERN + CHARGE_PATTERN +
                          CLASS_PATTERN + r'\]$')

AROMATIC_ATOMS = "B C N O P S Se As *".split()

ORBITAL_SIZES = [[2],  # 1s
                 [2,  6],  # 2s, 2p
                 [2,  6],  # 3s, 3p
                 [2,  10, 6],  # 4s, 3d, 4p
                 [2,  10, 6],  # 5s, 4d, 5p
                 [2,  14, 10, 6],]  # 6s, 4f, 5d, 6p


def parse_atom(atom):
    """
    Parses a SMILES atom token, and returns a dict with the information.

    Note
    ----
    Can not deal with stereochemical information yet. This gets discarded.

    Parameters
    ----------
    atom : str
        The atom string to interpret. Looks something like one of the
        following: "C", "c", "[13CH3-1:2]"

    Returns
    -------
    dict
        A dictionary containing at least 'element', 'aromatic', and 'charge'. If
        present, will also contain 'hcount', 'isotope', and 'class'.
    """
    defaults = {'charge': 0, 'hcount': 0, 'aromatic': False}
    if not atom.startswith('[') and not atom.endswith(']'):
        if atom != '*':
            # Don't specify hcount to signal we don't actually know anything
            # about it
            return {'element': atom.capitalize(), 'charge': 0,
                    'aromatic': atom.islower()}
        else:
            return defaults.copy()
    match = ATOM_PATTERN.match(atom)
    if match is None:
        raise ValueError('The atom {} is malformatted'.format(atom))
    out = defaults.copy()
    out.update({k: v for k, v in match.groupdict().items() if v is not None})

    if out.get('element', 'X').islower():
        out['aromatic'] = True

    parse_helpers = {
        'isotope': int,
        'element': str.capitalize,
        'stereo': lambda x: x,
        'hcount': parse_hcount,
        'charge': parse_charge,
        'class': int,
        'aromatic': lambda x: x,
    }

    for attr, val_str in out.items():
        out[attr] = parse_helpers[attr](val_str)

    if out['element'] == '*':
        del out['element']

    if out.get('element') == 'H' and out.get('hcount', 0):
        raise ValueError("A hydrogen atom can't have hydrogens")

    if 'stereo' in out:
        LOGGER.warning('Atom "%s" contains stereochemical information that will be discarded.', atom)

    return out


def format_atom(molecule, node_key, default_element='*'):
    """
    Formats a node following SMILES conventions. Uses the attributes `element`,
    `charge`, `hcount`, `stereo`, `isotope` and `class`.

    Parameters
    ----------
    molecule : nx.Graph
        The molecule containing the atom.
    node_key : hashable
        The node key of the atom in `molecule`.
    default_element : str
        The element to use if the attribute is not present in the node.

    Returns
    -------
    str
        The atom as SMILES string.
    """
    node = molecule.nodes[node_key]
    name = node.get('element', default_element)
    charge = node.get('charge', 0)
    hcount = node.get('hcount', 0)
    stereo = node.get('stereo', None)
    isotope = node.get('isotope', '')
    class_ = node.get('class', '')
    aromatic = node.get('aromatic', False)
    default_h = has_default_h_count(molecule, node_key)

    if stereo is not None:  # pragma: nocover
        raise NotImplementedError

    if aromatic:
        name = name.lower()

    if (stereo is None and isotope == '' and charge == 0 and default_h and class_ == '' and
            (name.lower() in 'b c n o p s *'.split() or name in 'F Cl Br I'.split())):
        return name

    if hcount:
        hcountstr = 'H'
        if hcount > 1:
            hcountstr += str(hcount)
    else:
        hcountstr = ''

    if charge > 0:
        chargestr = '+'
        if charge > 1:
            chargestr += str(charge)
    elif charge < 0:
        chargestr = '-'
        if charge < -1:
            chargestr += str(-charge)
    else:
        chargestr = ''

    if class_ != '':
        class_ = ':{}'.format(class_)
    fmt = '[{isotope}{name}{stereo}{hcount}{charge}{class_}]'
    return fmt.format(isotope=isotope, name=name, stereo='', hcount=hcountstr,
                      charge=chargestr, class_=class_)


def parse_hcount(hcount_str):
    """
    Parses a SMILES hydrogen count specifications.

    Parameters
    ----------
    hcount_str : str
        The hydrogen count specification to parse.

    Returns
    -------
    int
        The number of hydrogens specified.
    """
    if not hcount_str:
        return 0
    if hcount_str == 'H':
        return 1
    return int(hcount_str[1:])


def parse_charge(charge_str):
    """
    Parses a SMILES charge specification.

    Parameters
    ----------
    charge_str : str
        The charge specification to parse.

    Returns
    -------
    int
        The charge.
    """
    if not charge_str:
        return 0
    signs = {'-': -1, '+': 1}
    sign = signs[charge_str[0]]
    if len(charge_str) > 1 and charge_str[1].isdigit():
        charge = sign * int(charge_str[1:])
    else:
        charge = sign * charge_str.count(charge_str[0])
    return charge


def add_explicit_hydrogens(mol):
    """
    Adds explicit hydrogen nodes to `mol`, the amount is determined by the node
    attribute 'hcount'. Will remove the 'hcount' attribute.

    Parameters
    ----------
    mol : nx.Graph
        The molecule to which explicit hydrogens should be added. Is modified
        in-place.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    h_atom = parse_atom('[H]')
    if 'hcount' in h_atom:  # pragma: nocover; defensive
        del h_atom['hcount']
    for n_idx in list(mol.nodes):
        hcount = mol.nodes[n_idx].get('hcount', 0)
        idxs = range(max(mol) + 1, max(mol) + hcount + 1)
        # Get the defaults from parse_atom.
        mol.add_nodes_from(idxs, **h_atom.copy())
        mol.add_edges_from([(n_idx, jdx) for jdx in idxs], order=1)
        if 'hcount' in mol.nodes[n_idx]:
            del mol.nodes[n_idx]['hcount']


def remove_explicit_hydrogens(mol):
    """
    Removes all explicit, simple hydrogens from `mol`. Simple means it is
    identical to the SMILES string "[H]", and has exactly one bond. Increments
    'hcount' where appropriate.

    Parameters
    ----------
    mol : nx.Graph
        The molecule whose explicit hydrogens should be removed. Is modified
        in-place.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    to_remove = set()
#    defaults = parse_atom('[H]')
    for n_idx in mol.nodes:
        node = mol.nodes[n_idx]
        neighbors = list(mol[n_idx])
        # TODO: get these defaults from parsing [H]. But do something smart
        #       with the hcount attribute.
        if (node.get('charge', 0) == 0 and node.get('element', '') == 'H' and
                'isotope' not in node and node.get('class', 0) == 0 and
                len(neighbors) == 1):
            neighbor = neighbors[0]
            if (mol.nodes[neighbor].get('element', '') == 'H' or
                    mol.edges[n_idx, neighbor].get('order', 1) != 1):
                # The molecule is H2, or the bond order is not 1.
                continue
            to_remove.add(n_idx)
            mol.nodes[neighbor]['hcount'] = mol.nodes[neighbor].get('hcount', 0) + 1
    mol.remove_nodes_from(to_remove)
    for n_idx in mol.nodes:
        if 'hcount' not in mol.nodes[n_idx]:
            mol.nodes[n_idx]['hcount'] = 0


def fill_valence(mol, respect_hcount=True, respect_bond_order=True,
                 max_bond_order=3):
    """
    Sets the attribute 'hcount' on all nodes in `mol` that don't have it yet.
    The value to which it is set is based on the node's 'element', and the
    number of bonds it has. Default valences are as specified by the global
    variable VALENCES.

    Parameters
    ----------
    mol : nx.Graph
        The molecule whose nodes should get a 'hcount'. Is modified in-place.
    respect_hcount : bool
        If True, don't change the hcount on nodes that already have it set.
    respect_bond_order : bool
        If False, first try to fill the valence by increasing bond orders, and
        add hydrogens after.
    max_bond_order : number
        Only meaningful if respect_bond_order is False. This is the highest
        bond order that will be set.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    if not respect_bond_order:
        increment_bond_orders(mol, max_bond_order=max_bond_order)
    for n_idx in mol:
        node = mol.nodes[n_idx]
        if ('hcount' in node and respect_hcount) or node.get('element') == 'H':
            continue
        missing = max(bonds_missing(mol, n_idx), 0)
        node['hcount'] = node.get('hcount', 0) + missing


def bonds_missing(mol, node_idx, use_order=True):
    """
    Returns how much the specified node is under valence. If use_order is
    False, treat all bonds as if they are order 1.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    node_idx : hashable
        The node to look at. Should be in mol.
    use_order : bool
        If False, treat all bonds as single.

    Returns
    -------
    int
        The number of missing bonds.
    """
    bonds = _bonds(mol, node_idx, use_order)
    bonds += mol.nodes[node_idx].get('hcount', 0)

    val = valence(mol.nodes[node_idx])
    val = [v for v in val if v >= bonds] or [0]
    return int(max(val[0] - bonds, 0))


def valence(atom):
    """
    Returns the valence of the atom. Since some elements can have
    multiple valences. The valence is returned as list.

    Parameters
    ----------
    atom: dict

    Returns
    -------
    list[int]
        The valences for the given atom.
    """
    element = atom.get("element", '*')
    if element == '*':
        electrons = 0
    else:
        electrons = PTE[element.capitalize()]['AtomicNumber']
    electrons -= atom.get('charge', 0)

    if electrons < 0:
        raise ValueError(f"Atom {atom} has a negative number of electrons: {electrons}")

    # Let's start by filling complete shells:
    for idx, shell in enumerate(ORBITAL_SIZES):
        shell_size = sum(shell)
        if shell_size <= electrons:
            electrons -= shell_size
        else:
            break
    else:  # nobreak
        raise ValueError(f'Too many electrons for sanity for {atom}: {electrons+sum(map(sum, ORBITAL_SIZES))}')

    # Any electrons we have leftover we distribute over the orbitals. First 1
    # electron in each, then we start making pairs. The resulting valence will
    # be the number of unpaired electrons. Added bonus/complication: electrons
    # in pairs we can excite to higher shells, increasing the number of unpaired
    # electrons
    shell = ORBITAL_SIZES[idx]
    shell_size = sum(shell)
    single_electrons, paired_electrons = divmod(electrons, shell_size//2)
    if single_electrons == 0:
        single_electrons, paired_electrons, empty_orbitals = paired_electrons, 0, shell_size//2 - paired_electrons
    elif single_electrons == 1:
        single_electrons = shell_size//2 - paired_electrons
        single_electrons, paired_electrons, empty_orbitals = single_electrons, paired_electrons, 0
    else:  # pragma: nocover
        raise AssertionError(f"{single_electrons=}")

    if idx >= 2:
        # Excite any paired electrons to a higher orbital to deal with
        # multi-valency. Naturally, each electron pair consists of 2 electrons,
        # and results in 2 new bonding electrons
        val = [single_electrons + 2*n for n in range(paired_electrons+1)]
    else:
        val = [single_electrons]

    return val


def _bonds(mol, node_idx, use_order=True):
    """
    Returns how many explicit bonds the specified node has. If use_order is
    False, treat all bonds as if they are order 1.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    node_idx : hashable
        The node to look at. Should be in mol.
    use_order : bool
        If False, treat all bonds as single.

    Returns
    -------
    int
        The number of bonds.
    """
    if use_order:
        bond_orders = map(operator.itemgetter(2),
                          mol.edges(nbunch=node_idx, data='order', default=1))
        bonds = sum(bond_orders)
    else:
        bonds = len(mol[node_idx])
    return bonds


def has_default_h_count(mol, node_idx, use_order=True):
    """
    Returns whether the hydrogen count for this atom is non-standard.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    node_idx : hashable
        The node to look at. Should be in mol.
    use_order : bool
        If False, treat all bonds as single.

    Returns
    -------
    bool
    """
    bonds = _bonds(mol, node_idx, use_order)
    val = valence(mol.nodes[node_idx])
    val = [v for v in val if v >= bonds] or [0]
    hcount = mol.nodes[node_idx].get('hcount', 0)
    return max(val[0] - bonds, 0) == hcount


def _prune_nodes(nodes, mol):
    new_nodes = []
    for node in nodes:
        # all wild card nodes are eligible
        if mol.nodes[node].get('element', '*') == '*':
            new_nodes.append(node)
            continue
        missing = bonds_missing(mol, node, use_order=True)
        if missing > 0:
            new_nodes.append(node)
    return mol.subgraph(new_nodes)


def mark_aromatic_atoms(mol, atoms=None):
    """
    Properly kekeulizes molecules and sets the aromatic attribute.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.
    atoms: collections.abc.Iterable
        The atoms to act on; all other nodes are pruned

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    if atoms is None:
        atoms = set(mol.nodes)
    # prune all nodes from molecule that are eligible and have
    # full valency
    ds_graph = _prune_nodes(atoms, mol)

    # set the aromatic attribute to False for all nodes
    # as a precaution
    nx.set_node_attributes(mol, False, 'aromatic')

    for sub_ds in nx.connected_components(ds_graph):
        sub_ds_graph = mol.subgraph(sub_ds)
        max_match = nx.max_weight_matching(sub_ds_graph)
        # we check if a maximum matching exists and
        # if it is perfect. if it is not perfect,
        # this graph originates from a completely invalid
        # smiles and we raise an error
        if not nx.is_perfect_matching(sub_ds_graph, max_match):
            msg = "Your molecule is invalid and cannot be kekulized."
            raise SyntaxError(msg)

        # we consider a node aromatic if it can take part in DIME
        # to do so all nodes in a delocalized subgraph have to be
        # part of a cycle system
        cycles = nx.cycle_basis(sub_ds_graph)
        nodes_in_cycles = []
        for cycle in cycles:
            nodes_in_cycles += cycle
            nx.set_node_attributes(mol, {node: True for node in cycle}, 'aromatic')

        for node in sub_ds_graph.nodes:
            if node not in nodes_in_cycles:
                mol.nodes[node]['aromatic'] = False

        for edge in max_match:
            if not edge[0] in nodes_in_cycles and not edge[1] in nodes_in_cycles:
                mol.edges[edge]['order'] = 2

def mark_aromatic_edges(mol):
    """
    Set all bonds between aromatic atoms (attribute 'aromatic' is `True`) to
    1.5. Gives all other bonds that don't have an order yet an order of 1.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    for edge in mol.edges:
        if all(mol.nodes[node].get('aromatic', 'False') for node in edge):
            mol.edges[edge]['order'] = 1.5
        elif 'order' not in mol.edges[edge]:
            mol.edges[edge]['order'] = 1

def correct_aromatic_rings(mol):
    """
    Sets hcount for all atoms, marks aromaticity for all atoms, and the order of
    all aromatic bonds to 1.5.

    Parameters
    ----------
    mol : nx.Graph
        The molecule.

    Returns
    -------
    None
        `mol` is modified in-place.
    """
    fill_valence(mol)
    mark_aromatic_atoms(mol, atoms=mol.nodes)
    mark_aromatic_edges(mol)


def increment_bond_orders(molecule, max_bond_order=3):
    """
    Increments bond orders up to what the atom's valence allows.

    Parameters
    ----------
    molecule : nx.Graph
        The molecule to process.
    max_bond_order : number
        The highest bond order allowed to make.

    Returns
    -------
    None
        molecule is modified in-place.
    """
    # Gather the number of open spots for all atoms beforehand, since some
    # might have multiple oxidation states (e.g. S). We don't want to change
    # oxidation state halfway through for some funny reason. It shouldn't be
    # nescessary, but it can't hurt.
    missing_bonds = {}
    for idx in molecule:
        missing_bonds[idx] = max(bonds_missing(molecule, idx), 0)

    for idx, jdx in molecule.edges:
        missing_idx = missing_bonds[idx]
        missing_jdx = missing_bonds[jdx]
        edge_missing = min(missing_idx, missing_jdx)
        current_order = molecule.edges[idx, jdx].get("order", 1)
        if current_order == 1.5:
            continue
        new_order = edge_missing + current_order
        new_order = min(new_order, max_bond_order)
        molecule.edges[idx, jdx]['order'] = new_order
        missing_bonds[idx] -= edge_missing
        missing_bonds[jdx] -= edge_missing

1	# -- coding: utf-8 --
2	# Copyright 2018 Peter C Kroon
3
4	# Licensed under the Apache License, Version 2.0 (the "License");
5	# you may not use this file except in compliance with the License.
6	# You may obtain a copy of the License at
7	#
8	# http://www.apache.org/licenses/LICENSE-2.0
9	#
10	# Unless required by applicable law or agreed to in writing, software
11	# distributed under the License is distributed on an "AS IS" BASIS,
12	# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13	# See the License for the specific language governing permissions and
14	# limitations under the License.
15
16	"""
17	Contains helper functions for parsing and writing SMILES strings, as well as
18	some convenience functions for adding hydrogens, and detecting aromaticity.
19	"""
20
21	import logging	4✔
22	import re	4✔
23	import operator	4✔
24
25	import networkx as nx	4✔
26	from . import PTE	4✔
27
28	LOGGER = logging.getLogger(__name__)	4✔
29
30	ISOTOPE_PATTERN = r'(?P<isotope>[\d]+)?'	4✔
31	ELEMENT_PATTERN = r'(?P<element>b\|c\|n\|o\|s\|p\|as\|se\|\*\|[A-Z][a-z]{0,2})'	4✔
32	STEREO_PATTERN = r'(?P<stereo>@\|@@\|@TH[1-2]\|@AL[1-2]\|@SP[1-3]\|@OH[\d]{1,2}\|'\	4✔
33	r'@TB[\d]{1,2})?'
34	HCOUNT_PATTERN = r'(?P<hcount>H[\d]?)?'	4✔
35	CHARGE_PATTERN = r'(?P<charge>(-\|\+)(\++\|-+\|[\d]{1,2})?)?'	4✔
36	CLASS_PATTERN = r'(?::(?P<class>[\d]+))?'	4✔
37	ATOM_PATTERN = re.compile(r'^\[' + ISOTOPE_PATTERN + ELEMENT_PATTERN +	4✔
38	STEREO_PATTERN + HCOUNT_PATTERN + CHARGE_PATTERN +
39	CLASS_PATTERN + r'\]$')
40
41	AROMATIC_ATOMS = "B C N O P S Se As *".split()	4✔
42
43	ORBITAL_SIZES = [[2], # 1s	4✔
44	[2, 6], # 2s, 2p
45	[2, 6], # 3s, 3p
46	[2, 10, 6], # 4s, 3d, 4p
47	[2, 10, 6], # 5s, 4d, 5p
48	[2, 14, 10, 6],] # 6s, 4f, 5d, 6p
49
50
51	def parse_atom(atom):	4✔
52	"""
53	Parses a SMILES atom token, and returns a dict with the information.
54
55	Note
56	----
57	Can not deal with stereochemical information yet. This gets discarded.
58
59	Parameters
60	----------
61	atom : str
62	The atom string to interpret. Looks something like one of the
63	following: "C", "c", "[13CH3-1:2]"
64
65	Returns
66	-------
67	dict
68	A dictionary containing at least 'element', 'aromatic', and 'charge'. If
69	present, will also contain 'hcount', 'isotope', and 'class'.
70	"""
71	defaults = {'charge': 0, 'hcount': 0, 'aromatic': False}	4✔
72	if not atom.startswith('[') and not atom.endswith(']'):	4✔
73	if atom != '*':	4✔
74	# Don't specify hcount to signal we don't actually know anything
75	# about it
76	return {'element': atom.capitalize(), 'charge': 0,	4✔
77	'aromatic': atom.islower()}
78	else:
79	return defaults.copy()	4✔
80	match = ATOM_PATTERN.match(atom)	4✔
81	if match is None:	4✔
82	raise ValueError('The atom {} is malformatted'.format(atom))	4✔
83	out = defaults.copy()	4✔
84	out.update({k: v for k, v in match.groupdict().items() if v is not None})	4✔
85
86	if out.get('element', 'X').islower():	4✔
87	out['aromatic'] = True	4✔
88
89	parse_helpers = {	4✔
90	'isotope': int,
91	'element': str.capitalize,
92	'stereo': lambda x: x,
93	'hcount': parse_hcount,
94	'charge': parse_charge,
95	'class': int,
96	'aromatic': lambda x: x,
97	}
98
99	for attr, val_str in out.items():	4✔
100	out[attr] = parse_helpers[attr](val_str)	4✔
101
102	if out['element'] == '*':	4✔
103	del out['element']	4✔
104
105	if out.get('element') == 'H' and out.get('hcount', 0):	4✔
106	raise ValueError("A hydrogen atom can't have hydrogens")	4✔
107
108	if 'stereo' in out:	4✔
109	LOGGER.warning('Atom "%s" contains stereochemical information that will be discarded.', atom)	4✔
110
111	return out	4✔
112
113
114	def format_atom(molecule, node_key, default_element='*'):	4✔
115	"""
116	Formats a node following SMILES conventions. Uses the attributes `element`,
117	`charge`, `hcount`, `stereo`, `isotope` and `class`.
118
119	Parameters
120	----------
121	molecule : nx.Graph
122	The molecule containing the atom.
123	node_key : hashable
124	The node key of the atom in `molecule`.
125	default_element : str
126	The element to use if the attribute is not present in the node.
127
128	Returns
129	-------
130	str
131	The atom as SMILES string.
132	"""
133	node = molecule.nodes[node_key]	4✔
134	name = node.get('element', default_element)	4✔
135	charge = node.get('charge', 0)	4✔
136	hcount = node.get('hcount', 0)	4✔
137	stereo = node.get('stereo', None)	4✔
138	isotope = node.get('isotope', '')	4✔
139	class_ = node.get('class', '')	4✔
140	aromatic = node.get('aromatic', False)	4✔
141	default_h = has_default_h_count(molecule, node_key)	4✔
142
143	if stereo is not None: # pragma: nocover
144	raise NotImplementedError
145
146	if aromatic:	4✔
147	name = name.lower()	4✔
148
149	if (stereo is None and isotope == '' and charge == 0 and default_h and class_ == '' and	4✔
150	(name.lower() in 'b c n o p s *'.split() or name in 'F Cl Br I'.split())):
151	return name	4✔
152
153	if hcount:	4✔
154	hcountstr = 'H'	4✔
155	if hcount > 1:	4✔
156	hcountstr += str(hcount)	4✔
157	else:
158	hcountstr = ''	4✔
159
160	if charge > 0:	4✔
161	chargestr = '+'	4✔
162	if charge > 1:	4✔
163	chargestr += str(charge)	4✔
164	elif charge < 0:	4✔
165	chargestr = '-'	4✔
166	if charge < -1:	4✔
167	chargestr += str(-charge)	4✔
168	else:
169	chargestr = ''	4✔
170
171	if class_ != '':	4✔
172	class_ = ':{}'.format(class_)	4✔
173	fmt = '[{isotope}{name}{stereo}{hcount}{charge}{class_}]'	4✔
174	return fmt.format(isotope=isotope, name=name, stereo='', hcount=hcountstr,	4✔
175	charge=chargestr, class_=class_)
176
177
178	def parse_hcount(hcount_str):	4✔
179	"""
180	Parses a SMILES hydrogen count specifications.
181
182	Parameters
183	----------
184	hcount_str : str
185	The hydrogen count specification to parse.
186
187	Returns
188	-------
189	int
190	The number of hydrogens specified.
191	"""
192	if not hcount_str:	4✔
193	return 0	4✔
194	if hcount_str == 'H':	4✔
195	return 1	4✔
196	return int(hcount_str[1:])	4✔
197
198
199	def parse_charge(charge_str):	4✔
200	"""
201	Parses a SMILES charge specification.
202
203	Parameters
204	----------
205	charge_str : str
206	The charge specification to parse.
207
208	Returns
209	-------
210	int
211	The charge.
212	"""
213	if not charge_str:	4✔
214	return 0	4✔
215	signs = {'-': -1, '+': 1}	4✔
216	sign = signs[charge_str[0]]	4✔
217	if len(charge_str) > 1 and charge_str[1].isdigit():	4✔
218	charge = sign * int(charge_str[1:])	4✔
219	else:
220	charge = sign * charge_str.count(charge_str[0])	4✔
221	return charge	4✔
222
223
224	def add_explicit_hydrogens(mol):	4✔
225	"""
226	Adds explicit hydrogen nodes to `mol`, the amount is determined by the node
227	attribute 'hcount'. Will remove the 'hcount' attribute.
228
229	Parameters
230	----------
231	mol : nx.Graph
232	The molecule to which explicit hydrogens should be added. Is modified
233	in-place.
234
235	Returns
236	-------
237	None
238	`mol` is modified in-place.
239	"""
240	h_atom = parse_atom('[H]')	4✔
241	if 'hcount' in h_atom: # pragma: nocover; defensive
242	del h_atom['hcount']
243	for n_idx in list(mol.nodes):	4✔
244	hcount = mol.nodes[n_idx].get('hcount', 0)	4✔
245	idxs = range(max(mol) + 1, max(mol) + hcount + 1)	4✔
246	# Get the defaults from parse_atom.
247	mol.add_nodes_from(idxs, **h_atom.copy())	4✔
248	mol.add_edges_from([(n_idx, jdx) for jdx in idxs], order=1)	4✔
249	if 'hcount' in mol.nodes[n_idx]:	4✔
250	del mol.nodes[n_idx]['hcount']	4✔
251
252
253	def remove_explicit_hydrogens(mol):	4✔
254	"""
255	Removes all explicit, simple hydrogens from `mol`. Simple means it is
256	identical to the SMILES string "[H]", and has exactly one bond. Increments
257	'hcount' where appropriate.
258
259	Parameters
260	----------
261	mol : nx.Graph
262	The molecule whose explicit hydrogens should be removed. Is modified
263	in-place.
264
265	Returns
266	-------
267	None
268	`mol` is modified in-place.
269	"""
270	to_remove = set()	4✔
271	# defaults = parse_atom('[H]')
272	for n_idx in mol.nodes:	4✔
273	node = mol.nodes[n_idx]	4✔
274	neighbors = list(mol[n_idx])	4✔
275	# TODO: get these defaults from parsing [H]. But do something smart
276	# with the hcount attribute.
277	if (node.get('charge', 0) == 0 and node.get('element', '') == 'H' and	4✔
278	'isotope' not in node and node.get('class', 0) == 0 and
279	len(neighbors) == 1):
280	neighbor = neighbors[0]	4✔
281	if (mol.nodes[neighbor].get('element', '') == 'H' or	4✔
282	mol.edges[n_idx, neighbor].get('order', 1) != 1):
283	# The molecule is H2, or the bond order is not 1.
284	continue	4✔
285	to_remove.add(n_idx)	4✔
286	mol.nodes[neighbor]['hcount'] = mol.nodes[neighbor].get('hcount', 0) + 1	4✔
287	mol.remove_nodes_from(to_remove)	4✔
288	for n_idx in mol.nodes:	4✔
289	if 'hcount' not in mol.nodes[n_idx]:	4✔
290	mol.nodes[n_idx]['hcount'] = 0	4✔
291
292
293	def fill_valence(mol, respect_hcount=True, respect_bond_order=True,	4✔
294	max_bond_order=3):
295	"""
296	Sets the attribute 'hcount' on all nodes in `mol` that don't have it yet.
297	The value to which it is set is based on the node's 'element', and the
298	number of bonds it has. Default valences are as specified by the global
299	variable VALENCES.
300
301	Parameters
302	----------
303	mol : nx.Graph
304	The molecule whose nodes should get a 'hcount'. Is modified in-place.
305	respect_hcount : bool
306	If True, don't change the hcount on nodes that already have it set.
307	respect_bond_order : bool
308	If False, first try to fill the valence by increasing bond orders, and
309	add hydrogens after.
310	max_bond_order : number
311	Only meaningful if respect_bond_order is False. This is the highest
312	bond order that will be set.
313
314	Returns
315	-------
316	None
317	`mol` is modified in-place.
318	"""
319	if not respect_bond_order:	4✔
320	increment_bond_orders(mol, max_bond_order=max_bond_order)	4✔
321	for n_idx in mol:	4✔
322	node = mol.nodes[n_idx]	4✔
323	if ('hcount' in node and respect_hcount) or node.get('element') == 'H':	4✔
324	continue	4✔
325	missing = max(bonds_missing(mol, n_idx), 0)	4✔
326	node['hcount'] = node.get('hcount', 0) + missing	4✔
327
328
329	def bonds_missing(mol, node_idx, use_order=True):	4✔
330	"""
331	Returns how much the specified node is under valence. If use_order is
332	False, treat all bonds as if they are order 1.
333
334	Parameters
335	----------
336	mol : nx.Graph
337	The molecule.
338	node_idx : hashable
339	The node to look at. Should be in mol.
340	use_order : bool
341	If False, treat all bonds as single.
342
343	Returns
344	-------
345	int
346	The number of missing bonds.
347	"""
348	bonds = _bonds(mol, node_idx, use_order)	4✔
349	bonds += mol.nodes[node_idx].get('hcount', 0)	4✔
350
351	val = valence(mol.nodes[node_idx])	4✔
352	val = [v for v in val if v >= bonds] or [0]	4✔
353	return int(max(val[0] - bonds, 0))	4✔
354
355
356	def valence(atom):	4✔
357	"""
358	Returns the valence of the atom. Since some elements can have
359	multiple valences. The valence is returned as list.
360
361	Parameters
362	----------
363	atom: dict
364
365	Returns
366	-------
367	list[int]
368	The valences for the given atom.
369	"""
370	element = atom.get("element", '*')	4✔
371	if element == '*':	4✔
372	electrons = 0	4✔
373	else:
374	electrons = PTE[element.capitalize()]['AtomicNumber']	4✔
375	electrons -= atom.get('charge', 0)	4✔
376
377	if electrons < 0:	4✔
378	raise ValueError(f"Atom {atom} has a negative number of electrons: {electrons}")	4✔
379
380	# Let's start by filling complete shells:
381	for idx, shell in enumerate(ORBITAL_SIZES):	4✔
382	shell_size = sum(shell)	4✔
383	if shell_size <= electrons:	4✔
384	electrons -= shell_size	4✔
385	else:
386	break	4✔
387	else: # nobreak
388	raise ValueError(f'Too many electrons for sanity for {atom}: {electrons+sum(map(sum, ORBITAL_SIZES))}')	4✔
389
390	# Any electrons we have leftover we distribute over the orbitals. First 1
391	# electron in each, then we start making pairs. The resulting valence will
392	# be the number of unpaired electrons. Added bonus/complication: electrons
393	# in pairs we can excite to higher shells, increasing the number of unpaired
394	# electrons
395	shell = ORBITAL_SIZES[idx]	4✔
396	shell_size = sum(shell)	4✔
397	single_electrons, paired_electrons = divmod(electrons, shell_size//2)	4✔
398	if single_electrons == 0:	4✔
399	single_electrons, paired_electrons, empty_orbitals = paired_electrons, 0, shell_size//2 - paired_electrons	4✔
400	elif single_electrons == 1:	4✔
401	single_electrons = shell_size//2 - paired_electrons	4✔
402	single_electrons, paired_electrons, empty_orbitals = single_electrons, paired_electrons, 0	4✔
403	else: # pragma: nocover
404	raise AssertionError(f"{single_electrons=}")
405
406	if idx >= 2:	4✔
407	# Excite any paired electrons to a higher orbital to deal with
408	# multi-valency. Naturally, each electron pair consists of 2 electrons,
409	# and results in 2 new bonding electrons
410	val = [single_electrons + 2*n for n in range(paired_electrons+1)]	4✔
411	else:
412	val = [single_electrons]	4✔
413
414	return val	4✔
415
416
417	def _bonds(mol, node_idx, use_order=True):	4✔
418	"""
419	Returns how many explicit bonds the specified node has. If use_order is
420	False, treat all bonds as if they are order 1.
421
422	Parameters
423	----------
424	mol : nx.Graph
425	The molecule.
426	node_idx : hashable
427	The node to look at. Should be in mol.
428	use_order : bool
429	If False, treat all bonds as single.
430
431	Returns
432	-------
433	int
434	The number of bonds.
435	"""
436	if use_order:	4!
437	bond_orders = map(operator.itemgetter(2),	4✔
438	mol.edges(nbunch=node_idx, data='order', default=1))
439	bonds = sum(bond_orders)	4✔
440	else:
UNCOV 441	bonds = len(mol[node_idx])	×
442	return bonds	4✔
443
444
445	def has_default_h_count(mol, node_idx, use_order=True):	4✔
446	"""
447	Returns whether the hydrogen count for this atom is non-standard.
448
449	Parameters
450	----------
451	mol : nx.Graph
452	The molecule.
453	node_idx : hashable
454	The node to look at. Should be in mol.
455	use_order : bool
456	If False, treat all bonds as single.
457
458	Returns
459	-------
460	bool
461	"""
462	bonds = _bonds(mol, node_idx, use_order)	4✔
463	val = valence(mol.nodes[node_idx])	4✔
464	val = [v for v in val if v >= bonds] or [0]	4✔
465	hcount = mol.nodes[node_idx].get('hcount', 0)	4✔
466	return max(val[0] - bonds, 0) == hcount	4✔
467
468
469	def _prune_nodes(nodes, mol):	4✔
470	new_nodes = []	4✔
471	for node in nodes:	4✔
472	# all wild card nodes are eligible
473	if mol.nodes[node].get('element', '') == '':	4✔
474	new_nodes.append(node)	4✔
475	continue	4✔
476	missing = bonds_missing(mol, node, use_order=True)	4✔
477	if missing > 0:	4✔
478	new_nodes.append(node)	4✔
479	return mol.subgraph(new_nodes)	4✔
480
481
482	def mark_aromatic_atoms(mol, atoms=None):	4✔
483	"""
484	Properly kekeulizes molecules and sets the aromatic attribute.
485
486	Parameters
487	----------
488	mol : nx.Graph
489	The molecule.
490	atoms: collections.abc.Iterable
491	The atoms to act on; all other nodes are pruned
492
493	Returns
494	-------
495	None
496	`mol` is modified in-place.
497	"""
498	if atoms is None:	4✔
499	atoms = set(mol.nodes)	4✔
500	# prune all nodes from molecule that are eligible and have
501	# full valency
502	ds_graph = _prune_nodes(atoms, mol)	4✔
503
504	# set the aromatic attribute to False for all nodes
505	# as a precaution
506	nx.set_node_attributes(mol, False, 'aromatic')	4✔
507
508	for sub_ds in nx.connected_components(ds_graph):	4✔
509	sub_ds_graph = mol.subgraph(sub_ds)	4✔
510	max_match = nx.max_weight_matching(sub_ds_graph)	4✔
511	# we check if a maximum matching exists and
512	# if it is perfect. if it is not perfect,
513	# this graph originates from a completely invalid
514	# smiles and we raise an error
515	if not nx.is_perfect_matching(sub_ds_graph, max_match):	4✔
516	msg = "Your molecule is invalid and cannot be kekulized."	4✔
517	raise SyntaxError(msg)	4✔
518
519	# we consider a node aromatic if it can take part in DIME
520	# to do so all nodes in a delocalized subgraph have to be
521	# part of a cycle system
522	cycles = nx.cycle_basis(sub_ds_graph)	4✔
523	nodes_in_cycles = []	4✔
524	for cycle in cycles:	4✔
525	nodes_in_cycles += cycle	4✔
526	nx.set_node_attributes(mol, {node: True for node in cycle}, 'aromatic')	4✔
527
528	for node in sub_ds_graph.nodes:	4✔
529	if node not in nodes_in_cycles:	4✔
530	mol.nodes[node]['aromatic'] = False	4✔
531
532	for edge in max_match:	4✔
533	if not edge[0] in nodes_in_cycles and not edge[1] in nodes_in_cycles:	4✔
534	mol.edges[edge]['order'] = 2	4✔
535
536	def mark_aromatic_edges(mol):	4✔
537	"""
538	Set all bonds between aromatic atoms (attribute 'aromatic' is `True`) to
539	1.5. Gives all other bonds that don't have an order yet an order of 1.
540
541	Parameters
542	----------
543	mol : nx.Graph
544	The molecule.
545
546	Returns
547	-------
548	None
549	`mol` is modified in-place.
550	"""
551	for edge in mol.edges:	4✔
552	if all(mol.nodes[node].get('aromatic', 'False') for node in edge):	4✔
553	mol.edges[edge]['order'] = 1.5	4✔
554	elif 'order' not in mol.edges[edge]:	4✔
555	mol.edges[edge]['order'] = 1	4✔
556
557	def correct_aromatic_rings(mol):	4✔
558	"""
559	Sets hcount for all atoms, marks aromaticity for all atoms, and the order of
560	all aromatic bonds to 1.5.
561
562	Parameters
563	----------
564	mol : nx.Graph
565	The molecule.
566
567	Returns
568	-------
569	None
570	`mol` is modified in-place.
571	"""
572	fill_valence(mol)	4✔
573	mark_aromatic_atoms(mol, atoms=mol.nodes)	4✔
574	mark_aromatic_edges(mol)	4✔
575
576
577	def increment_bond_orders(molecule, max_bond_order=3):	4✔
578	"""
579	Increments bond orders up to what the atom's valence allows.
580
581	Parameters
582	----------
583	molecule : nx.Graph
584	The molecule to process.
585	max_bond_order : number
586	The highest bond order allowed to make.
587
588	Returns
589	-------
590	None
591	molecule is modified in-place.
592	"""
593	# Gather the number of open spots for all atoms beforehand, since some
594	# might have multiple oxidation states (e.g. S). We don't want to change
595	# oxidation state halfway through for some funny reason. It shouldn't be
596	# nescessary, but it can't hurt.
597	missing_bonds = {}	4✔
598	for idx in molecule:	4✔
599	missing_bonds[idx] = max(bonds_missing(molecule, idx), 0)	4✔
600
601	for idx, jdx in molecule.edges:	4✔
602	missing_idx = missing_bonds[idx]	4✔
603	missing_jdx = missing_bonds[jdx]	4✔
604	edge_missing = min(missing_idx, missing_jdx)	4✔
605	current_order = molecule.edges[idx, jdx].get("order", 1)	4✔
606	if current_order == 1.5:	4✔
607	continue	4✔
608	new_order = edge_missing + current_order	4✔
609	new_order = min(new_order, max_bond_order)	4✔
610	molecule.edges[idx, jdx]['order'] = new_order	4✔
611	missing_bonds[idx] -= edge_missing	4✔
612	missing_bonds[jdx] -= edge_missing	4✔

pckroon / pysmiles / 9265907197

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous