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