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