1 #
   2 # File: TorsionAlertsUtil.py
   3 # Author: Manish Sud <msud@san.rr.com>
   4 #
   5 # Copyright (C) 2026 Manish Sud. All rights reserved.
   6 #
   7 # This file is part of MayaChemTools.
   8 #
   9 # MayaChemTools is free software; you can redistribute it and/or modify it under
  10 # the terms of the GNU Lesser General Public License as published by the Free
  11 # Software Foundation; either version 3 of the License, or (at your option) any
  12 # later version.
  13 #
  14 # MayaChemTools is distributed in the hope that it will be useful, but without
  15 # any warranty; without even the implied warranty of merchantability of fitness
  16 # for a particular purpose.  See the GNU Lesser General Public License for more
  17 # details.
  18 #
  19 # You should have received a copy of the GNU Lesser General Public License
  20 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
  21 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
  22 # Boston, MA, 02111-1307, USA.
  23 #
  24 
  25 import sys
  26 import re
  27 import xml.etree.ElementTree as ET
  28 
  29 from rdkit import Chem
  30 
  31 __all__ = [
  32     "CalculateTorsionAngleDifference",
  33     "DoesSMARTSContainsMappedAtoms",
  34     "DoesSMARTSContainValidSubClassMappedAtoms",
  35     "DoesSMARTSContainValidTorsionRuleMappedAtoms",
  36     "FilterSubstructureMatchesByAtomMapNumbers",
  37     "GetAtomPositions",
  38     "GetGenericHierarchyClassElementNode",
  39     "GetHeavyAtomNeighbors",
  40     "IdentifyRotatableBondsForTorsionLibraryMatch",
  41     "IsSpecificHierarchyClass",
  42     "ListTorsionLibraryInfo",
  43     "RemoveLastHierarchyClassElementNodeFromTracking",
  44     "RemoveLastHierarchySubClassElementNodeFromTracking",
  45     "RetrieveTorsionLibraryInfo",
  46     "SetupHierarchyClassAndSubClassNamesForRotatableBond",
  47     "SetupHierarchySubClassElementPatternMol",
  48     "SetupTorsionRuleElementPatternMol",
  49     "SetupTorsionLibraryInfoForMatchingRotatableBonds",
  50     "TrackHierarchyClassElementNode",
  51     "TrackHierarchySubClassElementNode",
  52 ]
  53 
  54 
  55 def RetrieveTorsionLibraryInfo(TorsionLibraryFilePath, Quiet=True):
  56     """Retrieve torsion library information.
  57 
  58     Arguments:
  59         TorsionLibraryFilePath (str):  Torsion library XML file path.
  60 
  61     Returns:
  62         object: An object returned by xml.etree.ElementTree.parse function.
  63 
  64     Notes:
  65         The XML file is parsed using xml.etree.ElementTree.parse function and
  66         object created by the parse function is simply returned.
  67 
  68     """
  69 
  70     if not Quiet:
  71         _PrintInfoMsg("\nRetrieving data from torsion library file %s..." % TorsionLibraryFilePath)
  72 
  73     TorsionLibElementTree = ET.parse(TorsionLibraryFilePath)
  74 
  75     return TorsionLibElementTree
  76 
  77 
  78 def ListTorsionLibraryInfo(TorsionLibElementTree):
  79     """List torsion library information using XML tree object. The following
  80     information is listed:
  81 
  82         Summary:
  83 
  84             Total number of HierarchyClass nodes: <Number>
  85             Total number of HierarchyClassSubClass nodes: <Number
  86             Total number of TorsionRule nodes: <Number
  87 
  88         Details:
  89 
  90             HierarchyClass: <Name>; HierarchySubClass nodes: <Number>;
  91                 TorsionRule nodes: <SMARTS>
  92              ... ... ...
  93 
  94     Arguments:
  95         TorsionLibElementTree (object): XML tree object.
  96 
  97     Returns:
  98         Nothing.
  99 
 100     """
 101     HierarchyClassesInfo = {}
 102     HierarchyClassesInfo["HierarchyClassNames"] = []
 103     HierarchyClassesInfo["HierarchySubClassCount"] = {}
 104     HierarchyClassesInfo["TorsionRuleCount"] = {}
 105 
 106     HierarchyClassCount, HierarchySubClassCount, TorsionRuleCount = [0] * 3
 107 
 108     for HierarchyClassNode in TorsionLibElementTree.findall("hierarchyClass"):
 109         HierarchyClassCount += 1
 110         HierarchyClassName = HierarchyClassNode.get("name")
 111         if HierarchyClassName in HierarchyClassesInfo["HierarchyClassNames"]:
 112             _PrintWarningMsg("Hierarchy class name, %s, already exists..." % HierarchyClassName)
 113         HierarchyClassesInfo["HierarchyClassNames"].append(HierarchyClassName)
 114 
 115         SubClassCount = 0
 116         for HierarchySubClassNode in HierarchyClassNode.iter("hierarchySubClass"):
 117             SubClassCount += 1
 118         HierarchyClassesInfo["HierarchySubClassCount"][HierarchyClassName] = SubClassCount
 119         HierarchySubClassCount += SubClassCount
 120 
 121         RuleCount = 0
 122         for TorsionRuleNode in HierarchyClassNode.iter("torsionRule"):
 123             RuleCount += 1
 124 
 125         HierarchyClassesInfo["TorsionRuleCount"][HierarchyClassName] = RuleCount
 126         TorsionRuleCount += RuleCount
 127 
 128     _PrintInfoMsg("\nTotal number of HierarchyClass nodes: %s" % HierarchyClassCount)
 129     _PrintInfoMsg("Total number of HierarchyClassSubClass nodes: %s" % HierarchySubClassCount)
 130     _PrintInfoMsg("Total number of TorsionRule nodes: %s" % TorsionRuleCount)
 131 
 132     # List info for each hierarchyClass...
 133     _PrintInfoMsg("")
 134 
 135     # Generic class first...
 136     GenericClassName = "GG"
 137     if GenericClassName in HierarchyClassesInfo["HierarchyClassNames"]:
 138         _PrintInfoMsg(
 139             "HierarchyClass: %s; HierarchySubClass nodes: %s; TorsionRule nodes: %s"
 140             % (
 141                 GenericClassName,
 142                 HierarchyClassesInfo["HierarchySubClassCount"][GenericClassName],
 143                 HierarchyClassesInfo["TorsionRuleCount"][GenericClassName],
 144             )
 145         )
 146 
 147     for HierarchyClassName in sorted(HierarchyClassesInfo["HierarchyClassNames"]):
 148         if HierarchyClassName == GenericClassName:
 149             continue
 150         _PrintInfoMsg(
 151             "HierarchyClass: %s; HierarchySubClass nodes: %s; TorsionRule nodes: %s"
 152             % (
 153                 HierarchyClassName,
 154                 HierarchyClassesInfo["HierarchySubClassCount"][HierarchyClassName],
 155                 HierarchyClassesInfo["TorsionRuleCount"][HierarchyClassName],
 156             )
 157         )
 158 
 159 
 160 def SetupTorsionLibraryInfoForMatchingRotatableBonds(TorsionLibraryInfo):
 161     """Setup torsion  library information for matching rotatable bonds. The
 162     following information is initialized and updated in torsion library
 163     dictionary for matching rotatable bonds:
 164 
 165             TorsionLibraryInfo["GenericClass"] = None
 166             TorsionLibraryInfo["GenericClassElementNode"] = None
 167 
 168             TorsionLibraryInfo["SpecificClasses"] = {}
 169             TorsionLibraryInfo["SpecificClasses"]["Names"] = []
 170             TorsionLibraryInfo["SpecificClasses"]["ElementNode"] = {}
 171 
 172             TorsionLibraryInfo["HierarchyClassNodes"] = []
 173             TorsionLibraryInfo["HierarchySubClassNodes"] = []
 174 
 175             TorsionLibraryInfo["DataCache"] = {}
 176             TorsionLibraryInfo["DataCache"]["SubClassPatternMol"] = {}
 177 
 178             TorsionLibraryInfo["DataCache"]["TorsionRulePatternMol"] = {}
 179             TorsionLibraryInfo["DataCache"]["TorsionRuleAnglesInfo"] = {}
 180 
 181     Arguments:
 182         TorsionLibraryInfo (dict): A dictionary containing root node for
 183             torsion library element tree.
 184 
 185     Returns:
 186         Nonthing. The torsion library information dictionary is updated.
 187 
 188     """
 189     _SetupTorsionLibraryHierarchyClassesInfoForMatchingRotatableBonds(TorsionLibraryInfo)
 190     _SetupTorsionLibraryDataCacheInfoForMatchingRotatableBonds(TorsionLibraryInfo)
 191 
 192 
 193 def _SetupTorsionLibraryHierarchyClassesInfoForMatchingRotatableBonds(TorsionLibraryInfo):
 194     """Setup  hierarchy classes information for generic and specific classes."""
 195 
 196     RootElementNode = TorsionLibraryInfo["TorsionLibElementTree"]
 197 
 198     TorsionLibraryInfo["GenericClass"] = None
 199     TorsionLibraryInfo["GenericClassElementNode"] = None
 200 
 201     TorsionLibraryInfo["SpecificClasses"] = {}
 202     TorsionLibraryInfo["SpecificClasses"]["Names"] = []
 203     TorsionLibraryInfo["SpecificClasses"]["ElementNode"] = {}
 204 
 205     # Class name stacks for tracking names during processing of torsion rules..
 206     TorsionLibraryInfo["HierarchyClassNodes"] = []
 207     TorsionLibraryInfo["HierarchySubClassNodes"] = []
 208 
 209     ElementNames = []
 210     for ElementNode in RootElementNode.findall("hierarchyClass"):
 211         ElementName = ElementNode.get("name")
 212         if ElementName in ElementNames:
 213             _PrintWarningMsg("Hierarchy class name, %s, already exists. Ignoring duplicate name..." % ElementName)
 214             continue
 215 
 216         ElementNames.append(ElementName)
 217 
 218         if re.match("^GG$", ElementName, re.I):
 219             TorsionLibraryInfo["GenericClass"] = ElementName
 220             TorsionLibraryInfo["GenericClassElementNode"] = ElementNode
 221         else:
 222             TorsionLibraryInfo["SpecificClasses"]["Names"].append(ElementName)
 223             TorsionLibraryInfo["SpecificClasses"]["ElementNode"][ElementName] = ElementNode
 224 
 225 
 226 def _SetupTorsionLibraryDataCacheInfoForMatchingRotatableBonds(TorsionLibraryInfo):
 227     """Setup information for caching molecules for hierarchy subclass and torsion rule patterns."""
 228 
 229     TorsionLibElementTree = TorsionLibraryInfo["TorsionLibElementTree"]
 230 
 231     # Initialize data cache for pattern molecules corresponding to SMARTS patterns for
 232     # hierarchy subclasses and torsion rules. The pattern mols are generated and cached
 233     # later.
 234     TorsionLibraryInfo["DataCache"] = {}
 235     TorsionLibraryInfo["DataCache"]["SubClassPatternMol"] = {}
 236 
 237     TorsionLibraryInfo["DataCache"]["TorsionRulePatternMol"] = {}
 238     TorsionLibraryInfo["DataCache"]["TorsionRuleLonePairMapNumber"] = {}
 239     TorsionLibraryInfo["DataCache"]["TorsionRuleAnglesInfo"] = {}
 240 
 241     HierarchyClassID, HierarchySubClassID, TorsionRuleID = [0] * 3
 242 
 243     for HierarchyClassNode in TorsionLibElementTree.findall("hierarchyClass"):
 244         HierarchyClassID += 1
 245 
 246         for HierarchySubClassNode in HierarchyClassNode.iter("hierarchySubClass"):
 247             HierarchySubClassID += 1
 248             # Add unique ID to node...
 249             HierarchySubClassNode.set("NodeID", HierarchySubClassID)
 250 
 251         for TorsionRuleNode in HierarchyClassNode.iter("torsionRule"):
 252             TorsionRuleID += 1
 253             # Add unique ID to node...
 254             TorsionRuleNode.set("NodeID", TorsionRuleID)
 255 
 256 
 257 def IdentifyRotatableBondsForTorsionLibraryMatch(TorsionLibraryInfo, Mol, RotBondsPatternMol):
 258     """Identify rotatable bonds in a molecule for torsion library match.
 259 
 260     Arguments:
 261         TorsionLibraryInfo (dict): A dictionary containing information for
 262             matching rotatable bonds.
 263         Mol (object): RDKit molecule object.
 264         RotBondsPatternMol (object): RDKit molecule object for SMARTS pattern
 265             corresponding to rotatable bonds.
 266 
 267     Returns:
 268         bool: True - Rotatable bonds present in molecule; Otherwise, false.
 269         None or dict: None - For no rotatable bonds in molecule; otherwise, a
 270             dictionary containing the following informations for rotatable bonds
 271             matched to RotBondsPatternMol:
 272 
 273                 RotBondsInfo["IDs"] = []
 274                 RotBondsInfo["AtomIndices"] = {}
 275                 RotBondsInfo["HierarchyClass"] = {}
 276 
 277     """
 278 
 279     # Match rotatable bonds...
 280     RotBondsMatches = FilterSubstructureMatchesByAtomMapNumbers(
 281         Mol, RotBondsPatternMol, Mol.GetSubstructMatches(RotBondsPatternMol, useChirality=False)
 282     )
 283 
 284     #  Check and filter rotatable bond matches...
 285     RotBondsMatches = _FilterRotatableBondMatches(Mol, RotBondsMatches)
 286 
 287     if not len(RotBondsMatches):
 288         return False, None
 289 
 290     # Initialize rotatable bonds info...
 291     RotBondsInfo = {}
 292     RotBondsInfo["IDs"] = []
 293     RotBondsInfo["AtomIndices"] = {}
 294     RotBondsInfo["HierarchyClass"] = {}
 295 
 296     # Setup rotatable bonds info...
 297     ID = 0
 298     for RotBondAtomIndices in RotBondsMatches:
 299         ID += 1
 300 
 301         RotBondAtoms = [Mol.GetAtomWithIdx(RotBondAtomIndices[0]), Mol.GetAtomWithIdx(RotBondAtomIndices[1])]
 302         RotBondAtomSymbols = [RotBondAtoms[0].GetSymbol(), RotBondAtoms[1].GetSymbol()]
 303 
 304         ClassID = "%s%s" % (RotBondAtomSymbols[0], RotBondAtomSymbols[1])
 305         if ClassID not in TorsionLibraryInfo["SpecificClasses"]["Names"]:
 306             ReverseClassID = "%s%s" % (RotBondAtomSymbols[1], RotBondAtomSymbols[0])
 307             if ReverseClassID in TorsionLibraryInfo["SpecificClasses"]["Names"]:
 308                 ClassID = ReverseClassID
 309                 # Reverse atom indices and related information...
 310                 RotBondAtomIndices = list(reversed(RotBondAtomIndices))
 311                 RotBondAtoms = list(reversed(RotBondAtoms))
 312                 RotBondAtomSymbols = list(reversed(RotBondAtomSymbols))
 313 
 314         # Track information...
 315         RotBondsInfo["IDs"].append(ID)
 316         RotBondsInfo["AtomIndices"][ID] = RotBondAtomIndices
 317         RotBondsInfo["HierarchyClass"][ID] = ClassID
 318 
 319     return True, RotBondsInfo
 320 
 321 
 322 def FilterSubstructureMatchesByAtomMapNumbers(Mol, PatternMol, AtomIndicesList):
 323     """Filter a list of lists containing matched atom indices by map atom numbers
 324     present in a pattern molecule. The list of atom indices correspond to a list retrieved by
 325     RDKit function GetSubstructureMatches using SMILES/SMARTS pattern. The
 326     atom map numbers are mapped to appropriate atom indices during the generation
 327     of molecules. For example: [O:1]=[S:2](=[O])[C:3][C:4].
 328 
 329     Arguments:
 330         Mol (object): RDKit molecule object.
 331         PatternMol (object): RDKit molecule object for a SMILES/SMARTS pattern.
 332         AtomIndicesList (list): A list of lists containing atom indices.
 333 
 334     Returns:
 335         list : A list of lists containing filtered atom indices.
 336 
 337     """
 338 
 339     AtomMapIndices = _GetAtomMapIndices(PatternMol)
 340 
 341     MatchedAtomIndicesList = []
 342     for AtomIndices in AtomIndicesList:
 343         MatchedAtomIndicesList.append(
 344             _FilterSubstructureMatchByAtomMapNumbers(Mol, PatternMol, AtomIndices, AtomMapIndices)
 345         )
 346 
 347     return MatchedAtomIndicesList
 348 
 349 
 350 def _GetAtomMapIndices(Mol):
 351     """Get a list of available atom indices corresponding to sorted atom map
 352     numbers present in a SMILES/SMARTS pattern used for creating a molecule.
 353     """
 354 
 355     AtomMapIndices, AtomMapNumbers = _GetAtomMapIndicesAndMapNumbers(Mol)
 356 
 357     return AtomMapIndices
 358 
 359 
 360 def _GetAtomMapIndicesAndMapNumbers(Mol):
 361     """Get a list of available atom indices and atom map numbers present
 362     in  a SMILES/SMARTS pattern used for creating a molecule. Both lists
 363     are sorted in ascending order by atom map numbers.
 364     """
 365 
 366     # Setup a atom map number to atom indices map..
 367     AtomMapNumToIndices = {}
 368     for Atom in Mol.GetAtoms():
 369         AtomMapNum = Atom.GetAtomMapNum()
 370 
 371         if AtomMapNum:
 372             AtomMapNumToIndices[AtomMapNum] = Atom.GetIdx()
 373 
 374     # Setup atom indices corresponding to sorted atom map numbers...
 375     AtomMapIndices = None
 376     AtomMapNumbers = None
 377     if len(AtomMapNumToIndices):
 378         AtomMapNumbers = sorted(AtomMapNumToIndices)
 379         AtomMapIndices = [AtomMapNumToIndices[AtomMapNum] for AtomMapNum in AtomMapNumbers]
 380 
 381     return (AtomMapIndices, AtomMapNumbers)
 382 
 383 
 384 def _FilterSubstructureMatchByAtomMapNumbers(Mol, PatternMol, AtomIndices, AtomMapIndices):
 385     """Filter substructure match atom indices by atom map indices corresponding to
 386     atom map numbers.
 387     """
 388 
 389     if AtomMapIndices is None:
 390         return list(AtomIndices)
 391 
 392     return [AtomIndices[Index] for Index in AtomMapIndices]
 393 
 394 
 395 def _FilterRotatableBondMatches(Mol, RotBondsMatches):
 396     """Filter rotatable bond matches to ensure that each rotatable bond atom
 397     is attached to at least two heavy atoms. Otherwise, the torsion rules might match
 398     hydrogens."""
 399 
 400     FilteredRotBondMatches = []
 401 
 402     # Go over rotatable bonds...
 403     for RotBondMatch in RotBondsMatches:
 404         SkipRotBondMatch = False
 405         for AtomIndex in RotBondMatch:
 406             Atom = Mol.GetAtomWithIdx(AtomIndex)
 407 
 408             HeavyAtomNbrCount = len(GetHeavyAtomNeighbors(Atom))
 409             if HeavyAtomNbrCount <= 1:
 410                 SkipRotBondMatch = True
 411                 break
 412 
 413         if not SkipRotBondMatch:
 414             FilteredRotBondMatches.append(RotBondMatch)
 415 
 416     return FilteredRotBondMatches
 417 
 418 
 419 def SetupHierarchySubClassElementPatternMol(TorsionLibraryInfo, ElementNode):
 420     """Setup pattern molecule for SMARTS pattern in hierarchy subclass element.
 421 
 422     Arguments:
 423         TorsionLibraryInfo (dict): A dictionary containing information for
 424             matching rotatable bonds.
 425         ElementNode (object): A hierarchy sub class element node being matched
 426            in torsion library XML tree.
 427 
 428     Returns:
 429         object: RDKit molecule object corresponding to SMARTS pattern for
 430             hierarchy sub class element node.
 431 
 432     """
 433 
 434     # Check data cache...
 435     SubClassNodeID = ElementNode.get("NodeID")
 436     if SubClassNodeID in TorsionLibraryInfo["DataCache"]["SubClassPatternMol"]:
 437         return TorsionLibraryInfo["DataCache"]["SubClassPatternMol"][SubClassNodeID]
 438 
 439     # Setup and track pattern mol...
 440     SubClassSMARTSPattern = ElementNode.get("smarts")
 441     SubClassPatternMol = Chem.MolFromSmarts(SubClassSMARTSPattern)
 442 
 443     if SubClassPatternMol is None:
 444         _PrintWarningMsg(
 445             "Ignoring hierachical subclass, %s, containing invalid SMARTS pattern %s"
 446             % (ElementNode.get("name"), SubClassSMARTSPattern)
 447         )
 448 
 449     if not DoesSMARTSContainValidSubClassMappedAtoms(SubClassSMARTSPattern):
 450         SubClassPatternMol = None
 451         _PrintWarningMsg(
 452             "Ignoring hierachical subclass, %s, containing invalid map atom numbers in SMARTS pattern %s"
 453             % (ElementNode.get("name"), SubClassSMARTSPattern)
 454         )
 455 
 456     TorsionLibraryInfo["DataCache"]["SubClassPatternMol"][SubClassNodeID] = SubClassPatternMol
 457 
 458     return SubClassPatternMol
 459 
 460 
 461 def SetupTorsionRuleElementPatternMol(TorsionLibraryInfo, ElementNode, TorsionRuleNodeID, TorsionSMARTSPattern):
 462     """Setup pattern molecule for SMARTS pattern in torsion rule element.
 463 
 464     Arguments:
 465         TorsionLibraryInfo (dict): A dictionary containing information for
 466             matching rotatable bonds.
 467         ElementNode (object): A torsion rule element node being matched in
 468            torsion library XML tree.
 469         TorsionRuleNodeID (int): Torsion rule element node ID.
 470         TorsionSMARTSPattern (str): SMARTS pattern for torsion rule element node.
 471 
 472     Returns:
 473         object: RDKit molecule object corresponding to SMARTS pattern for
 474             torsion rule element node.
 475 
 476     """
 477 
 478     if TorsionRuleNodeID in TorsionLibraryInfo["DataCache"]["TorsionRulePatternMol"]:
 479         return TorsionLibraryInfo["DataCache"]["TorsionRulePatternMol"][TorsionRuleNodeID]
 480 
 481     TorsionPatternMol = Chem.MolFromSmarts(TorsionSMARTSPattern)
 482     if TorsionPatternMol is None:
 483         _PrintWarningMsg("Ignoring torsion rule element containing invalid SMARTS pattern %s" % TorsionSMARTSPattern)
 484 
 485     if not DoesSMARTSContainValidTorsionRuleMappedAtoms(TorsionSMARTSPattern):
 486         TorsionPatternMol = None
 487         _PrintWarningMsg(
 488             "Ignoring torsion rule element containing invalid map atoms numbers in SMARTS pattern %s"
 489             % TorsionSMARTSPattern
 490         )
 491 
 492     TorsionLibraryInfo["DataCache"]["TorsionRulePatternMol"][TorsionRuleNodeID] = TorsionPatternMol
 493 
 494     return TorsionPatternMol
 495 
 496 
 497 def SetupHierarchyClassAndSubClassNamesForRotatableBond(TorsionLibraryInfo):
 498     """Setup hierarchy class and subclass names for a rotatable bond matched to
 499     a torsion rule element node.
 500 
 501     Returns:
 502         TorsionLibraryInfo (dict): A dictionary containing information for
 503             matching rotatable bonds.
 504 
 505     Returns:
 506         str: A back slash delimited string containing hierarchy class names at
 507             the level of torsion rule element node.
 508         str: A back slash delimited string containing hierarchy sub class names
 509           at the level of torsion rule element node.
 510 
 511     """
 512 
 513     HierarchyClassName, HierarchyClassSubName = ["None"] * 2
 514 
 515     # Setup hierarchy class name...
 516     if len(TorsionLibraryInfo["HierarchyClassNodes"]):
 517         HierarchyClassElementNode = TorsionLibraryInfo["HierarchyClassNodes"][-1]
 518         HierarchyClassName = HierarchyClassElementNode.get("name")
 519         if len(HierarchyClassName) == 0:
 520             HierarchyClassName = "None"
 521 
 522     # Setup hierarchy class name...
 523     if len(TorsionLibraryInfo["HierarchySubClassNodes"]):
 524         HierarchySubClassNames = []
 525         for ElementNode in TorsionLibraryInfo["HierarchySubClassNodes"]:
 526             Name = ElementNode.get("name")
 527             if len(Name) == 0:
 528                 Name = "None"
 529             HierarchySubClassNames.append(Name)
 530 
 531         HierarchyClassSubName = "/".join(HierarchySubClassNames)
 532 
 533     # Replace spaces by underscores in class and subclass names...
 534     if HierarchyClassName is not None:
 535         if " " in HierarchyClassName:
 536             HierarchyClassName = HierarchyClassName.replace(" ", "_")
 537 
 538     if HierarchyClassSubName is not None:
 539         if " " in HierarchyClassSubName:
 540             HierarchyClassSubName = HierarchyClassSubName.replace(" ", "_")
 541 
 542     return (HierarchyClassName, HierarchyClassSubName)
 543 
 544 
 545 def SetupTorsionRuleAnglesInfo(TorsionLibraryInfo, TorsionRuleElementNode):
 546     """Setup torsion angles and energy info for matching a torsion rule.
 547 
 548     Arguments:
 549         TorsionLibraryInfo (dict): A dictionary containing information for
 550             matching rotatable bonds.
 551         TorsionRuleElementNode (object): A torsion rule element node being
 552            matched in torsion library XML tree.
 553 
 554     Returns:
 555         dict: A dictionary containing the following information for torsion rule
 556             being matched to a rotatable bond:
 557 
 558             RuleAnglesInfo = {}
 559 
 560             RuleAnglesInfo["IDs"] = []
 561             RuleAnglesInfo["Value"] = {}
 562             RuleAnglesInfo["Score"] = {}
 563             RuleAnglesInfo["Tolerance1"] = {}
 564             RuleAnglesInfo["Tolerance2"] = {}
 565 
 566             RuleAnglesInfo["ValuesList"] = []
 567             RuleAnglesInfo["ValuesIn360RangeList"] = []
 568             RuleAnglesInfo["Tolerances1List"] = []
 569             RuleAnglesInfo["Tolerances2List"] = []
 570 
 571             # Strain energy calculations...
 572             RuleAnglesInfo["EnergyMethod"] = None
 573             RuleAnglesInfo["EnergyMethodExact"] = None
 574             RuleAnglesInfo["EnergyMethodApproximate"] = None
 575 
 576             # For approximate strain energy calculation...
 577             RuleAnglesInfo["Beta1"] = {}
 578             RuleAnglesInfo["Beta2"] = {}
 579             RuleAnglesInfo["Theta0"] = {}
 580 
 581             # For exact strain energy calculation...
 582             RuleAnglesInfo["HistogramEnergy"] = []
 583             RuleAnglesInfo["HistogramEnergyLowerBound"] = []
 584             RuleAnglesInfo["HistogramEnergyUpperBound"] = []
 585 
 586     """
 587 
 588     # Check data cache...
 589     TorsionRuleNodeID = TorsionRuleElementNode.get("NodeID")
 590     if TorsionRuleNodeID in TorsionLibraryInfo["DataCache"]["TorsionRuleAnglesInfo"]:
 591         return TorsionLibraryInfo["DataCache"]["TorsionRuleAnglesInfo"][TorsionRuleNodeID]
 592 
 593     # Initialize rule angles info...
 594     RuleAnglesInfo = {}
 595 
 596     RuleAnglesInfo["IDs"] = []
 597     RuleAnglesInfo["Value"] = {}
 598     RuleAnglesInfo["Score"] = {}
 599     RuleAnglesInfo["Tolerance1"] = {}
 600     RuleAnglesInfo["Tolerance2"] = {}
 601 
 602     RuleAnglesInfo["ValuesList"] = []
 603     RuleAnglesInfo["ValuesIn360RangeList"] = []
 604     RuleAnglesInfo["Tolerances1List"] = []
 605     RuleAnglesInfo["Tolerances2List"] = []
 606 
 607     # Strain energy calculations...
 608     RuleAnglesInfo["EnergyMethod"] = None
 609     RuleAnglesInfo["EnergyMethodExact"] = None
 610     RuleAnglesInfo["EnergyMethodApproximate"] = None
 611 
 612     # For approximate strain energy calculation....
 613     RuleAnglesInfo["Beta1"] = {}
 614     RuleAnglesInfo["Beta2"] = {}
 615     RuleAnglesInfo["Theta0"] = {}
 616 
 617     # For exact strain energy calculation...
 618     RuleAnglesInfo["HistogramEnergy"] = []
 619     RuleAnglesInfo["HistogramEnergyLowerBound"] = []
 620     RuleAnglesInfo["HistogramEnergyUpperBound"] = []
 621 
 622     # Setup strain energy calculation information...
 623     EnergyMethod, EnergyMethodExact, EnergyMethodApproximate = [None] * 3
 624     EnergyMethod = TorsionRuleElementNode.get("method")
 625     if EnergyMethod is not None:
 626         EnergyMethodExact = True if re.match("^exact$", EnergyMethod, re.I) else False
 627         EnergyMethodApproximate = True if re.match("^approximate$", EnergyMethod, re.I) else False
 628 
 629     RuleAnglesInfo["EnergyMethod"] = EnergyMethod
 630     RuleAnglesInfo["EnergyMethodExact"] = EnergyMethodExact
 631     RuleAnglesInfo["EnergyMethodApproximate"] = EnergyMethodApproximate
 632 
 633     # Setup angles information....
 634     AngleID = 0
 635     for AngleListElementNode in TorsionRuleElementNode.findall("angleList"):
 636         for AngleNode in AngleListElementNode.iter("angle"):
 637             AngleID += 1
 638             Value = float(AngleNode.get("value"))
 639             Tolerance1 = float(AngleNode.get("tolerance1"))
 640             Tolerance2 = float(AngleNode.get("tolerance2"))
 641             Score = float(AngleNode.get("score"))
 642 
 643             # Track values...
 644             RuleAnglesInfo["IDs"].append(AngleID)
 645             RuleAnglesInfo["Value"][AngleID] = Value
 646             RuleAnglesInfo["Score"][AngleID] = Score
 647             RuleAnglesInfo["Tolerance1"][AngleID] = Tolerance1
 648             RuleAnglesInfo["Tolerance2"][AngleID] = Tolerance2
 649 
 650             RuleAnglesInfo["ValuesList"].append(Value)
 651             RuleAnglesInfo["Tolerances1List"].append(Tolerance1)
 652             RuleAnglesInfo["Tolerances2List"].append(Tolerance2)
 653 
 654             # Map value to 0 to 360 range...
 655             MappedValue = Value + 360 if Value < 0 else Value
 656             RuleAnglesInfo["ValuesIn360RangeList"].append(MappedValue)
 657 
 658             # Approximate strain energy calculation information...
 659             if EnergyMethodApproximate:
 660                 Beta1 = float(AngleNode.get("beta_1"))
 661                 Beta2 = float(AngleNode.get("beta_2"))
 662                 Theta0 = float(AngleNode.get("theta_0"))
 663 
 664                 RuleAnglesInfo["Beta1"][AngleID] = Beta1
 665                 RuleAnglesInfo["Beta2"][AngleID] = Beta2
 666                 RuleAnglesInfo["Theta0"][AngleID] = Theta0
 667 
 668     #  Exact energy method  information...
 669     if EnergyMethodExact:
 670         for HistogramBinNode in TorsionRuleElementNode.find("histogram_converted").iter("bin"):
 671             Energy = float(HistogramBinNode.get("energy"))
 672             Lower = float(HistogramBinNode.get("lower"))
 673             Upper = float(HistogramBinNode.get("upper"))
 674 
 675             RuleAnglesInfo["HistogramEnergy"].append(Energy)
 676             RuleAnglesInfo["HistogramEnergyLowerBound"].append(Lower)
 677             RuleAnglesInfo["HistogramEnergyUpperBound"].append(Upper)
 678 
 679     if len(RuleAnglesInfo["IDs"]) == 0:
 680         RuleAnglesInfo = None
 681 
 682     # Cache data...
 683     TorsionLibraryInfo["DataCache"]["TorsionRuleAnglesInfo"][TorsionRuleNodeID] = RuleAnglesInfo
 684 
 685     return RuleAnglesInfo
 686 
 687 
 688 def DoesSMARTSContainValidSubClassMappedAtoms(SMARTS):
 689     """Check for the presence of two central mapped atoms in SMARTS pattern.
 690     A valid SMARTS pattern must contain only two mapped atoms corresponding
 691     to map atom numbers ':2' and ':3'.
 692 
 693     Arguments:
 694         SMARTS (str): SMARTS pattern for sub class in torsion library XML tree.
 695 
 696     Returns:
 697         bool: True - A valid pattern; Otherwise, false.
 698 
 699     """
 700 
 701     MatchedMappedAtoms = re.findall(":[0-9]", SMARTS, re.I)
 702     if len(MatchedMappedAtoms) < 2 or len(MatchedMappedAtoms) > 4:
 703         return False
 704 
 705     # Check for the presence of two central atom map numbers for a torsion...
 706     for MapAtomNum in [":2", ":3"]:
 707         if MapAtomNum not in MatchedMappedAtoms:
 708             return False
 709 
 710     return True
 711 
 712 
 713 def DoesSMARTSContainValidTorsionRuleMappedAtoms(SMARTS):
 714     """Check for the presence of four mapped atoms in a SMARTS pattern.
 715     A valid SMARTS pattern must contain only four mapped atoms corresponding
 716     to map atom numbers ':1', ':2', ':3' and ':4'.
 717 
 718     Arguments:
 719         SMARTS (str): SMARTS pattern for torsion rule in torsion library XML
 720             tree.
 721 
 722     Returns:
 723         bool: True - A valid pattern; Otherwise, false.
 724 
 725     """
 726 
 727     MatchedMappedAtoms = re.findall(":[0-9]", SMARTS, re.I)
 728     if len(MatchedMappedAtoms) != 4:
 729         return False
 730 
 731     # Check for the presence of four atom map numbers for a torsion...
 732     for MapAtomNum in [":1", ":2", ":3", ":4"]:
 733         if MapAtomNum not in MatchedMappedAtoms:
 734             return False
 735 
 736     return True
 737 
 738 
 739 def DoesSMARTSContainsMappedAtoms(SMARTS, MappedAtomNumsList):
 740     """Check for the presence of specified mapped atoms in SMARTS pattern.
 741     The mapped atom numbers in the list are specified as ':1', ':2', ':3' etc.
 742 
 743     Arguments:
 744         SMARTS (str): SMARTS pattern in torsion library XML tree.
 745         MappedAtoms (list): Mapped atom numbers as ":1", ":2" etc.
 746 
 747     Returns:
 748         bool: True - All mapped atoms present in pattern; Otherwise, false.
 749 
 750     """
 751 
 752     MatchedMappedAtoms = re.findall(":[0-9]", SMARTS, re.I)
 753     if len(MatchedMappedAtoms) == 0:
 754         return False
 755 
 756     # Check for the presence of specified mapped atoms in pattern...
 757     for MapAtomNum in MappedAtomNumsList:
 758         if MapAtomNum not in MatchedMappedAtoms:
 759             return False
 760 
 761     return True
 762 
 763 
 764 def IsSpecificHierarchyClass(TorsionLibraryInfo, HierarchyClass):
 765     """Check whether it's a specific hierarchy class.
 766 
 767     Arguments:
 768         TorsionLibraryInfo (dict): A dictionary containing information for
 769             matching rotatable bonds.
 770         HierarchyClass (str): Hierarchy class name.
 771 
 772     Returns:
 773         bool: True - A valid hierarchy class name; Otherwise, false.
 774 
 775     """
 776     return True if HierarchyClass in TorsionLibraryInfo["SpecificClasses"]["ElementNode"] else False
 777 
 778 
 779 def GetGenericHierarchyClassElementNode(TorsionLibraryInfo):
 780     """Get generic hierarchy class element node.
 781 
 782     Arguments:
 783         TorsionLibraryInfo (dict): A dictionary containing information for
 784             matching rotatable bonds.
 785 
 786     Returns:
 787         object: Generic hierarchy class element node in torsion library XML
 788             tree.
 789 
 790     """
 791     return TorsionLibraryInfo["GenericClassElementNode"]
 792 
 793 
 794 def TrackHierarchyClassElementNode(TorsionLibraryInfo, ElementNode):
 795     """Track hierarchy class element node using a stack.
 796 
 797     Arguments:
 798         TorsionLibraryInfo (dict): A dictionary containing information for
 799             matching rotatable bonds.
 800         ElementNode (object): Hierarchy class element node in torsion library
 801             XML tree.
 802 
 803     Returns:
 804         Nothing. The torsion library info is updated.
 805 
 806     """
 807     TorsionLibraryInfo["HierarchyClassNodes"].append(ElementNode)
 808 
 809 
 810 def RemoveLastHierarchyClassElementNodeFromTracking(TorsionLibraryInfo):
 811     """Remove last hierarchy class element node from tracking by removing it
 812     from a stack.
 813 
 814     Arguments:
 815         TorsionLibraryInfo (dict): A dictionary containing information for
 816             matching rotatable bonds.
 817 
 818     Returns:
 819         Nothing. The torsion library info is updated.
 820 
 821     """
 822     TorsionLibraryInfo["HierarchyClassNodes"].pop()
 823 
 824 
 825 def TrackHierarchySubClassElementNode(TorsionLibraryInfo, ElementNode):
 826     """Track hierarchy sub class element node using a stack.
 827 
 828     Arguments:
 829         TorsionLibraryInfo (dict): A dictionary containing information for
 830             matching rotatable bonds.
 831         ElementNode (object): Hierarchy sub class element node in torsion
 832             library XML tree.
 833 
 834     Returns:
 835         Nothing. The torsion library info is updated.
 836 
 837     """
 838     TorsionLibraryInfo["HierarchySubClassNodes"].append(ElementNode)
 839 
 840 
 841 def RemoveLastHierarchySubClassElementNodeFromTracking(TorsionLibraryInfo):
 842     """Remove last hierarchy sub class element node from tracking by removing it
 843     from a stack.
 844 
 845     Arguments:
 846         TorsionLibraryInfo (dict): A dictionary containing information for
 847             matching rotatable bonds.
 848 
 849     Returns:
 850         Nothing. The torsion library info is updated.
 851 
 852     """
 853     TorsionLibraryInfo["HierarchySubClassNodes"].pop()
 854 
 855 
 856 def CalculateTorsionAngleDifference(TorsionAngle1, TorsionAngle2):
 857     """Calculate torsion angle difference in the range from 0 to 180.
 858 
 859     Arguments:
 860         TorsionAngle1 (float): First torsion angle.
 861         TorsionAngle2 (float): Second torsion angle.
 862 
 863     Returns:
 864         float: Difference between first and second torsion angle.
 865 
 866     """
 867 
 868     # Map angles to 0 to 360 range...
 869     if TorsionAngle1 < 0:
 870         TorsionAngle1 = TorsionAngle1 + 360
 871     if TorsionAngle2 < 0:
 872         TorsionAngle2 = TorsionAngle2 + 360
 873 
 874     # Calculate and map angle difference in the range from 0 to 180 range...
 875     TorsionAngleDiff = abs(TorsionAngle1 - TorsionAngle2)
 876     if TorsionAngleDiff > 180.0:
 877         TorsionAngleDiff = abs(TorsionAngleDiff - 360)
 878 
 879     return TorsionAngleDiff
 880 
 881 
 882 def _PrintInfoMsg(Msg=""):
 883     """Print message to stderr along with flushing stderr."""
 884 
 885     print(Msg, sep=" ", end="\n", file=sys.stderr)
 886     sys.stderr.flush()
 887 
 888 
 889 def _PrintWarningMsg(Msg):
 890     """Print message to stderr along with flushing stderr. An `Warning` prefix
 891     is placed before the message."""
 892 
 893     _PrintInfoMsg("Warning: %s" % Msg)
 894 
 895 
 896 def GetHeavyAtomNeighbors(Atom):
 897     """Get a list of heavy atom neighbors.
 898 
 899     Arguments:
 900         Atom (object): RDKit atom object.
 901 
 902     Returns:
 903         list : List of heavy atom neighbors.
 904 
 905     """
 906 
 907     AtomNeighbors = []
 908     for AtomNbr in Atom.GetNeighbors():
 909         if AtomNbr.GetAtomicNum() > 1:
 910             AtomNeighbors.append(AtomNbr)
 911 
 912     return AtomNeighbors
 913 
 914 
 915 def GetAtomPositions(Mol, ConfID=-1):
 916     """Retrieve a list of lists containing coordinates of all atoms in a
 917     molecule.
 918 
 919     Arguments:
 920         Mol (object): RDKit molecule object.
 921         ConfID (int): Conformer number.
 922 
 923     Returns:
 924         list : List of lists containing atom positions.
 925 
 926     Examples:
 927 
 928         for AtomPosition in TorsionAlertsUtil.GetAtomPositions(Mol):
 929             print("X: %s; Y: %s; Z: %s" % (AtomPosition[0], AtomPosition[1], AtomPosition[2]))
 930 
 931     """
 932 
 933     return Mol.GetConformer(id=ConfID).GetPositions().tolist()