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