1 #
2 # File: TorsionStrainEnergyAlerts.py
3 # Author: Manish Sud <msud@san.rr.com>
4 #
5 # Collaborator: Pat Walters
6 #
7 # Copyright (C) 2026 Manish Sud. All rights reserved.
8 #
9 # This module uses the torsion strain energy library developed by Gu, S.;
10 # Smith, M. S.; Yang, Y.; Irwin, J. J.; Shoichet, B. K. [ Ref 153 ].
11 #
12 # The torsion strain enegy library is based on the Torsion Library jointly
13 # developed by the University of Hamburg, Center for Bioinformatics,
14 # Hamburg, Germany and F. Hoffmann-La-Roche Ltd., Basel, Switzerland.
15 #
16 # This file is part of MayaChemTools.
17 #
18 # MayaChemTools is free software; you can redistribute it and/or modify it under
19 # the terms of the GNU Lesser General Public License as published by the Free
20 # Software Foundation; either version 3 of the License, or (at your option) any
21 # later version.
22 #
23 # MayaChemTools is distributed in the hope that it will be useful, but without
24 # any warranty; without even the implied warranty of merchantability of fitness
25 # for a particular purpose. See the GNU Lesser General Public License for more
26 # details.
27 #
28 # You should have received a copy of the GNU Lesser General Public License
29 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
30 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
31 # Boston, MA, 02111-1307, USA.
32 #
33
34 import os
35 import re
36 import math
37
38 from rdkit import Chem
39 from rdkit.Chem import rdMolTransforms
40
41 from . import TorsionAlertsUtil
42
43
44 class TorsionStrainEnergyAlerts:
45 def __init__(
46 self,
47 AlertsMode="TotalEnergy",
48 TotalEnergyCutoff=6.0,
49 MaxSingleEnergyCutoff=1.8,
50 RotBondsSMARTSMode="SemiStrict",
51 RotBondsSMARTSPattern=None,
52 TorsionLibraryFilePath="auto",
53 AlertTorsionsNotObserved=False,
54 ):
55 """Identify strained molecules based on torsion strain energy library [ Ref 153 ]
56 alerts by matching rotatable bonds against SMARTS patterns specified for
57 torsion rules in a torsion energy library file. The molecules must have 3D coordinates.
58 The default torsion strain energy library file, TorsionStrainEnergyLibrary.xml, is
59 ia available in the directory containing this file.
60
61 The data in torsion strain energy library file is organized in a hierarchical
62 manner. It consists of one generic class and six specific classes at the highest
63 level. Each class contains multiple subclasses corresponding to named functional
64 groups or substructure patterns. The subclasses consist of torsion rules sorted
65 from specific to generic torsion patterns. The torsion rule, in turn, contains a
66 list of peak values for torsion angles and two tolerance values. A pair of tolerance
67 values define torsion bins around a torsion peak value.
68
69 A strain energy calculation method, 'exact' or 'approximate' [ Ref 153 ], is
70 associated with each torsion rule for calculating torsion strain energy. The 'exact'
71 stain energy calculation relies on the energy bins available under the energy histogram
72 consisting of 36 bins covering angles from -180 to 180. The width of each bin is 10
73 degree. The energy bins are are defined at the right end points. The first and the
74 last energy bins correspond to -170 and 180 respectively. The torsion angle is mapped
75 to a energy bin. An angle offset is calculated for the torsion angle from the the right
76 end point angle of the bin. The strain energy is estimated for the angle offset based
77 on the energy difference between the current and previous bins. The torsion strain
78 energy, in terms of torsion energy units (TEUs), corresponds to the sum of bin strain
79 energy and the angle offset strain energy.
80
81 Energy = BinEnergyDiff/10.0 * BinAngleOffset + BinEnergy[BinNum]
82
83 Where:
84
85 BinEnergyDiff = BinEnergy[BinNum] - BinEnergy[PreviousBinNum]
86 BinAngleOffset = TorsionAngle - BinAngleRightSide
87
88 The 'approximate' strain energy calculation relies on the angle difference between a
89 torsion angle and the torsion peaks observed for the torsion rules in the torsion
90 energy library. The torsion angle is matched to a torsion peak based on the value of
91 torsion angle difference. It must be less than or equal to the value for the second
92 tolerance 'tolerance2'. Otherwise, the torsion angle is not observed in the torsion
93 energy library and a value of 'NA' is assigned for torsion energy along with the lower
94 and upper bounds on energy at 95% confidence interval. The 'approximate' torsion
95 energy (TEUs) for observed torsion angle is calculated using the following formula:
96
97 Energy = beta_1 * (AngleDiff ** 2) + beta_2 * (AngleDiff ** 4)
98
99 The coefficients 'beta_1' and 'beta_2' are available for the observed angles in
100 the torsion strain energy library. The 'AngleDiff' is the difference between the
101 torsion angle and the matched torsion peak.
102
103 For example:
104
105 <library>
106 <hierarchyClass id1="G" id2="G" name="GG">
107 ...
108 </hierarchyClass>
109 <hierarchyClass id1="C" id2="O" name="CO">
110 <hierarchySubClass name="Ester bond I" smarts="O=[C:2][O:3]">
111 <torsionRule method="exact" smarts=
112 "[O:1]=[C:2]!@[O:3]~[CH0:4]">
113 <angleList>
114 <angle score="56.52" tolerance1="20.00"
115 tolerance2="25.00" value="0.0"/>
116 </angleList>
117 <histogram>
118 <bin count="1"/>
119 ...
120 </histogram>
121 <histogram_shifted>
122 <bin count="0"/>
123 ...
124 </histogram_shifted>
125 <histogram_converted>
126 <bin energy="4.67... lower="2.14..." upper="Inf"/>
127 ...
128 <bin energy="1.86..." lower="1.58..." upper="2.40..."/>
129 ...
130 </histogram_converted>
131 </torsionRule>
132 <torsionRule method="approximate" smarts=
133 "[cH0:1][c:2]([cH0])!@[O:3][p:4]">
134 <angleList>
135 <angle beta_1="0.002..." beta_2="-7.843...e-07"
136 score="27.14" theta_0="-90.0" tolerance1="30.00"
137 tolerance2="45.00" value="-90.0"/>
138 ...
139 </angleList>
140 <histogram>
141 <bin count="0"/>
142 ...
143 </histogram>
144 <histogram_shifted>
145 <bin count="0"/>
146 ...
147 </histogram_shifted>
148 </torsionRule>
149 ...
150 ...
151 </hierarchyClass>
152 <hierarchyClass id1="N" id2="C" name="NC">
153 ...
154 </hierarchyClass>
155 <hierarchyClass id1="S" id2="N" name="SN">
156 ...
157 </hierarchyClass>
158 <hierarchyClass id1="C" id2="S" name="CS">
159 ...
160 </hierarchyClass>
161 <hierarchyClass id1="C" id2="C" name="CC">
162 ...
163 </hierarchyClass>
164 <hierarchyClass id1="S" id2="S" name="SS">
165 ...
166 </hierarchyClass>
167 </library>
168
169 The rotatable bonds in a 3D molecule are identified using a default SMARTS pattern.
170 A custom SMARTS pattern may be optionally specified to detect rotatable bonds.
171 Each rotatable bond is matched to a torsion rule in the torsion strain energy library.
172 The strain energy is calculated for each rotatable bond using the calculation
173 method, 'exact' or 'approximate', associated with the matched torsion rule.
174
175 The total strain energy (TEUs) of a molecule corresponds to the sum of 'exact' and
176 'approximate' strain energies calculated for all matched rotatable bonds in the
177 molecule. The total strain energy is set to 'NA' for molecules containing a 'approximate'
178 energy estimate for a torsion angle not observed in the torsion energy library. In
179 addition, the lower and upper bounds on energy at 95% confidence interval are
180 set to 'NA'.
181
182 Arguments:
183 AlertsMode (str): Torsion strain energy library alert types to use
184 for issuing alerts. Possible values: TotalEnergy,
185 MaxSingleEnergy, or TotalOrMaxSingleEnergy.
186 TotalEnergyCutoff (float): Total strain strain energy (TEUs) cutoff.
187 MaxSingleEnergyCutoff (float): Maximum single strain energy (TEUs)
188 cutoff.
189 RotBondsSMARTSMode (str): SMARTS pattern to use for identifying
190 rotatable bonds in a molecule. Possible values: NonStrict,
191 SemiStrict, Strict or Specify.
192 RotBondsSMARTSPattern (str):SMARTS pattern for identifying rotatable
193 bonds. This paramater is only valid for 'Specify' value
194 'RotBondsSMARTSMode'.
195 TorsionLibraryFilePath (str): A XML file name containing data for
196 torsion starin energy library.
197 AlertTorsionsNotObserved (bool): Issue alerts about molecules
198 containing torsion angles not observed in torsion strain energy
199 library.
200
201 Returns:
202 object: An instantiated class object.
203
204 Examples:
205
206 StrainEnergyAlertsHandle = TorsionStrainEnergyAlerts()
207
208 StrainEnergyAlertsHandle = TorsionStrainEnergyAlerts(AlertsMode = "TotalEnergy",
209 TotalEnergyCutoff = 6.0, MaxSingleEnergyCutoff = 1.8, RotBondsSMARTSMode = "SemiStrict",
210 RotBondsSMARTSPattern = None, TorsionLibraryFilePath = "auto", AlertTorsionsNotObserved = False)
211
212 Notes:
213 The following sections provide additional details for the parameters.
214
215 AlertsMode: Torsion strain energy library alert types to use for issuing
216 alerts about molecules containing rotatable bonds based on the calculated
217 values for the total torsion strain energy of a molecule and the maximum
218 single strain energy of a rotatable bond in a molecule.
219
220 Possible values: TotalEnergy, MaxSingleEnergy, or TotalOrMaxSingleEnergy
221
222 The strain energy cutoff values in terms of torsion energy units (TEUs) are
223 used to filter molecules as shown below:
224
225 AlertsMode AlertsEnergyCutoffs (TEUs)
226
227 TotalEnergy >= TotalEnergyCutoff
228
229 MaxSingleEnergy >= MaxSingleEnergyCutoff
230
231 TotalOrMaxSingleEnergy >= TotalEnergyCutoff
232 or >= MaxSingleEnergyCutoff
233
234 TotalEnergyCutoff: Total strain strain energy (TEUs) cutoff [ Ref 153 ] for
235 issuing alerts based on total strain energy for all rotatable bonds in a
236 molecule. This option is used during 'TotalEnergy' or 'TotalOrMaxSingleEnergy'
237 values of 'AlertsMode' parameter.
238
239 The total strain energy must be greater than or equal to the specified
240 cutoff value to identify a molecule containing strained torsions.
241
242 MaxSingleEnergyCutoff: Maximum single strain energy (TEUs) cutoff [ Ref 153 ]
243 for issuing alerts based on the maximum value of a single strain energy of a
244 rotatable bond in a molecule. This option is used during 'MaxSingleEnergy' or
245 'TotalOrMaxSingleEnergy' values of 'AlertsMode' parameter.
246
247 The maximum single strain energy must be greater than or equal to the
248 specified cutoff valueo to identify a molecule containing strained torsions.
249
250 RotBondsSMARTSMode: SMARTS pattern to use for identifying rotatable bonds in
251 a molecule for matching against torsion rules in the torsion library. Possible values:
252 NonStrict, SemiStrict, Strict or Specify. The rotatable bond SMARTS matches
253 are filtered to ensure that each atom in the rotatable bond is attached to
254 at least two heavy atoms.
255
256 The following SMARTS patterns are used to identify rotatable bonds for
257 different modes:
258
259 NonStrict: [!$(*#*)&!D1]-&!@[!$(*#*)&!D1]
260
261 SemiStrict:
262 [!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)
263 &!$(C([CH3])([CH3])[CH3])]-!@[!$(*#*)&!D1&!$(C(F)(F)F)
264 &!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]
265
266 Strict:
267 [!$(*#*)&!D1&!$(C(F)(F)F)&!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)
268 &!$(C([CH3])([CH3])[CH3])&!$([CD3](=[N,O,S])-!@[#7,O,S!D1])
269 &!$([#7,O,S!D1]-!@[CD3]=[N,O,S])&!$([CD3](=[N+])-!@[#7!D1])
270 &!$([#7!D1]-!@[CD3]=[N+])]-!@[!$(*#*)&!D1&!$(C(F)(F)F)
271 &!$(C(Cl)(Cl)Cl)&!$(C(Br)(Br)Br)&!$(C([CH3])([CH3])[CH3])]
272
273 The 'NonStrict' and 'Strict' SMARTS patterns are available in RDKit. The
274 'NonStrict' SMARTS pattern corresponds to original Daylight SMARTS
275 specification for rotatable bonds. The 'SemiStrict' SMARTS pattern is
276 derived from 'Strict' SMARTS pattern.
277
278 TorsionLibraryFilePath (str): Specify a XML file name containing data for
279 torsion starin energy library hierarchy or use default file,
280 TorsionEnergyLibrary.xml, available in the directory containing this file.
281
282 AlertTorsionsNotObserved: Issue alerts abpout molecules con ntaining torsion
283 angles not observed in torsion strain energy library. It's not possible to
284 calculate torsion strain energies for these torsions during 'approximate'
285 match to a specified torsion in the library.
286
287 The 'approximate' strain energy calculation relies on the angle difference
288 between a torsion angle and the torsion peaks observed for the torsion
289 rules in the torsion energy library. The torsion angle is matched to a
290 torsion peak based on the value of torsion angle difference. It must be
291 less than or equal to the value for the second tolerance 'tolerance2'.
292 Otherwise, the torsion angle is not observed in the torsion energy library
293 and a value of 'NA' is assigned for torsion energy along with the lower and
294 upper bounds on energy at 95% confidence interval.
295
296 """
297
298 self._ProcessTorsionStrainEnergyAlertsParameters(
299 AlertsMode,
300 TotalEnergyCutoff,
301 MaxSingleEnergyCutoff,
302 RotBondsSMARTSMode,
303 RotBondsSMARTSPattern,
304 TorsionLibraryFilePath,
305 AlertTorsionsNotObserved,
306 )
307
308 self._InitializeTorsionStrainEnergyAlerts()
309
310 def _ProcessTorsionStrainEnergyAlertsParameters(
311 self,
312 AlertsMode,
313 TotalEnergyCutoff,
314 MaxSingleEnergyCutoff,
315 RotBondsSMARTSMode,
316 RotBondsSMARTSPattern,
317 TorsionLibraryFilePath,
318 AlertTorsionsNotObserved,
319 ):
320 """Process torsion strain energy alerts paramaters."""
321
322 # Process AltertsMode paramater...
323 self._ProcessAlertsModeParameter(AlertsMode)
324
325 # Process TotalEnergyCutoff parameter...
326 if TotalEnergyCutoff <= 0.0:
327 raise ValueError(
328 "The value, %s, specified for TotalEnergyCutoff parameter is not valid. Supported value: > 0.0"
329 % TotalEnergyCutoff
330 )
331 self.TotalEnergyCutoff = TotalEnergyCutoff
332
333 # Process MaxSingleEnergyCutoff parameter...
334 if MaxSingleEnergyCutoff <= 0.0:
335 raise ValueError(
336 "The value, %s, specified for MaxSingleEnergyCutoff parameter is not valid. Supported value: > 0.0"
337 % MaxSingleEnergyCutoff
338 )
339 self.MaxSingleEnergyCutoff = MaxSingleEnergyCutoff
340
341 # Process RotBondsSMARTSMode and RotBondsSMARTSPattern parameters...
342 self._ProcessRotBondsParameters(RotBondsSMARTSMode, RotBondsSMARTSPattern)
343
344 # Process TorsionLibraryFilePath parameter...
345 self._ProcessTorsionLibraryFilePathParameter(TorsionLibraryFilePath)
346
347 # AlertTorsionsNotObserved parameter...
348 self.AlertTorsionsNotObserved = AlertTorsionsNotObserved
349
350 def _InitializeTorsionStrainEnergyAlerts(self):
351 """Initialize tosion strain energy alerts."""
352
353 # Retrieve and setup torsion library info for matching rotatable bonds...
354 TorsionLibraryInfo = {}
355
356 TorsionLibElementTree = TorsionAlertsUtil.RetrieveTorsionLibraryInfo(self.TorsionLibraryFilePath)
357 TorsionLibraryInfo["TorsionLibElementTree"] = TorsionLibElementTree
358
359 TorsionAlertsUtil.SetupTorsionLibraryInfoForMatchingRotatableBonds(TorsionLibraryInfo)
360
361 self.TorsionLibraryInfo = TorsionLibraryInfo
362
363 def _ProcessAlertsModeParameter(self, AlertsMode):
364 """Process AlertsMode parameter."""
365
366 TotalEnergyMode, MaxSingleEnergyMode, TotalOrMaxSingleEnergyMode = [False] * 3
367 if re.match("^TotalEnergy$", AlertsMode, re.I):
368 TotalEnergyMode = True
369 elif re.match("^MaxSingleEnergy$", AlertsMode, re.I):
370 MaxSingleEnergyMode = True
371 elif re.match("^TotalOrMaxSingleEnergy$", AlertsMode, re.I):
372 TotalOrMaxSingleEnergyMode = True
373 else:
374 raise ValueError(
375 "Invalid value, %s, specified for AlertsMode parameter. Valid values: TotalEnergy, MaxSingleEnergy, or TotalOrMaxSingleEnergy"
376 % (AlertsMode)
377 )
378
379 self.AltersMode = AlertsMode
380 self.TotalEnergyMode = TotalEnergyMode
381 self.MaxSingleEnergyMode = MaxSingleEnergyMode
382 self.TotalOrMaxSingleEnergyMode = TotalOrMaxSingleEnergyMode
383
384 def _ProcessRotBondsParameters(self, RotBondsSMARTSMode, RotBondsSMARTSPattern):
385 """Process RotBondsSMARTSMode and RotBondsSMARTSPattern parameters."""
386
387 if re.match("^NonStrict$", RotBondsSMARTSMode, re.I):
388 RotBondsSMARTSPattern = "[!$(*#*)&!D1]-&!@[!$(*#*)&!D1]"
389 elif re.match("^SemiStrict$", RotBondsSMARTSMode, re.I):
390 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])]"
391 elif re.match("^Strict$", RotBondsSMARTSMode, re.I):
392 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])]"
393 elif re.match("^Specify$", RotBondsSMARTSMode, re.I):
394 if RotBondsSMARTSPattern is None:
395 raise ValueError(
396 "The RotBondsSMARTSPattern parameter value must be specified during Specify value for RotBondsSMARTSMode parameter."
397 )
398 RotBondsSMARTSPattern = RotBondsSMARTSPattern.strip()
399 if not len(RotBondsSMARTSPattern):
400 raise ValueError("Empty value specified for RotBondsSMARTSPattern parameter.")
401 else:
402 raise ValueError(
403 "Invalid value, %s, specified for RotBondsSMARTSMode parameter. Valid values: NonStrict, SemiStrict, Strict or Specify"
404 % (RotBondsSMARTSMode)
405 )
406
407 RotBondsPatternMol = Chem.MolFromSmarts(RotBondsSMARTSPattern)
408 if RotBondsPatternMol is None:
409 if re.match("Specify", RotBondsSMARTSMode, re.I):
410 raise ValueError(
411 'Failed to create rotatable bonds pattern molecule. The rotatable bonds SMARTS pattern, "%s", specified using RotBondsSMARTSPattern parameter is not valid.'
412 % (RotBondsSMARTSPattern)
413 )
414 else:
415 raise ValueError(
416 'Failed to create rotatable bonds pattern molecule. The default rotatable bonds SMARTS pattern, "%s", used for %s value of RotBondsSMARTSMode parameter is not valid.'
417 % (RotBondsSMARTSPattern, RotBondsSMARTSMode)
418 )
419
420 self.RotBondsSMARTSMode = RotBondsSMARTSMode
421 self.RotBondsSMARTSPattern = RotBondsSMARTSPattern
422 self.RotBondsPatternMol = RotBondsPatternMol
423
424 def _ProcessTorsionLibraryFilePathParameter(self, TorsionLibraryFilePath):
425 """Process torsion library path."""
426
427 # TorsionLibraryFilePath parameter...
428 if re.match("^auto$", TorsionLibraryFilePath):
429 TorsionLibraryFile = "TorsionStrainEnergyLibrary.xml"
430 TorsionLibraryFilePath = os.path.join(os.path.dirname(os.path.abspath(__file__)), TorsionLibraryFile)
431 if not os.path.isfile(TorsionLibraryFilePath):
432 raise ValueError("The default torsion alerts library file %s doesn't exist." % TorsionLibraryFilePath)
433 else:
434 if not os.path.isfile(TorsionLibraryFilePath):
435 raise ValueError(
436 "The file specified, %s, for parameter TorsionLibraryFilePath doesn't exist."
437 % TorsionLibraryFilePath
438 )
439 TorsionLibraryFilePath = os.path.abspath(TorsionLibraryFilePath)
440
441 self.TorsionLibraryFilePath = TorsionLibraryFilePath
442
443 def GetTorsionLibraryFilePath(self):
444 """Get torsion strain library file path.
445
446 Arguments:
447 None
448
449 Returns:
450 FilePath (str): Torsion strain library path.
451
452 """
453
454 return self.TorsionLibraryFilePath
455
456 def ListTorsionLibraryInfo(self):
457 """List torsion strain library information.
458
459 Arguments:
460 None
461
462 Returns:
463 None
464
465 """
466
467 return TorsionAlertsUtil.ListTorsionLibraryInfo(self.TorsionLibraryInfo["TorsionLibElementTree"])
468
469 def IdentifyTorsionLibraryAlertsForRotatableBonds(self, Mol):
470 """Identify torsion strain library alerts for a molecule by matching rotatable
471 bonds against SMARTS patterns specified for torsion rules in torsion energy
472 library file.
473
474 Arguments:
475 Mol (object): RDKit molecule object.
476
477 Returns:
478 bool: True - Molecule contains strained torsions; False - Molecule
479 contains no strained torsions or rotatable bonds.
480 dict or None: Torsion alerts information regarding matching of
481 rotatable bonds to torsion strain library.
482
483 Examples:
484
485 from TorsionAlerts.TorsionStrainEnergyAlerts import
486 TorsionStrainEnergyAlerts
487
488 StrainEnergyAlerts = TorsionStrainEnergyAlerts()
489 AlertsStatus, AlertsInfo = StrainEnergyAlerts.
490 IdentifyTorsionLibraryAlertsForRotatableBonds(RDKitMol)
491
492 RotBondsAlertsStatus = AlertsInfo["RotBondsAlertsStatus"]
493 TotalEnergy = AlertsInfo["TotalEnergy"]
494 TotalEnergyLowerBound = AlertsInfo["TotalEnergyLowerBound"]
495 TotalEnergyUpperBound = AlertsInfo["TotalEnergyUpperBound"]
496 AnglesNotObservedCount = AlertsInfo["AnglesNotObservedCount"]
497 MaxSingleEnergy = AlertsInfo["MaxSingleEnergy"])
498 MaxSingleEnergyAlertsCount = AlertsInfo[
499 "MaxSingleEnergyAlertsCount"]
500
501 # List of rotatable bond IDs...
502 RotatableBondIDs = AlertsInfo["IDs"]
503
504 # Dictionaries containing information for rotatable bonds by using
505 # bond ID as key...
506 for ID in AlertsInfo["IDs"]:
507 MatchStatus = AlertsInfo["MatchStatus"][ID]
508 MaxSingleEnergyAlertStatus = AlertsInfo[
509 "MaxSingleEnergyAlertStatus"][ID]
510 AtomIndices = AlertsInfo["AtomIndices"][ID]
511 TorsionAtomIndices = AlertsInfo["TorsionAtomIndices"][ID]
512 TorsionAngle = AlertsInfo["TorsionAngle"][ID]
513 HierarchyClassName = AlertsInfo["HierarchyClassName"][ID]
514 HierarchySubClassName = AlertsInfo["HierarchySubClassName"][ID]
515 TorsionRuleNodeID = AlertsInfo["TorsionRuleNodeID"][ID]
516 TorsionRuleSMARTS = AlertsInfo["TorsionRuleSMARTS"][ID]
517 EnergyMethod = AlertsInfo["EnergyMethod"][ID]
518 AngleNotObserved = AlertsInfo["AngleNotObserved"][ID]
519 Energy = AlertsInfo["Energy"][ID]
520 EnergyLowerBound = AlertsInfo["EnergyLowerBound"][ID]
521 EnergyUpperBound = AlertsInfo["EnergyUpperBound"][ID]
522
523 """
524
525 # Identify rotatable bonds...
526 RotBondsStatus, RotBondsInfo = TorsionAlertsUtil.IdentifyRotatableBondsForTorsionLibraryMatch(
527 self.TorsionLibraryInfo, Mol, self.RotBondsPatternMol
528 )
529
530 if not RotBondsStatus:
531 return (False, None)
532
533 # Identify alerts for rotatable bonds...
534 RotBondsAlertsStatus, RotBondsAlertsInfo = self._MatchRotatableBondsToTorsionLibrary(Mol, RotBondsInfo)
535
536 return (RotBondsAlertsStatus, RotBondsAlertsInfo)
537
538 def _MatchRotatableBondsToTorsionLibrary(self, Mol, RotBondsInfo):
539 """Match rotatable bond to torsion library."""
540
541 # Initialize...
542 RotBondsAlertsInfo = self._InitializeRotatableBondsAlertsInfo()
543
544 # Match rotatable bonds to torsion library...
545 for ID in RotBondsInfo["IDs"]:
546 AtomIndices = RotBondsInfo["AtomIndices"][ID]
547 HierarchyClass = RotBondsInfo["HierarchyClass"][ID]
548
549 MatchStatus, MatchInfo = self._MatchRotatableBondToTorsionLibrary(Mol, AtomIndices, HierarchyClass)
550 self._TrackRotatableBondsAlertsInfo(RotBondsAlertsInfo, ID, AtomIndices, MatchStatus, MatchInfo)
551
552 RotBondsAlertsStatus = self._SetupRotatableBondsAlertStatusTotalStrainEnergiesInfo(RotBondsAlertsInfo)
553
554 return (RotBondsAlertsStatus, RotBondsAlertsInfo)
555
556 def _InitializeRotatableBondsAlertsInfo(self):
557 """Initialize alerts information for rotatable bonds."""
558
559 RotBondsAlertsInfo = {}
560 RotBondsAlertsInfo["IDs"] = []
561
562 for DataLabel in [
563 "RotBondsAlertsStatus",
564 "TotalEnergy",
565 "TotalEnergyLowerBound",
566 "TotalEnergyUpperBound",
567 "AnglesNotObservedCount",
568 "MaxSingleEnergy",
569 "MaxSingleEnergyAlertsCount",
570 ]:
571 RotBondsAlertsInfo[DataLabel] = None
572
573 for DataLabel in [
574 "MatchStatus",
575 "MaxSingleEnergyAlertStatus",
576 "AtomIndices",
577 "TorsionAtomIndices",
578 "TorsionAngle",
579 "HierarchyClassName",
580 "HierarchySubClassName",
581 "TorsionRuleNodeID",
582 "TorsionRuleSMARTS",
583 "EnergyMethod",
584 "AngleNotObserved",
585 "Energy",
586 "EnergyLowerBound",
587 "EnergyUpperBound",
588 ]:
589 RotBondsAlertsInfo[DataLabel] = {}
590
591 return RotBondsAlertsInfo
592
593 def _TrackRotatableBondsAlertsInfo(self, RotBondsAlertsInfo, ID, AtomIndices, MatchStatus, MatchInfo):
594 """Track alerts information for rotatable bonds."""
595
596 if MatchInfo is None or len(MatchInfo) == 0:
597 (
598 TorsionAtomIndices,
599 TorsionAngle,
600 HierarchyClassName,
601 HierarchySubClassName,
602 TorsionRuleNodeID,
603 TorsionRuleSMARTS,
604 EnergyMethod,
605 AngleNotObserved,
606 Energy,
607 EnergyLowerBound,
608 EnergyUpperBound,
609 ) = [None] * 11
610 else:
611 (
612 TorsionAtomIndices,
613 TorsionAngle,
614 HierarchyClassName,
615 HierarchySubClassName,
616 TorsionRuleNodeID,
617 TorsionRuleSMARTS,
618 EnergyMethod,
619 AngleNotObserved,
620 Energy,
621 EnergyLowerBound,
622 EnergyUpperBound,
623 ) = MatchInfo
624
625 # Track torsion match information...
626 RotBondsAlertsInfo["IDs"].append(ID)
627 RotBondsAlertsInfo["MatchStatus"][ID] = MatchStatus
628 RotBondsAlertsInfo["AtomIndices"][ID] = AtomIndices
629 RotBondsAlertsInfo["TorsionAtomIndices"][ID] = TorsionAtomIndices
630 RotBondsAlertsInfo["TorsionAngle"][ID] = TorsionAngle
631 RotBondsAlertsInfo["HierarchyClassName"][ID] = HierarchyClassName
632 RotBondsAlertsInfo["HierarchySubClassName"][ID] = HierarchySubClassName
633 RotBondsAlertsInfo["TorsionRuleNodeID"][ID] = TorsionRuleNodeID
634 RotBondsAlertsInfo["TorsionRuleSMARTS"][ID] = TorsionRuleSMARTS
635 RotBondsAlertsInfo["EnergyMethod"][ID] = EnergyMethod
636 RotBondsAlertsInfo["AngleNotObserved"][ID] = AngleNotObserved
637 RotBondsAlertsInfo["Energy"][ID] = Energy
638 RotBondsAlertsInfo["EnergyLowerBound"][ID] = EnergyLowerBound
639 RotBondsAlertsInfo["EnergyUpperBound"][ID] = EnergyUpperBound
640
641 def _SetupRotatableBondsAlertStatusTotalStrainEnergiesInfo(self, RotBondsAlertsInfo):
642 """Setup rotatable bonds alert status along with total strain energies."""
643
644 # Initialize...
645 RotBondsAlertsStatus = False
646 TotalEnergy, TotalEnergyLowerBound, TotalEnergyUpperBound, AnglesNotObservedCount = [None, None, None, None]
647 MaxSingleEnergy, MaxSingleEnergyAlertsCount = [None, None]
648
649 # Initialize max single energy alert status...
650 for ID in RotBondsAlertsInfo["IDs"]:
651 RotBondsAlertsInfo["MaxSingleEnergyAlertStatus"][ID] = None
652
653 # Check for torsion angles not obervered in the strain library...
654 AnglesNotObservedCount = 0
655 for ID in RotBondsAlertsInfo["IDs"]:
656 AngleNotObserved = RotBondsAlertsInfo["AngleNotObserved"][ID]
657 if AngleNotObserved is not None and AngleNotObserved:
658 AnglesNotObservedCount += 1
659
660 # Setup alert status for rotable bonds...
661 if AnglesNotObservedCount > 0:
662 if self.AlertTorsionsNotObserved:
663 RotBondsAlertsStatus = True
664 else:
665 TotalEnergy = 0.0
666 for ID in RotBondsAlertsInfo["IDs"]:
667 Energy = RotBondsAlertsInfo["Energy"][ID]
668 TotalEnergy += Energy
669 if self.TotalEnergyMode:
670 if TotalEnergy > self.TotalEnergyCutoff:
671 RotBondsAlertsStatus = True
672 break
673 elif self.MaxSingleEnergyMode:
674 if Energy > self.MaxSingleEnergyCutoff:
675 RotBondsAlertsStatus = True
676 break
677 elif self.TotalOrMaxSingleEnergyMode:
678 if TotalEnergy > self.TotalEnergyCutoff or Energy > self.MaxSingleEnergyCutoff:
679 RotBondsAlertsStatus = True
680 break
681
682 # Setup energy infomation...
683 TotalEnergy, TotalEnergyLowerBound, TotalEnergyUpperBound = [0.0, 0.0, 0.0]
684 if self.MaxSingleEnergyMode or self.TotalOrMaxSingleEnergyMode:
685 MaxSingleEnergy, MaxSingleEnergyAlertsCount = [0.0, 0]
686
687 for ID in RotBondsAlertsInfo["IDs"]:
688 Energy = RotBondsAlertsInfo["Energy"][ID]
689
690 # Setup total energy along with the lower and upper bounds...
691 TotalEnergy += Energy
692 TotalEnergyLowerBound += RotBondsAlertsInfo["EnergyLowerBound"][ID]
693 TotalEnergyUpperBound += RotBondsAlertsInfo["EnergyUpperBound"][ID]
694
695 # Setup max single energy and max single energy alerts count...
696 if self.MaxSingleEnergyMode or self.TotalOrMaxSingleEnergyMode:
697 MaxSingleEnergyAlertStatus = False
698
699 if Energy > MaxSingleEnergy:
700 MaxSingleEnergy = Energy
701 if Energy > self.MaxSingleEnergyCutoff:
702 MaxSingleEnergyAlertStatus = True
703 MaxSingleEnergyAlertsCount += 1
704
705 RotBondsAlertsInfo["MaxSingleEnergyAlertStatus"][ID] = MaxSingleEnergyAlertStatus
706
707 RotBondsAlertsInfo["RotBondsAlertsStatus"] = RotBondsAlertsStatus
708
709 RotBondsAlertsInfo["TotalEnergy"] = TotalEnergy
710 RotBondsAlertsInfo["TotalEnergyLowerBound"] = TotalEnergyLowerBound
711 RotBondsAlertsInfo["TotalEnergyUpperBound"] = TotalEnergyUpperBound
712
713 RotBondsAlertsInfo["AnglesNotObservedCount"] = AnglesNotObservedCount
714
715 RotBondsAlertsInfo["MaxSingleEnergy"] = MaxSingleEnergy
716 RotBondsAlertsInfo["MaxSingleEnergyAlertsCount"] = MaxSingleEnergyAlertsCount
717
718 return RotBondsAlertsStatus
719
720 def _MatchRotatableBondToTorsionLibrary(self, Mol, RotBondAtomIndices, RotBondHierarchyClass):
721 """Match rotatable bond to torsion library."""
722
723 if TorsionAlertsUtil.IsSpecificHierarchyClass(self.TorsionLibraryInfo, RotBondHierarchyClass):
724 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstSpecificHierarchyClass(
725 Mol, RotBondAtomIndices, RotBondHierarchyClass
726 )
727 if not MatchStatus:
728 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstGenericHierarchyClass(
729 Mol, RotBondAtomIndices, RotBondHierarchyClass
730 )
731 else:
732 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstGenericHierarchyClass(
733 Mol, RotBondAtomIndices, RotBondHierarchyClass
734 )
735
736 return (MatchStatus, MatchInfo)
737
738 def _MatchRotatableBondAgainstSpecificHierarchyClass(self, Mol, RotBondAtomIndices, RotBondHierarchyClass):
739 """Match rotatable bond against a specific hierarchy class."""
740
741 TorsionLibraryInfo = self.TorsionLibraryInfo
742
743 HierarchyClassElementNode = None
744 if RotBondHierarchyClass in TorsionLibraryInfo["SpecificClasses"]["ElementNode"]:
745 HierarchyClassElementNode = TorsionLibraryInfo["SpecificClasses"]["ElementNode"][RotBondHierarchyClass]
746
747 if HierarchyClassElementNode is None:
748 return (False, None, None, None)
749
750 TorsionAlertsUtil.TrackHierarchyClassElementNode(TorsionLibraryInfo, HierarchyClassElementNode)
751 MatchStatus, MatchInfo = self._ProcessElementForRotatableBondMatch(
752 Mol, RotBondAtomIndices, HierarchyClassElementNode
753 )
754 TorsionAlertsUtil.RemoveLastHierarchyClassElementNodeFromTracking(TorsionLibraryInfo)
755
756 return (MatchStatus, MatchInfo)
757
758 def _MatchRotatableBondAgainstGenericHierarchyClass(self, Mol, RotBondAtomIndices, RotBondHierarchyClass):
759 """Match rotatable bond against a generic hierarchy class."""
760
761 TorsionLibraryInfo = self.TorsionLibraryInfo
762
763 HierarchyClassElementNode = TorsionAlertsUtil.GetGenericHierarchyClassElementNode(TorsionLibraryInfo)
764 if HierarchyClassElementNode is None:
765 return (False, None)
766
767 TorsionAlertsUtil.TrackHierarchyClassElementNode(TorsionLibraryInfo, HierarchyClassElementNode)
768
769 # Match hierarchy subclasses before matching torsion rules...
770 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstGenericHierarchySubClasses(
771 Mol, RotBondAtomIndices, HierarchyClassElementNode
772 )
773
774 if not MatchStatus:
775 MatchStatus, MatchInfo = self._MatchRotatableBondAgainstGenericHierarchyTorsionRules(
776 Mol, RotBondAtomIndices, HierarchyClassElementNode
777 )
778
779 TorsionAlertsUtil.RemoveLastHierarchyClassElementNodeFromTracking(TorsionLibraryInfo)
780
781 return (MatchStatus, MatchInfo)
782
783 def _MatchRotatableBondAgainstGenericHierarchySubClasses(self, Mol, RotBondAtomIndices, HierarchyClassElementNode):
784 """Match rotatable bond againat generic hierarchy subclasses."""
785
786 for ElementChildNode in HierarchyClassElementNode:
787 if ElementChildNode.tag != "hierarchySubClass":
788 continue
789
790 SubClassMatchStatus = self._ProcessHierarchySubClassElementForRotatableBondMatch(
791 Mol, RotBondAtomIndices, ElementChildNode
792 )
793
794 if SubClassMatchStatus:
795 MatchStatus, MatchInfo = self._ProcessElementForRotatableBondMatch(
796 Mol, RotBondAtomIndices, ElementChildNode
797 )
798
799 if MatchStatus:
800 return (MatchStatus, MatchInfo)
801
802 return (False, None)
803
804 def _MatchRotatableBondAgainstGenericHierarchyTorsionRules(
805 self, Mol, RotBondAtomIndices, HierarchyClassElementNode
806 ):
807 """Match rotatable bond againat torsion rules generic hierarchy class."""
808
809 for ElementChildNode in HierarchyClassElementNode:
810 if ElementChildNode.tag != "torsionRule":
811 continue
812
813 MatchStatus, MatchInfo = self._ProcessTorsionRuleElementForRotatableBondMatch(
814 Mol, RotBondAtomIndices, ElementChildNode
815 )
816
817 if MatchStatus:
818 return (MatchStatus, MatchInfo)
819
820 return (False, None)
821
822 def _ProcessElementForRotatableBondMatch(self, Mol, RotBondAtomIndices, ElementNode):
823 """Process element node to recursively match rotatable bond against hierarchy
824 subclasses and torsion rules."""
825
826 TorsionLibraryInfo = self.TorsionLibraryInfo
827
828 for ElementChildNode in ElementNode:
829 if ElementChildNode.tag == "hierarchySubClass":
830 SubClassMatchStatus = self._ProcessHierarchySubClassElementForRotatableBondMatch(
831 Mol, RotBondAtomIndices, ElementChildNode
832 )
833
834 if SubClassMatchStatus:
835 TorsionAlertsUtil.TrackHierarchySubClassElementNode(TorsionLibraryInfo, ElementChildNode)
836
837 MatchStatus, MatchInfo = self._ProcessElementForRotatableBondMatch(
838 Mol, RotBondAtomIndices, ElementChildNode
839 )
840 if MatchStatus:
841 TorsionAlertsUtil.RemoveLastHierarchySubClassElementNodeFromTracking(TorsionLibraryInfo)
842 return (MatchStatus, MatchInfo)
843
844 TorsionAlertsUtil.RemoveLastHierarchySubClassElementNodeFromTracking(TorsionLibraryInfo)
845
846 elif ElementChildNode.tag == "torsionRule":
847 MatchStatus, MatchInfo = self._ProcessTorsionRuleElementForRotatableBondMatch(
848 Mol, RotBondAtomIndices, ElementChildNode
849 )
850
851 if MatchStatus:
852 return (MatchStatus, MatchInfo)
853
854 return (False, None)
855
856 def _ProcessHierarchySubClassElementForRotatableBondMatch(self, Mol, RotBondAtomIndices, ElementNode):
857 """Process hierarchy subclass element to match rotatable bond."""
858
859 # Setup subclass SMARTS pattern mol...
860 SubClassPatternMol = TorsionAlertsUtil.SetupHierarchySubClassElementPatternMol(
861 self.TorsionLibraryInfo, ElementNode
862 )
863 if SubClassPatternMol is None:
864 return False
865
866 # Match SMARTS pattern...
867 SubClassPatternMatches = TorsionAlertsUtil.FilterSubstructureMatchesByAtomMapNumbers(
868 Mol, SubClassPatternMol, Mol.GetSubstructMatches(SubClassPatternMol, useChirality=False)
869 )
870 if len(SubClassPatternMatches) == 0:
871 return False
872
873 # Match rotatable bond indices...
874 RotBondAtomIndex1, RotBondAtomIndex2 = RotBondAtomIndices
875 MatchStatus = False
876 for SubClassPatternMatch in SubClassPatternMatches:
877 if len(SubClassPatternMatch) == 2:
878 # Matched to pattern containing map atom numbers ":2" and ":3"...
879 CentralAtomsIndex1, CentralAtomsIndex2 = SubClassPatternMatch
880 elif len(SubClassPatternMatch) == 4:
881 # Matched to pattern containing map atom numbers ":1", ":2", ":3" and ":4"...
882 CentralAtomsIndex1 = SubClassPatternMatch[1]
883 CentralAtomsIndex2 = SubClassPatternMatch[2]
884 elif len(SubClassPatternMatch) == 3:
885 SubClassSMARTSPattern = ElementNode.get("smarts")
886 if TorsionAlertsUtil.DoesSMARTSContainsMappedAtoms(SubClassSMARTSPattern, [":2", ":3", ":4"]):
887 # Matched to pattern containing map atom numbers ":2", ":3" and ":4"...
888 CentralAtomsIndex1 = SubClassPatternMatch[0]
889 CentralAtomsIndex2 = SubClassPatternMatch[1]
890 else:
891 # Matched to pattern containing map atom numbers ":1", ":2" and ":3"...
892 CentralAtomsIndex1 = SubClassPatternMatch[1]
893 CentralAtomsIndex2 = SubClassPatternMatch[2]
894 else:
895 continue
896
897 if CentralAtomsIndex1 != CentralAtomsIndex2:
898 if (CentralAtomsIndex1 == RotBondAtomIndex1 and CentralAtomsIndex2 == RotBondAtomIndex2) or (
899 CentralAtomsIndex1 == RotBondAtomIndex2 and CentralAtomsIndex2 == RotBondAtomIndex1
900 ):
901 MatchStatus = True
902 break
903
904 return MatchStatus
905
906 def _ProcessTorsionRuleElementForRotatableBondMatch(self, Mol, RotBondAtomIndices, ElementNode):
907 """Process torsion rule element to match rotatable bond."""
908
909 TorsionLibraryInfo = self.TorsionLibraryInfo
910
911 # Retrieve torsions matched to rotatable bond...
912 TorsionAtomIndicesList, TorsionAnglesList = self._MatchTorsionRuleToRotatableBond(
913 Mol, RotBondAtomIndices, ElementNode
914 )
915 if TorsionAtomIndicesList is None:
916 return (False, None)
917
918 # Setup torsion angles and enery bin information for matched torsion rule...
919 TorsionRuleAnglesInfo = TorsionAlertsUtil.SetupTorsionRuleAnglesInfo(TorsionLibraryInfo, ElementNode)
920 if TorsionRuleAnglesInfo is None:
921 return (False, None)
922
923 # Setup highest strain energy for matched torsions...
924 TorsionAtomIndices, TorsionAngle, EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound = (
925 self._SelectHighestStrainEnergyTorsionForRotatableBond(
926 TorsionRuleAnglesInfo, TorsionAtomIndicesList, TorsionAnglesList
927 )
928 )
929
930 # Setup hierarchy class and subclass names...
931 HierarchyClassName, HierarchySubClassName = (
932 TorsionAlertsUtil.SetupHierarchyClassAndSubClassNamesForRotatableBond(TorsionLibraryInfo)
933 )
934
935 # Setup rule node ID...
936 TorsionRuleNodeID = ElementNode.get("NodeID")
937
938 # Setup SMARTS...
939 TorsionRuleSMARTS = ElementNode.get("smarts")
940 if " " in TorsionRuleSMARTS:
941 TorsionRuleSMARTS = TorsionRuleSMARTS.replace(" ", "")
942
943 # Setup match info...
944 MatchInfo = [
945 TorsionAtomIndices,
946 TorsionAngle,
947 HierarchyClassName,
948 HierarchySubClassName,
949 TorsionRuleNodeID,
950 TorsionRuleSMARTS,
951 EnergyMethod,
952 AngleNotObserved,
953 Energy,
954 EnergyLowerBound,
955 EnergyUpperBound,
956 ]
957
958 # Setup match status...
959 MatchStatus = True
960
961 return (MatchStatus, MatchInfo)
962
963 def _SelectHighestStrainEnergyTorsionForRotatableBond(
964 self, TorsionRuleAnglesInfo, TorsionAtomIndicesList, TorsionAnglesList
965 ):
966 """Select highest strain energy torsion matched to a rotatable bond."""
967
968 TorsionAtomIndices, TorsionAngle, EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound = [
969 None
970 ] * 7
971 ValidEnergyValue, ValidCurrentEnergyValue = [False] * 2
972
973 FirstTorsion = True
974 for Index in range(0, len(TorsionAtomIndicesList)):
975 CurrentTorsionAtomIndices = TorsionAtomIndicesList[Index]
976 CurrentTorsionAngle = TorsionAnglesList[Index]
977
978 if FirstTorsion:
979 FirstTorsion = False
980 EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound = (
981 self._SetupStrainEnergyForRotatableBond(TorsionRuleAnglesInfo, CurrentTorsionAngle)
982 )
983 TorsionAtomIndices = CurrentTorsionAtomIndices
984 TorsionAngle = CurrentTorsionAngle
985 ValidEnergyValue = self._IsEnergyValueValid(Energy)
986 continue
987
988 # Select highest strain energy...
989 (
990 CurrentEnergyMethod,
991 CurrentAngleNotObserved,
992 CurrentEnergy,
993 CurrentEnergyLowerBound,
994 CurrentEnergyUpperBound,
995 ) = self._SetupStrainEnergyForRotatableBond(TorsionRuleAnglesInfo, CurrentTorsionAngle)
996 ValidCurrentEnergyValue = self._IsEnergyValueValid(CurrentEnergy)
997
998 UpdateValues = False
999 if ValidEnergyValue and ValidCurrentEnergyValue:
1000 if CurrentEnergy > Energy:
1001 UpdateValues = True
1002 elif ValidCurrentEnergyValue:
1003 if not ValidEnergyValue:
1004 UpdateValues = True
1005
1006 if UpdateValues:
1007 EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound = [
1008 CurrentEnergyMethod,
1009 CurrentAngleNotObserved,
1010 CurrentEnergy,
1011 CurrentEnergyLowerBound,
1012 CurrentEnergyUpperBound,
1013 ]
1014 TorsionAtomIndices = CurrentTorsionAtomIndices
1015 TorsionAngle = CurrentTorsionAngle
1016
1017 return (
1018 TorsionAtomIndices,
1019 TorsionAngle,
1020 EnergyMethod,
1021 AngleNotObserved,
1022 Energy,
1023 EnergyLowerBound,
1024 EnergyUpperBound,
1025 )
1026
1027 def _IsEnergyValueValid(self, Value):
1028 """Check for valid energy value."""
1029
1030 return False if (Value is None or math.isnan(Value) or math.isinf(Value)) else True
1031
1032 def _SetupStrainEnergyForRotatableBond(self, TorsionRuleAnglesInfo, TorsionAngle):
1033 """Setup strain energy for rotatable bond."""
1034
1035 if TorsionRuleAnglesInfo["EnergyMethodExact"]:
1036 return self._SetupStrainEnergyForRotatableBondByExactMethod(TorsionRuleAnglesInfo, TorsionAngle)
1037 elif TorsionRuleAnglesInfo["EnergyMethodApproximate"]:
1038 return self._SetupStrainEnergyForRotatableBondByApproximateMethod(TorsionRuleAnglesInfo, TorsionAngle)
1039 else:
1040 EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound = [None, None, None, None, None]
1041 return (EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound)
1042
1043 def _SetupStrainEnergyForRotatableBondByExactMethod(self, TorsionRuleAnglesInfo, TorsionAngle):
1044 """Setup strain energy for rotatable bond by exact method."""
1045
1046 EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound = ["Exact", None, None, None, None]
1047
1048 # Map angle to energy bin numbers...
1049 BinNum = math.ceil(TorsionAngle / 10) + 17
1050 PreviousBinNum = (BinNum + 35) % 36
1051
1052 # Bin angle from -170 to 180 by the right end points...
1053 BinAngleRightSide = (BinNum - 17) * 10
1054
1055 # Angle offset towards the left of the bin from the right end point...
1056 AngleOffset = TorsionAngle - BinAngleRightSide
1057
1058 BinEnergy = TorsionRuleAnglesInfo["HistogramEnergy"][BinNum]
1059 PreviousBinEnergy = TorsionRuleAnglesInfo["HistogramEnergy"][PreviousBinNum]
1060 Energy = BinEnergy + (BinEnergy - PreviousBinEnergy) / 10.0 * AngleOffset
1061
1062 BinEnergyLowerBound = TorsionRuleAnglesInfo["HistogramEnergyLowerBound"][BinNum]
1063 PreviousBinEnergyLowerBound = TorsionRuleAnglesInfo["HistogramEnergyLowerBound"][PreviousBinNum]
1064 EnergyLowerBound = (
1065 BinEnergyLowerBound + (BinEnergyLowerBound - PreviousBinEnergyLowerBound) / 10.0 * AngleOffset
1066 )
1067
1068 BinEnergyUpperBound = TorsionRuleAnglesInfo["HistogramEnergyUpperBound"][BinNum]
1069 PreviousBinEnergyUpperBound = TorsionRuleAnglesInfo["HistogramEnergyUpperBound"][PreviousBinNum]
1070 EnergyUpperBound = (
1071 BinEnergyUpperBound + (BinEnergyUpperBound - PreviousBinEnergyUpperBound) / 10.0 * AngleOffset
1072 )
1073
1074 return (EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound)
1075
1076 def _SetupStrainEnergyForRotatableBondByApproximateMethod(self, TorsionRuleAnglesInfo, TorsionAngle):
1077 """Setup strain energy for rotatable bond by approximate method."""
1078
1079 EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound = [
1080 "Approximate",
1081 True,
1082 None,
1083 None,
1084 None,
1085 ]
1086
1087 for AngleID in TorsionRuleAnglesInfo["IDs"]:
1088 Tolerance2 = TorsionRuleAnglesInfo["Tolerance2"][AngleID]
1089 Theta0 = TorsionRuleAnglesInfo["Theta0"][AngleID]
1090
1091 AngleDiff = TorsionAlertsUtil.CalculateTorsionAngleDifference(TorsionAngle, Theta0)
1092 if abs(AngleDiff) <= Tolerance2:
1093 Beta1 = TorsionRuleAnglesInfo["Beta1"][AngleID]
1094 Beta2 = TorsionRuleAnglesInfo["Beta2"][AngleID]
1095
1096 Energy = Beta1 * (AngleDiff**2) + Beta2 * (AngleDiff**4)
1097
1098 # Estimates of lower and upper bound are not available for
1099 # approximate method...
1100 EnergyLowerBound = Energy
1101 EnergyUpperBound = Energy
1102
1103 AngleNotObserved = False
1104
1105 break
1106
1107 return (EnergyMethod, AngleNotObserved, Energy, EnergyLowerBound, EnergyUpperBound)
1108
1109 def _MatchTorsionRuleToRotatableBond(self, Mol, RotBondAtomIndices, ElementNode):
1110 """Retrieve matched torsions for torsion rule matched to rotatable bond."""
1111
1112 # Get torsion matches...
1113 TorsionMatches = self._GetMatchesForTorsionRule(Mol, ElementNode)
1114 if TorsionMatches is None or len(TorsionMatches) == 0:
1115 return (None, None)
1116
1117 # Identify all torsion matches corresponding to central atoms in RotBondAtomIndices...
1118 RotBondAtomIndex1, RotBondAtomIndex2 = RotBondAtomIndices
1119 RotBondTorsionMatches, RotBondTorsionAngles = [None] * 2
1120
1121 for TorsionMatch in TorsionMatches:
1122 CentralAtomIndex1 = TorsionMatch[1]
1123 CentralAtomIndex2 = TorsionMatch[2]
1124
1125 if (CentralAtomIndex1 == RotBondAtomIndex1 and CentralAtomIndex2 == RotBondAtomIndex2) or (
1126 CentralAtomIndex1 == RotBondAtomIndex2 and CentralAtomIndex2 == RotBondAtomIndex1
1127 ):
1128 TorsionAngle = self._CalculateTorsionAngle(Mol, TorsionMatch)
1129 if RotBondTorsionMatches is None:
1130 RotBondTorsionMatches = []
1131 RotBondTorsionAngles = []
1132 RotBondTorsionMatches.append(TorsionMatch)
1133 RotBondTorsionAngles.append(TorsionAngle)
1134
1135 return (RotBondTorsionMatches, RotBondTorsionAngles)
1136
1137 def _CalculateTorsionAngle(self, Mol, TorsionMatch):
1138 """Calculate torsion angle."""
1139
1140 # Calculate torsion angle using torsion atom indices..
1141 MolConf = Mol.GetConformer(0)
1142 TorsionAngle = rdMolTransforms.GetDihedralDeg(
1143 MolConf, TorsionMatch[0], TorsionMatch[1], TorsionMatch[2], TorsionMatch[3]
1144 )
1145 TorsionAngle = round(TorsionAngle, 2)
1146
1147 return TorsionAngle
1148
1149 def _GetMatchesForTorsionRule(self, Mol, ElementNode):
1150 """Get matches for torsion rule."""
1151
1152 # Match torsions...
1153 TorsionMatches = self._GetSubstructureMatchesForTorsionRule(Mol, ElementNode)
1154
1155 if TorsionMatches is None or len(TorsionMatches) == 0:
1156 return TorsionMatches
1157
1158 # Filter torsion matches...
1159 FiltertedTorsionMatches = []
1160 for TorsionMatch in TorsionMatches:
1161 if len(TorsionMatch) != 4:
1162 continue
1163
1164 # Ignore matches containing hydrogen atoms as first or last atom...
1165 if Mol.GetAtomWithIdx(TorsionMatch[0]).GetAtomicNum() == 1:
1166 continue
1167 if Mol.GetAtomWithIdx(TorsionMatch[3]).GetAtomicNum() == 1:
1168 continue
1169
1170 FiltertedTorsionMatches.append(TorsionMatch)
1171
1172 return FiltertedTorsionMatches
1173
1174 def _GetSubstructureMatchesForTorsionRule(self, Mol, ElementNode):
1175 """Get substructure matches for a torsion rule."""
1176
1177 # Setup torsion rule SMARTS pattern mol....
1178 TorsionRuleNodeID = ElementNode.get("NodeID")
1179 TorsionSMARTSPattern = ElementNode.get("smarts")
1180 TorsionPatternMol = TorsionAlertsUtil.SetupTorsionRuleElementPatternMol(
1181 self.TorsionLibraryInfo, ElementNode, TorsionRuleNodeID, TorsionSMARTSPattern
1182 )
1183 if TorsionPatternMol is None:
1184 return None
1185
1186 # Match torsions...
1187 TorsionMatches = TorsionAlertsUtil.FilterSubstructureMatchesByAtomMapNumbers(
1188 Mol, TorsionPatternMol, Mol.GetSubstructMatches(TorsionPatternMol, useChirality=False)
1189 )
1190
1191 return TorsionMatches
