1 #!/bin/env python
2 #
3 # File: RDKitPerformTorsionScan.py
4 # Author: Manish Sud <msud@san.rr.com>
5 #
6 # Copyright (C) 2024 Manish Sud. All rights reserved.
7 #
8 # The functionality available in this script is implemented using RDKit, an
9 # open source toolkit for cheminformatics developed by Greg Landrum.
10 #
11 # This file is part of MayaChemTools.
12 #
13 # MayaChemTools is free software; you can redistribute it and/or modify it under
14 # the terms of the GNU Lesser General Public License as published by the Free
15 # Software Foundation; either version 3 of the License, or (at your option) any
16 # later version.
17 #
18 # MayaChemTools is distributed in the hope that it will be useful, but without
19 # any warranty; without even the implied warranty of merchantability of fitness
20 # for a particular purpose. See the GNU Lesser General Public License for more
21 # details.
22 #
23 # You should have received a copy of the GNU Lesser General Public License
24 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
25 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
26 # Boston, MA, 02111-1307, USA.
27 #
28
29 from __future__ import print_function
30
31 # Add local python path to the global path and import standard library modules...
32 import os
33 import sys; sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
34 import time
35 import re
36 import glob
37 import multiprocessing as mp
38
39 import matplotlib.pyplot as plt
40 import seaborn as sns
41
42 # RDKit imports...
43 try:
44 from rdkit import rdBase
45 from rdkit import Chem
46 from rdkit.Chem import AllChem
47 from rdkit.Chem import rdMolTransforms
48 except ImportError as ErrMsg:
49 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
50 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
51 sys.exit(1)
52
53 # MayaChemTools imports...
54 try:
55 from docopt import docopt
56 import MiscUtil
57 import RDKitUtil
58 except ImportError as ErrMsg:
59 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
60 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
61 sys.exit(1)
62
63 ScriptName = os.path.basename(sys.argv[0])
64 Options = {}
65 OptionsInfo = {}
66
67 def main():
68 """Start execution of the script."""
69
70 MiscUtil.PrintInfo("\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n" % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime()))
71
72 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
73
74 # Retrieve command line arguments and options...
75 RetrieveOptions()
76
77 # Process and validate command line arguments and options...
78 ProcessOptions()
79
80 # Perform actions required by the script...
81 PerformTorsionScan()
82
83 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
84 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
85
86 def PerformTorsionScan():
87 """Perform torsion scan."""
88
89 # Setup a molecule reader for input file...
90 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
91 OptionsInfo["InfileParams"]["AllowEmptyMols"] = True
92 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
93
94 PlotExt = OptionsInfo["OutPlotParams"]["OutExt"]
95 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
96 MiscUtil.PrintInfo("Generating output files %s_*.sdf, %s_*Torsion*Match*.sdf, %s_*Torsion*Match*Energies.csv, %s_*Torsion*Match*Plot.%s..." % (FileName, FileName, FileName, FileName, PlotExt))
97
98 MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount = ProcessMolecules(Mols)
99
100 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
101 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
102 MiscUtil.PrintInfo("Number of molecules failed during initial minimization: %d" % MinimizationFailedCount)
103 MiscUtil.PrintInfo("Number of molecules without any matched torsions: %d" % TorsionsMissingCount)
104 MiscUtil.PrintInfo("Number of molecules failed during torsion scan: %d" % TorsionsScanFailedCount)
105 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + TorsionsMissingCount + MinimizationFailedCount + TorsionsScanFailedCount))
106
107 def ProcessMolecules(Mols):
108 """Process molecules to perform torsion scan."""
109
110 if OptionsInfo["MPMode"]:
111 return ProcessMoleculesUsingMultipleProcesses(Mols)
112 else:
113 return ProcessMoleculesUsingSingleProcess( Mols)
114
115 def ProcessMoleculesUsingSingleProcess(Mols):
116 """Process molecules to perform torsion scan using a single process."""
117
118 MolInfoText = "first molecule"
119 if not OptionsInfo["FirstMolMode"]:
120 MolInfoText = "all molecules"
121
122 if OptionsInfo["TorsionMinimize"]:
123 MiscUtil.PrintInfo("\nPerforming torsion scan on %s by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
124 else:
125 MiscUtil.PrintInfo("\nPerforming torsion scan on %s by skipping generation of conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
126
127 SetupTorsionsPatternsInfo()
128
129 (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
130
131 for Mol in Mols:
132 MolCount += 1
133
134 if OptionsInfo["FirstMolMode"] and MolCount > 1:
135 MolCount -= 1
136 break
137
138 if Mol is None:
139 continue
140
141 if RDKitUtil.IsMolEmpty(Mol):
142 if not OptionsInfo["QuietMode"]:
143 MolName = RDKitUtil.GetMolName(Mol, MolCount)
144 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
145 continue
146 ValidMolCount += 1
147
148 Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(Mol, MolCount)
149
150 if not MinimizationCalcStatus:
151 MinimizationFailedCount += 1
152 continue
153
154 if not TorsionsMatchStatus:
155 TorsionsMissingCount += 1
156 continue
157
158 if not TorsionsScanCalcStatus:
159 TorsionsScanFailedCount += 1
160 continue
161
162 return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
163
164 def ProcessMoleculesUsingMultipleProcesses(Mols):
165 """Process and minimize molecules using multiprocessing."""
166
167 if OptionsInfo["MPLevelTorsionAnglesMode"]:
168 return ProcessMoleculesUsingMultipleProcessesAtTorsionAnglesLevel(Mols)
169 elif OptionsInfo["MPLevelMoleculesMode"]:
170 return ProcessMoleculesUsingMultipleProcessesAtMoleculesLevel(Mols)
171 else:
172 MiscUtil.PrintError("The value, %s, option \"--mpLevel\" is not supported." % (OptionsInfo["MPLevel"]))
173
174 def ProcessMoleculesUsingMultipleProcessesAtMoleculesLevel(Mols):
175 """Process molecules to perform torsion scan using multiprocessing at molecules level."""
176
177 MolInfoText = "first molecule"
178 if not OptionsInfo["FirstMolMode"]:
179 MolInfoText = "all molecules"
180
181 if OptionsInfo["TorsionMinimize"]:
182 MiscUtil.PrintInfo("\nPerforming torsion scan on %s using multiprocessing at molecules level by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
183 else:
184 MiscUtil.PrintInfo("\nPerforming torsion scan %s using multiprocessing at molecules level by skipping generation of conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
185
186 MPParams = OptionsInfo["MPParams"]
187
188 # Setup data for initializing a worker process...
189 MiscUtil.PrintInfo("\nEncoding options info...")
190
191 InitializeWorkerProcessArgs = (MiscUtil.ObjectToBase64EncodedString(Options), MiscUtil.ObjectToBase64EncodedString(OptionsInfo))
192
193 if OptionsInfo["FirstMolMode"]:
194 Mol = Mols[0]
195 Mols = [Mol]
196
197 # Setup a encoded mols data iterable for a worker process...
198 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
199
200 # Setup process pool along with data initialization for each process...
201 MiscUtil.PrintInfo("\nConfiguring multiprocessing using %s method..." % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()"))
202 MiscUtil.PrintInfo("NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" % (MPParams["NumProcesses"], MPParams["InputDataMode"], ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"])))
203
204 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
205
206 # Start processing...
207 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
208 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
209 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
210 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
211 else:
212 MiscUtil.PrintError("The value, %s, specified for \"--inputDataMode\" is not supported." % (MPParams["InputDataMode"]))
213
214 (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
215
216 for Result in Results:
217 MolCount += 1
218
219 MolIndex, EncodedMol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = Result
220
221 if EncodedMol is None:
222 continue
223 ValidMolCount += 1
224
225 if not MinimizationCalcStatus:
226 MinimizationFailedCount += 1
227 continue
228
229 if not TorsionsMatchStatus:
230 TorsionsMissingCount += 1
231 continue
232
233 if not TorsionsScanCalcStatus:
234 TorsionsScanFailedCount += 1
235 continue
236
237 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
238
239 return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
240
241 def InitializeWorkerProcess(*EncodedArgs):
242 """Initialize data for a worker process."""
243
244 global Options, OptionsInfo
245
246 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
247
248 # Decode Options and OptionInfo...
249 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
250 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
251
252 # Initialize torsion patterns info...
253 SetupTorsionsPatternsInfo()
254
255 def WorkerProcess(EncodedMolInfo):
256 """Process data for a worker process."""
257
258 MolIndex, EncodedMol = EncodedMolInfo
259
260 (MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus) = [False] * 3
261
262 if EncodedMol is None:
263 return [MolIndex, None, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
264
265 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
266 if RDKitUtil.IsMolEmpty(Mol):
267 if not OptionsInfo["QuietMode"]:
268 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1))
269 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
270 return [MolIndex, None, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
271
272 Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(Mol, (MolIndex + 1))
273
274 return [MolIndex, RDKitUtil.MolToBase64EncodedMolString(Mol, PropertyPickleFlags = Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps), MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
275
276 def ProcessMoleculesUsingMultipleProcessesAtTorsionAnglesLevel(Mols):
277 """Process molecules to perform torsion scan using multiprocessing at torsion angles level."""
278
279 MolInfoText = "first molecule"
280 if not OptionsInfo["FirstMolMode"]:
281 MolInfoText = "all molecules"
282
283 if OptionsInfo["TorsionMinimize"]:
284 MiscUtil.PrintInfo("\nPerforming torsion scan on %s using multiprocessing at torsion angles level by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
285 else:
286 MiscUtil.PrintInfo("\nPerforming torsion scan %s using multiprocessing at torsion angles level by skipping generation of conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
287
288 SetupTorsionsPatternsInfo()
289
290 (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
291
292 for Mol in Mols:
293 MolCount += 1
294
295 if OptionsInfo["FirstMolMode"] and MolCount > 1:
296 MolCount -= 1
297 break
298
299 if Mol is None:
300 continue
301
302 if RDKitUtil.IsMolEmpty(Mol):
303 if not OptionsInfo["QuietMode"]:
304 MolName = RDKitUtil.GetMolName(Mol, MolCount)
305 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
306 continue
307 ValidMolCount += 1
308
309 Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(Mol, MolCount, UseMultiProcessingAtTorsionAnglesLevel = True)
310
311 if not MinimizationCalcStatus:
312 MinimizationFailedCount += 1
313 continue
314
315 if not TorsionsMatchStatus:
316 TorsionsMissingCount += 1
317 continue
318
319 if not TorsionsScanCalcStatus:
320 TorsionsScanFailedCount += 1
321 continue
322
323 return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
324
325 def ScanSingleTorsionInMolUsingMultipleProcessesAtTorsionAnglesLevel(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum):
326 """Perform torsion scan for a molecule using multiple processses at torsion angles
327 level along with constrained energy minimization.
328 """
329
330 if OptionsInfo["MPLevelMoleculesMode"]:
331 MiscUtil.PrintError("Single torison scanning for a molecule is not allowed in multiprocessing mode at molecules level.\n")
332
333 Mols, Angles = SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum)
334
335 MPParams = OptionsInfo["MPParams"]
336
337 # Setup data for initializing a worker process...
338 MiscUtil.PrintInfo("\nEncoding options info...")
339
340 # Track and avoid encoding TorsionsPatternsInfo as it contains RDKit molecule object...
341 TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
342 OptionsInfo["TorsionsPatternsInfo"] = None
343
344 InitializeWorkerProcessArgs = (MiscUtil.ObjectToBase64EncodedString(Options), MiscUtil.ObjectToBase64EncodedString(OptionsInfo))
345
346 # Restore TorsionsPatternsInfo...
347 OptionsInfo["TorsionsPatternsInfo"] = TorsionsPatternsInfo
348
349 # Setup a encoded mols data iterable for a worker process...
350 WorkerProcessDataIterable = GenerateBase64EncodedMolStringsWithTorsionScanInfo(Mol, (MolNum -1), Mols, Angles, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches)
351
352 # Setup process pool along with data initialization for each process...
353 MiscUtil.PrintInfo("\nConfiguring multiprocessing using %s method..." % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()"))
354 MiscUtil.PrintInfo("NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" % (MPParams["NumProcesses"], MPParams["InputDataMode"], ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"])))
355
356 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeTorsionAngleWorkerProcess, InitializeWorkerProcessArgs)
357
358 # Start processing...
359 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
360 Results = ProcessPool.imap(TorsionAngleWorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
361 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
362 Results = ProcessPool.map(TorsionAngleWorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
363 else:
364 MiscUtil.PrintError("The value, %s, specified for \"--inputDataMode\" is not supported." % (MPParams["InputDataMode"]))
365
366 TorsionMols = []
367 TorsionEnergies = []
368 TorsionAngles = []
369
370 for Result in Results:
371 EncodedTorsionMol, CalcStatus, Angle, Energy = Result
372
373 if not CalcStatus:
374 return (Mol, False, None, None, None)
375
376 if EncodedTorsionMol is None:
377 return (Mol, False, None, None, None)
378 TorsionMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionMol)
379
380 if OptionsInfo["RemoveHydrogens"]:
381 TorsionMol = Chem.RemoveHs(TorsionMol)
382
383 TorsionMols.append(TorsionMol)
384 TorsionEnergies.append(Energy)
385 TorsionAngles.append(Angle)
386
387 return (Mol, True, TorsionMols, TorsionEnergies, TorsionAngles)
388
389 def InitializeTorsionAngleWorkerProcess(*EncodedArgs):
390 """Initialize data for a worker process."""
391
392 global Options, OptionsInfo
393
394 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
395
396 # Decode Options and OptionInfo...
397 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
398 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
399
400 # Initialize torsion patterns info...
401 SetupTorsionsPatternsInfo()
402
403 def TorsionAngleWorkerProcess(EncodedMolInfo):
404 """Process data for a worker process."""
405
406 MolIndex, EncodedMol, EncodedTorsionMol, TorsionAngle, TorsionID, TorsionPattern, EncodedTorsionPatternMol, TorsionMatches = EncodedMolInfo
407
408 if EncodedMol is None or EncodedTorsionMol is None or EncodedTorsionPatternMol is None:
409 return (None, False, None, None)
410
411 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
412 TorsionMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionMol)
413 TorsionPatternMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionPatternMol)
414
415 TorsionMol, CalcStatus, Energy = MinimizeCalculateEnergyForTorsionMol(Mol, TorsionMol, TorsionAngle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, (MolIndex + 1))
416
417 return (RDKitUtil.MolToBase64EncodedMolString(TorsionMol, PropertyPickleFlags = Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps), CalcStatus, TorsionAngle, Energy)
418
419 def GenerateBase64EncodedMolStringsWithTorsionScanInfo(Mol, MolIndex, TorsionMols, TorsionAngles, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, PropertyPickleFlags = Chem.PropertyPickleOptions.AllProps):
420 """Set up an iterator for generating base64 encoded molecule string for
421 a torsion in a molecule along with appropriate trosion scan information.
422 """
423
424 for Index, TorsionMol in enumerate(TorsionMols):
425 yield [MolIndex, None, None, TorsionAngles[Index], TorsionID, TorsionPattern, None, TorsionMatches] if (Mol is None or TorsionMol is None) else [MolIndex, RDKitUtil.MolToBase64EncodedMolString(Mol, PropertyPickleFlags), RDKitUtil.MolToBase64EncodedMolString(TorsionMol, PropertyPickleFlags), TorsionAngles[Index], TorsionID, TorsionPattern, RDKitUtil.MolToBase64EncodedMolString(TorsionPatternMol, PropertyPickleFlags), TorsionMatches]
426
427 def PerformMinimizationAndTorsionScan(Mol, MolNum, UseMultiProcessingAtTorsionAnglesLevel = False):
428 """Perform minimization and torsions scan."""
429
430 if not OptionsInfo["QuietMode"]:
431 MiscUtil.PrintInfo("\nProcessing molecule %s..." % (RDKitUtil.GetMolName(Mol, MolNum)))
432
433 Mol = AddHydrogens(Mol)
434
435 if not OptionsInfo["Infile3D"]:
436 Mol, MinimizationCalcStatus = MinimizeMolecule(Mol, MolNum)
437 if not MinimizationCalcStatus:
438 return (Mol, False, False, False)
439
440 TorsionsMolInfo = SetupTorsionsMolInfo(Mol, MolNum)
441 if TorsionsMolInfo["NumOfMatches"] == 0:
442 return (Mol, True, False, False)
443
444 Mol, ScanCalcStatus = ScanAllTorsionsInMol(Mol, TorsionsMolInfo, MolNum, UseMultiProcessingAtTorsionAnglesLevel)
445 if not ScanCalcStatus:
446 return (Mol, True, True, False)
447
448 return (Mol, True, True, True)
449
450 def ScanAllTorsionsInMol(Mol, TorsionsMolInfo, MolNum, UseMultiProcessingAtTorsionAnglesLevel = False):
451 """Perform scans on all torsions in a molecule."""
452
453 if TorsionsMolInfo["NumOfMatches"] == 0:
454 return Mol, True
455
456 MolName = RDKitUtil.GetMolName(Mol, MolNum)
457
458 FirstTorsionMode = OptionsInfo["FirstTorsionMode"]
459 TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
460
461 TorsionPatternCount, TorsionScanCount, TorsionMatchCount = [0] * 3
462 TorsionMaxMatches = OptionsInfo["TorsionMaxMatches"]
463
464 for TorsionID in TorsionsPatternsInfo["IDs"]:
465 TorsionPatternCount += 1
466 TorsionPattern = TorsionsPatternsInfo["Pattern"][TorsionID]
467 TorsionPatternMol = TorsionsPatternsInfo["Mol"][TorsionID]
468
469 TorsionsMatches = TorsionsMolInfo["Matches"][TorsionID]
470
471 if TorsionsMatches is None:
472 continue
473
474 if FirstTorsionMode and TorsionPatternCount > 1:
475 if not OptionsInfo["QuietMode"]:
476 MiscUtil.PrintWarning("Already scaned first torsion pattern, \"%s\" for molecule %s during \"%s\" value of \"--modeTorsions\" option . Abandoning torsion scan...\n" % (TorsionPattern, MolName, OptionsInfo["ModeTorsions"]))
477 break
478
479 for Index, TorsionMatches in enumerate(TorsionsMatches):
480 TorsionMatchNum = Index + 1
481 TorsionMatchCount += 1
482
483 if TorsionMatchCount > TorsionMaxMatches:
484 if not OptionsInfo["QuietMode"]:
485 MiscUtil.PrintWarning("Already scaned a maximum of %s torsion matches for molecule %s specified by \"--torsionMaxMatches\" option. Abandoning torsion scan...\n" % (TorsionMaxMatches, MolName))
486 break
487
488 TmpMol, TorsionScanStatus, TorsionMols, TorsionEnergies, TorsionAngles = ScanSingleTorsionInMol(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, TorsionMatchNum, MolNum, UseMultiProcessingAtTorsionAnglesLevel)
489 if not TorsionScanStatus:
490 continue
491
492 TorsionScanCount += 1
493 GenerateOutputFiles(Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
494
495 if TorsionMatchCount > TorsionMaxMatches:
496 break
497
498 if OptionsInfo["RemoveHydrogens"]:
499 Mol = Chem.RemoveHs(Mol)
500
501 if TorsionScanCount:
502 GenerateStartingTorsionScanStructureOutfile(Mol, MolNum)
503
504 Status = True if TorsionScanCount else False
505
506 return (Mol, Status)
507
508 def ScanSingleTorsionInMol(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, TorsionMatchNum, MolNum, UseMultiProcessingAtTorsionAnglesLevel):
509 """Perform torsion scan for a molecule along with constrained energy minimization."""
510
511 if not OptionsInfo["QuietMode"]:
512 MiscUtil.PrintInfo("Processing torsion pattern, %s, match number, %s, in molecule %s..." % (TorsionPattern, TorsionMatchNum, RDKitUtil.GetMolName(Mol, MolNum)))
513
514 if UseMultiProcessingAtTorsionAnglesLevel:
515 return ScanSingleTorsionInMolUsingMultipleProcessesAtTorsionAnglesLevel(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum)
516 else:
517 return ScanSingleTorsionInMolUsingSingleProcess(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum)
518
519 def ScanSingleTorsionInMolUsingSingleProcess(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum):
520 """Perform torsion scan for a molecule using single processs along with constrained
521 energy minimization."""
522
523 TorsionMols = []
524 TorsionEnergies = []
525 TorsionAngles = []
526
527 Mols, Angles = SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum)
528
529 for Index, Angle in enumerate(Angles):
530 TorsionMol = Mols[Index]
531 TorsionMol, CalcStatus, Energy = MinimizeCalculateEnergyForTorsionMol(Mol, TorsionMol, Angle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum)
532
533 if not CalcStatus:
534 return (Mol, False, None, None, None)
535
536 if OptionsInfo["RemoveHydrogens"]:
537 TorsionMol = Chem.RemoveHs(TorsionMol)
538
539 TorsionMols.append(TorsionMol)
540 TorsionEnergies.append(Energy)
541 TorsionAngles.append(Angle)
542
543 return (Mol, True, TorsionMols, TorsionEnergies, TorsionAngles)
544
545 def SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum = None):
546 """Setup molecules corresponding to all torsion angles in a molecule."""
547
548 AtomIndex1, AtomIndex2, AtomIndex3, AtomIndex4 = TorsionMatches
549
550 TorsionMols = []
551 TorsionAngles = OptionsInfo["TorsionAngles"]
552
553 for Angle in TorsionAngles:
554 TorsionMol = Chem.Mol(Mol)
555 TorsionMolConf = TorsionMol.GetConformer(0)
556
557 rdMolTransforms.SetDihedralDeg(TorsionMolConf, AtomIndex1, AtomIndex2, AtomIndex3, AtomIndex4, Angle)
558 TorsionMols.append(TorsionMol)
559
560 return (TorsionMols, TorsionAngles)
561
562 def MinimizeCalculateEnergyForTorsionMol(Mol, TorsionMol, TorsionAngle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum):
563 """"Calculate energy of a torsion molecule by performing an optional constrained
564 energy minimzation.
565 """
566
567 if OptionsInfo["TorsionMinimize"]:
568 # Perform constrained minimization...
569 TorsionMatchesMol = RDKitUtil.MolFromSubstructureMatch(TorsionMol, TorsionPatternMol, TorsionMatches)
570 TorsionMol, CalcStatus, Energy = ConstrainAndMinimizeMolecule(TorsionMol, TorsionMatchesMol, TorsionMatches, MolNum)
571
572 if not CalcStatus:
573 if not OptionsInfo["QuietMode"]:
574 MolName = RDKitUtil.GetMolName(Mol, MolNum)
575 MiscUtil.PrintWarning("Failed to perform constrained minimization for molecule %s with torsion angle set to %s during torsion scan for torsion pattern %s. Abandoning torsion scan..." % (MolName, TorsionAngle, TorsionPattern))
576 return (TorsionMol, False, None)
577 else:
578 # Calculate energy...
579 CalcStatus, Energy = GetEnergy(TorsionMol)
580 if not CalcStatus:
581 if not OptionsInfo["QuietMode"]:
582 MolName = RDKitUtil.GetMolName(Mol, MolNum)
583 MiscUtil.PrintWarning("Failed to retrieve calculated energy for molecule %s with torsion angle set to %s during torsion scan for torsion pattern %s. Abandoning torsion scan..." % (MolName, TorsionAngle, TorsionPattern))
584 return (TorsionMol, False, None)
585
586 return (TorsionMol, CalcStatus, Energy)
587
588 def SetupTorsionsMolInfo(Mol, MolNum = None):
589 """Setup torsions info for a molecule."""
590
591 TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
592
593 # Initialize...
594 TorsionsMolInfo = {}
595 TorsionsMolInfo["IDs"] = []
596 TorsionsMolInfo["NumOfMatches"] = 0
597 TorsionsMolInfo["Matches"] = {}
598 for TorsionID in TorsionsPatternsInfo["IDs"]:
599 TorsionsMolInfo["IDs"].append(TorsionID)
600 TorsionsMolInfo["Matches"][TorsionID] = None
601
602 MolName = RDKitUtil.GetMolName(Mol, MolNum)
603 UseChirality = OptionsInfo["UseChirality"]
604
605 for TorsionID in TorsionsPatternsInfo["IDs"]:
606 # Match torsions..
607 TorsionPattern = TorsionsPatternsInfo["Pattern"][TorsionID]
608 TorsionPatternMol = TorsionsPatternsInfo["Mol"][TorsionID]
609 TorsionsMatches = RDKitUtil.FilterSubstructureMatchesByAtomMapNumbers(Mol, TorsionPatternMol, Mol.GetSubstructMatches(TorsionPatternMol, useChirality = UseChirality))
610
611 # Validate tosion matches...
612 ValidTorsionsMatches = []
613 for Index, TorsionMatch in enumerate(TorsionsMatches):
614 if len(TorsionMatch) != 4:
615 if not OptionsInfo["QuietMode"]:
616 MiscUtil.PrintWarning("Ignoring invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: It must match exactly 4 atoms." % (TorsionMatch, TorsionPattern, MolName))
617 continue
618
619 if not RDKitUtil.AreAtomIndicesSequentiallyConnected(Mol, TorsionMatch):
620 if not OptionsInfo["QuietMode"]:
621 MiscUtil.PrintInfo("")
622 MiscUtil.PrintWarning("Invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices must be sequentially connected." % (TorsionMatch, TorsionPattern, MolName))
623 MiscUtil.PrintWarning("Reordering matched atom indices in a sequentially connected manner...")
624
625 Status, ReorderdTorsionMatch = RDKitUtil.ReorderAtomIndicesInSequentiallyConnectedManner(Mol, TorsionMatch)
626 if Status:
627 TorsionMatch = ReorderdTorsionMatch
628 if not OptionsInfo["QuietMode"]:
629 MiscUtil.PrintWarning("Successfully reordered torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices are now sequentially connected." % (TorsionMatch, TorsionPattern, MolName))
630 else:
631 if not OptionsInfo["QuietMode"]:
632 MiscUtil.PrintWarning("Ignoring torsion match. Failed to reorder torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices are not sequentially connected." % (TorsionMatch, TorsionPattern, MolName))
633 continue
634
635 Bond = Mol.GetBondBetweenAtoms(TorsionMatch[1], TorsionMatch[2])
636 if Bond.IsInRing():
637 if not OptionsInfo["QuietMode"]:
638 MiscUtil.PrintWarning("Ignoring invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices, %s and %s, are not allowed to be in a ring." % (TorsionMatch, TorsionPattern, MolName, TorsionMatch[1], TorsionMatch[2]))
639 continue
640
641 # Filter matched torsions...
642 if OptionsInfo["FilterTorsionsByAtomIndicesMode"]:
643 InvalidAtomIndices = []
644 for AtomIndex in TorsionMatch:
645 if AtomIndex not in OptionsInfo["TorsionsFilterByAtomIndicesList"]:
646 InvalidAtomIndices.append(AtomIndex)
647 if len(InvalidAtomIndices):
648 if not OptionsInfo["QuietMode"]:
649 MiscUtil.PrintWarning("Ignoring invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices, %s, must be present in the list, %s, specified using option \"--torsionsFilterbyAtomIndices\"." % (TorsionMatch, TorsionPattern, MolName, InvalidAtomIndices, OptionsInfo["TorsionsFilterByAtomIndicesList"]))
650 continue
651
652 ValidTorsionsMatches.append(TorsionMatch)
653
654 # Track valid matches...
655 if len(ValidTorsionsMatches):
656 TorsionsMolInfo["NumOfMatches"] += len(ValidTorsionsMatches)
657 TorsionsMolInfo["Matches"][TorsionID] = ValidTorsionsMatches
658
659 if TorsionsMolInfo["NumOfMatches"] == 0:
660 if not OptionsInfo["QuietMode"]:
661 MiscUtil.PrintWarning("Failed to match any torsions in molecule %s" % (MolName))
662
663 return TorsionsMolInfo
664
665 def SetupTorsionsPatternsInfo():
666 """Setup torsions patterns info."""
667
668 TorsionsPatternsInfo = {}
669 TorsionsPatternsInfo["IDs"] = []
670 TorsionsPatternsInfo["Pattern"] = {}
671 TorsionsPatternsInfo["Mol"] = {}
672
673 TorsionID = 0
674 for TorsionPattern in OptionsInfo["TorsionPatternsList"]:
675 TorsionID += 1
676
677 TorsionMol = Chem.MolFromSmarts(TorsionPattern)
678 if TorsionMol is None:
679 MiscUtil.PrintError("Failed to create torsion pattern molecule. The torsion SMILES/SMARTS pattern, \"%s\", specified using \"-t, --torsions\" option is not valid." % (TorsionPattern))
680
681 TorsionsPatternsInfo["IDs"].append(TorsionID)
682 TorsionsPatternsInfo["Pattern"][TorsionID] = TorsionPattern
683 TorsionsPatternsInfo["Mol"][TorsionID] = TorsionMol
684
685 OptionsInfo["TorsionsPatternsInfo"] = TorsionsPatternsInfo
686
687 def MinimizeMolecule(Mol, MolNum = None):
688 """Generate and minimize conformers for a molecule to get the lowest energy conformer."""
689
690 if not OptionsInfo["QuietMode"]:
691 MiscUtil.PrintInfo("Minimizing molecule %s..." % (RDKitUtil.GetMolName(Mol, MolNum)))
692
693 ConfIDs = EmbedMolecule(Mol, MolNum)
694 if not len(ConfIDs):
695 if not OptionsInfo["QuietMode"]:
696 MolName = RDKitUtil.GetMolName(Mol, MolNum)
697 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s: Embedding failed...\n" % MolName)
698 return (Mol, False)
699
700 CalcEnergyMap = {}
701 for ConfID in ConfIDs:
702 try:
703 if OptionsInfo["UseUFF"]:
704 Status = AllChem.UFFOptimizeMolecule(Mol, confId = ConfID, maxIters = OptionsInfo["MaxIters"])
705 elif OptionsInfo["UseMMFF"]:
706 Status = AllChem.MMFFOptimizeMolecule(Mol, confId = ConfID, maxIters = OptionsInfo["MaxIters"], mmffVariant = OptionsInfo["MMFFVariant"])
707 else:
708 MiscUtil.PrintError("Minimization couldn't be performed: Specified forcefield, %s, is not supported" % OptionsInfo["ForceField"])
709 except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
710 if not OptionsInfo["QuietMode"]:
711 MolName = RDKitUtil.GetMolName(Mol, MolNum)
712 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s:\n%s\n" % (MolName, ErrMsg))
713 return (Mol, False)
714
715 EnergyStatus, Energy = GetEnergy(Mol, ConfID)
716 if not EnergyStatus:
717 if not OptionsInfo["QuietMode"]:
718 MolName = RDKitUtil.GetMolName(Mol, MolNum)
719 MiscUtil.PrintWarning("Failed to retrieve calculated energy for conformation number %d of molecule %s. Try again after removing any salts or cleaing up the molecule...\n" % (ConfID, MolName))
720 return (Mol, False)
721
722 if Status != 0:
723 if not OptionsInfo["QuietMode"]:
724 MolName = RDKitUtil.GetMolName(Mol, MolNum)
725 MiscUtil.PrintWarning("Minimization failed to converge for conformation number %d of molecule %s in %d steps. Try using higher value for \"--maxIters\" option...\n" % (ConfID, MolName, OptionsInfo["MaxIters"]))
726
727 CalcEnergyMap[ConfID] = Energy
728
729 SortedConfIDs = sorted(ConfIDs, key = lambda ConfID: CalcEnergyMap[ConfID])
730 MinEnergyConfID = SortedConfIDs[0]
731
732 for ConfID in [Conf.GetId() for Conf in Mol.GetConformers()]:
733 if ConfID == MinEnergyConfID:
734 continue
735 Mol.RemoveConformer(ConfID)
736
737 # Set ConfID to 0 for MinEnergyConf...
738 Mol.GetConformer(MinEnergyConfID).SetId(0)
739
740 return (Mol, True)
741
742 def ConstrainAndMinimizeMolecule(Mol, RefMolCore, RefMolMatches = None, MolNum = None):
743 """Constrain and Minimize molecule."""
744
745 # Setup forcefield function to use for constrained minimization...
746 ForceFieldFunction = None
747 ForceFieldName = None
748 if OptionsInfo["UseUFF"]:
749 ForceFieldFunction = lambda mol, confId = -1 : AllChem.UFFGetMoleculeForceField(mol, confId = confId)
750 ForeceFieldName = "UFF"
751 else:
752 ForceFieldFunction = lambda mol, confId = -1 : AllChem.MMFFGetMoleculeForceField(mol, AllChem.MMFFGetMoleculeProperties(mol, mmffVariant = OptionsInfo["MMFFVariant"]) , confId = confId)
753 ForeceFieldName = "MMFF"
754
755 if ForceFieldFunction is None:
756 if not OptionsInfo["QuietMode"]:
757 MiscUtil.PrintWarning("Failed to setup forcefield %s for molecule: %s\n" % (ForceFieldName, RDKitUtil.GetMolName(Mol, MolNum)))
758 return (None, False, None)
759
760 MaxConfs = OptionsInfo["MaxConfsTorsion"]
761 EnforceChirality = OptionsInfo["EnforceChirality"]
762 UseExpTorsionAnglePrefs = OptionsInfo["UseExpTorsionAnglePrefs"]
763 ETVersion = OptionsInfo["ETVersion"]
764 UseBasicKnowledge = OptionsInfo["UseBasicKnowledge"]
765 UseTethers = OptionsInfo["UseTethers"]
766
767 CalcEnergyMap = {}
768 MolConfsMap = {}
769 ConfIDs = [ConfID for ConfID in range(0, MaxConfs)]
770
771 for ConfID in ConfIDs:
772 try:
773 MolConf = Chem.Mol(Mol)
774 RDKitUtil.ConstrainAndEmbed(MolConf, RefMolCore, coreMatchesMol = RefMolMatches, useTethers = UseTethers, coreConfId = -1, randomseed = ConfID, getForceField = ForceFieldFunction, enforceChirality = EnforceChirality, useExpTorsionAnglePrefs = UseExpTorsionAnglePrefs, useBasicKnowledge = UseBasicKnowledge, ETversion = ETVersion)
775 except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
776 if not OptionsInfo["QuietMode"]:
777 MolName = RDKitUtil.GetMolName(Mol, MolNum)
778 MiscUtil.PrintWarning("Constrained embedding couldn't be performed for molecule %s:\n%s\n" % (RDKitUtil.GetMolName(Mol, MolNum), ErrMsg))
779 return (None, False, None)
780
781 EnergyStatus, Energy = GetEnergy(MolConf)
782
783 if not EnergyStatus:
784 if not OptionsInfo["QuietMode"]:
785 MolName = RDKitUtil.GetMolName(Mol, MolNum)
786 MiscUtil.PrintWarning("Failed to retrieve calculated energy for conformation number %d of molecule %s. Try again after removing any salts or cleaing up the molecule...\n" % (ConfID, MolName))
787 return (None, False, None)
788
789 CalcEnergyMap[ConfID] = Energy
790 MolConfsMap[ConfID] = MolConf
791
792 SortedConfIDs = sorted(ConfIDs, key = lambda ConfID: CalcEnergyMap[ConfID])
793 MinEnergyConfID = SortedConfIDs[0]
794
795 MinEnergy = CalcEnergyMap[MinEnergyConfID]
796 MinEnergyMolConf = MolConfsMap[MinEnergyConfID]
797
798 MinEnergyMolConf.ClearProp('EmbedRMS')
799
800 return (MinEnergyMolConf, True, MinEnergy)
801
802 def GetEnergy(Mol, ConfID = None):
803 """Calculate energy."""
804
805 Status = True
806 Energy = None
807
808 if ConfID is None:
809 ConfID = -1
810
811 if OptionsInfo["UseUFF"]:
812 UFFMoleculeForcefield = AllChem.UFFGetMoleculeForceField(Mol, confId = ConfID)
813 if UFFMoleculeForcefield is None:
814 Status = False
815 else:
816 Energy = UFFMoleculeForcefield.CalcEnergy()
817 elif OptionsInfo["UseMMFF"]:
818 MMFFMoleculeProperties = AllChem.MMFFGetMoleculeProperties(Mol, mmffVariant = OptionsInfo["MMFFVariant"])
819 MMFFMoleculeForcefield = AllChem.MMFFGetMoleculeForceField(Mol, MMFFMoleculeProperties, confId = ConfID)
820 if MMFFMoleculeForcefield is None:
821 Status = False
822 else:
823 Energy = MMFFMoleculeForcefield.CalcEnergy()
824 else:
825 MiscUtil.PrintError("Couldn't retrieve conformer energy: Specified forcefield, %s, is not supported" % OptionsInfo["ForceField"])
826
827 return (Status, Energy)
828
829 def EmbedMolecule(Mol, MolNum = None):
830 """Embed conformations."""
831
832 ConfIDs = []
833
834 MaxConfs = OptionsInfo["MaxConfs"]
835 RandomSeed = OptionsInfo["RandomSeed"]
836 EnforceChirality = OptionsInfo["EnforceChirality"]
837 UseExpTorsionAnglePrefs = OptionsInfo["UseExpTorsionAnglePrefs"]
838 ETVersion = OptionsInfo["ETVersion"]
839 UseBasicKnowledge = OptionsInfo["UseBasicKnowledge"]
840
841 try:
842 ConfIDs = AllChem.EmbedMultipleConfs(Mol, numConfs = MaxConfs, randomSeed = RandomSeed, enforceChirality = EnforceChirality, useExpTorsionAnglePrefs = UseExpTorsionAnglePrefs, useBasicKnowledge = UseBasicKnowledge, ETversion = ETVersion)
843 except ValueError as ErrMsg:
844 if not OptionsInfo["QuietMode"]:
845 MolName = RDKitUtil.GetMolName(Mol, MolNum)
846 MiscUtil.PrintWarning("Embedding failed for molecule %s:\n%s\n" % (MolName, ErrMsg))
847 ConfIDs = []
848
849 return ConfIDs
850
851 def GenerateStartingTorsionScanStructureOutfile(Mol, MolNum):
852 """Write out the structure of molecule used for starting tosion scan."""
853
854 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
855 MolName = GetOutputFileMolName(Mol, MolNum)
856
857 Outfile = "%s_%s.%s" % (FileName, MolName, FileExt)
858
859 # Set up a molecule writer...
860 Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
861 if Writer is None:
862 MiscUtil.PrintWarning("Failed to setup a writer for output fie %s " % Outfile)
863 return
864
865 Writer.write(Mol)
866
867 if Writer is not None:
868 Writer.close()
869
870 def GenerateOutputFiles(Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
871 """Generate output files."""
872
873 StructureOutfile, EnergyTextOutfile, PlotOutfile = SetupOutputFileNames(Mol, MolNum, TorsionID, TorsionMatchNum)
874
875 GenerateScannedTorsionsStructureOutfile(StructureOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
876 GenerateEnergyTextOutfile(EnergyTextOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
877 GeneratePlotOutfile(PlotOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
878
879 def GenerateScannedTorsionsStructureOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
880 """Write out structures generated after torsion scan along with associated data."""
881
882 # Set up a molecule writer...
883 Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
884 if Writer is None:
885 MiscUtil.PrintWarning("Failed to setup a writer for output fie %s " % Outfile)
886 return
887
888 MolName = RDKitUtil.GetMolName(Mol, MolNum)
889
890 RelativeTorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
891 for Index, TorsionMol in enumerate(TorsionMols):
892 TorsionAngle = "%s" % TorsionAngles[Index]
893 TorsionMol.SetProp("Torsion_Angle", TorsionAngle)
894
895 TorsionEnergy = "%.2f" % TorsionEnergies[Index]
896 TorsionMol.SetProp(OptionsInfo["EnergyLabel"], TorsionEnergy)
897
898 RelativeTorsionEnergy = "%.2f" % RelativeTorsionEnergies[Index]
899 TorsionMol.SetProp(OptionsInfo["RelativeEnergyLabel"], RelativeTorsionEnergy)
900
901 TorsionMolName = "%s_Deg%s" % (MolName, TorsionAngle)
902 TorsionMol.SetProp("_Name", TorsionMolName)
903
904 Writer.write(TorsionMol)
905
906 if Writer is not None:
907 Writer.close()
908
909 def GenerateEnergyTextOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
910 """Write out torsion angles and energies."""
911
912 # Setup a writer...
913 Writer = open(Outfile, "w")
914 if Writer is None:
915 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
916
917 # Write headers...
918 Writer.write("TorsionAngle,%s,%s\n" % (OptionsInfo["EnergyLabel"], OptionsInfo["RelativeEnergyLabel"]))
919
920 RelativeTorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
921 for Index, TorsionAngle in enumerate(TorsionAngles):
922 Writer.write("%d,%.2f,%.2f\n" % (TorsionAngle, TorsionEnergies[Index], RelativeTorsionEnergies[Index]))
923
924 if Writer is not None:
925 Writer.close()
926
927 def GeneratePlotOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
928 """Generate a plot corresponding to torsion angles and energies."""
929
930 OutPlotParams = OptionsInfo["OutPlotParams"]
931
932 # Initialize seaborn and matplotlib paramaters...
933 if not OptionsInfo["OutPlotInitialized"]:
934 OptionsInfo["OutPlotInitialized"] = True
935 RCParams = {"figure.figsize":(OutPlotParams["Width"], OutPlotParams["Height"]),
936 "axes.titleweight": OutPlotParams["TitleWeight"],
937 "axes.labelweight": OutPlotParams["LabelWeight"]}
938 sns.set(context = OutPlotParams["Context"], style = OutPlotParams["Style"], palette = OutPlotParams["Palette"], font = OutPlotParams["Font"], font_scale = OutPlotParams["FontScale"], rc = RCParams)
939
940 # Create a new figure...
941 plt.figure()
942
943 if OptionsInfo["OutPlotRelativeEnergy"]:
944 TorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
945
946 # Draw plot...
947 PlotType = OutPlotParams["Type"]
948 if re.match("linepoint", PlotType, re.I):
949 Axis = sns.lineplot(x = TorsionAngles, y = TorsionEnergies, marker = "o", legend = False)
950 elif re.match("scatter", PlotType, re.I):
951 Axis = sns.scatterplot(x = TorsionAngles, y = TorsionEnergies, legend = False)
952 elif re.match("line", PlotType, re.I):
953 Axis = sns.lineplot(x = TorsionAngles, y = TorsionEnergies, legend = False)
954 else:
955 MiscUtil.PrintError("The value, %s, specified for \"type\" using option \"--outPlotParams\" is not supported. Valid plot types: linepoint, scatter or line" % (PlotType))
956
957 # Setup title and labels...
958 Title = OutPlotParams["Title"]
959 if OptionsInfo["OutPlotTitleTorsionSpec"]:
960 TorsionPattern = OptionsInfo["TorsionsPatternsInfo"]["Pattern"][TorsionID]
961 Title = "%s: %s" % (OutPlotParams["Title"], TorsionPattern)
962
963 # Set labels and title...
964 Axis.set(xlabel = OutPlotParams["XLabel"], ylabel = OutPlotParams["YLabel"], title = Title)
965
966 # Save figure...
967 plt.savefig(Outfile)
968
969 # Close the plot...
970 plt.close()
971
972 def SetupRelativeEnergies(Energies):
973 """Set up a list of relative energies."""
974
975 SortedEnergies = sorted(Energies)
976 MinEnergy = SortedEnergies[0]
977 RelativeEnergies = [(Energy - MinEnergy) for Energy in Energies]
978
979 return RelativeEnergies
980
981 def SetupOutputFileNames(Mol, MolNum, TorsionID, TorsionMatchNum):
982 """Setup names of output files."""
983
984 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
985 MolName = GetOutputFileMolName(Mol, MolNum)
986
987 OutfileRoot = "%s_%s_Torsion%s_Match%s" % (FileName, MolName, TorsionID, TorsionMatchNum)
988
989 StructureOutfile = "%s.%s" % (OutfileRoot, FileExt)
990 EnergyTextOutfile = "%s_Energies.csv" % (OutfileRoot)
991 PlotExt = OptionsInfo["OutPlotParams"]["OutExt"]
992 PlotOutfile = "%s_Plot.%s" % (OutfileRoot, PlotExt)
993
994 return (StructureOutfile, EnergyTextOutfile, PlotOutfile)
995
996 def GetOutputFileMolName(Mol, MolNum):
997 """Get output file prefix."""
998
999 MolName = "Mol%s" % MolNum
1000 if OptionsInfo["OutfileMolName"]:
1001 MolName = re.sub("[^a-zA-Z0-9]", "_", RDKitUtil.GetMolName(Mol, MolNum), re.I)
1002
1003 return MolName
1004
1005 def AddHydrogens(Mol, AddCoords = True):
1006 """Check and add hydrogens."""
1007
1008 if not OptionsInfo["AddHydrogens"]:
1009 return Mol
1010
1011 return Chem.AddHs(Mol, addCoords = AddCoords)
1012
1013 def ProcessTorsionRangeOptions():
1014 """Process tosion range options. """
1015
1016 TosionRangeMode = Options["--torsionRangeMode"]
1017 OptionsInfo["TosionRangeMode"] = TosionRangeMode
1018
1019 if re.match("^Range$", TosionRangeMode, re.I):
1020 ProcessTorsionRangeValues()
1021 elif re.match("^Angles$", TosionRangeMode, re.I):
1022 ProcessTorsionAnglesValues()
1023 else:
1024 MiscUtil.PrintError("The value, %s, option \"--torsionRangeMode\" is not supported." % TosionRangeMode)
1025
1026 def ProcessTorsionRangeValues():
1027 """Process tosion range values. """
1028
1029 TorsionRange = Options["--torsionRange"]
1030 if re.match("^auto$", TorsionRange, re.I):
1031 TorsionRange = "0,360,5"
1032 TorsionRangeWords = TorsionRange.split(",")
1033
1034 TorsionStart = int(TorsionRangeWords[0])
1035 TorsionStop = int(TorsionRangeWords[1])
1036 TorsionStep = int(TorsionRangeWords[2])
1037
1038 if TorsionStart >= TorsionStop:
1039 MiscUtil.PrintError("The start value, %d, specified for option \"--torsionRange\" in string \"%s\" must be less than stop value, %s." % (TorsionStart, Options["--torsionRange"], TorsionStop))
1040 if TorsionStep == 0:
1041 MiscUtil.PrintError("The step value, %d, specified for option \"--torsonRange\" in string \"%s\" must be > 0." % (TorsionStep, Options["--torsionRange"]))
1042 if TorsionStep >= (TorsionStop - TorsionStart):
1043 MiscUtil.PrintError("The step value, %d, specified for option \"--torsonRange\" in string \"%s\" must be less than, %s." % (TorsionStep, Options["--torsionRange"], (TorsionStop - TorsionStart)))
1044
1045 if TorsionStart < 0:
1046 if TorsionStart < -180:
1047 MiscUtil.PrintError("The start value, %d, specified for option \"--torsionRange\" in string \"%s\" must be >= -180 to use scan range from -180 to 180." % (TorsionStart, Options["--torsionRange"]))
1048 if TorsionStop > 180:
1049 MiscUtil.PrintError("The stop value, %d, specified for option \"--torsionRange\" in string \"%s\" must be <= 180 to use scan range from -180 to 180." % (TorsionStop, Options["--torsionRange"]))
1050 else:
1051 if TorsionStop > 360:
1052 MiscUtil.PrintError("The stop value, %d, specified for option \"--torsionRange\" in string \"%s\" must be <= 360 to use scan range from 0 to 360." % (TorsionStop, Options["--torsionRange"]))
1053
1054 TorsionAngles = [Angle for Angle in range(TorsionStart, TorsionStop, TorsionStep)]
1055 TorsionAngles.append(TorsionStop)
1056
1057 OptionsInfo["TorsionRange"] = TorsionRange
1058 OptionsInfo["TorsionStart"] = TorsionStart
1059 OptionsInfo["TorsionStop"] = TorsionStop
1060 OptionsInfo["TorsionStep"] = TorsionStep
1061
1062 OptionsInfo["TorsionAngles"] = TorsionAngles
1063
1064 def ProcessTorsionAnglesValues():
1065 """Process tosion angle values. """
1066
1067 TorsionRange = Options["--torsionRange"]
1068 if re.match("^auto$", TorsionRange, re.I):
1069 MiscUtil.PrintError("The value specified, %s, for option \"--torsionRange\" is not valid." % (TorsionRange))
1070
1071 TorsionAngles = []
1072
1073 for TorsionAngle in TorsionRange.split(","):
1074 TorsionAngle = int(TorsionAngle)
1075
1076 if TorsionAngle < -180:
1077 MiscUtil.PrintError("The stop value, %d, specified for option \"--torsionRange\" in string \"%s\" must be >= -180." % (TorsionAngle, TorsionRange))
1078
1079 if TorsionAngle > 360:
1080 MiscUtil.PrintError("The stop value, %d, specified for option \"--torsionRange\" in string \"%s\" must be <= 360." % (TorsionAngle, TorsionRange))
1081
1082 TorsionAngles.append(TorsionAngle)
1083
1084 OptionsInfo["TorsionRange"] = TorsionRange
1085 OptionsInfo["TorsionStart"] = None
1086 OptionsInfo["TorsionStop"] = None
1087 OptionsInfo["TorsionStep"] = None
1088
1089 OptionsInfo["TorsionAngles"] = sorted(TorsionAngles)
1090
1091 def ProcesssConformerGeneratorOption():
1092 """Process comformer generator option. """
1093
1094 ConfGenParams = MiscUtil.ProcessOptionConformerGenerator("--conformerGenerator", Options["--conformerGenerator"])
1095
1096 OptionsInfo["ConformerGenerator"] = ConfGenParams["ConformerGenerator"]
1097 OptionsInfo["UseBasicKnowledge"] = ConfGenParams["UseBasicKnowledge"]
1098 OptionsInfo["UseExpTorsionAnglePrefs"] = ConfGenParams["UseExpTorsionAnglePrefs"]
1099 OptionsInfo["ETVersion"] = ConfGenParams["ETVersion"]
1100
1101 def ProcessOptions():
1102 """Process and validate command line arguments and options."""
1103
1104 MiscUtil.PrintInfo("Processing options...")
1105
1106 # Validate options...
1107 ValidateOptions()
1108
1109 OptionsInfo["ModeMols"] = Options["--modeMols"]
1110 OptionsInfo["FirstMolMode"] = True if re.match("^First$", Options["--modeMols"], re.I) else False
1111
1112 OptionsInfo["ModeTorsions"] = Options["--modeTorsions"]
1113 OptionsInfo["FirstTorsionMode"] = True if re.match("^First$", Options["--modeTorsions"], re.I) else False
1114
1115 OptionsInfo["Infile"] = Options["--infile"]
1116 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"], Options["--infile"])
1117 OptionsInfo["Infile3D"] = True if re.match("^yes$", Options["--infile3D"], re.I) else False
1118
1119 OptionsInfo["Outfile"] = Options["--outfile"]
1120 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters("--outfileParams", Options["--outfileParams"])
1121
1122 OptionsInfo["OutfileMolName"] = True if re.match("^yes$", Options["--outfileMolName"], re.I) else False
1123
1124 OptionsInfo["OutPlotRelativeEnergy"] = True if re.match("^yes$", Options["--outPlotRelativeEnergy"], re.I) else False
1125 OptionsInfo["OutPlotTitleTorsionSpec"] = True if re.match("^yes$", Options["--outPlotTitleTorsionSpec"], re.I) else False
1126
1127 # The default width and height, 10.0 and 7.5, map to aspect raito of 16/9 (1.778)...
1128 DefaultValues = {'Type': 'linepoint', 'Width': 10.0, 'Height': 5.6, 'Title': 'RDKit Torsion Scan', 'XLabel': 'Torsion Angle (degrees)', 'YLabel': 'Energy (kcal/mol)'}
1129 if OptionsInfo["OutPlotRelativeEnergy"]:
1130 DefaultValues["YLabel"] = "Relative Energy (kcal/mol)"
1131 OptionsInfo["OutPlotParams"] = MiscUtil.ProcessOptionSeabornPlotParameters("--outPlotParams", Options["--outPlotParams"], DefaultValues)
1132 if not re.match("^(linepoint|scatter|Line)$", OptionsInfo["OutPlotParams"]["Type"], re.I):
1133 MiscUtil.PrintError("The value, %s, specified for \"type\" using option \"--outPlotParams\" is not supported. Valid plot types: linepoint, scatter or line" % (OptionsInfo["OutPlotParams"]["Type"]))
1134
1135 OptionsInfo["OutPlotInitialized"] = False
1136
1137 # Procsss and validate specified SMILES/SMARTS torsion patterns...
1138 TorsionPatterns = Options["--torsions"]
1139 TorsionPatternsList = []
1140 for TorsionPattern in TorsionPatterns.split(","):
1141 TorsionPattern = TorsionPattern.strip()
1142 if not len(TorsionPattern):
1143 MiscUtil.PrintError("Empty value specified for SMILES/SMARTS pattern in \"-t, --torsions\" option: %s" % TorsionPatterns)
1144
1145 TorsionMol = Chem.MolFromSmarts(TorsionPattern)
1146 if TorsionMol is None:
1147 MiscUtil.PrintError("Failed to create torsion pattern molecule. The torsion SMILES/SMARTS pattern, \"%s\", specified using \"-t, --torsions\" option, \"%s\", is not valid." % (TorsionPattern, TorsionPatterns))
1148 TorsionPatternsList.append(TorsionPattern)
1149
1150 OptionsInfo["TorsionPatterns"] = TorsionPatterns
1151 OptionsInfo["TorsionPatternsList"] = TorsionPatternsList
1152
1153 # Process and validate any specified torsion atom indices for filtering torsion matches...
1154 TorsionsFilterByAtomIndices = Options["--torsionsFilterbyAtomIndices"]
1155 TorsionsFilterByAtomIndicesList = []
1156 if not re.match("^None$", TorsionsFilterByAtomIndices, re.I):
1157 for AtomIndex in TorsionsFilterByAtomIndices.split(","):
1158 AtomIndex = AtomIndex.strip()
1159 if not MiscUtil.IsInteger(AtomIndex):
1160 MiscUtil.PrintError("The value specified, %s, for option \"--torsionsFilterbyAtomIndices\" must be an integer." % AtomIndex)
1161 AtomIndex = int(AtomIndex)
1162 if AtomIndex < 0:
1163 MiscUtil.PrintError("The value specified, %s, for option \"--torsionsFilterbyAtomIndices\" must be >= 0." % AtomIdex)
1164 TorsionsFilterByAtomIndicesList.append(AtomIndex)
1165
1166 if len(TorsionsFilterByAtomIndicesList) < 4:
1167 MiscUtil.PrintError("The number of values, %s, specified, %s, for option \"--torsionsFilterbyAtomIndices\" must be >=4." % (len(TorsionsFilterByAtomIndicesList), TorsionsFilterByAtomIndices))
1168
1169 OptionsInfo["TorsionsFilterByAtomIndices"] = TorsionsFilterByAtomIndices
1170 OptionsInfo["TorsionsFilterByAtomIndicesList"] = TorsionsFilterByAtomIndicesList
1171 OptionsInfo["FilterTorsionsByAtomIndicesMode"] = True if len(TorsionsFilterByAtomIndicesList) > 0 else False
1172
1173 OptionsInfo["Overwrite"] = Options["--overwrite"]
1174
1175 OptionsInfo["AddHydrogens"] = True if re.match("^yes$", Options["--addHydrogens"], re.I) else False
1176
1177 ProcesssConformerGeneratorOption()
1178
1179 if re.match("^UFF$", Options["--forceField"], re.I):
1180 ForceField = "UFF"
1181 UseUFF = True
1182 UseMMFF = False
1183 elif re.match("^MMFF$", Options["--forceField"], re.I):
1184 ForceField = "MMFF"
1185 UseUFF = False
1186 UseMMFF = True
1187 else:
1188 MiscUtil.PrintError("The value, %s, specified for \"--forceField\" is not supported." % (Options["--forceField"],))
1189
1190 MMFFVariant = "MMFF94" if re.match("^MMFF94$", Options["--forceFieldMMFFVariant"], re.I) else "MMFF94s"
1191
1192 OptionsInfo["ForceField"] = ForceField
1193 OptionsInfo["MMFFVariant"] = MMFFVariant
1194 OptionsInfo["UseMMFF"] = UseMMFF
1195 OptionsInfo["UseUFF"] = UseUFF
1196
1197 EnergyUnits = "kcal/mol"
1198 if UseMMFF:
1199 OptionsInfo["EnergyLabel"] = "%s_Energy" % MMFFVariant
1200 OptionsInfo["RelativeEnergyLabel"] = "%s_Relative_Energy" % MMFFVariant
1201 else:
1202 OptionsInfo["EnergyLabel"] = "%s_Energy" % ForceField
1203 OptionsInfo["RelativeEnergyLabel"] = "%s_Relative_Energy" % ForceField
1204
1205 OptionsInfo["EnforceChirality"] = True if re.match("^yes$", Options["--enforceChirality"], re.I) else False
1206
1207 OptionsInfo["MaxConfs"] = int(Options["--maxConfs"])
1208 OptionsInfo["MaxConfsTorsion"] = int(Options["--maxConfsTorsion"])
1209 OptionsInfo["MaxIters"] = int(Options["--maxIters"])
1210
1211 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
1212 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
1213
1214 # Multiprocessing level...
1215 MPLevelMoleculesMode = False
1216 MPLevelTorsionAnglesMode = False
1217 MPLevel = Options["--mpLevel"]
1218 if re.match("^Molecules$", MPLevel, re.I):
1219 MPLevelMoleculesMode = True
1220 elif re.match("^TorsionAngles$", MPLevel, re.I):
1221 MPLevelTorsionAnglesMode = True
1222 else:
1223 MiscUtil.PrintError("The value, %s, specified for option \"--mpLevel\" is not valid. " % MPLevel)
1224 OptionsInfo["MPLevel"] = MPLevel
1225 OptionsInfo["MPLevelMoleculesMode"] = MPLevelMoleculesMode
1226 OptionsInfo["MPLevelTorsionAnglesMode"] = MPLevelTorsionAnglesMode
1227
1228 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
1229
1230 RandomSeed = -1
1231 if not re.match("^auto$", Options["--randomSeed"], re.I):
1232 RandomSeed = int(Options["--randomSeed"])
1233 OptionsInfo["RandomSeed"] = RandomSeed
1234
1235 OptionsInfo["RemoveHydrogens"] = True if re.match("^yes$", Options["--removeHydrogens"], re.I) else False
1236
1237 OptionsInfo["TorsionMaxMatches"] = int(Options["--torsionMaxMatches"])
1238 OptionsInfo["TorsionMinimize"] = True if re.match("^yes$", Options["--torsionMinimize"], re.I) else False
1239
1240 ProcessTorsionRangeOptions()
1241
1242 OptionsInfo["UseTethers"] = True if re.match("^yes$", Options["--useTethers"], re.I) else False
1243 OptionsInfo["UseChirality"] = True if re.match("^yes$", Options["--useChirality"], re.I) else False
1244
1245 def RetrieveOptions():
1246 """Retrieve command line arguments and options."""
1247
1248 # Get options...
1249 global Options
1250 Options = docopt(_docoptUsage_)
1251
1252 # Set current working directory to the specified directory...
1253 WorkingDir = Options["--workingdir"]
1254 if WorkingDir:
1255 os.chdir(WorkingDir)
1256
1257 # Handle examples option...
1258 if "--examples" in Options and Options["--examples"]:
1259 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
1260 sys.exit(0)
1261
1262 def ValidateOptions():
1263 """Validate option values."""
1264
1265 MiscUtil.ValidateOptionTextValue("-a, --addHydrogens", Options["--addHydrogens"], "yes no")
1266 MiscUtil.ValidateOptionTextValue("-c, --conformerGenerator", Options["--conformerGenerator"], "SDG KDG ETDG ETKDG ETKDGv2")
1267
1268 MiscUtil.ValidateOptionTextValue("-f, --forceField", Options["--forceField"], "UFF MMFF")
1269 MiscUtil.ValidateOptionTextValue(" --forceFieldMMFFVariant", Options["--forceFieldMMFFVariant"], "MMFF94 MMFF94s")
1270
1271 MiscUtil.ValidateOptionTextValue("--enforceChirality ", Options["--enforceChirality"], "yes no")
1272
1273 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
1274 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi txt csv tsv")
1275 MiscUtil.ValidateOptionTextValue("--infile3D", Options["--infile3D"], "yes no")
1276
1277 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
1278 MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
1279 MiscUtil.ValidateOptionsDistinctFileNames("-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"])
1280
1281 if not Options["--overwrite"]:
1282 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
1283 FileNames = glob.glob("%s_*" % FileName)
1284 if len(FileNames):
1285 MiscUtil.PrintError("The outfile names, %s_*, generated from file specified, %s, for option \"-o, --outfile\" already exist. Use option \"--overwrite\" or \"--ov\" and try again.\n" % (FileName, Options["--outfile"]))
1286
1287 MiscUtil.ValidateOptionTextValue("--outPlotRelativeEnergy", Options["--outPlotRelativeEnergy"], "yes no")
1288 MiscUtil.ValidateOptionTextValue("--outPlotTitleTorsionSpec", Options["--outPlotTitleTorsionSpec"], "yes no")
1289
1290 MiscUtil.ValidateOptionTextValue("--outfileMolName ", Options["--outfileMolName"], "yes no")
1291
1292 MiscUtil.ValidateOptionTextValue("--modeMols", Options["--modeMols"], "First All")
1293 MiscUtil.ValidateOptionTextValue("--modeTorsions", Options["--modeTorsions"], "First All")
1294
1295 MiscUtil.ValidateOptionIntegerValue("--maxConfs", Options["--maxConfs"], {">": 0})
1296 MiscUtil.ValidateOptionIntegerValue("--maxConfsTorsion", Options["--maxConfsTorsion"], {">": 0})
1297 MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0})
1298
1299 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
1300 MiscUtil.ValidateOptionTextValue("--mpLevel", Options["--mpLevel"], "Molecules TorsionAngles")
1301 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
1302
1303 if not re.match("^auto$", Options["--randomSeed"], re.I):
1304 MiscUtil.ValidateOptionIntegerValue("--randomSeed", Options["--randomSeed"], {})
1305
1306 MiscUtil.ValidateOptionTextValue("-r, --removeHydrogens", Options["--removeHydrogens"], "yes no")
1307
1308 MiscUtil.ValidateOptionIntegerValue("--torsionMaxMatches", Options["--torsionMaxMatches"], {">": 0})
1309 MiscUtil.ValidateOptionTextValue("--torsionMinimize", Options["--torsionMinimize"], "yes no")
1310
1311 MiscUtil.ValidateOptionTextValue("--torsionRangeMode", Options["--torsionRangeMode"], "Range or Angles")
1312 TorsionRange = Options["--torsionRange"]
1313 if re.match("^Range$", Options["--torsionRangeMode"], re.I):
1314 if not re.match("^auto$", TorsionRange, re.I):
1315 MiscUtil.ValidateOptionNumberValues("--torsionRange", TorsionRange, 3, ",", "integer", {})
1316 else:
1317 if re.match("^auto$", TorsionRange, re.I):
1318 MiscUtil.PrintError("The value, %s, specified for option \"-torsionRange\" is not valid for \"%s\" value of \"--torsionRangeMode\" option. You must specify a torsion angle or a comma delimited list of torsion angles." % (TorsionRange, Options["--torsionRangeMode"]))
1319 TorsionAngles = []
1320 for TorsionAngle in TorsionRange.split(","):
1321 TorsionAngle = TorsionAngle.strip()
1322 if not MiscUtil.IsInteger(TorsionAngle):
1323 MiscUtil.PrintError("The value specified, %s, for option \"--torsionRange\" in string \"%s\" must be an integer." % (TorsionAngle, TorsionRange))
1324 if TorsionAngle in TorsionAngles:
1325 MiscUtil.PrintError("The value specified, %s, for option \"--torsionRange\" in string \"%s\" is a duplicate value." % (TorsionAngle, TorsionRange))
1326 TorsionAngles.append(TorsionAngle)
1327
1328 MiscUtil.ValidateOptionTextValue("--useChirality", Options["--useChirality"], "yes no")
1329 MiscUtil.ValidateOptionTextValue("--useTethers", Options["--useTethers"], "yes no")
1330
1331 # Setup a usage string for docopt...
1332 _docoptUsage_ = """
1333 RDKitPerformTorsionScan.py - Perform torsion scan
1334
1335 Usage:
1336 RDKitPerformTorsionScan.py [--addHydrogens <yes or no>] [--conformerGenerator <SDG, KDG, ETDG, ETKDG, ETKDGv2>]
1337 [--forceField <UFF, or MMFF>] [--forceFieldMMFFVariant <MMFF94 or MMFF94s>]
1338 [--enforceChirality <yes or no>] [--infile3D <yes or no>] [--infileParams <Name,Value,...>]
1339 [--modeMols <First or All>] [--modeTorsions <First or All> ] [--maxConfs <number>]
1340 [--maxConfsTorsion <number>] [--maxIters <number>] [--mp <yes or no>]
1341 [--mpLevel <Molecules or TorsionAngles>] [--mpParams <Name,Value,...>]
1342 [--outfileMolName <yes or no>] [--outfileParams <Name,Value,...>] [--outPlotParams <Name,Value,...>]
1343 [--outPlotRelativeEnergy <yes or no>] [--outPlotTitleTorsionSpec <yes or no>] [--overwrite] [--quiet <yes or no>]
1344 [--removeHydrogens <yes or no>] [--randomSeed <number>] [--torsionsFilterbyAtomIndices <Index1, Index2, ...>]
1345 [--torsionMaxMatches <number>] [--torsionMinimize <yes or no>] [--torsionRangeMode <Range or Angles>]
1346 [--torsionRange <Start,Stop,Step or Angle1,Angle2,...>] [--useChirality <yes or no>] [--useTethers <yes or no>]
1347 [-w <dir>] -t <torsions> -i <infile> -o <outfile>
1348 RDKitPerformTorsionScan.py -h | --help | -e | --examples
1349
1350 Description:
1351 Perform torsion scan for molecules around torsion angles specified using
1352 SMILES/SMARTS patterns. A molecule is optionally minimized before performing
1353 a torsion scan. A set of initial 3D structures are generated for a molecule
1354 by scanning the torsion angle across the specified range and updating the 3D
1355 coordinates of the molecule. A conformation ensemble is optionally generated
1356 for each 3D structure representing a specific torsion angle. The conformation
1357 with the lowest energy is selected to represent the torsion angle. An option
1358 is available to skip the generation of the conformation ensemble and simply
1359 calculate the energy for the initial 3D structure for a specific torsion angle.
1360
1361 The torsions are specified using SMILES or SMARTS patterns. A substructure match
1362 is performed to select torsion atoms in a molecule. The SMILES pattern match must
1363 correspond to four torsion atoms. The SMARTS patterns containing atom map numbers
1364 may match more than four atoms. The atom map numbers, however, must match
1365 exactly four torsion atoms. For example: [s:1][c:2]([aX2,cH1])!@[CX3:3](O)=[O:4] for
1366 thiophene esters and carboxylates as specified in Torsion Library (TorLib) [Ref 146].
1367
1368 A set of four output files is generated for each torsion match in each
1369 molecule. The names of the output files are generated using the root of
1370 the specified output file. They may either contain sequential molecule
1371 numbers or molecule names as shown below:
1372
1373 <OutfileRoot>_Mol<Num>.sdf
1374 <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>.sdf
1375 <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Energies.csv
1376 <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Plot.<ImgExt>
1377
1378 or
1379
1380 <OutfileRoot>_<MolName>.sdf
1381 <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>.sdf
1382 <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Energies.csv
1383 <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Plot.<ImgExt>
1384
1385 The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi,
1386 .csv, .tsv, .txt)
1387
1388 The supported output file formats are: SD (.sdf, .sd)
1389
1390 Options:
1391 -a, --addHydrogens <yes or no> [default: yes]
1392 Add hydrogens before minimization.
1393 -c, --conformerGenerator <text> [default: ETKDGv2]
1394 Conformation generation methodology for generating initial 3D structure
1395 of a molecule and conformation ensemble representing a specific torsion
1396 angle. No conformation ensemble is generated for 'No' value of
1397 '--torsionMinimize' option.
1398
1399 The possible values along with a brief description are shown below:
1400
1401 SDG: Standard Distance Geometry
1402 KDG: basic Knowledge-terms with Distance Geometry
1403 ETDG: Experimental Torsion-angle preference with Distance Geometry
1404 ETKDG: Experimental Torsion-angle preference along with basic
1405 Knowledge-terms and Distance Geometry [Ref 129]
1406 ETKDGv2: Experimental Torsion-angle preference along with basic
1407 Knowledge-terms and Distance Geometry [Ref 167]
1408 -f, --forceField <UFF, MMFF> [default: MMFF]
1409 Forcefield method to use for energy minimization of initial 3D structure
1410 of a molecule and conformation ensemble representing a specific torsion.
1411 No conformation ensemble is generated during for 'No' value of '--torsionMinimze'
1412 option and constrained energy minimization is not performed. Possible values:
1413 Universal Force Field (UFF) [ Ref 81 ] or Merck Molecular Mechanics Force
1414 Field [ Ref 83-87 ] .
1415 --forceFieldMMFFVariant <MMFF94 or MMFF94s> [default: MMFF94]
1416 Variant of MMFF forcefield to use for energy minimization.
1417 --enforceChirality <yes or no> [default: Yes]
1418 Enforce chirality for defined chiral centers during generation of conformers.
1419 -e, --examples
1420 Print examples.
1421 -h, --help
1422 Print this help message.
1423 -i, --infile <infile>
1424 Input file name.
1425 --infile3D <yes or no> [default: no]
1426 Skip generation and minimization of initial 3D structures for molecules in
1427 input file containing 3D coordinates.
1428 --infileParams <Name,Value,...> [default: auto]
1429 A comma delimited list of parameter name and value pairs for reading
1430 molecules from files. The supported parameter names for different file
1431 formats, along with their default values, are shown below:
1432
1433 SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes
1434
1435 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
1436 smilesTitleLine,auto,sanitize,yes
1437
1438 Possible values for smilesDelimiter: space, comma or tab.
1439 --modeMols <First or All> [default: First]
1440 Perform torsion scan for the first molecule or all molecules in input
1441 file.
1442 --modeTorsions <First or All> [default: First]
1443 Perform torsion scan for the first or all specified torsion pattern in
1444 molecules up to a maximum number of matches for each torsion
1445 specification as indicated by '--torsionMaxMatches' option.
1446 --maxConfs <number> [default: 250]
1447 Maximum number of conformations to generate for initial 3D structure of a
1448 molecule. The lowest energy conformation is written to the output file.
1449 --maxConfsTorsion <number> [default: 50]
1450 Maximum number of conformations to generate for conformation ensemble
1451 representing a specific torsion. A constrained minimization is performed
1452 using the coordinates of the specified torsion and the lowest energy
1453 conformation is written to the output file.
1454 --maxIters <number> [default: 500]
1455 Maximum number of iterations to perform for a molecule during minimization
1456 to generation initial 3D structures. This option is ignored during 'yes' value
1457 of '--infile3D' option.
1458 --mp <yes or no> [default: no]
1459 Use multiprocessing.
1460
1461 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
1462 function employing lazy RDKit data iterable. This allows processing of
1463 arbitrary large data sets without any additional requirements memory.
1464
1465 All input data may be optionally loaded into memory by mp.Pool.map()
1466 before starting worker processes in a process pool by setting the value
1467 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
1468
1469 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
1470 data mode may adversely impact the performance. The '--mpParams' section
1471 provides additional information to tune the value of 'chunkSize'.
1472 --mpLevel <Molecules or TorsionAngles> [default: Molecules]
1473 Perform multiprocessing at molecules or torsion angles level. Possible values:
1474 Molecules or TorsionAngles. The 'Molecules' value starts a process pool at the
1475 molecules level. All torsion angles of a molecule are processed in a single
1476 process. The 'TorsionAngles' value, however, starts a process pool at the
1477 torsion angles level. Each torsion angle in a torsion match for a molecule is
1478 processed in an individual process in the process pool.
1479 --mpParams <Name,Value,...> [default: auto]
1480 A comma delimited list of parameter name and value pairs to configure
1481 multiprocessing.
1482
1483 The supported parameter names along with their default and possible
1484 values are shown below:
1485
1486 chunkSize, auto
1487 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
1488 numProcesses, auto [ Default: mp.cpu_count() ]
1489
1490 These parameters are used by the following functions to configure and
1491 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
1492 mp.Pool.imap().
1493
1494 The chunkSize determines chunks of input data passed to each worker
1495 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
1496 The default value of chunkSize is dependent on the value of 'inputDataMode'.
1497
1498 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
1499 automatically converts RDKit data iterable into a list, loads all data into
1500 memory, and calculates the default chunkSize using the following method
1501 as shown in its code:
1502
1503 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
1504 if extra: chunkSize += 1
1505
1506 For example, the default chunkSize will be 7 for a pool of 4 worker processes
1507 and 100 data items.
1508
1509 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
1510 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
1511 data into memory. Consequently, the size of input data is not known a priori.
1512 It's not possible to estimate an optimal value for the chunkSize. The default
1513 chunkSize is set to 1.
1514
1515 The default value for the chunkSize during 'Lazy' data mode may adversely
1516 impact the performance due to the overhead associated with exchanging
1517 small chunks of data. It is generally a good idea to explicitly set chunkSize to
1518 a larger value during 'Lazy' input data mode, based on the size of your input
1519 data and number of processes in the process pool.
1520
1521 The mp.Pool.map() function waits for all worker processes to process all
1522 the data and return the results. The mp.Pool.imap() function, however,
1523 returns the the results obtained from worker processes as soon as the
1524 results become available for specified chunks of data.
1525
1526 The order of data in the results returned by both mp.Pool.map() and
1527 mp.Pool.imap() functions always corresponds to the input data.
1528 -o, --outfile <outfile>
1529 Output file name. The output file root is used for generating the names
1530 of the output files corresponding to structures, energies, and plots during
1531 the torsion scan.
1532 --outfileMolName <yes or no> [default: no]
1533 Append molecule name to output file root during the generation of the names
1534 for output files. The default is to use <MolNum>. The non alphabetical
1535 characters in molecule names are replaced by underscores.
1536 --outfileParams <Name,Value,...> [default: auto]
1537 A comma delimited list of parameter name and value pairs for writing
1538 molecules to files. The supported parameter names for different file
1539 formats, along with their default values, are shown below:
1540
1541 SD: kekulize,yes,forceV3000,no
1542
1543 --outPlotParams <Name,Value,...> [default: auto]
1544 A comma delimited list of parameter name and value pairs for generating
1545 plots using Seaborn module. The supported parameter names along with their
1546 default values are shown below:
1547
1548 type,linepoint,outExt,svg,width,10,height,5.6,
1549 title,auto,xlabel,auto,ylabel,auto,titleWeight,bold,labelWeight,bold
1550 style,darkgrid,palette,deep,font,sans-serif,fontScale,1,
1551 context,notebook
1552
1553 Possible values:
1554
1555 type: linepoint, scatter, or line. Both points and lines are drawn
1556 for linepoint plot type.
1557 outExt: Any valid format supported by Python module Matplotlib.
1558 For example: PDF (.pdf), PNG (.png), PS (.ps), SVG (.svg)
1559 titleWeight, labelWeight: Font weight for title and axes labels.
1560 Any valid value.
1561 style: darkgrid, whitegrid, dark, white, ticks
1562 palette: deep, muted, pastel, dark, bright, colorblind
1563 font: Any valid font name
1564 context: paper, notebook, talk, poster, or any valid name
1565
1566 --outPlotRelativeEnergy <yes or no> [default: yes]
1567 Plot relative energies in the torsion plot. The minimum energy value is
1568 subtracted from energy values to calculate relative energies.
1569 --outPlotTitleTorsionSpec <yes or no> [default: yes]
1570 Append torsion specification to the title of the torsion plot.
1571 --overwrite
1572 Overwrite existing files.
1573 -q, --quiet <yes or no> [default: no]
1574 Use quiet mode. The warning and information messages will not be printed.
1575 --randomSeed <number> [default: auto]
1576 Seed for the random number generator for generating initial 3D coordinates.
1577 Default is to use a random seed.
1578 --removeHydrogens <yes or no> [default: Yes]
1579 Remove hydrogens after minimization.
1580 -t, --torsions <SMILES/SMARTS,...,...>
1581 SMILES/SMARTS patterns corresponding to torsion specifications. It's a
1582 comma delimited list of valid SMILES/SMART patterns.
1583
1584 A substructure match is performed to select torsion atoms in a molecule.
1585 The SMILES pattern match must correspond to four torsion atoms. The
1586 SMARTS patterns containing atom map numbers may match more than four
1587 atoms. The atom map numbers, however, must match exactly four torsion
1588 atoms. For example: [s:1][c:2]([aX2,cH1])!@[CX3:3](O)=[O:4] for thiophene
1589 esters and carboxylates as specified in Torsion Library (TorLib) [Ref 146].
1590 --torsionsFilterbyAtomIndices <Index1, Index2, ...> [default: none]
1591 Comma delimited list of atom indices for filtering torsion matches
1592 corresponding to torsion specifications "-t, --torsions". The atom indices
1593 must be valid. No explicit validation is performed. The list must contain at
1594 least 4 atom indices.
1595
1596 The torsion atom indices, matched by "-t, --torsions" specifications, must be
1597 present in the list. Otherwise, the torsion matches are ignored.
1598 --torsionMaxMatches <number> [default: 5]
1599 Maximum number of torsions to match for each torsion specification in a
1600 molecule.
1601 --torsionMinimize <yes or no> [default: no]
1602 Perform constrained energy minimization on a conformation ensemble
1603 for a specific torsion angle and select the lowest energy conformation
1604 representing the torsion angle.
1605 --torsionRangeMode <Range or Angles> [default: Range]
1606 Perform torsion scan using torsion angles corresponding to a torsion range
1607 or an explicit list of torsion angles. Possible values: Range or Angles. You
1608 may use '--torsionRange' option to specify values for torsion angle or
1609 torsion angles.
1610 --torsionRange <Start,Stop,Step or Angle1,Angle2...> [default: auto]
1611 Start, stop, and step size angles or a comma delimited list of angles in
1612 degrees for a torsion scan.
1613
1614 This value is '--torsionRangeMode' specific. It must be a triplet corresponding
1615 to 'start,Stop,Step' for 'Range' value of '--torsionRange' option. Otherwise, it
1616 is comma delimited list of one or more torsion angles for 'Angles' value of
1617 '--torsionRange' option.
1618
1619 The default values, based on '--torsionRangeMode' option, are shown below:
1620
1621 TorsionRangeMode Default value
1622 Range 0,360,5
1623 Angles None
1624
1625 You must explicitly provide a list of torsion angle(s) for 'Angles' of
1626 '--torsionRangeMode' option.
1627 --useChirality <yes or no> [default: no]
1628 Use chirrality during substructure matches for identification of torsions.
1629 --useTethers <yes or no> [default: yes]
1630 Use tethers to optimize the final conformation by applying a series of extra forces
1631 to align matching atoms to the positions of the core atoms. Otherwise, use simple
1632 distance constraints during the optimization.
1633 -w, --workingdir <dir>
1634 Location of working directory which defaults to the current directory.
1635
1636 Examples:
1637 To perform a torsion scan from 0 to 360 degrees with a stepsize of 5 on the
1638 first molecule in a SMILES file using a minimum energy structure of the molecule
1639 selected from an ensemble of conformations, skipping generation of conformation
1640 ensembles for specific torsion angles and constrained energy minimization of the
1641 ensemble, generate output files corresponding to structure, energy and torsion
1642 plot, type:
1643
1644 % RDKitPerformTorsionScan.py -t "O=CNC" -i SampleSeriesD3R.smi
1645 -o SampleOut.sdf
1646
1647 To run the previous example for performing a torsion scan using a specific list
1648 of torsion angles, type:
1649
1650 % RDKitPerformTorsionScan.py -t "O=CNC" -i SampleSeriesD3R.smi
1651 -o SampleOut.sdf --torsionRangleMode Angles
1652 --torsionRange "160,220,280"
1653
1654 To run the previous example on all molecules in a SD file, type:
1655
1656 % RDKitPerformTorsionScan.py -t "O=CNC" --modeMols All
1657 -i SampleSeriesD3R.sdf -o SampleOut.sdf
1658
1659 To perform a torsion scan on the first molecule in a SMILES file using a minimum
1660 energy structure of the molecule selected from an ensemble of conformations,
1661 generation of conformation ensembles for specific torsion angles and constrained
1662 energy minimization of the ensemble, generate output files corresponding to
1663 structure, energy and torsion plot, type:
1664
1665 % RDKitPerformTorsionScan.py -t "O=CNC" --torsionMinimize Yes
1666 -i SampleSeriesD3R.smi -o SampleOut.sdf
1667
1668 To run the previous example on all molecules in a SD file, type:
1669
1670 % RDKitPerformTorsionScan.py -t "O=CNC" --modeMols All
1671 --torsionMinimize Yes -i SampleSeriesD3R.sdf -o SampleOut.sdf
1672
1673 To run the previous example in multiprocessing mode at molecules level
1674 on all available CPUs without loading all data into memory and write out
1675 a SD file, type:
1676
1677 % RDKitPerformTorsionScan.py -t "O=CNC" -i SampleSeriesD3R.smi
1678 -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
1679
1680 To run the previous example in multiprocessing mode at torsion angles level
1681 on all available CPUs without loading all data into memory and write out
1682 a SD file, type:
1683
1684 % RDKitPerformTorsionScan.py -t "O=CNC" -i SampleSeriesD3R.smi
1685 -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
1686 --mpLevel TorsionAngles
1687
1688 To run the previous example in multiprocessing mode on all available CPUs
1689 by loading all data into memory and write out a SD file, type:
1690
1691 % RDKitPerformTorsionScan.py -t "O=CNC" -i SampleSeriesD3R.smi
1692 -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
1693 --mpParams "inputDataMode,InMemory"
1694
1695 To run the previous example in multiprocessing mode on specific number of
1696 CPUs and chunk size without loading all data into memory and write out a SD file,
1697 type:
1698
1699 % RDKitPerformTorsionScan.py -t "O=CNC" -i SampleSeriesD3R.smi
1700 -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
1701 --mpParams "inputDataMode,Lazy,numProcesses,4,chunkSize,8"
1702
1703 To perform a torsion scan on first molecule in a SD file containing 3D coordinates,
1704 skipping generation of conformation ensembles for specific torsion angles and
1705 constrained energy minimization of the ensemble, generate output files
1706 corresponding to structure, energy and torsion plot, type:
1707
1708 % RDKitPerformTorsionScan.py -t "O=CNC" --infile3D yes
1709 -i SampleSeriesD3R3D.sdf -o SampleOut.sdf
1710
1711 To perform a torsion scan using multiple torsion specifications on all molecules in
1712 a SD file containing 3D coordinates, generation of conformation ensembles for specific
1713 torsion angles and constrained energy minimization of the ensemble, generate output files
1714 corresponding to structure, energy and torsion plot, type:
1715
1716 % RDKitPerformTorsionScan.py -t "O=CNC,[O:1]=[C:2](c)[N:3][C:4]"
1717 --infile3D yes --modeMols All --modeTorsions All
1718 --torsionMinimize Yes -i SampleSeriesD3R3D.sdf -o SampleOut.sdf
1719
1720 To run the previous example using a specific torsion scan range, type:
1721
1722 % RDKitPerformTorsionScan.py -t "O=CNC,[O:1]=[C:2](c)[N:3][C:4]"
1723 --infile3D yes --modeMols All --modeTorsions All --torsionMinimize
1724 Yes --torsionRange 0,360,10 -i SampleSeriesD3R.smi -o SampleOut.sdf
1725
1726 Author:
1727 Manish Sud(msud@san.rr.com)
1728
1729 See also:
1730 RDKitCalculateRMSD.py, RDKitCalculateMolecularDescriptors.py, RDKitCompareMoleculeShapes.py,
1731 RDKitConvertFileFormat.py, RDKitPerformConstrainedMinimization.py
1732
1733 Copyright:
1734 Copyright (C) 2024 Manish Sud. All rights reserved.
1735
1736 The functionality available in this script is implemented using RDKit, an
1737 open source toolkit for cheminformatics developed by Greg Landrum.
1738
1739 This file is part of MayaChemTools.
1740
1741 MayaChemTools is free software; you can redistribute it and/or modify it under
1742 the terms of the GNU Lesser General Public License as published by the Free
1743 Software Foundation; either version 3 of the License, or (at your option) any
1744 later version.
1745
1746 """
1747
1748 if __name__ == "__main__":
1749 main()
