9109869248

Committed 16 May 2024 09:26AM UTC coverage: 97.525% (+1.2%) from 96.296%

Build # 9109869248

Build Type

Pull #42

travis-ci

Committed by

web-flow

Commit Message

Merge e36eebe71 into a1684afb9

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

Run Details

235 of 247 branches covered (95.14%)

Branch coverage included in aggregate %.

30 of 31 new or added lines in 1 file covered. (96.77%)

2 existing lines in 1 file now uncovered.

435 of 440 relevant lines covered (98.86%)

3.95 hits per line

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

97.99

/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

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

ELECTRON_COUNTS = {"B": 5, "C": 6, "N": 7, "O": 8, "F": 9, "P": 15,
                   "S": 16, "Cl": 17, "As": 33, "Se": 34, "Br": 35, "I": 53}


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:
        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:
        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.
    """
    electrons = ELECTRON_COUNTS.get(atom.get('element', '*').capitalize(), 0) - atom.get('charge', 0)
    # 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)
    to_assign = min(electrons, shell_size // 2)
    single_electrons = to_assign
    electrons -= to_assign

    to_assign = min(electrons, shell_size // 2)
    electron_pairs = to_assign
    single_electrons -= to_assign
    electrons -= to_assign
    assert electrons == 0

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

pckroon / pysmiles / 9109869248

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