# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- https://www.mdanalysis.org # Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors # (see the file AUTHORS for the full list of names) # # Released under the GNU Lesser General Public License, v2 or higher # # Please cite your use of MDAnalysis in published work: # # R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, # D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. # MDAnalysis: A Python package for the rapid analysis of molecular dynamics # simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th # Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. # doi: 10.25080/majora-629e541a-00e # # N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # """RDKit bond order inferring --- :mod:`MDAnalysis.converters.RDKitInferring` ============================================================================= Bond order and formal charge inferring for RDKit molecules. Because most MD file formats don't preserve bond order information directly (or formal charges to some extent), the classes provided here give users different options to either provide or infer this information to the RDKit molecule constructed from the topology. Having bond orders and formal charges properly defined is a requirement for almost all cheminformatics-related task, hence the different approaches proposed here to cover most use cases. You can also defined your own function if need be, see the :mod:`~MDAnalysis.converters.RDKit` module for an example. These classes are meant to be passed directly to the RDKit converter:: >>> import MDAnalysis as mda >>> from rdkit import Chem >>> u = mda.Universe("aspirin.pdb") >>> template = Chem.MolFromSmiles("CC(=O)Oc1ccccc1C(=O)O") >>> inferrer = mda.converters.RDKitInferring.TemplateInferrer(template) >>> rdkit_mol = u.atoms.convert_to.rdkit(inferrer=inferrer) Classes ------- .. autoclass:: MDAnalysisInferrer :members: :private-members: .. autoclass:: TemplateInferrer :members: .. autoclass:: RDKitInferrer :members: .. autofunction:: sanitize_mol .. autofunction:: reorder_atoms """ import warnings from collections import defaultdict from contextlib import suppress from dataclasses import dataclass from typing import ClassVar, Dict, List, Optional, Tuple import numpy as np RDBONDORDER = {} with suppress(ImportError): from rdkit import Chem from rdkit.Chem.AllChem import AssignBondOrdersFromTemplate from rdkit.Chem.rdChemReactions import ChemicalReaction, ReactionFromSmarts RDBONDORDER = { 1: Chem.BondType.SINGLE, 1.5: Chem.BondType.AROMATIC, "ar": Chem.BondType.AROMATIC, 2: Chem.BondType.DOUBLE, 3: Chem.BondType.TRIPLE, } # add string version of the key for each bond RDBONDORDER.update({str(key): value for key, value in RDBONDORDER.items()}) PERIODIC_TABLE = Chem.GetPeriodicTable() with suppress(ImportError): from rdkit.Chem.rdDetermineBonds import ( DetermineBondOrders, ) # available since 2022.09.1 def reorder_atoms( mol: "Chem.Mol", field: str = "_MDAnalysis_index" ) -> "Chem.Mol": """Reorder atoms based on the given field. Defaults to sorting in the same order as the input AtomGroup. Notes ----- If the field is not found on the molecule, skips reordering. """ with suppress(KeyError): order = np.argsort([atom.GetIntProp(field) for atom in mol.GetAtoms()]) return Chem.RenumberAtoms(mol, order.astype(int).tolist()) # not a molecule converted by MDAnalysis warnings.warn( f"{field!r} not available on the input mol atoms, skipping reordering of atoms." ) return mol def sanitize_mol(mol: "Chem.Mol") -> None: """Sanitizes the molecule.""" Chem.SanitizeMol(mol) @dataclass(frozen=True) class MDAnalysisInferrer: """Bond order and formal charge inferring as originally implemented for the RDKit converter. This algorithm only relies on the topology with explicit hydrogens to assign bond orders and formal charges. Attributes ---------- max_iter : int Maximum number of iterations to standardize conjugated systems. See :meth:`~MDAnalysisInferrer._rebuild_conjugated_bonds` sanitize : bool Whether to sanitize the molecule or not. MONATOMIC_CATION_CHARGES : ClassVar[Dict[int, int]] Charges that should be assigned to monatomic cations. Maps atomic number to their formal charge. Anion charges are directly handled by the code using the typical valence of the atom. STANDARDIZATION_REACTIONS : ClassVar[List[str]] Reactions uses by :meth:`~MDAnalysisInferrer._standardize_patterns` to fix challenging cases must have single reactant and product, and cannot add any atom. Notes ----- There are some molecules containing specific substructures that this inferrer cannot currently tackle correctly. See `Issue #3339 `__ for more info. .. versionadded:: 2.10.0 """ MONATOMIC_CATION_CHARGES: ClassVar[Dict[int, int]] = { 3: 1, 11: 1, 19: 1, 37: 1, 47: 1, 55: 1, 12: 2, 20: 2, 29: 2, 30: 2, 38: 2, 56: 2, 26: 2, # Fe could also be +3 13: 3, } STANDARDIZATION_REACTIONS: ClassVar[List[str]] = [ "[C-;X2;H0:1]=[O:2]>>[C+0:1]=[O:2]", # Cterm "[N-;X2;H1;$(N-[*^3]):1]>>[N+0:1]", # Nterm "[#6-:1]-[#6:2]=[O:3]>>[#6+0:1]=[#6:2]-[O-:3]", # keto-enolate "[C-;v3:1]-[#7+0;v3;H2:2]>>[#6+0:1]=[#7+:2]", # ARG "[#6+0;H0:1]=[#6+0:2]-[#7;X3:3]-[#6-;X3:4]" ">>[#6:1]=[#6:2]-[#7+:3]=[#6+0:4]", # HIP "[S;D4;!v6:1]-[*-:2]>>[S;v6:1]=[*+0:2]", # sulfone "[#7+0;X3:1]-[*-:2]>>[#7+:1]=[*+0:2]", # charged-N ] max_iter: int = 200 sanitize: bool = True def __call__(self, mol: "Chem.Mol") -> "Chem.Mol": """Infer bond orders and formal charges in the molecule. Parameters ---------- mol : Chem.Mol The molecule to infer bond orders and charges for. Returns ------- A new molecule with proper bond orders and charges. """ self._infer_bo_and_charges(mol) mol = self._standardize_patterns(mol) mol = reorder_atoms(mol) if self.sanitize: sanitize_mol(mol) return mol @classmethod def _atom_sorter(cls, atom: "Chem.Atom") -> Tuple[int, int]: """Sorts atoms in the molecule in a way that makes it easy for the bond order and charge infering code to get the correct state on the first try. Currently sorts by number of unpaired electrons, then by number of heavy atom neighbors (i.e. atoms at the edge first).""" num_heavy_neighbors = len( [ neighbor for neighbor in atom.GetNeighbors() if neighbor.GetAtomicNum() > 1 ] ) return (-cls._get_nb_unpaired_electrons(atom)[0], num_heavy_neighbors) @classmethod def _infer_bo_and_charges(cls, mol: "Chem.Mol") -> None: """Infer bond orders and formal charges from a molecule. Since most MD topology files don't explicitly retain information on bond orders or charges, it has to be guessed from the topology. This is done by looping over each atom and comparing its expected valence to the current valence to get the Number of Unpaired Electrons (NUE). If an atom has a negative NUE, it needs a positive formal charge (-NUE). If two neighbouring atoms have UEs, the bond between them most likely has to be increased by the value of the smallest NUE. If after this process, an atom still has UEs, it needs a negative formal charge of -NUE. Parameters ---------- mol : rdkit.Chem.rdchem.RWMol The molecule is modified inplace and must have all hydrogens added Notes ----- This algorithm is order dependant. For example, for a carboxylate group ``R-C(-O)-O`` the first oxygen read will receive a double bond and the other one will be charged. It will also affect more complex conjugated systems. """ # heavy atoms sorted by number of heavy atom neighbors (lower first) # then NUE (higher first) atoms = sorted( [atom for atom in mol.GetAtoms() if atom.GetAtomicNum() > 1], key=cls._atom_sorter, ) for atom in atoms: # monatomic ions if atom.GetDegree() == 0: atom.SetFormalCharge( cls.MONATOMIC_CATION_CHARGES.get( atom.GetAtomicNum(), -cls._get_nb_unpaired_electrons(atom)[0], ) ) mol.UpdatePropertyCache(strict=False) continue # get NUE for each possible valence nue = cls._get_nb_unpaired_electrons(atom) # if there's only one possible valence state and the corresponding # NUE is negative, it means we can only add a positive charge to # the atom if (len(nue) == 1) and (nue[0] < 0): atom.SetFormalCharge(-nue[0]) mol.UpdatePropertyCache(strict=False) # go to next atom if above case or atom has no unpaired electron if (len(nue) == 1) and (nue[0] <= 0): continue neighbors = sorted( atom.GetNeighbors(), reverse=True, key=lambda a: cls._get_nb_unpaired_electrons(a)[0], ) # check if one of the neighbors has a common NUE for na in neighbors: # get NUE for the neighbor na_nue = cls._get_nb_unpaired_electrons(na) # smallest common NUE common_nue = min( [i for i in [*nue, *na_nue] if i >= 0], default=0, ) # a common NUE of 0 means we don't need to do anything if common_nue != 0: # increase bond order bond = mol.GetBondBetweenAtoms(atom.GetIdx(), na.GetIdx()) order = common_nue + 1 bond.SetBondType(RDBONDORDER[order]) mol.UpdatePropertyCache(strict=False) # go to next atom if one of the valences is complete nue = cls._get_nb_unpaired_electrons(atom) if any(n == 0 for n in nue): break # if atom valence is still not filled nue = cls._get_nb_unpaired_electrons(atom) if not any(n == 0 for n in nue): # transform nue to charge atom.SetFormalCharge(-nue[0]) atom.SetNumRadicalElectrons(0) mol.UpdatePropertyCache(strict=False) @staticmethod def _get_nb_unpaired_electrons(atom: "Chem.Atom") -> List[int]: """Calculate the number of unpaired electrons (NUE) of an atom Parameters ---------- atom: rdkit.Chem.rdchem.Atom The atom for which the NUE will be computed Returns ------- nue : List[int] The NUE for each possible valence of the atom """ expected_vs = PERIODIC_TABLE.GetValenceList(atom.GetAtomicNum()) current_v = atom.GetTotalValence() - atom.GetFormalCharge() return [v - current_v for v in expected_vs] def _standardize_patterns( self, mol: "Chem.Mol", max_iter: Optional[int] = None ) -> "Chem.Mol": """Standardizes functional groups Uses :meth:`~MDAnalysisInferrer._rebuild_conjugated_bonds` to standardize conjugated systems, and SMARTS reactions for other functional groups. Due to the way reactions work, we first have to split the molecule by fragments. Then, for each fragment, we apply the standardization reactions. Finally, the fragments are recombined. Parameters ---------- mol : rdkit.Chem.rdchem.RWMol The molecule to standardize max_iter : Optional[int] Maximum number of iterations to standardize conjugated systems. .. deprecated:: 2.10.0 Will be removed in 3.0, use ``inferrer=MDAnalysisInferrer(max_iter=...)`` instead. Returns ------- mol : rdkit.Chem.rdchem.Mol The standardized molecule Notes ----- The following functional groups are transformed in this order: +---------------+------------------------------------------------------------------------------+ | Name | Reaction | +===============+==============================================================================+ | conjugated | ``[*-:1]-[*:2]=[*:3]-[*-:4]>>[*+0:1]=[*:2]-[*:3]=[*+0:4]`` | +---------------+------------------------------------------------------------------------------+ | conjugated N+ | ``[N;X3;v3:1]-[*:2]=[*:3]-[*-:4]>>[N+:1]=[*:2]-[*:3]=[*+0:4]`` | +---------------+------------------------------------------------------------------------------+ | conjugated O- | ``[O:1]=[#6+0,#7+:2]-[*:3]=[*:4]-[*-:5]>>[O-:1]-[*:2]=[*:3]-[*:4]=[*+0:5]`` | +---------------+------------------------------------------------------------------------------+ | conjug. S=O | ``[O-:1]-[S;D4;v4:2]-[*:3]=[*:4]-[*-:5]>>[O+0:1]=[*:2]=[*:3]-[*:4]=[*+0:5]`` | +---------------+------------------------------------------------------------------------------+ | Cterm | ``[C-;X2;H0:1]=[O:2]>>[C+0:1]=[O:2]`` | +---------------+------------------------------------------------------------------------------+ | Nterm | ``[N-;X2;H1;$(N-[*^3]):1]>>[N+0:1]`` | +---------------+------------------------------------------------------------------------------+ | keto-enolate | ``[#6-:1]-[#6:2]=[O:3]>>[#6+0:1]=[#6:2]-[O-:3]`` | +---------------+------------------------------------------------------------------------------+ | arginine | ``[C-;v3:1]-[#7+0;v3;H2:2]>>[#6+0:1]=[#7+:2]`` | +---------------+------------------------------------------------------------------------------+ | histidine | ``[#6+0;H0:1]=[#6+0:2]-[#7;X3:3]-[#6-;X3:4]>>[#6:1]=[#6:2]-[#7+:3]=[#6+0:4]``| +---------------+------------------------------------------------------------------------------+ | sulfone | ``[S;D4;!v6:1]-[*-:2]>>[S;v6:1]=[*+0:2]`` | +---------------+------------------------------------------------------------------------------+ | charged N | ``[#7+0;X3:1]-[*-:2]>>[#7+:1]=[*+0:2]`` | +---------------+------------------------------------------------------------------------------+ """ if max_iter is None: max_iter = self.max_iter else: warnings.warn( "Specifying `max_iter` is deprecated and will be removed in a " "future update. Directly create an instance of " "`MDAnalysisInferrer` with `MDAnalysisInferrer(max_iter=...)` " "instead.", DeprecationWarning, ) # standardize conjugated systems self._rebuild_conjugated_bonds(mol, max_iter) # list of sanitized reactions reactions = [ ReactionFromSmarts(rxn) for rxn in self.STANDARDIZATION_REACTIONS ] # fragment mol (reactions must have single reactant and product) fragments = list( Chem.GetMolFrags(mol, asMols=True, sanitizeFrags=self.sanitize) ) for reactant in fragments: self._apply_reactions(reactions, reactant) # recombine fragments newmol = fragments.pop(0) for fragment in fragments: newmol = Chem.CombineMols(newmol, fragment) return newmol @staticmethod def _apply_reactions( reactions: List["ChemicalReaction"], reactant: "Chem.Mol" ) -> None: """Applies a series of unimolecular reactions to a molecule. The reactions must not add any atom to the product. The molecule is modified inplace. Parameters ---------- reactions : List[rdkit.Chem.rdChemReactions.ChemicalReaction] Reactions from SMARTS. Each reaction is applied until no more transformations can be made. reactant : rdkit.Chem.rdchem.Mol The molecule to transform inplace """ reactant.UpdatePropertyCache(strict=False) Chem.KekulizeIfPossible(reactant) for reaction in reactions: while reaction.RunReactantInPlace(reactant): reactant.UpdatePropertyCache(strict=False) reactant.UpdatePropertyCache(strict=False) Chem.KekulizeIfPossible(reactant) def _rebuild_conjugated_bonds( self, mol: "Chem.Mol", max_iter: Optional[int] = None ) -> None: """Rebuild conjugated bonds without negatively charged atoms at the beginning and end of the conjugated system Depending on the order in which atoms are read during the conversion, the :meth:`~MDAnalysisInferrer._infer_bo_and_charges` function might write conjugated systems with a double bond less and both edges of the system as anions instead of the usual alternating single and double bonds. This function corrects this behaviour by using an iterative procedure. The problematic molecules always follow the same pattern: ``anion[-*=*]n-anion`` instead of ``*=[*-*=]n*``, where ``n`` is the number of successive single and double bonds. The goal of the iterative procedure is to make ``n`` as small as possible by consecutively transforming ``anion-*=*`` into ``*=*-anion`` until it reaches the smallest pattern with ``n=1``. This last pattern is then transformed ``anion-*=*-anion`` to ``*=*-*=*``. Since ``anion-*=*`` is the same as ``*=*-anion`` in terms of SMARTS, we can control that we don't transform the same triplet of atoms back and forth by adding their indices to a list. This function also handles conjugated systems with triple bonds. The molecule needs to be kekulized first to also cover systems with aromatic rings. Parameters ---------- mol : rdkit.Chem.rdchem.RWMol The molecule to transform, modified inplace max_iter : Optional[int] Maximum number of iterations to standardize conjugated systems. .. deprecated:: 2.10.0 Will be removed in 3.0, use ``inferrer=MDAnalysisInferrer(max_iter=...)`` instead. Notes ----- The molecule is modified inplace """ max_iter = self.max_iter if max_iter is None else max_iter mol.UpdatePropertyCache(strict=False) Chem.KekulizeIfPossible(mol) # pattern used to find problematic conjugated bonds # there's usually an even number of matches for this pattern = Chem.MolFromSmarts("[*-{1-2}]-,=[*+0]=,#[*+0]") # pattern used to finish fixing a series of conjugated bonds base_end_pattern = Chem.MolFromSmarts( "[*-{1-2}]-,=[*+0]=,#[*+0]-,=[*-{1-2}]" ) # used when there's an odd number of matches for `pattern` odd_end_pattern = Chem.MolFromSmarts( "[*-]-[*+0]=[*+0]-[*-,$([#7;X3;v3]),$([#6+0,#7+1]=O),$([S;D4;v4]-[O-])]" ) # number of unique matches with the pattern n_matches = len( {match[0] for match in mol.GetSubstructMatches(pattern)} ) # nothing to standardize if n_matches == 0: return # single match (unusual) if n_matches == 1: # as a last resort, the only way to standardize is to find a # nitrogen that can accept a double bond and a positive charge # or a carbonyl that will become an enolate end_pattern = odd_end_pattern # at least 2 matches else: # give priority to base end pattern and then deal with the odd one # if necessary end_pattern = base_end_pattern backtrack = [] backtrack_cycles = 0 for _ in range(max_iter): # check for ending pattern end_match = mol.GetSubstructMatch(end_pattern) if end_match: # index of each atom anion1, a1, a2, anion2 = end_match term_atom = mol.GetAtomWithIdx(anion2) # edge-case 1: C-[O-] or [N+]-[O-] # [*-]-*=*-[C,N+]=O --> *=*-*=[C,N+]-[O-] # transform the =O to -[O-] if ( term_atom.GetAtomicNum() == 6 and term_atom.GetFormalCharge() == 0 ) or ( term_atom.GetAtomicNum() == 7 and term_atom.GetFormalCharge() == 1 ): for neighbor in term_atom.GetNeighbors(): bond = mol.GetBondBetweenAtoms( anion2, neighbor.GetIdx() ) if ( neighbor.GetAtomicNum() == 8 and bond.GetBondTypeAsDouble() == 2 ): bond.SetBondType(Chem.BondType.SINGLE) neighbor.SetFormalCharge(-1) break # edge-case 2: S=O # [*-]-*=*-[Sv4]-[O-] --> *=*-*=[Sv6]=O # transform -[O-] to =O elif ( term_atom.GetAtomicNum() == 16 and term_atom.GetFormalCharge() == 0 ): for neighbor in term_atom.GetNeighbors(): bond = mol.GetBondBetweenAtoms( anion2, neighbor.GetIdx() ) if ( neighbor.GetAtomicNum() == 8 and neighbor.GetFormalCharge() == -1 and bond.GetBondTypeAsDouble() == 1 ): bond.SetBondType(Chem.BondType.DOUBLE) neighbor.SetFormalCharge(0) break # not an edge case: # increment the charge else: term_atom.SetFormalCharge(term_atom.GetFormalCharge() + 1) # common to all cases: # [*-]-*=*-[R] --> *=*-*=[R] # increment the charge and switch single and double bonds a = mol.GetAtomWithIdx(anion1) a.SetFormalCharge(a.GetFormalCharge() + 1) b = mol.GetBondBetweenAtoms(anion1, a1) b.SetBondType(RDBONDORDER[b.GetBondTypeAsDouble() + 1]) b = mol.GetBondBetweenAtoms(a1, a2) b.SetBondType(RDBONDORDER[b.GetBondTypeAsDouble() - 1]) b = mol.GetBondBetweenAtoms(a2, anion2) b.SetBondType(RDBONDORDER[b.GetBondTypeAsDouble() + 1]) mol.UpdatePropertyCache(strict=False) continue # switch the position of the charge and the double bond matches = mol.GetSubstructMatches(pattern) if matches: # check if we haven't already transformed this triplet for match in matches: # store order-independent atom indices g = set(match) # already transformed --> try the next one if g in backtrack: continue # add to backtracking and start the switch anion, a1, a2 = match backtrack.append(g) break # already performed all changes else: if backtrack_cycles == 1: # might be stuck because there were 2 "odd" end # patterns misqualified as a single base one end_pattern = odd_end_pattern elif backtrack_cycles > 1: # likely stuck changing the same bonds over and over so # let's stop here mol.UpdatePropertyCache(strict=False) return match = matches[0] anion, a1, a2 = match backtrack = [set(match)] backtrack_cycles += 1 # switch charges a = mol.GetAtomWithIdx(anion) a.SetFormalCharge(a.GetFormalCharge() + 1) a = mol.GetAtomWithIdx(a2) a.SetFormalCharge(a.GetFormalCharge() - 1) # switch bond orders b = mol.GetBondBetweenAtoms(anion, a1) b.SetBondType(RDBONDORDER[b.GetBondTypeAsDouble() + 1]) b = mol.GetBondBetweenAtoms(a1, a2) b.SetBondType(RDBONDORDER[b.GetBondTypeAsDouble() - 1]) mol.UpdatePropertyCache(strict=False) # update number of matches for the end pattern n_matches = len({match[0] for match in matches}) if n_matches == 1: end_pattern = odd_end_pattern # start new iteration continue # no more changes to apply mol.UpdatePropertyCache(strict=False) return # reached max_iter warnings.warn( "The standardization could not be completed within a " "reasonable number of iterations" ) @dataclass(frozen=True) class TemplateInferrer: """Infer bond orders and charges by matching the molecule with a template molecule containing bond orders and charges. Attributes ---------- template : rdkit.Chem.rdchem.Mol Molecule containing the bond orders and charges. adjust_hydrogens: bool, default = True If ``True``, temporarily removes hydrogens from the molecule to assign bond orders and charges from the template, then adds them back. Useful to avoid adding explicit hydrogens on the template which can prevent RDKit from finding a match between the template and the molecule. Setting this to ``False`` can be useful to speed things up for inorganic molecules that don't have any hydrogens. .. versionadded:: 2.10.0 """ template: "Chem.Mol" adjust_hydrogens: bool = True def __call__(self, mol: "Chem.Mol") -> "Chem.Mol": if self.adjust_hydrogens: return self.assign_from_template_with_adjusted_hydrogens(mol) return AssignBondOrdersFromTemplate(self.template, mol) def assign_from_template_with_adjusted_hydrogens( self, mol: "Chem.Mol", index_field: str = "_MDAnalysis_index" ) -> "Chem.Mol": """Temporarily removes hydrogens from the molecule before assigning bond orders and charges from the template, then adds them back. Parameters ---------- mol : rdkit.Chem.rdchem.Mol Molecule to assign bond orders and charges to. index_field : str, default = "_MDAnalysis_index" Atom property corresponding to a unique integer that can used to put back the hydrogen atoms that were removed before matching, and sort them back to the original order. Must be accessible through ``atom.GetIntProp``. """ # remove explicit hydrogens and assign BO from template mol_no_h = Chem.RemoveAllHs(mol, sanitize=False) new = AssignBondOrdersFromTemplate(self.template, mol_no_h) # mapping between AtomGroup index of atom bearing H --> H atom(s) atoms_with_hydrogens = defaultdict(list) for atom in mol.GetAtoms(): if atom.GetAtomicNum() == 1: atoms_with_hydrogens[ atom.GetNeighbors()[0].GetIntProp(index_field) ].append(atom) # mapping between atom that should be bearing a H in RWMol and # corresponding H(s) reverse_mapping = {} for atom in new.GetAtoms(): if (idx := atom.GetIntProp(index_field)) in atoms_with_hydrogens: reverse_mapping[atom.GetIdx()] = atoms_with_hydrogens[idx] # Set ExplicitH to 0 to avoid valence errors when adding Hs atom.SetNumExplicitHs(0) # let UpdatePropertyCache recalculate the number of radicals later atom.SetNumRadicalElectrons(0) # add missing Hs rwmol = Chem.RWMol(new) for atom_idx, hydrogens in reverse_mapping.items(): for hydrogen in hydrogens: h_idx = rwmol.AddAtom(hydrogen) rwmol.AddBond(atom_idx, h_idx, Chem.BondType.SINGLE) new = rwmol.GetMol() # sanitize new.UpdatePropertyCache() sanitize_mol(new) # reorder atoms as input atomgroup (through _MDAnalysis_index) return reorder_atoms(new, field=index_field) @dataclass(frozen=True) class RDKitInferrer: """Uses RDKit's :func:`~rdkit.Chem.rdDetermineBonds.DetermineBondOrders` to infer bond orders and formal charges. This is the same algorithm used by the :ref:`xyz2mol ` package. .. versionadded:: 2.10.0 """ charge: int = 0 def __call__(self, mol: "Chem.Mol") -> "Chem.Mol": new = Chem.Mol(mol) DetermineBondOrders(new, charge=self.charge) return new