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