1 #
2 # File: TorsionLibraryAlerts.py
3 # Author: Manish Sud <msud@san.rr.com>
4 #
5 # Collaborator: Pat Walters
6 #
7 # Acknowledgments: Wolfgang Guba, Patrick Penner, and Levi Pierce
8 #
9 # Copyright (C) 2026 Manish Sud. All rights reserved.
10 #
11 # This module uses the Torsion Library jointly developed by the University
12 # of Hamburg, Center for Bioinformatics, Hamburg, Germany and
13 # F. Hoffmann-La-Roche Ltd., Basel, Switzerland.
14 #
15 # This file is part of MayaChemTools.
16 #
17 # MayaChemTools is free software; you can redistribute it and/or modify it under
18 # the terms of the GNU Lesser General Public License as published by the Free
19 # Software Foundation; either version 3 of the License, or (at your option) any
20 # later version.
21 #
22 # MayaChemTools is distributed in the hope that it will be useful, but without
23 # any warranty; without even the implied warranty of merchantability of fitness
24 # for a particular purpose. See the GNU Lesser General Public License for more
25 # details.
26 #
27 # You should have received a copy of the GNU Lesser General Public License
28 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
29 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
30 # Boston, MA, 02111-1307, USA.
31 #
32
33 import os
34 import re
35 import math
36 import numpy as np
37
38 from rdkit import Chem
39 from rdkit.Chem import rdMolTransforms
40
41 from . import TorsionAlertsUtil
42
43
44 class TorsionLibraryAlerts:
45 def __init__(
46 self,
47 AlertsMode="Red",
48 MinAlertsCount=1,
49 NitrogenLonePairAllowHydrogenNbrs=True,
50 NitrogenLonePairPlanarityTolerance=1.0,
51 RotBondsSMARTSMode="SemiStrict",
52 RotBondsSMARTSPattern=None,
53 TorsionLibraryFilePath="auto",
54 ):
55 """Identify strained molecules from an input file for torsion library [ Ref 146, 152, 159 ]
56 alerts by matching rotatable bonds against SMARTS patterns specified for torsion
57 rules in a torsion library file, The molecules must have 3D coordinates. The default
58 torsion library file, TorsionLibrary.xml, is available in the directory containing this file.
59
60 The data in torsion library file is organized in a hierarchical manner. It consists
61 of one generic class and six specific classes at the highest level. Each class
62 contains multiple subclasses corresponding to named functional groups or
63 substructure patterns. The subclasses consist of torsion rules sorted from
64 specific to generic torsion patterns. The torsion rule, in turn, contains a list
65 of peak values for torsion angles and two tolerance values. A pair of tolerance
66 values define torsion bins around a torsion peak value. For example:
67
68 <library>
69 <hierarchyClass name="GG" id1="G" id2="G">
70 ...
71 </hierarchyClass>
72 <hierarchyClass name="CO" id1="C" id2="O">
73 <hierarchySubClass name="Ester bond I" smarts="O=[C:2][O:3]">
74 <torsionRule smarts="[O:1]=[C:2]!@[O:3]~[CH0:4]">
75 <angleList>
76 <angle value="0.0" tolerance1="20.00"
77 tolerance2="25.00" score="56.52"/>
78 </angleList>
79 </torsionRule>
80 ...
81 ...
82 ...
83 </hierarchyClass>
84 <hierarchyClass name="NC" id1="N" id2="C">
85 ...
86 </hierarchyClass>
87 <hierarchyClass name="SN" id1="S" id2="N">
88 ...
89 </hierarchyClass>
90 <hierarchyClass name="CS" id1="C" id2="S">
91 ...
92 </hierarchyClass>
93 <hierarchyClass name="CC" id1="C" id2="C">
94 ...
95 </hierarchyClass>
96 <hierarchyClass name="SS" id1="S" id2="S">
97 ...
98 </hierarchyClass>
99 </library>
100
101 The rotatable bonds in a 3D molecule are identified using a default SMARTS pattern.
102 A custom SMARTS pattern may be optionally specified to detect rotatable bonds.
103 Each rotatable bond is matched to a torsion rule in the torsion library and
104 assigned one of the following three alert categories: Green, Orange or Red. The
105 rotatable bond is marked Green or Orange for the measured angle of the torsion
106 pattern within the first or second tolerance bins around a torsion peak.
107 Otherwise, it's marked Red implying that the measured angle is not observed in
108 the structure databases employed to generate the torsion library.
109
110 Arguments:
111 AlertsMode(str): Torsion library alert types to use for issuing
112 alerts about molecules.
113 MinAlertsCount(int): Minimum number of alerts allowed in a molecule.
114 NitrogenLonePairAllowHydrogenNbrs (bool): Use hydrogen neighbors
115 attached to nitrogen during the determination of its planarity.
116 NitrogenLonePairPlanarityTolerance (float): Angle tolerance in
117 degrees allowed for nitrogen to be considered coplanar with its
118 three neighbors.
119 RotBondsSMARTSMode (str): SMARTS pattern to use for identifying
120 rotatable bonds in a molecule. Possible values: NonStrict,
121 SemiStrict, Strict or Specify.
122 RotBondsSMARTSPattern (str):SMARTS pattern for identifying rotatable
123 bonds. This paramater is only valid for 'Specify' value
124 'RotBondsSMARTSMode'.
125 TorsionLibraryFilePath (str): A XML file name containing data for
126 torsion library.
127
128 Returns:
129 object: An instantiated class object.
130
131 Notes:
132 The following sections provide additional details for the parameters.
133
134 AlertsMode: Torsion library alert types to use for issuing alerts about
135 molecules containing rotatable bonds marked with Green, Orange, or
136 Red alerts. Possible values: Red or RedAndOrange.
137
138 MinAlertsCount: Minimum number of rotatable bond alerts allowed in a
139 molecule.
140
141 NitrogenLonePairAllowHydrogenNbrs, NitrogenLonePairPlanarityTolerance:
142 These parameters are used during the matching of torsion rules containing
143 'N_lp' in their SMARTS patterns. The 'allowHydrogensNbrs' allows the use
144 hydrogen neighbors attached to nitrogen during the determination of its
145 planarity. The 'planarityTolerance' in degrees represents the tolerance
146 allowed for nitrogen to be considered coplanar with its three neighbors.
147
148 The torsion rules containing 'N_lp' in their SMARTS patterns are categorized
149 into the following two types of rules:
150
151 TypeOne:
152
153 [CX4:1][CX4H2:2]!@[NX3;"N_lp":3][CX4:4]
154 [C:1][CX4H2:2]!@[NX3;"N_lp":3][C:4]
155 ... ... ...
156
157 TypeTwo:
158
159 [!#1:1][CX4:2]!@[NX3;"N_lp":3]
160 [C:1][$(S(=O)=O):2]!@["N_lp":3]
161 ... ... ...
162
163 The torsions are matched to torsion rules containing 'N_lp' using specified
164 SMARTS patterns without the 'N_lp' along with additional constraints using
165 the following methodology:
166
167 TypeOne:
168
169 . SMARTS pattern must contain four mapped atoms and the third
170 mapped atom must be a nitrogen matched with 'NX3:3'
171 . Nitrogen atom must have 3 neighbors. The 'allowHydrogens'
172 parameter controls inclusion of hydrogens as its neighbors.
173 . Nitrogen atom and its 3 neighbors must be coplanar.
174 'planarityTolerance' parameter provides tolerance in degrees
175 for nitrogen to be considered coplanar with its 3 neighbors.
176
177 TypeTwo:
178
179 . SMARTS pattern must contain three mapped atoms and the third
180 mapped atom must be a nitrogen matched with 'NX3:3'. The
181 third mapped atom may contain only 'N_lp:3' The missing 'NX3'
182 is automatically detected.
183 . Nitrogen atom must have 3 neighbors. 'allowHydrogens'
184 parameter controls inclusion of hydrogens as neighbors.
185 . Nitrogen atom and its 3 neighbors must not be coplanar.
186 'planarityTolerance' parameter provides tolerance in degrees
187 for nitrogen to be considered coplanar with its 3 neighbors.
188 . Nitrogen lone pair position equivalent to VSEPR theory is
189 determined based on the position of nitrogen and its neighbors.
190 A vector normal to 3 nitrogen neighbors is calculated and added
191 to the coordinates of nitrogen atom to determine the approximate
192 position of the lone pair. It is used as the fourth position to
193 calculate the torsion angle.
194
195 RotBondsSMARTSMode: SMARTS pattern to use for identifying rotatable bonds in
196 a molecule for matching against torsion rules in the torsion library. Possible
197 values: NonStrict, SemiStrict, Strict or Specify. The rotatable bond SMARTS
198 matches are filtered to ensure that each atom in the rotatable bond is attached
199 to at least two heavy atoms.
200
201 The following SMARTS patterns are used to identify rotatable bonds for
202 different modes:
203
204 NonStrict: [!$(*#*)&!D1]-&!@[!$(*#*)&!D1]
205
206 SemiStrict:
207 [!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)
208 &!$(C([CH3])([CH3])[CH3])]-!@[!$(*#*)&!D1&!$(C(F)(F)F)
209 &!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]
210
211 Strict:
212 [!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)
213 &!$(C([CH3])([CH3])[CH3])&!$([CD3](=[N,O,S])-!@[#7,O,S!D1])
214 &!$([#7,O,S!D1]-!@[CD3]=[N,O,S])&!$([CD3](=[N+])-!@[#7!D1])
215 &!$([#7!D1]-!@[CD3]=[N+])]-!@[!$(*#*)&!D1&!$(C(F)(F)F)
216 &!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]
217
218 The 'NonStrict' and 'Strict' SMARTS patterns are available in RDKit. The
219 'NonStrict' SMARTS pattern corresponds to original Daylight SMARTS
220 specification for rotatable bonds. The 'SemiStrict' SMARTS pattern is
221 derived from 'Strict' SMARTS pattern.
222
223 TorsionLibraryFilePath: Specify a XML file name containing data for
224 torsion library hierarchy or use default file, TorsionLibrary.xml, available in
225 the directory containing this file.
226
227 """
228
229 self._ProcessTorsionLibraryAlertsParameters(
230 AlertsMode,
231 MinAlertsCount,
232 NitrogenLonePairAllowHydrogenNbrs,
233 NitrogenLonePairPlanarityTolerance,
234 RotBondsSMARTSMode,
235 RotBondsSMARTSPattern,
236 TorsionLibraryFilePath,
237 )
238
239 self._InitializeTorsionLibraryAlerts()
240
241 def _ProcessTorsionLibraryAlertsParameters(
242 self,
243 AlertsMode,
244 MinAlertsCount,
245 NitrogenLonePairAllowHydrogenNbrs,
246 NitrogenLonePairPlanarityTolerance,
247 RotBondsSMARTSMode,
248 RotBondsSMARTSPattern,
249 TorsionLibraryFilePath,
250 ):
251 """Process torsion library alerts paramaters."""
252
253 # Process AltertsMode paramater...
254 self._ProcessAlertsModeParameter(AlertsMode)
255
256 # Process MinAlertsCount...
257 if MinAlertsCount < 1:
258 raise ValueError(
259 "The value, %s, specified for AltertsMinCount parameter is not valid. Supported value: >=1"
260 % MinAlertsCount
261 )
262 self.MinAlertsCount = MinAlertsCount
263
264 # Process nitrogen lone pair paramaters...
265 self.NitrogenLonePairAllowHydrogenNbrs = NitrogenLonePairAllowHydrogenNbrs
266 if NitrogenLonePairPlanarityTolerance < 0:
267 raise ValueError(
268 "The value, %s, specified for NitrogenLonePairPlanarityTolerance parameter is not valid. Supported value: >= 0"
269 % NitrogenLonePairPlanarityTolerance
270 )
271 self.NitrogenLonePairPlanarityTolerance = NitrogenLonePairPlanarityTolerance
272
273 # Process RotBondsSMARTSMode and RotBondsSMARTSPattern parameters...
274 self._ProcessRotBondsParameters(RotBondsSMARTSMode, RotBondsSMARTSPattern)
275
276 # Process TorsionLibraryFilePath parameter...
277 self._ProcessTorsionLibraryFilePathParameter(TorsionLibraryFilePath)
278
279 def _InitializeTorsionLibraryAlerts(self):
280 """Initialize tosion library alerts."""
281
282 # Retrieve and setup torsion library info for matching rotatable bonds...
283 TorsionLibraryInfo = {}
284
285 TorsionLibElementTree = TorsionAlertsUtil.RetrieveTorsionLibraryInfo(self.TorsionLibraryFilePath)
286 TorsionLibraryInfo["TorsionLibElementTree"] = TorsionLibElementTree
287
288 TorsionAlertsUtil.SetupTorsionLibraryInfoForMatchingRotatableBonds(TorsionLibraryInfo)
289
290 self.TorsionLibraryInfo = TorsionLibraryInfo
291
292 def _ProcessAlertsModeParameter(self, AlertsMode):
293 """Process AlertsMode parameter."""
294
295 SpecifiedAlertsModeList = []
296 if re.match("^Red$", AlertsMode, re.I):
297 SpecifiedAlertsModeList.append("Red")
298 elif re.match("^RedAndOrange$", AlertsMode, re.I):
299 SpecifiedAlertsModeList.append("Red")
300 SpecifiedAlertsModeList.append("Orange")
301 else:
302 raise ValueError(
303 "Invalid value, %s, specified for AlertsMode parameter. Valid values: Red, or RedAndOrange"
304 % (AlertsMode)
305 )
306
307 self.AltersMode = AlertsMode
308 self.SpecifiedAlertsModeList = SpecifiedAlertsModeList
309
310 def _ProcessRotBondsParameters(self, RotBondsSMARTSMode, RotBondsSMARTSPattern):
311 """Process RotBondsSMARTSMode and RotBondsSMARTSPattern parameters."""
312
313 if re.match("^NonStrict$", RotBondsSMARTSMode, re.I):
314 RotBondsSMARTSPattern = "[!$(*#*)&!D1]-&!@[!$(*#*)&!D1]"
315 elif re.match("^SemiStrict$", RotBondsSMARTSMode, re.I):
316 RotBondsSMARTSPattern = "[!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]-!@[!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]"
317 elif re.match("^Strict$", RotBondsSMARTSMode, re.I):
318 RotBondsSMARTSPattern = "[!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])&!$([CD3](=[N,O,S])-!@[#7,O,S!D1])&!$([#7,O,S!D1]-!@[CD3]=[N,O,S])&!$([CD3](=[N+])-!@[#7!D1])&!$([#7!D1]-!@[CD3]=[N+])]-!@[!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]"
319 elif re.match("^Specify$", RotBondsSMARTSMode, re.I):
320 if RotBondsSMARTSPattern is None:
321 raise ValueError(
322 "The RotBondsSMARTSPattern parameter value must be specified during Specify value for RotBondsSMARTSMode parameter."
323 )
324 RotBondsSMARTSPattern = RotBondsSMARTSPattern.strip()
325 if not len(RotBondsSMARTSPattern):
326 raise ValueError("Empty value specified for RotBondsSMARTSPattern parameter.")
327 else:
328 raise ValueError(
329 "Invalid value, %s, specified for RotBondsSMARTSMode parameter. Valid values: NonStrict, SemiStrict, Strict or Specify"
330 % (RotBondsSMARTSMode)
331 )
332
333 RotBondsPatternMol = Chem.MolFromSmarts(RotBondsSMARTSPattern)
334 if RotBondsPatternMol is None:
335 if re.match("Specify", RotBondsSMARTSMode, re.I):
336 raise ValueError(
337 'Failed to create rotatable bonds pattern molecule. The rotatable bonds SMARTS pattern, "%s", specified using RotBondsSMARTSPattern parameter is not valid.'
338 % (RotBondsSMARTSPattern)
339 )
340 else:
341 raise ValueError(
342 'Failed to create rotatable bonds pattern molecule. The default rotatable bonds SMARTS pattern, "%s", used for %s value of RotBondsSMARTSMode parameter is not valid.'
343 % (RotBondsSMARTSPattern, RotBondsSMARTSMode)
344 )
345
346 self.RotBondsSMARTSMode = RotBondsSMARTSMode
347 self.RotBondsSMARTSPattern = RotBondsSMARTSPattern
348 self.RotBondsPatternMol = RotBondsPatternMol
349
350 def _ProcessTorsionLibraryFilePathParameter(self, TorsionLibraryFilePath):
351 """Process torsion library path parameter."""
352
353 # TorsionLibraryFilePath parameter...
354 if re.match("^auto$", TorsionLibraryFilePath):
355 TorsionLibraryFile = "TorsionLibrary.xml"
356 TorsionLibraryFilePath = os.path.join(os.path.dirname(os.path.abspath(__file__)), TorsionLibraryFile)
357 if not os.path.isfile(TorsionLibraryFilePath):
358 raise ValueError("The default torsion alerts library file %s doesn't exist." % TorsionLibraryFilePath)
359 else:
360 if not os.path.isfile(TorsionLibraryFilePath):
361 raise ValueError(
362 "The file specified, %s, for parameter TorsionLibraryFilePath doesn't exist."
363 % TorsionLibraryFilePath
364 )
365 TorsionLibraryFilePath = os.path.abspath(TorsionLibraryFilePath)
366
367 self.TorsionLibraryFilePath = TorsionLibraryFilePath
368
369 def GetTorsionLibraryFilePath(self):
370 """Get torsion strain library file path.
371
372 Arguments:
373 Nothing.
374
375 Returns:
376 FilePath (str): Torsion strain library path.
377
378 """
379
380 return self.TorsionLibraryFilePath
381
382 def ListTorsionLibraryInfo(self):
383 """List torsion strain library information.
384
385 Arguments:
386 Nothing.
387
388 Returns:
389 Nothing. The torsion library information is printed.
390
391 """
392
393 return TorsionAlertsUtil.ListTorsionLibraryInfo(self.TorsionLibraryInfo["TorsionLibElementTree"])
394
395 def IdentifyTorsionLibraryAlertsForRotatableBonds(self, Mol):
396 """Identify torsion library alerts for a molecule by matching rotatable bonds
397 against SMARTS patterns specified for torsion rules in torsion energy library
398 file.
399
400 Arguments:
401 Mol (object): RDKit molecule object.
402
403 Returns:
404 bool: True - Molecule contains strained torsions; False - Molecule
405 contains no strained torsions or rotatable bonds.
406 dict or None: Torsion alerts information regarding matching of
407 rotatable bonds to torsion strain library.
408
409 Examples:
410
411 from TorsionAlerts.TorsionLibraryAlerts import TorsionLibraryAlerts
412
413 LibraryAlerts = TorsionLibraryAlerts()
414 AlertsStatus, AlertsInfo = LibraryAlerts.
415 IdentifyTorsionLibraryAlertsForRotatableBonds(RDKitMol)
416
417 # List of rotatable bond IDs...
418 RotatableBondIDs = AlertsInfo["IDs"]
419
420 # Dictionaries containing information for rotatable bonds by using
421 # bond ID as key...
422 for ID in AlertsInfo["IDs"]:
423 MatchStatus = AlertsInfo["MatchStatus"][ID]
424 AlertTypes = AlertsInfo["AlertTypes"][ID]
425 AtomIndices = AlertsInfo["AtomIndices"][ID]
426 TorsionAtomIndices = AlertsInfo["TorsionAtomIndices"][ID]
427 TorsionAngles = AlertsInfo["TorsionAngles"][ID]
428 TorsionAngleViolations = AlertsInfo["TorsionAngleViolations"][ID]
429 HierarchyClassNames = AlertsInfo["HierarchyClassNames"][ID]
430 HierarchySubClassNames = AlertsInfo["HierarchySubClassNames"][ID]
431 TorsionRuleNodeID = AlertsInfo["TorsionRuleNodeID"][ID]
432 TorsionRulePeaks = AlertsInfo["TorsionRulePeaks"][ID]
433 TorsionRuleTolerances1 = AlertsInfo["TorsionRuleTolerances1"][ID]
434 TorsionRuleTolerances2 = AlertsInfo["TorsionRuleTolerances2"][ID]
435 TorsionRuleSMARTS = AlertsInfo["TorsionRuleSMARTS"][ID]
436
437 """
438
439 # Identify rotatable bonds...
440 RotBondsStatus, RotBondsInfo = TorsionAlertsUtil.IdentifyRotatableBondsForTorsionLibraryMatch(
441 self.TorsionLibraryInfo, Mol, self.RotBondsPatternMol
442 )
443 if not RotBondsStatus:
444 return (False, None)
445
446 # Identify alerts for rotatable bonds...
447 RotBondsAlertsStatus, RotBondsAlertsInfo = self._MatchRotatableBondsToTorsionLibrary(Mol, RotBondsInfo)
448
449 return (RotBondsAlertsStatus, RotBondsAlertsInfo)
450
451 def _MatchRotatableBondsToTorsionLibrary(self, Mol, RotBondsInfo):
452 """Match rotatable bond to torsion library."""
453
454 # Initialize...
455 RotBondsAlertsInfo = self._InitializeRotatableBondsAlertsInfo()
456
457 # Match rotatable bonds to torsion library...
458 for ID in RotBondsInfo["IDs"]:
459 AtomIndices = RotBondsInfo["AtomIndices"][ID]
460 HierarchyClass = RotBondsInfo["HierarchyClass"][ID]
461
462 MatchStatus, MatchInfo = self._MatchRotatableBondToTorsionLibrary(Mol, AtomIndices, HierarchyClass)
463
464 if MatchInfo is None or len(MatchInfo) == 0:
465 (
466 AlertType,
467 TorsionAtomIndices,
468 TorsionAngle,
469 TorsionAngleViolation,
470 HierarchyClassName,
471 HierarchySubClassName,
472 TorsionRuleNodeID,
473 TorsionRulePeaks,
474 TorsionRuleTolerances1,
475 TorsionRuleTolerances2,
476 TorsionRuleSMARTS,
477 ) = [None] * 11
478 else:
479 (
480 AlertType,
481 TorsionAtomIndices,
482 TorsionAngle,
483 TorsionAngleViolation,
484 HierarchyClassName,
485 HierarchySubClassName,
486 TorsionRuleNodeID,
487 TorsionRulePeaks,
488 TorsionRuleTolerances1,
489 TorsionRuleTolerances2,
490 TorsionRuleSMARTS,
491 ) = MatchInfo
492
493 # Track alerts information...
494 RotBondsAlertsInfo["IDs"].append(ID)
495 RotBondsAlertsInfo["MatchStatus"][ID] = MatchStatus
496 RotBondsAlertsInfo["AlertTypes"][ID] = AlertType
497 RotBondsAlertsInfo["AtomIndices"][ID] = AtomIndices
498 RotBondsAlertsInfo["TorsionAtomIndices"][ID] = TorsionAtomIndices
499 RotBondsAlertsInfo["TorsionAngles"][ID] = TorsionAngle
500 RotBondsAlertsInfo["TorsionAngleViolations"][ID] = TorsionAngleViolation
501 RotBondsAlertsInfo["HierarchyClassNames"][ID] = HierarchyClassName
502 RotBondsAlertsInfo["HierarchySubClassNames"][ID] = HierarchySubClassName
503 RotBondsAlertsInfo["TorsionRuleNodeID"][ID] = TorsionRuleNodeID
504 RotBondsAlertsInfo["TorsionRulePeaks"][ID] = TorsionRulePeaks
505 RotBondsAlertsInfo["TorsionRuleTolerances1"][ID] = TorsionRuleTolerances1
506 RotBondsAlertsInfo["TorsionRuleTolerances2"][ID] = TorsionRuleTolerances2
507 RotBondsAlertsInfo["TorsionRuleSMARTS"][ID] = TorsionRuleSMARTS
508
509 # Count alert types...
510 if AlertType is not None:
511 if AlertType not in RotBondsAlertsInfo["Count"]:
512 RotBondsAlertsInfo["Count"][AlertType] = 0
513 RotBondsAlertsInfo["Count"][AlertType] += 1
514
515 # Setup alert status for rotatable bonds...
516 RotBondsAlertsStatus = False
517 AlertsCount = 0
518 for ID in RotBondsInfo["IDs"]:
519 if RotBondsAlertsInfo["AlertTypes"][ID] in self.SpecifiedAlertsModeList:
520 AlertsCount += 1
521 if AlertsCount >= self.MinAlertsCount:
522 RotBondsAlertsStatus = True
523 break
524
525 return (RotBondsAlertsStatus, RotBondsAlertsInfo)
526
527 def _InitializeRotatableBondsAlertsInfo(
528 self,
529 ):
530 """Initialize alerts information for rotatable bonds."""
531
532 RotBondsAlertsInfo = {}
533 RotBondsAlertsInfo["IDs"] = []
534
535 for DataLabel in [
536 "MatchStatus",
537 "AlertTypes",
538 "AtomIndices",
539 "TorsionAtomIndices",
540 "TorsionAngles",
541 "TorsionAngleViolations",
542 "HierarchyClassNames",
543 "HierarchySubClassNames",
544 "TorsionRuleNodeID",
545 "TorsionRulePeaks",
546 "TorsionRuleTolerances1",
547 "TorsionRuleTolerances2",
548 "TorsionRuleSMARTS",
549 "Count",
550 ]:
551 RotBondsAlertsInfo[DataLabel] = {}
552
553 return RotBondsAlertsInfo
554
555 def _MatchRotatableBondToTorsionLibrary(self, Mol, RotBondAtomIndices, RotBondHierarchyClass):
556 """Match rotatable bond to torsion library."""
557
558 if TorsionAlertsUtil.IsSpecificHierarchyClass(self.TorsionLibraryInfo, RotBondHierarchyClass):
559 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstSpecificHierarchyClass(
560 Mol, RotBondAtomIndices, RotBondHierarchyClass
561 )
562 if not MatchStatus:
563 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstGenericHierarchyClass(
564 Mol, RotBondAtomIndices, RotBondHierarchyClass
565 )
566 else:
567 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstGenericHierarchyClass(
568 Mol, RotBondAtomIndices, RotBondHierarchyClass
569 )
570
571 return (MatchStatus, MatchInfo)
572
573 def _MatchRotatableBondAgainstSpecificHierarchyClass(self, Mol, RotBondAtomIndices, RotBondHierarchyClass):
574 """Match rotatable bond against a specific hierarchy class."""
575
576 TorsionLibraryInfo = self.TorsionLibraryInfo
577
578 HierarchyClassElementNode = None
579 if RotBondHierarchyClass in TorsionLibraryInfo["SpecificClasses"]["ElementNode"]:
580 HierarchyClassElementNode = TorsionLibraryInfo["SpecificClasses"]["ElementNode"][RotBondHierarchyClass]
581
582 if HierarchyClassElementNode is None:
583 return (False, None, None, None)
584
585 TorsionAlertsUtil.TrackHierarchyClassElementNode(TorsionLibraryInfo, HierarchyClassElementNode)
586 MatchStatus, MatchInfo = self._ProcessElementForRotatableBondMatch(
587 Mol, RotBondAtomIndices, HierarchyClassElementNode
588 )
589 TorsionAlertsUtil.RemoveLastHierarchyClassElementNodeFromTracking(TorsionLibraryInfo)
590
591 return (MatchStatus, MatchInfo)
592
593 def _MatchRotatableBondAgainstGenericHierarchyClass(self, Mol, RotBondAtomIndices, RotBondHierarchyClass):
594 """Match rotatable bond against a generic hierarchy class."""
595
596 TorsionLibraryInfo = self.TorsionLibraryInfo
597
598 HierarchyClassElementNode = TorsionAlertsUtil.GetGenericHierarchyClassElementNode(TorsionLibraryInfo)
599 if HierarchyClassElementNode is None:
600 return (False, None)
601
602 TorsionAlertsUtil.TrackHierarchyClassElementNode(TorsionLibraryInfo, HierarchyClassElementNode)
603
604 # Match hierarchy subclasses before matching torsion rules...
605 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstGenericHierarchySubClasses(
606 Mol, RotBondAtomIndices, HierarchyClassElementNode
607 )
608
609 if not MatchStatus:
610 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstGenericHierarchyTorsionRules(
611 Mol, RotBondAtomIndices, HierarchyClassElementNode
612 )
613
614 TorsionAlertsUtil.RemoveLastHierarchyClassElementNodeFromTracking(TorsionLibraryInfo)
615
616 return (MatchStatus, MatchInfo)
617
618 def _MatchRotatableBondAgainstGenericHierarchySubClasses(self, Mol, RotBondAtomIndices, HierarchyClassElementNode):
619 """Match rotatable bond againat generic hierarchy subclasses."""
620
621 for ElementChildNode in HierarchyClassElementNode:
622 if ElementChildNode.tag != "hierarchySubClass":
623 continue
624
625 SubClassMatchStatus = self._ProcessHierarchySubClassElementForRotatableBondMatch(
626 Mol, RotBondAtomIndices, ElementChildNode
627 )
628
629 if SubClassMatchStatus:
630 MatchStatus, MatchInfo = self._ProcessElementForRotatableBondMatch(
631 Mol, RotBondAtomIndices, ElementChildNode
632 )
633
634 if MatchStatus:
635 return (MatchStatus, MatchInfo)
636
637 return (False, None)
638
639 def _MatchRotatableBondAgainstGenericHierarchyTorsionRules(
640 self, Mol, RotBondAtomIndices, HierarchyClassElementNode
641 ):
642 """Match rotatable bond againat torsion rules generic hierarchy class."""
643
644 for ElementChildNode in HierarchyClassElementNode:
645 if ElementChildNode.tag != "torsionRule":
646 continue
647
648 MatchStatus, MatchInfo = self._ProcessTorsionRuleElementForRotatableBondMatch(
649 Mol, RotBondAtomIndices, ElementChildNode
650 )
651
652 if MatchStatus:
653 return (MatchStatus, MatchInfo)
654
655 return (False, None)
656
657 def _ProcessElementForRotatableBondMatch(self, Mol, RotBondAtomIndices, ElementNode):
658 """Process element node to recursively match rotatable bond against hierarchy
659 subclasses and torsion rules."""
660
661 TorsionLibraryInfo = self.TorsionLibraryInfo
662
663 for ElementChildNode in ElementNode:
664 if ElementChildNode.tag == "hierarchySubClass":
665 SubClassMatchStatus = self._ProcessHierarchySubClassElementForRotatableBondMatch(
666 Mol, RotBondAtomIndices, ElementChildNode
667 )
668
669 if SubClassMatchStatus:
670 TorsionAlertsUtil.TrackHierarchySubClassElementNode(TorsionLibraryInfo, ElementChildNode)
671
672 MatchStatus, MatchInfo = self._ProcessElementForRotatableBondMatch(
673 Mol, RotBondAtomIndices, ElementChildNode
674 )
675 if MatchStatus:
676 TorsionAlertsUtil.RemoveLastHierarchySubClassElementNodeFromTracking(TorsionLibraryInfo)
677 return (MatchStatus, MatchInfo)
678
679 TorsionAlertsUtil.RemoveLastHierarchySubClassElementNodeFromTracking(TorsionLibraryInfo)
680
681 elif ElementChildNode.tag == "torsionRule":
682 MatchStatus, MatchInfo = self._ProcessTorsionRuleElementForRotatableBondMatch(
683 Mol, RotBondAtomIndices, ElementChildNode
684 )
685
686 if MatchStatus:
687 return (MatchStatus, MatchInfo)
688
689 return (False, None)
690
691 def _ProcessHierarchySubClassElementForRotatableBondMatch(self, Mol, RotBondAtomIndices, ElementNode):
692 """Process hierarchy subclass element to match rotatable bond."""
693
694 # Setup subclass SMARTS pattern mol...
695 SubClassPatternMol = TorsionAlertsUtil.SetupHierarchySubClassElementPatternMol(
696 self.TorsionLibraryInfo, ElementNode
697 )
698 if SubClassPatternMol is None:
699 return False
700
701 # Match SMARTS pattern...
702 SubClassPatternMatches = TorsionAlertsUtil.FilterSubstructureMatchesByAtomMapNumbers(
703 Mol, SubClassPatternMol, Mol.GetSubstructMatches(SubClassPatternMol, useChirality=False)
704 )
705 if len(SubClassPatternMatches) == 0:
706 return False
707
708 # Match rotatable bond indices...
709 RotBondAtomIndex1, RotBondAtomIndex2 = RotBondAtomIndices
710 MatchStatus = False
711 for SubClassPatternMatch in SubClassPatternMatches:
712 if len(SubClassPatternMatch) == 2:
713 # Matched to pattern containing map atom numbers ":2" and ":3"...
714 CentralAtomsIndex1, CentralAtomsIndex2 = SubClassPatternMatch
715 elif len(SubClassPatternMatch) == 4:
716 # Matched to pattern containing map atom numbers ":1", ":2", ":3" and ":4"...
717 CentralAtomsIndex1 = SubClassPatternMatch[1]
718 CentralAtomsIndex2 = SubClassPatternMatch[2]
719 elif len(SubClassPatternMatch) == 3:
720 SubClassSMARTSPattern = ElementNode.get("smarts")
721 if TorsionAlertsUtil.DoesSMARTSContainsMappedAtoms(SubClassSMARTSPattern, [":2", ":3", ":4"]):
722 # Matched to pattern containing map atom numbers ":2", ":3" and ":4"...
723 CentralAtomsIndex1 = SubClassPatternMatch[0]
724 CentralAtomsIndex2 = SubClassPatternMatch[1]
725 else:
726 # Matched to pattern containing map atom numbers ":1", ":2" and ":3"...
727 CentralAtomsIndex1 = SubClassPatternMatch[1]
728 CentralAtomsIndex2 = SubClassPatternMatch[2]
729 else:
730 continue
731
732 if CentralAtomsIndex1 != CentralAtomsIndex2:
733 if (CentralAtomsIndex1 == RotBondAtomIndex1 and CentralAtomsIndex2 == RotBondAtomIndex2) or (
734 CentralAtomsIndex1 == RotBondAtomIndex2 and CentralAtomsIndex2 == RotBondAtomIndex1
735 ):
736 MatchStatus = True
737 break
738
739 return MatchStatus
740
741 def _ProcessTorsionRuleElementForRotatableBondMatch(self, Mol, RotBondAtomIndices, ElementNode):
742 """Process torsion rule element to match rotatable bond."""
743
744 TorsionLibraryInfo = self.TorsionLibraryInfo
745
746 # Retrieve torsion matched to rotatable bond...
747 TorsionAtomIndices, TorsionAngle = self._MatchTorsionRuleToRotatableBond(Mol, RotBondAtomIndices, ElementNode)
748 if TorsionAtomIndices is None:
749 return (False, None)
750
751 # Setup torsion angles info for matched torsion rule...
752 TorsionAnglesInfo = TorsionAlertsUtil.SetupTorsionRuleAnglesInfo(TorsionLibraryInfo, ElementNode)
753 if TorsionAnglesInfo is None:
754 return (False, None)
755
756 # Setup torsion alert type and angle violation...
757 AlertType, TorsionAngleViolation = self._SetupTorsionAlertTypeForRotatableBond(TorsionAnglesInfo, TorsionAngle)
758
759 # Setup hierarchy class and subclass names...
760 HierarchyClassName, HierarchySubClassName = (
761 TorsionAlertsUtil.SetupHierarchyClassAndSubClassNamesForRotatableBond(TorsionLibraryInfo)
762 )
763
764 # Setup rule node ID...
765 TorsionRuleNodeID = ElementNode.get("NodeID")
766
767 # Setup SMARTS...
768 TorsionRuleSMARTS = ElementNode.get("smarts")
769 if " " in TorsionRuleSMARTS:
770 TorsionRuleSMARTS = TorsionRuleSMARTS.replace(" ", "")
771
772 # Setup torsion peaks and tolerances...
773 TorsionRulePeaks = TorsionAnglesInfo["ValuesList"]
774 TorsionRuleTolerances1 = TorsionAnglesInfo["Tolerances1List"]
775 TorsionRuleTolerances2 = TorsionAnglesInfo["Tolerances2List"]
776
777 MatchInfo = [
778 AlertType,
779 TorsionAtomIndices,
780 TorsionAngle,
781 TorsionAngleViolation,
782 HierarchyClassName,
783 HierarchySubClassName,
784 TorsionRuleNodeID,
785 TorsionRulePeaks,
786 TorsionRuleTolerances1,
787 TorsionRuleTolerances2,
788 TorsionRuleSMARTS,
789 ]
790
791 # Setup match status...
792 MatchStatus = True
793
794 return (MatchStatus, MatchInfo)
795
796 def _MatchTorsionRuleToRotatableBond(self, Mol, RotBondAtomIndices, ElementNode):
797 """Retrieve matched torsion for torsion rule matched to rotatable bond."""
798
799 # Get torsion matches...
800 TorsionMatches = self._GetMatchesForTorsionRule(Mol, ElementNode)
801 if TorsionMatches is None or len(TorsionMatches) == 0:
802 return (None, None)
803
804 # Identify the first torsion match corresponding to central atoms in RotBondAtomIndices...
805 RotBondAtomIndex1, RotBondAtomIndex2 = RotBondAtomIndices
806 for TorsionMatch in TorsionMatches:
807 CentralAtomIndex1 = TorsionMatch[1]
808 CentralAtomIndex2 = TorsionMatch[2]
809
810 if (CentralAtomIndex1 == RotBondAtomIndex1 and CentralAtomIndex2 == RotBondAtomIndex2) or (
811 CentralAtomIndex1 == RotBondAtomIndex2 and CentralAtomIndex2 == RotBondAtomIndex1
812 ):
813 TorsionAngle = self._CalculateTorsionAngle(Mol, TorsionMatch)
814
815 return (TorsionMatch, TorsionAngle)
816
817 return (None, None)
818
819 def _CalculateTorsionAngle(self, Mol, TorsionMatch):
820 """Calculate torsion angle."""
821
822 if type(TorsionMatch[3]) is list:
823 return self._CalculateTorsionAngleUsingNitrogenLonePairPosition(Mol, TorsionMatch)
824
825 # Calculate torsion angle using torsion atom indices..
826 MolConf = Mol.GetConformer(0)
827 TorsionAngle = rdMolTransforms.GetDihedralDeg(
828 MolConf, TorsionMatch[0], TorsionMatch[1], TorsionMatch[2], TorsionMatch[3]
829 )
830 TorsionAngle = round(TorsionAngle, 2)
831
832 return TorsionAngle
833
834 def _CalculateTorsionAngleUsingNitrogenLonePairPosition(self, Mol, TorsionMatch):
835 """Calculate torsion angle using nitrogen lone pair positon."""
836
837 # Setup a carbon atom as position holder for lone pair position...
838 TmpMol = Chem.RWMol(Mol)
839 LonePairAtomIndex = TmpMol.AddAtom(Chem.Atom(6))
840
841 TmpMolConf = TmpMol.GetConformer(0)
842 TmpMolConf.SetAtomPosition(LonePairAtomIndex, TorsionMatch[3])
843
844 TorsionAngle = rdMolTransforms.GetDihedralDeg(
845 TmpMolConf, TorsionMatch[0], TorsionMatch[1], TorsionMatch[2], LonePairAtomIndex
846 )
847 TorsionAngle = round(TorsionAngle, 2)
848
849 return TorsionAngle
850
851 def _GetMatchesForTorsionRule(self, Mol, ElementNode):
852 """Get matches for torsion rule."""
853
854 # Match torsions...
855 TorsionMatches = None
856 if self._IsNitogenLonePairTorsionRule(ElementNode):
857 TorsionMatches = self._GetSubstructureMatchesForNitrogenLonePairTorsionRule(Mol, ElementNode)
858 else:
859 TorsionMatches = self._GetSubstructureMatchesForTorsionRule(Mol, ElementNode)
860
861 if TorsionMatches is None or len(TorsionMatches) == 0:
862 return TorsionMatches
863
864 # Filter torsion matches...
865 FiltertedTorsionMatches = []
866 for TorsionMatch in TorsionMatches:
867 if len(TorsionMatch) != 4:
868 continue
869
870 # Ignore matches containing hydrogen atoms as first or last atom...
871 if Mol.GetAtomWithIdx(TorsionMatch[0]).GetAtomicNum() == 1:
872 continue
873 if type(TorsionMatch[3]) is int:
874 # May contains a list for type two nitrogen lone pair match...
875 if Mol.GetAtomWithIdx(TorsionMatch[3]).GetAtomicNum() == 1:
876 continue
877 FiltertedTorsionMatches.append(TorsionMatch)
878
879 return FiltertedTorsionMatches
880
881 def _GetSubstructureMatchesForTorsionRule(self, Mol, ElementNode):
882 """Get substructure matches for a torsion rule."""
883
884 # Setup torsion rule SMARTS pattern mol....
885 TorsionRuleNodeID = ElementNode.get("NodeID")
886 TorsionSMARTSPattern = ElementNode.get("smarts")
887 TorsionPatternMol = TorsionAlertsUtil.SetupTorsionRuleElementPatternMol(
888 self.TorsionLibraryInfo, ElementNode, TorsionRuleNodeID, TorsionSMARTSPattern
889 )
890 if TorsionPatternMol is None:
891 return None
892
893 # Match torsions...
894 TorsionMatches = TorsionAlertsUtil.FilterSubstructureMatchesByAtomMapNumbers(
895 Mol, TorsionPatternMol, Mol.GetSubstructMatches(TorsionPatternMol, useChirality=False)
896 )
897
898 return TorsionMatches
899
900 def _GetSubstructureMatchesForNitrogenLonePairTorsionRule(self, Mol, ElementNode):
901 """Get substructure matches for a torsion rule containing N_lp."""
902
903 if self._IsTypeOneNitogenLonePairTorsionRule(ElementNode):
904 return self._GetSubstructureMatchesForTypeOneNitrogenLonePairTorsionRule(Mol, ElementNode)
905 elif self._IsTypeTwoNitogenLonePairTorsionRule(ElementNode):
906 return self._GetSubstructureMatchesForTypeTwoNitrogenLonePairTorsionRule(Mol, ElementNode)
907
908 return None
909
910 def _GetSubstructureMatchesForTypeOneNitrogenLonePairTorsionRule(self, Mol, ElementNode):
911 """Get substructure matches for a torsion rule containing N_lp and four mapped atoms."""
912
913 # For example:
914 # [CX4:1][CX4H2:2]!@[NX3;"N_lp":3][CX4:4]
915 # [C:1][CX4H2:2]!@[NX3;"N_lp":3][C:4]
916 # ... ... ...
917
918 TorsionRuleNodeID = ElementNode.get("NodeID")
919 TorsionPatternMol, LonePairMapNumber = self._SetupNitrogenLonePairTorsionRuleElementInfo(
920 ElementNode, TorsionRuleNodeID
921 )
922
923 if TorsionPatternMol is None:
924 return None
925
926 if LonePairMapNumber is None:
927 return None
928
929 # Match torsions...
930 TorsionMatches = TorsionAlertsUtil.FilterSubstructureMatchesByAtomMapNumbers(
931 Mol, TorsionPatternMol, Mol.GetSubstructMatches(TorsionPatternMol, useChirality=False)
932 )
933
934 # Filter matches...
935 FiltertedTorsionMatches = []
936 for TorsionMatch in TorsionMatches:
937 if len(TorsionMatch) != 4:
938 continue
939
940 # Check for Nitogen atom at LonePairMapNumber...
941 LonePairNitrogenAtom = Mol.GetAtomWithIdx(TorsionMatch[LonePairMapNumber - 1])
942 if LonePairNitrogenAtom.GetSymbol() != "N":
943 continue
944
945 # Make sure LonePairNitrogenAtom is planar...
946 # test
947 PlanarityStatus = self._IsLonePairNitrogenAtomPlanar(Mol, LonePairNitrogenAtom)
948 if PlanarityStatus is None:
949 continue
950
951 if not PlanarityStatus:
952 continue
953
954 FiltertedTorsionMatches.append(TorsionMatch)
955
956 return FiltertedTorsionMatches
957
958 def _GetSubstructureMatchesForTypeTwoNitrogenLonePairTorsionRule(self, Mol, ElementNode):
959 """Get substructure matches for a torsion rule containing N_lp and three mapped atoms."""
960
961 # For example:
962 # [!#1:1][CX4:2]!@[NX3;"N_lp":3]
963 # [!#1:1][$(S(=O)=O):2]!@["N_lp":3]@[Cr3]
964 # [C:1][$(S(=O)=O):2]!@["N_lp":3]
965 # [c:1][$(S(=O)=O):2]!@["N_lp":3]
966 # [!#1:1][$(S(=O)=O):2]!@["N_lp":3]
967 # ... ... ...
968
969 TorsionRuleNodeID = ElementNode.get("NodeID")
970 TorsionPatternMol, LonePairMapNumber = self._SetupNitrogenLonePairTorsionRuleElementInfo(
971 ElementNode, TorsionRuleNodeID
972 )
973
974 if TorsionPatternMol is None:
975 return None
976
977 if not self._IsValidTypeTwoNitrogenLonePairMapNumber(LonePairMapNumber):
978 return None
979
980 # Match torsions...
981 TorsionMatches = TorsionAlertsUtil.FilterSubstructureMatchesByAtomMapNumbers(
982 Mol, TorsionPatternMol, Mol.GetSubstructMatches(TorsionPatternMol, useChirality=False)
983 )
984
985 # Filter matches...
986 FiltertedTorsionMatches = []
987 for TorsionMatch in TorsionMatches:
988 if len(TorsionMatch) != 3:
989 continue
990
991 # Check for Nitogen atom at LonePairMapNumber...
992 LonePairNitrogenAtom = Mol.GetAtomWithIdx(TorsionMatch[LonePairMapNumber - 1])
993 if LonePairNitrogenAtom.GetSymbol() != "N":
994 continue
995
996 # Make sure LonePairNitrogenAtom is not planar...
997 PlanarityStatus = self._IsLonePairNitrogenAtomPlanar(Mol, LonePairNitrogenAtom)
998
999 if PlanarityStatus is None:
1000 continue
1001
1002 if PlanarityStatus:
1003 continue
1004
1005 # Calculate lone pair coordinates for a non-planar nitrogen...
1006 LonePairPosition = self._CalculateLonePairCoordinatesForNitrogenAtom(Mol, LonePairNitrogenAtom)
1007 if LonePairPosition is None:
1008 continue
1009
1010 # Append lone pair coodinate list to list of torsion match containing atom indices...
1011 TorsionMatch.append(LonePairPosition)
1012
1013 # Track torsion matches...
1014 FiltertedTorsionMatches.append(TorsionMatch)
1015
1016 return FiltertedTorsionMatches
1017
1018 def _SetupNitrogenLonePairTorsionRuleElementInfo(self, ElementNode, TorsionRuleNodeID):
1019 """Setup pattern molecule and lone pair map number for type one and type
1020 two nitrogen lone pair rules."""
1021
1022 TorsionLibraryInfo = self.TorsionLibraryInfo
1023 TorsionPatternMol, LonePairMapNumber = [None] * 2
1024
1025 if TorsionRuleNodeID in TorsionLibraryInfo["DataCache"]["TorsionRulePatternMol"]:
1026 TorsionPatternMol = TorsionLibraryInfo["DataCache"]["TorsionRulePatternMol"][TorsionRuleNodeID]
1027 LonePairMapNumber = TorsionLibraryInfo["DataCache"]["TorsionRuleLonePairMapNumber"][TorsionRuleNodeID]
1028 else:
1029 # Setup torsion pattern...
1030 TorsionSMARTSPattern = ElementNode.get("smarts")
1031 TorsionSMARTSPattern, LonePairMapNumber = self._ProcessSMARTSForNitrogenLonePairTorsionRule(
1032 TorsionSMARTSPattern
1033 )
1034
1035 # Setup torsion pattern mol...
1036 TorsionPatternMol = Chem.MolFromSmarts(TorsionSMARTSPattern)
1037 if TorsionPatternMol is None:
1038 print(
1039 "Warning: Ignoring torsion rule element containing invalid map atoms numbers in SMARTS pattern %s"
1040 % TorsionSMARTSPattern
1041 )
1042
1043 # Cache data...
1044 TorsionLibraryInfo["DataCache"]["TorsionRulePatternMol"][TorsionRuleNodeID] = TorsionPatternMol
1045 TorsionLibraryInfo["DataCache"]["TorsionRuleLonePairMapNumber"][TorsionRuleNodeID] = LonePairMapNumber
1046
1047 return (TorsionPatternMol, LonePairMapNumber)
1048
1049 def _IsLonePairNitrogenAtomPlanar(self, Mol, NitrogenAtom):
1050 """Check for the planarity of nitrogen atom and its three neighbors."""
1051
1052 AllowHydrogenNbrs = self.NitrogenLonePairAllowHydrogenNbrs
1053 Tolerance = self.NitrogenLonePairPlanarityTolerance
1054
1055 # Get neighbors...
1056 if AllowHydrogenNbrs:
1057 AtomNeighbors = NitrogenAtom.GetNeighbors()
1058 else:
1059 AtomNeighbors = TorsionAlertsUtil.GetHeavyAtomNeighbors(NitrogenAtom)
1060
1061 if len(AtomNeighbors) != 3:
1062 return None
1063
1064 # Setup atom positions...
1065 AtomPositions = []
1066 MolAtomsPositions = TorsionAlertsUtil.GetAtomPositions(Mol)
1067
1068 # Neighbor positions...
1069 for AtomNbr in AtomNeighbors:
1070 AtomNbrIndex = AtomNbr.GetIdx()
1071 AtomPositions.append(MolAtomsPositions[AtomNbrIndex])
1072
1073 # Nitrogen position...
1074 NitrogenAtomIndex = NitrogenAtom.GetIdx()
1075 AtomPositions.append(MolAtomsPositions[NitrogenAtomIndex])
1076
1077 Status = self._AreFourPointsCoplanar(
1078 AtomPositions[0], AtomPositions[1], AtomPositions[2], AtomPositions[3], Tolerance
1079 )
1080
1081 return Status
1082
1083 def _AreFourPointsCoplanar(self, Point1, Point2, Point3, Point4, Tolerance=1.0):
1084 """Check whether four points are coplanar with in the threshold of 1 degree."""
1085
1086 # Setup normalized direction vectors...
1087 VectorP2P1 = self._NormalizeVector(np.subtract(Point2, Point1))
1088 VectorP3P1 = self._NormalizeVector(np.subtract(Point3, Point1))
1089 VectorP1P4 = self._NormalizeVector(np.subtract(Point1, Point4))
1090
1091 # Calculate angle between VectorP1P4 and normal to vectors VectorP2P1 and VectorP3P1...
1092 PlaneP1P2P3Normal = self._NormalizeVector(np.cross(VectorP2P1, VectorP3P1))
1093 PlanarityAngle = np.arccos(np.clip(np.dot(PlaneP1P2P3Normal, VectorP1P4), -1.0, 1.0))
1094
1095 Status = math.isclose(PlanarityAngle, math.radians(90), abs_tol=math.radians(Tolerance))
1096
1097 return Status
1098
1099 def _NormalizeVector(self, Vector):
1100 """Normalize vector."""
1101
1102 Norm = np.linalg.norm(Vector)
1103
1104 return Vector if math.isclose(Norm, 0.0, abs_tol=1e-08) else Vector / Norm
1105
1106 def _CalculateLonePairCoordinatesForNitrogenAtom(self, Mol, NitrogenAtom):
1107 """Calculate approximate lone pair coordinates for non-plannar nitrogen atom."""
1108
1109 AllowHydrogenNbrs = self.NitrogenLonePairAllowHydrogenNbrs
1110
1111 # Get neighbors...
1112 if AllowHydrogenNbrs:
1113 AtomNeighbors = NitrogenAtom.GetNeighbors()
1114 else:
1115 AtomNeighbors = TorsionAlertsUtil.GetHeavyAtomNeighbors(NitrogenAtom)
1116
1117 if len(AtomNeighbors) != 3:
1118 return None
1119
1120 # Setup positions for nitrogen and its neghbors...
1121 MolAtomsPositions = TorsionAlertsUtil.GetAtomPositions(Mol)
1122
1123 NitrogenPosition = MolAtomsPositions[NitrogenAtom.GetIdx()]
1124 NbrPositions = []
1125 for AtomNbr in AtomNeighbors:
1126 NbrPositions.append(MolAtomsPositions[AtomNbr.GetIdx()])
1127 Nbr1Position, Nbr2Position, Nbr3Position = NbrPositions
1128
1129 # Setup normalized direction vectors...
1130 VectorP2P1 = self._NormalizeVector(np.subtract(Nbr2Position, Nbr1Position))
1131 VectorP3P1 = self._NormalizeVector(np.subtract(Nbr3Position, Nbr1Position))
1132 VectorP1P4 = self._NormalizeVector(np.subtract(Nbr1Position, NitrogenPosition))
1133
1134 # Calculate angle between VectorP1P4 and normal to vectors VectorP2P1 and VectorP3P1...
1135 PlaneP1P2P3Normal = self._NormalizeVector(np.cross(VectorP2P1, VectorP3P1))
1136 PlanarityAngle = np.arccos(np.clip(np.dot(PlaneP1P2P3Normal, VectorP1P4), -1.0, 1.0))
1137
1138 # Check for reversing the direction of the normal...
1139 if PlanarityAngle < math.radians(90):
1140 PlaneP1P2P3Normal = PlaneP1P2P3Normal * -1
1141
1142 # Add normal to nitrogen cooridnates for the approximate coordinates of the
1143 # one pair. The exact VSEPR coordinates of the lone pair are not necessary to
1144 # calculate the torsion angle...
1145 LonePairPosition = NitrogenPosition + PlaneP1P2P3Normal
1146
1147 return list(LonePairPosition)
1148
1149 def _ProcessSMARTSForNitrogenLonePairTorsionRule(self, SMARTSPattern):
1150 """Process SMARTS pattern for a torion rule containing N_lp."""
1151
1152 LonePairMapNumber = self._GetNitrogenLonePairMapNumber(SMARTSPattern)
1153
1154 # Remove double quotes around N_lp..
1155 SMARTSPattern = re.sub('"N_lp"', "N_lp", SMARTSPattern, flags=re.I)
1156
1157 # Remove N_lp specification from SMARTS pattern for torsion rule...
1158 if re.search(r"\[N_lp", SMARTSPattern, re.I):
1159 # Handle missing NX3...
1160 SMARTSPattern = re.sub(r"\[N_lp", "[NX3", SMARTSPattern)
1161 else:
1162 SMARTSPattern = re.sub(";N_lp", "", SMARTSPattern)
1163
1164 return (SMARTSPattern, LonePairMapNumber)
1165
1166 def _GetNitrogenLonePairMapNumber(self, SMARTSPattern):
1167 """Get atom map number for nitrogen involved in N_lp."""
1168
1169 LonePairMapNumber = None
1170
1171 SMARTSPattern = re.sub('"N_lp"', "N_lp", SMARTSPattern, flags=re.I)
1172 MatchedMappedAtoms = re.findall("N_lp:[0-9]", SMARTSPattern, re.I)
1173
1174 if len(MatchedMappedAtoms) == 1:
1175 LonePairMapNumber = int(re.sub("N_lp:", "", MatchedMappedAtoms[0]))
1176
1177 return LonePairMapNumber
1178
1179 def _IsNitogenLonePairTorsionRule(self, ElementNode):
1180 """Check for the presence of N_lp in SMARTS pattern for a torsion rule."""
1181
1182 if "N_lp" not in ElementNode.get("smarts"):
1183 return False
1184
1185 LonePairMatches = re.findall("N_lp", ElementNode.get("smarts"), re.I)
1186
1187 return True if len(LonePairMatches) == 1 else False
1188
1189 def _IsTypeOneNitogenLonePairTorsionRule(self, ElementNode):
1190 """Check for the presence four mapped atoms in a SMARTS pattern containing
1191 N_lp for a torsion rule."""
1192
1193 # For example:
1194 # [CX4:1][CX4H2:2]!@[NX3;"N_lp":3][CX4:4]
1195 # [C:1][CX4H2:2]!@[NX3;"N_lp":3][C:4]
1196 # ... ... ...
1197
1198 MatchedMappedAtoms = re.findall(":[0-9]", ElementNode.get("smarts"), re.I)
1199
1200 return True if len(MatchedMappedAtoms) == 4 else False
1201
1202 def _IsTypeTwoNitogenLonePairTorsionRule(self, ElementNode):
1203 """Check for the presence three mapped atoms in a SMARTS pattern containing
1204 N_lp for a torsion rule."""
1205
1206 # For example:
1207 # [!#1:1][CX4:2]!@[NX3;"N_lp":3]
1208 # [!#1:1][$(S(=O)=O):2]!@["N_lp":3]@[Cr3]
1209 # [C:1][$(S(=O)=O):2]!@["N_lp":3]
1210 # [c:1][$(S(=O)=O):2]!@["N_lp":3]
1211 # [!#1:1][$(S(=O)=O):2]!@["N_lp":3]
1212 #
1213
1214 MatchedMappedAtoms = re.findall(":[0-9]", ElementNode.get("smarts"), re.I)
1215
1216 return True if len(MatchedMappedAtoms) == 3 else False
1217
1218 def _IsValidTypeTwoNitogenLonePairTorsionRule(self, ElementNode):
1219 """Validate atom map number for nitrogen involved in N_lp for type two nitrogen
1220 lone pair torsion rule."""
1221
1222 LonePairMapNumber = self._GetNitrogenLonePairMapNumber(ElementNode.get("smarts"))
1223
1224 return self._IsValidTypeTwoNitrogenLonePairMapNumber(LonePairMapNumber)
1225
1226 def _IsValidTypeTwoNitrogenLonePairMapNumber(self, LonePairMapNumber):
1227 """Check that the atom map number is 3."""
1228
1229 return True if LonePairMapNumber is not None and LonePairMapNumber == 3 else False
1230
1231 def _SetupTorsionAlertTypeForRotatableBond(self, TorsionAnglesInfo, TorsionAngle):
1232 """Setup torsion alert type and angle violation for a rotatable bond."""
1233
1234 TorsionCategory, TorsionAngleViolation = [None, None]
1235
1236 for ID in TorsionAnglesInfo["IDs"]:
1237 if self._IsTorsionAngleInWithinTolerance(
1238 TorsionAngle, TorsionAnglesInfo["Value"][ID], TorsionAnglesInfo["Tolerance1"][ID]
1239 ):
1240 TorsionCategory = "Green"
1241 TorsionAngleViolation = 0.0
1242 break
1243
1244 if self._IsTorsionAngleInWithinTolerance(
1245 TorsionAngle, TorsionAnglesInfo["Value"][ID], TorsionAnglesInfo["Tolerance2"][ID]
1246 ):
1247 TorsionCategory = "Orange"
1248 TorsionAngleViolation = self._CalculateTorsionAngleViolation(
1249 TorsionAngle, TorsionAnglesInfo["ValuesIn360RangeList"], TorsionAnglesInfo["Tolerances1List"]
1250 )
1251 break
1252
1253 if TorsionCategory is None:
1254 TorsionCategory = "Red"
1255 TorsionAngleViolation = self._CalculateTorsionAngleViolation(
1256 TorsionAngle, TorsionAnglesInfo["ValuesIn360RangeList"], TorsionAnglesInfo["Tolerances2List"]
1257 )
1258
1259 return (TorsionCategory, TorsionAngleViolation)
1260
1261 def _IsTorsionAngleInWithinTolerance(self, TorsionAngle, TorsionPeak, TorsionTolerance):
1262 """Check torsion angle against torsion tolerance."""
1263
1264 TorsionAngleDiff = TorsionAlertsUtil.CalculateTorsionAngleDifference(TorsionPeak, TorsionAngle)
1265
1266 return True if (abs(TorsionAngleDiff) <= TorsionTolerance) else False
1267
1268 def _CalculateTorsionAngleViolation(self, TorsionAngle, TorsionPeaks, TorsionTolerances):
1269 """Calculate torsion angle violation."""
1270
1271 TorsionAngleViolation = None
1272
1273 # Map angle to 0 to 360 range. TorsionPeaks values must be in this range...
1274 if TorsionAngle < 0:
1275 TorsionAngle = TorsionAngle + 360
1276
1277 # Identify the closest torsion peak index...
1278 if len(TorsionPeaks) == 1:
1279 NearestPeakIndex = 0
1280 else:
1281 NearestPeakIndex = min(range(len(TorsionPeaks)), key=lambda Index: abs(TorsionPeaks[Index] - TorsionAngle))
1282
1283 # Calculate torsion angle violation from the nearest peak and its tolerance value...
1284 TorsionAngleDiff = TorsionAlertsUtil.CalculateTorsionAngleDifference(
1285 TorsionPeaks[NearestPeakIndex], TorsionAngle
1286 )
1287 TorsionAngleViolation = abs(abs(TorsionAngleDiff) - TorsionTolerances[NearestPeakIndex])
1288
1289 return TorsionAngleViolation
