1 #!/bin/env python
2 #
3 # File: Psi4PerformTorsionScan.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 Psi4, an
11 # open source quantum chemistry software package, and RDKit, an open
12 # source toolkit for cheminformatics developed by Greg Landrum.
13 #
14 # This file is part of MayaChemTools.
15 #
16 # MayaChemTools is free software; you can redistribute it and/or modify it under
17 # the terms of the GNU Lesser General Public License as published by the Free
18 # Software Foundation; either version 3 of the License, or (at your option) any
19 # later version.
20 #
21 # MayaChemTools is distributed in the hope that it will be useful, but without
22 # any warranty; without even the implied warranty of merchantability of fitness
23 # for a particular purpose. See the GNU Lesser General Public License for more
24 # details.
25 #
26 # You should have received a copy of the GNU Lesser General Public License
27 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
28 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
29 # Boston, MA, 02111-1307, USA.
30 #
31 #
32
33 from __future__ import print_function
34
35 import os
36 import sys
37 import time
38 import re
39 import glob
40 import shutil
41 import multiprocessing as mp
42
43 import matplotlib.pyplot as plt
44 import seaborn as sns
45
46 # Psi4 imports...
47 if hasattr(shutil, "which") and shutil.which("psi4") is None:
48 sys.stderr.write("\nWarning: Failed to find 'psi4' in your PATH indicating potential issues with your\n")
49 sys.stderr.write("Psi4 environment. The 'import psi4' directive in the global scope of the script\n")
50 sys.stderr.write("interferes with the multiprocessing functionality. It is imported later in the\n")
51 sys.stderr.write("local scope during the execution of the script and may fail. Check/update your\n")
52 sys.stderr.write("Psi4 environment and try again.\n\n")
53
54 # RDKit imports...
55 try:
56 from rdkit import rdBase
57 from rdkit import Chem
58 from rdkit.Chem import AllChem
59 from rdkit.Chem import rdMolAlign
60 from rdkit.Chem import rdMolTransforms
61 except ImportError as ErrMsg:
62 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
63 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
64 sys.exit(1)
65
66 # MayaChemTools imports...
67 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
68 try:
69 from docopt import docopt
70 import MiscUtil
71 import Psi4Util
72 import RDKitUtil
73 except ImportError as ErrMsg:
74 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
75 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
76 sys.exit(1)
77
78 ScriptName = os.path.basename(sys.argv[0])
79 Options = {}
80 OptionsInfo = {}
81
82
83 def main():
84 """Start execution of the script."""
85
86 MiscUtil.PrintInfo(
87 "\n%s (Psi4: Imported later; RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
88 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime())
89 )
90
91 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
92
93 # Retrieve command line arguments and options...
94 RetrieveOptions()
95
96 # Process and validate command line arguments and options...
97 ProcessOptions()
98
99 # Perform actions required by the script...
100 PerformTorsionScan()
101
102 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
103 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
104
105
106 def PerformTorsionScan():
107 """Perform torsion scan."""
108
109 # Setup a molecule reader for input file...
110 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
111 OptionsInfo["InfileParams"]["AllowEmptyMols"] = True
112 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
113
114 PlotExt = OptionsInfo["OutPlotParams"]["OutExt"]
115 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
116 MiscUtil.PrintInfo(
117 "Generating output files %s_*.sdf, %s_*Torsion*Match*.sdf, %s_*Torsion*Match*Energies.csv, %s_*Torsion*Match*Plot.%s, %s_*Torsion*Match*Viewer.html..."
118 % (FileName, FileName, FileName, FileName, PlotExt, FileName)
119 )
120
121 MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount = ProcessMolecules(
122 Mols
123 )
124
125 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
126 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
127 MiscUtil.PrintInfo("Number of molecules failed during initial minimization: %d" % MinimizationFailedCount)
128 MiscUtil.PrintInfo("Number of molecules without any matched torsions: %d" % TorsionsMissingCount)
129 MiscUtil.PrintInfo("Number of molecules failed during torsion scan: %d" % TorsionsScanFailedCount)
130 MiscUtil.PrintInfo(
131 "Number of ignored molecules: %d"
132 % (MolCount - ValidMolCount + TorsionsMissingCount + MinimizationFailedCount + TorsionsScanFailedCount)
133 )
134
135
136 def ProcessMolecules(Mols):
137 """Process molecules to perform torsion scan."""
138
139 if OptionsInfo["MPMode"]:
140 return ProcessMoleculesUsingMultipleProcesses(Mols)
141 else:
142 return ProcessMoleculesUsingSingleProcess(Mols)
143
144
145 def ProcessMoleculesUsingSingleProcess(Mols):
146 """Process molecules to perform torsion scan using a single process."""
147
148 # Intialize Psi4...
149 MiscUtil.PrintInfo("\nInitializing Psi4...")
150 Psi4Handle = Psi4Util.InitializePsi4(
151 Psi4RunParams=OptionsInfo["Psi4RunParams"],
152 Psi4OptionsParams=OptionsInfo["Psi4OptionsParams"],
153 PrintVersion=True,
154 PrintHeader=True,
155 )
156 OptionsInfo["psi4"] = Psi4Handle
157
158 # Setup max iterations global variable...
159 Psi4Util.UpdatePsi4OptionsParameters(Psi4Handle, {"GEOM_MAXITER": OptionsInfo["MaxIters"]})
160
161 # Setup conversion factor for energy units...
162 SetupEnergyConversionFactor(Psi4Handle)
163
164 MolInfoText = "first molecule"
165 if not OptionsInfo["FirstMolMode"]:
166 MolInfoText = "all molecules"
167
168 if OptionsInfo["TorsionMinimize"]:
169 MiscUtil.PrintInfo(
170 "\nPerforming torsion scan on %s by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..."
171 % (MolInfoText)
172 )
173 else:
174 MiscUtil.PrintInfo(
175 "\nPerforming torsion scan on %s by skipping generation of conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..."
176 % (MolInfoText)
177 )
178
179 SetupTorsionsPatternsInfo()
180
181 (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
182
183 for Mol in Mols:
184 MolCount += 1
185
186 if OptionsInfo["FirstMolMode"] and MolCount > 1:
187 MolCount -= 1
188 break
189
190 if not CheckAndValidateMolecule(Mol, MolCount):
191 continue
192
193 # Setup 2D coordinates for SMILES input file...
194 if OptionsInfo["SMILESInfileStatus"]:
195 AllChem.Compute2DCoords(Mol)
196
197 ValidMolCount += 1
198
199 Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(
200 Mol, MolCount
201 )
202
203 if not MinimizationCalcStatus:
204 MinimizationFailedCount += 1
205 continue
206
207 if not TorsionsMatchStatus:
208 TorsionsMissingCount += 1
209 continue
210
211 if not TorsionsScanCalcStatus:
212 TorsionsScanFailedCount += 1
213 continue
214
215 return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
216
217
218 def ProcessMoleculesUsingMultipleProcesses(Mols):
219 """Process and minimize molecules using multiprocessing."""
220
221 if OptionsInfo["MPLevelTorsionAnglesMode"]:
222 return ProcessMoleculesUsingMultipleProcessesAtTorsionAnglesLevel(Mols)
223 elif OptionsInfo["MPLevelMoleculesMode"]:
224 return ProcessMoleculesUsingMultipleProcessesAtMoleculesLevel(Mols)
225 else:
226 MiscUtil.PrintError('The value, %s, option "--mpLevel" is not supported.' % (OptionsInfo["MPLevel"]))
227
228
229 def ProcessMoleculesUsingMultipleProcessesAtMoleculesLevel(Mols):
230 """Process molecules to perform torsion scan using multiprocessing at molecules level."""
231
232 MolInfoText = "first molecule"
233 if not OptionsInfo["FirstMolMode"]:
234 MolInfoText = "all molecules"
235
236 if OptionsInfo["TorsionMinimize"]:
237 MiscUtil.PrintInfo(
238 "\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..."
239 % (MolInfoText)
240 )
241 else:
242 MiscUtil.PrintInfo(
243 "\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..."
244 % (MolInfoText)
245 )
246
247 MPParams = OptionsInfo["MPParams"]
248
249 # Setup data for initializing a worker process...
250 InitializeWorkerProcessArgs = (
251 MiscUtil.ObjectToBase64EncodedString(Options),
252 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
253 )
254
255 if OptionsInfo["FirstMolMode"]:
256 Mol = Mols[0]
257 Mols = [Mol]
258
259 # Setup a encoded mols data iterable for a worker process...
260 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
261
262 # Setup process pool along with data initialization for each process...
263 MiscUtil.PrintInfo(
264 "\nConfiguring multiprocessing using %s method..."
265 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
266 )
267 MiscUtil.PrintInfo(
268 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
269 % (
270 MPParams["NumProcesses"],
271 MPParams["InputDataMode"],
272 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
273 )
274 )
275
276 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
277
278 # Start processing...
279 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
280 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
281 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
282 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
283 else:
284 MiscUtil.PrintError(
285 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
286 )
287
288 (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
289
290 for Result in Results:
291 MolCount += 1
292
293 MolIndex, EncodedMol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = Result
294
295 if EncodedMol is None:
296 continue
297 ValidMolCount += 1
298
299 if not MinimizationCalcStatus:
300 MinimizationFailedCount += 1
301 continue
302
303 if not TorsionsMatchStatus:
304 TorsionsMissingCount += 1
305 continue
306
307 if not TorsionsScanCalcStatus:
308 TorsionsScanFailedCount += 1
309 continue
310
311 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
312
313 return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
314
315
316 def InitializeWorkerProcess(*EncodedArgs):
317 """Initialize data for a worker process."""
318
319 global Options, OptionsInfo
320
321 if not OptionsInfo["QuietMode"]:
322 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
323
324 # Decode Options and OptionInfo...
325 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
326 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
327
328 # Initialize torsion patterns info...
329 SetupTorsionsPatternsInfo()
330
331 # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
332 # output files for each process...
333 OptionsInfo["Psi4Initialized"] = False
334
335
336 def WorkerProcess(EncodedMolInfo):
337 """Process data for a worker process."""
338
339 if not OptionsInfo["Psi4Initialized"]:
340 InitializePsi4ForWorkerProcess()
341
342 MolIndex, EncodedMol = EncodedMolInfo
343
344 MolNum = MolIndex + 1
345 (MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus) = [False] * 3
346
347 if EncodedMol is None:
348 return [MolIndex, None, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
349
350 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
351 if not CheckAndValidateMolecule(Mol, MolNum):
352 return [MolIndex, None, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
353
354 # Setup 2D coordinates for SMILES input file...
355 if OptionsInfo["SMILESInfileStatus"]:
356 AllChem.Compute2DCoords(Mol)
357
358 Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(
359 Mol, MolNum
360 )
361
362 return [
363 MolIndex,
364 RDKitUtil.MolToBase64EncodedMolString(
365 Mol, PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps
366 ),
367 MinimizationCalcStatus,
368 TorsionsMatchStatus,
369 TorsionsScanCalcStatus,
370 ]
371
372
373 def ProcessMoleculesUsingMultipleProcessesAtTorsionAnglesLevel(Mols):
374 """Process molecules to perform torsion scan using multiprocessing at torsion angles level."""
375
376 MolInfoText = "first molecule"
377 if not OptionsInfo["FirstMolMode"]:
378 MolInfoText = "all molecules"
379
380 if OptionsInfo["TorsionMinimize"]:
381 MiscUtil.PrintInfo(
382 "\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..."
383 % (MolInfoText)
384 )
385 else:
386 MiscUtil.PrintInfo(
387 "\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..."
388 % (MolInfoText)
389 )
390
391 SetupTorsionsPatternsInfo()
392
393 (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
394
395 for Mol in Mols:
396 MolCount += 1
397
398 if OptionsInfo["FirstMolMode"] and MolCount > 1:
399 MolCount -= 1
400 break
401
402 if not CheckAndValidateMolecule(Mol, MolCount):
403 continue
404
405 # Setup 2D coordinates for SMILES input file...
406 if OptionsInfo["SMILESInfileStatus"]:
407 AllChem.Compute2DCoords(Mol)
408
409 ValidMolCount += 1
410
411 Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(
412 Mol, MolCount, UseMultiProcessingAtTorsionAnglesLevel=True
413 )
414
415 if not MinimizationCalcStatus:
416 MinimizationFailedCount += 1
417 continue
418
419 if not TorsionsMatchStatus:
420 TorsionsMissingCount += 1
421 continue
422
423 if not TorsionsScanCalcStatus:
424 TorsionsScanFailedCount += 1
425 continue
426
427 return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
428
429
430 def ScanSingleTorsionInMolUsingMultipleProcessesAtTorsionAnglesLevel(
431 Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum
432 ):
433 """Perform torsion scan for a molecule using multiple processses at torsion angles
434 level along with constrained energy minimization.
435 """
436
437 if OptionsInfo["MPLevelMoleculesMode"]:
438 MiscUtil.PrintError(
439 "Single torison scanning for a molecule is not allowed in multiprocessing mode at molecules level.\n"
440 )
441
442 Mols, Angles = SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum)
443
444 MPParams = OptionsInfo["MPParams"]
445
446 # Setup data for initializing a worker process...
447
448 # Track and avoid encoding TorsionsPatternsInfo as it contains RDKit molecule object...
449 TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
450 OptionsInfo["TorsionsPatternsInfo"] = None
451
452 InitializeWorkerProcessArgs = (
453 MiscUtil.ObjectToBase64EncodedString(Options),
454 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
455 )
456
457 # Restore TorsionsPatternsInfo...
458 OptionsInfo["TorsionsPatternsInfo"] = TorsionsPatternsInfo
459
460 # Setup a encoded mols data iterable for a worker process...
461 WorkerProcessDataIterable = GenerateBase64EncodedMolStringsWithTorsionScanInfo(
462 Mol, (MolNum - 1), Mols, Angles, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches
463 )
464
465 # Setup process pool along with data initialization for each process...
466 MiscUtil.PrintInfo(
467 "\nConfiguring multiprocessing using %s method..."
468 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
469 )
470 MiscUtil.PrintInfo(
471 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
472 % (
473 MPParams["NumProcesses"],
474 MPParams["InputDataMode"],
475 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
476 )
477 )
478
479 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeTorsionAngleWorkerProcess, InitializeWorkerProcessArgs)
480
481 # Start processing...
482 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
483 Results = ProcessPool.imap(TorsionAngleWorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
484 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
485 Results = ProcessPool.map(TorsionAngleWorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
486 else:
487 MiscUtil.PrintError(
488 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
489 )
490
491 TorsionMols = []
492 TorsionEnergies = []
493 TorsionAngles = []
494
495 for Result in Results:
496 EncodedTorsionMol, CalcStatus, Angle, Energy = Result
497
498 if not CalcStatus:
499 return (Mol, False, None, None, None)
500
501 if EncodedTorsionMol is None:
502 return (Mol, False, None, None, None)
503 TorsionMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionMol)
504
505 TorsionMols.append(TorsionMol)
506 TorsionEnergies.append(Energy)
507 TorsionAngles.append(Angle)
508
509 return (Mol, True, TorsionMols, TorsionEnergies, TorsionAngles)
510
511
512 def InitializeTorsionAngleWorkerProcess(*EncodedArgs):
513 """Initialize data for a worker process."""
514
515 global Options, OptionsInfo
516
517 if not OptionsInfo["QuietMode"]:
518 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
519
520 # Decode Options and OptionInfo...
521 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
522 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
523
524 # Initialize torsion patterns info...
525 SetupTorsionsPatternsInfo()
526
527 # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
528 # output files for each process...
529 OptionsInfo["Psi4Initialized"] = False
530
531
532 def TorsionAngleWorkerProcess(EncodedMolInfo):
533 """Process data for a worker process."""
534
535 if not OptionsInfo["Psi4Initialized"]:
536 InitializePsi4ForWorkerProcess()
537
538 (
539 MolIndex,
540 EncodedMol,
541 EncodedTorsionMol,
542 TorsionAngle,
543 TorsionID,
544 TorsionPattern,
545 EncodedTorsionPatternMol,
546 TorsionMatches,
547 ) = EncodedMolInfo
548
549 if EncodedMol is None or EncodedTorsionMol is None or EncodedTorsionPatternMol is None:
550 return (None, False, None, None)
551
552 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
553 TorsionMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionMol)
554 TorsionPatternMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionPatternMol)
555
556 TorsionMol, CalcStatus, Energy = MinimizeCalculateEnergyForTorsionMol(
557 Mol, TorsionMol, TorsionAngle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, (MolIndex + 1)
558 )
559
560 return (
561 RDKitUtil.MolToBase64EncodedMolString(
562 TorsionMol,
563 PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps,
564 ),
565 CalcStatus,
566 TorsionAngle,
567 Energy,
568 )
569
570
571 def GenerateBase64EncodedMolStringsWithTorsionScanInfo(
572 Mol,
573 MolIndex,
574 TorsionMols,
575 TorsionAngles,
576 TorsionID,
577 TorsionPattern,
578 TorsionPatternMol,
579 TorsionMatches,
580 PropertyPickleFlags=Chem.PropertyPickleOptions.AllProps,
581 ):
582 """Set up an iterator for generating base64 encoded molecule string for
583 a torsion in a molecule along with appropriate trosion scan information.
584 """
585
586 for Index, TorsionMol in enumerate(TorsionMols):
587 yield (
588 [MolIndex, None, None, TorsionAngles[Index], TorsionID, TorsionPattern, None, TorsionMatches]
589 if (Mol is None or TorsionMol is None)
590 else [
591 MolIndex,
592 RDKitUtil.MolToBase64EncodedMolString(Mol, PropertyPickleFlags),
593 RDKitUtil.MolToBase64EncodedMolString(TorsionMol, PropertyPickleFlags),
594 TorsionAngles[Index],
595 TorsionID,
596 TorsionPattern,
597 RDKitUtil.MolToBase64EncodedMolString(TorsionPatternMol, PropertyPickleFlags),
598 TorsionMatches,
599 ]
600 )
601
602
603 def InitializePsi4ForWorkerProcess():
604 """Initialize Psi4 for a worker process."""
605
606 if OptionsInfo["Psi4Initialized"]:
607 return
608
609 OptionsInfo["Psi4Initialized"] = True
610
611 if OptionsInfo["MPLevelTorsionAnglesMode"] and re.match(
612 "auto", OptionsInfo["Psi4RunParams"]["OutputFileSpecified"], re.I
613 ):
614 # Run Psi4 in quiet mode during multiprocessing at Torsions level for 'auto' OutputFile...
615 OptionsInfo["Psi4RunParams"]["OutputFile"] = "quiet"
616 else:
617 # Update output file...
618 OptionsInfo["Psi4RunParams"]["OutputFile"] = Psi4Util.UpdatePsi4OutputFileUsingPID(
619 OptionsInfo["Psi4RunParams"]["OutputFile"], os.getpid()
620 )
621
622 # Intialize Psi4...
623 OptionsInfo["psi4"] = Psi4Util.InitializePsi4(
624 Psi4RunParams=OptionsInfo["Psi4RunParams"],
625 Psi4OptionsParams=OptionsInfo["Psi4OptionsParams"],
626 PrintVersion=False,
627 PrintHeader=True,
628 )
629
630 # Setup max iterations global variable...
631 Psi4Util.UpdatePsi4OptionsParameters(OptionsInfo["psi4"], {"GEOM_MAXITER": OptionsInfo["MaxIters"]})
632
633 # Setup conversion factor for energy units...
634 SetupEnergyConversionFactor(OptionsInfo["psi4"])
635
636
637 def PerformMinimizationAndTorsionScan(Mol, MolNum, UseMultiProcessingAtTorsionAnglesLevel=False):
638 """Perform minimization and torsions scan."""
639
640 if not OptionsInfo["Infile3D"]:
641 # Add hydrogens...
642 Mol = Chem.AddHs(Mol, addCoords=True)
643
644 Mol, MinimizationCalcStatus = MinimizeMolecule(Mol, MolNum)
645 if not MinimizationCalcStatus:
646 return (Mol, False, False, False)
647
648 TorsionsMolInfo = SetupTorsionsMolInfo(Mol, MolNum)
649 if TorsionsMolInfo["NumOfMatches"] == 0:
650 return (Mol, True, False, False)
651
652 Mol, ScanCalcStatus = ScanAllTorsionsInMol(Mol, TorsionsMolInfo, MolNum, UseMultiProcessingAtTorsionAnglesLevel)
653 if not ScanCalcStatus:
654 return (Mol, True, True, False)
655
656 return (Mol, True, True, True)
657
658
659 def ScanAllTorsionsInMol(Mol, TorsionsMolInfo, MolNum, UseMultiProcessingAtTorsionAnglesLevel=False):
660 """Perform scans on all torsions in a molecule."""
661
662 if TorsionsMolInfo["NumOfMatches"] == 0:
663 return Mol, True
664
665 MolName = RDKitUtil.GetMolName(Mol, MolNum)
666
667 FirstTorsionMode = OptionsInfo["FirstTorsionMode"]
668 TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
669
670 TorsionPatternCount, TorsionScanCount, TorsionMatchCount = [0] * 3
671 TorsionMaxMatches = OptionsInfo["TorsionMaxMatches"]
672
673 for TorsionID in TorsionsPatternsInfo["IDs"]:
674 TorsionPatternCount += 1
675 TorsionPattern = TorsionsPatternsInfo["Pattern"][TorsionID]
676 TorsionPatternMol = TorsionsPatternsInfo["Mol"][TorsionID]
677
678 TorsionsMatches = TorsionsMolInfo["Matches"][TorsionID]
679
680 if TorsionsMatches is None:
681 continue
682
683 if FirstTorsionMode and TorsionPatternCount > 1:
684 if not OptionsInfo["QuietMode"]:
685 MiscUtil.PrintWarning(
686 'Already scaned first torsion pattern, "%s" for molecule %s during "%s" value of "--modeTorsions" option . Abandoning torsion scan...\n'
687 % (TorsionPattern, MolName, OptionsInfo["ModeTorsions"])
688 )
689 break
690
691 for Index, TorsionMatches in enumerate(TorsionsMatches):
692 TorsionMatchNum = Index + 1
693 TorsionMatchCount += 1
694
695 if TorsionMatchCount > TorsionMaxMatches:
696 if not OptionsInfo["QuietMode"]:
697 MiscUtil.PrintWarning(
698 'Already scaned a maximum of %s torsion matches for molecule %s specified by "--torsionMaxMatches" option. Abandoning torsion scan...\n'
699 % (TorsionMaxMatches, MolName)
700 )
701 break
702
703 TmpMol, TorsionScanStatus, TorsionMols, TorsionEnergies, TorsionAngles = ScanSingleTorsionInMol(
704 Mol,
705 TorsionID,
706 TorsionPattern,
707 TorsionPatternMol,
708 TorsionMatches,
709 TorsionMatchNum,
710 MolNum,
711 UseMultiProcessingAtTorsionAnglesLevel,
712 )
713 if not TorsionScanStatus:
714 continue
715
716 TorsionScanCount += 1
717 GenerateOutputFiles(Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
718
719 if TorsionMatchCount > TorsionMaxMatches:
720 break
721
722 if TorsionScanCount:
723 GenerateStartingTorsionScanStructureOutfile(Mol, MolNum)
724
725 Status = True if TorsionScanCount else False
726
727 return (Mol, Status)
728
729
730 def ScanSingleTorsionInMol(
731 Mol,
732 TorsionID,
733 TorsionPattern,
734 TorsionPatternMol,
735 TorsionMatches,
736 TorsionMatchNum,
737 MolNum,
738 UseMultiProcessingAtTorsionAnglesLevel,
739 ):
740 """Perform torsion scan for a molecule along with constrained energy minimization."""
741
742 if not OptionsInfo["QuietMode"]:
743 MiscUtil.PrintInfo(
744 "\nProcessing torsion pattern, %s, match number, %s, in molecule %s..."
745 % (TorsionPattern, TorsionMatchNum, RDKitUtil.GetMolName(Mol, MolNum))
746 )
747
748 if OptionsInfo["TorsionMinimize"]:
749 MiscUtil.PrintInfo(
750 "Generating initial ensemble of constrained conformations by distance geometry and forcefield followed by Psi4 constraned minimization to select the lowest energy structure at specific torsion angles for molecule %s..."
751 % (RDKitUtil.GetMolName(Mol, MolNum))
752 )
753 else:
754 MiscUtil.PrintInfo(
755 "Calculating single point energy using Psi4 for molecule, %s, at specific torsion angles..."
756 % (RDKitUtil.GetMolName(Mol, MolNum))
757 )
758
759 if UseMultiProcessingAtTorsionAnglesLevel:
760 return ScanSingleTorsionInMolUsingMultipleProcessesAtTorsionAnglesLevel(
761 Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum
762 )
763 else:
764 return ScanSingleTorsionInMolUsingSingleProcess(
765 Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum
766 )
767
768
769 def ScanSingleTorsionInMolUsingSingleProcess(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum):
770 """Perform torsion scan for a molecule using single processs along with constrained
771 energy minimization."""
772
773 TorsionMols = []
774 TorsionEnergies = []
775 TorsionAngles = []
776
777 Mols, Angles = SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum)
778
779 for Index, Angle in enumerate(Angles):
780 TorsionMol = Mols[Index]
781 TorsionMol, CalcStatus, Energy = MinimizeCalculateEnergyForTorsionMol(
782 Mol, TorsionMol, Angle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum
783 )
784
785 if not CalcStatus:
786 return (Mol, False, None, None, None)
787
788 TorsionMols.append(TorsionMol)
789 TorsionEnergies.append(Energy)
790 TorsionAngles.append(Angle)
791
792 return (Mol, True, TorsionMols, TorsionEnergies, TorsionAngles)
793
794
795 def SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum=None):
796 """Setup molecules corresponding to all torsion angles in a molecule."""
797
798 AtomIndex1, AtomIndex2, AtomIndex3, AtomIndex4 = TorsionMatches
799
800 TorsionMols = []
801 TorsionAngles = OptionsInfo["TorsionAngles"]
802
803 for Angle in TorsionAngles:
804 TorsionMol = Chem.Mol(Mol)
805 TorsionMolConf = TorsionMol.GetConformer(0)
806
807 rdMolTransforms.SetDihedralDeg(TorsionMolConf, AtomIndex1, AtomIndex2, AtomIndex3, AtomIndex4, Angle)
808 TorsionMols.append(TorsionMol)
809
810 return (TorsionMols, TorsionAngles)
811
812
813 def MinimizeCalculateEnergyForTorsionMol(
814 Mol, TorsionMol, TorsionAngle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum
815 ):
816 """ "Calculate energy of a torsion molecule by performing an optional constrained
817 energy minimzation.
818 """
819
820 if OptionsInfo["TorsionMinimize"]:
821 if not OptionsInfo["QuietMode"]:
822 MolName = RDKitUtil.GetMolName(Mol, MolNum)
823 MiscUtil.PrintInfo("\nProcessing torsion angle %s for molecule %s..." % (TorsionAngle, MolName))
824
825 # Perform constrained minimization...
826 TorsionMatchesMol = RDKitUtil.MolFromSubstructureMatch(TorsionMol, TorsionPatternMol, TorsionMatches)
827 TorsionMol, CalcStatus, Energy = ConstrainAndMinimizeMolecule(
828 TorsionMol, TorsionAngle, TorsionMatchesMol, TorsionMatches, MolNum
829 )
830
831 if not CalcStatus:
832 if not OptionsInfo["QuietMode"]:
833 MolName = RDKitUtil.GetMolName(Mol, MolNum)
834 MiscUtil.PrintWarning(
835 "Failed to perform constrained minimization for molecule %s with torsion angle set to %s during torsion scan for torsion pattern %s. Abandoning torsion scan..."
836 % (MolName, TorsionAngle, TorsionPattern)
837 )
838 return (TorsionMol, False, None)
839 else:
840 # Setup a Psi4 molecule...
841 Psi4Mol = SetupPsi4Mol(OptionsInfo["psi4"], TorsionMol, MolNum)
842 if Psi4Mol is None:
843 return (TorsionMol, False, None)
844
845 # Calculate single point Psi4 energy...
846 CalcStatus, Energy = CalculateEnergyUsingPsi4(OptionsInfo["psi4"], Psi4Mol, TorsionMol, MolNum)
847 if not CalcStatus:
848 if not OptionsInfo["QuietMode"]:
849 MolName = RDKitUtil.GetMolName(Mol, MolNum)
850 MiscUtil.PrintWarning(
851 "Failed to calculate Psi4 energy for molecule %s with torsion angle set to %s during torsion scan for torsion pattern %s. Abandoning torsion scan..."
852 % (MolName, TorsionAngle, TorsionPattern)
853 )
854 return (TorsionMol, False, None)
855
856 return (TorsionMol, CalcStatus, Energy)
857
858
859 def SetupTorsionsMolInfo(Mol, MolNum=None):
860 """Setup torsions info for a molecule."""
861
862 TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
863
864 # Initialize...
865 TorsionsMolInfo = {}
866 TorsionsMolInfo["IDs"] = []
867 TorsionsMolInfo["NumOfMatches"] = 0
868 TorsionsMolInfo["Matches"] = {}
869 for TorsionID in TorsionsPatternsInfo["IDs"]:
870 TorsionsMolInfo["IDs"].append(TorsionID)
871 TorsionsMolInfo["Matches"][TorsionID] = None
872
873 MolName = RDKitUtil.GetMolName(Mol, MolNum)
874 UseChirality = OptionsInfo["UseChirality"]
875
876 for TorsionID in TorsionsPatternsInfo["IDs"]:
877 # Match torsions..
878 TorsionPattern = TorsionsPatternsInfo["Pattern"][TorsionID]
879 TorsionPatternMol = TorsionsPatternsInfo["Mol"][TorsionID]
880 TorsionsMatches = RDKitUtil.FilterSubstructureMatchesByAtomMapNumbers(
881 Mol, TorsionPatternMol, Mol.GetSubstructMatches(TorsionPatternMol, useChirality=UseChirality)
882 )
883
884 # Validate tosion matches...
885 ValidTorsionsMatches = []
886 for Index, TorsionMatch in enumerate(TorsionsMatches):
887 if len(TorsionMatch) != 4:
888 if not OptionsInfo["QuietMode"]:
889 MiscUtil.PrintWarning(
890 "Ignoring invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: It must match exactly 4 atoms."
891 % (TorsionMatch, TorsionPattern, MolName)
892 )
893 continue
894
895 if not RDKitUtil.AreAtomIndicesSequentiallyConnected(Mol, TorsionMatch):
896 if not OptionsInfo["QuietMode"]:
897 MiscUtil.PrintInfo("")
898 MiscUtil.PrintWarning(
899 "Invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices must be sequentially connected."
900 % (TorsionMatch, TorsionPattern, MolName)
901 )
902 MiscUtil.PrintWarning("Reordering matched atom indices in a sequentially connected manner...")
903
904 Status, ReorderdTorsionMatch = RDKitUtil.ReorderAtomIndicesInSequentiallyConnectedManner(
905 Mol, TorsionMatch
906 )
907 if Status:
908 TorsionMatch = ReorderdTorsionMatch
909 if not OptionsInfo["QuietMode"]:
910 MiscUtil.PrintWarning(
911 "Successfully reordered torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices are now sequentially connected."
912 % (TorsionMatch, TorsionPattern, MolName)
913 )
914 else:
915 if not OptionsInfo["QuietMode"]:
916 MiscUtil.PrintWarning(
917 "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."
918 % (TorsionMatch, TorsionPattern, MolName)
919 )
920 continue
921
922 Bond = Mol.GetBondBetweenAtoms(TorsionMatch[1], TorsionMatch[2])
923 if Bond.IsInRing():
924 if not OptionsInfo["QuietMode"]:
925 MiscUtil.PrintWarning(
926 "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."
927 % (TorsionMatch, TorsionPattern, MolName, TorsionMatch[1], TorsionMatch[2])
928 )
929 continue
930
931 # Filter matched torsions...
932 if OptionsInfo["FilterTorsionsByAtomIndicesMode"]:
933 InvalidAtomIndices = []
934 for AtomIndex in TorsionMatch:
935 if AtomIndex not in OptionsInfo["TorsionsFilterByAtomIndicesList"]:
936 InvalidAtomIndices.append(AtomIndex)
937 if len(InvalidAtomIndices):
938 if not OptionsInfo["QuietMode"]:
939 MiscUtil.PrintWarning(
940 '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".'
941 % (
942 TorsionMatch,
943 TorsionPattern,
944 MolName,
945 InvalidAtomIndices,
946 OptionsInfo["TorsionsFilterByAtomIndicesList"],
947 )
948 )
949 continue
950
951 ValidTorsionsMatches.append(TorsionMatch)
952
953 # Track valid matches...
954 if len(ValidTorsionsMatches):
955 TorsionsMolInfo["NumOfMatches"] += len(ValidTorsionsMatches)
956 TorsionsMolInfo["Matches"][TorsionID] = ValidTorsionsMatches
957
958 if TorsionsMolInfo["NumOfMatches"] == 0:
959 if not OptionsInfo["QuietMode"]:
960 MiscUtil.PrintWarning("Failed to match any torsions in molecule %s" % (MolName))
961
962 return TorsionsMolInfo
963
964
965 def SetupTorsionsPatternsInfo():
966 """Setup torsions patterns info."""
967
968 TorsionsPatternsInfo = {}
969 TorsionsPatternsInfo["IDs"] = []
970 TorsionsPatternsInfo["Pattern"] = {}
971 TorsionsPatternsInfo["Mol"] = {}
972
973 TorsionID = 0
974 for TorsionPattern in OptionsInfo["TorsionPatternsList"]:
975 TorsionID += 1
976
977 TorsionMol = Chem.MolFromSmarts(TorsionPattern)
978 if TorsionMol is None:
979 MiscUtil.PrintError(
980 'Failed to create torsion pattern molecule. The torsion SMILES/SMARTS pattern, "%s", specified using "-t, --torsions" option is not valid.'
981 % (TorsionPattern)
982 )
983
984 TorsionsPatternsInfo["IDs"].append(TorsionID)
985 TorsionsPatternsInfo["Pattern"][TorsionID] = TorsionPattern
986 TorsionsPatternsInfo["Mol"][TorsionID] = TorsionMol
987
988 OptionsInfo["TorsionsPatternsInfo"] = TorsionsPatternsInfo
989
990
991 def MinimizeMolecule(Mol, MolNum=None):
992 """Minimize molecule."""
993
994 return GenerateAndMinimizeConformersUsingForceField(Mol, MolNum)
995
996
997 def GenerateAndMinimizeConformersUsingForceField(Mol, MolNum=None):
998 """Generate and minimize conformers for a molecule to get the lowest energy conformer
999 as the minimized structure."""
1000
1001 MolName = RDKitUtil.GetMolName(Mol, MolNum)
1002
1003 # Setup conformers...
1004 ConfIDs = EmbedMolecule(Mol, MolNum)
1005 if not len(ConfIDs):
1006 if not OptionsInfo["QuietMode"]:
1007 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s: Embedding failed...\n" % MolName)
1008 return (Mol, False)
1009
1010 if not OptionsInfo["QuietMode"]:
1011 MiscUtil.PrintInfo(
1012 "Performing initial minimization of molecule, %s, using forcefield by generating a conformation ensemble and selecting the lowest energy conformer - EmbedRMSDCutoff: %s; Size: %s; Size after RMSD filtering: %s"
1013 % (
1014 MolName,
1015 OptionsInfo["ConfGenerationParams"]["EmbedRMSDCutoff"],
1016 OptionsInfo["ConfGenerationParams"]["MaxConfs"],
1017 len(ConfIDs),
1018 )
1019 )
1020
1021 # Minimize conformers...
1022 CalcEnergyMap = {}
1023 for ConfID in ConfIDs:
1024 # Perform forcefield minimization...
1025 Status, ConvergeStatus = MinimizeMoleculeUsingForceField(Mol, MolNum, ConfID)
1026 if not Status:
1027 return (Mol, False)
1028
1029 EnergyStatus, Energy = CalculateEnergyUsingForceField(Mol, ConfID)
1030 if not EnergyStatus:
1031 if not OptionsInfo["QuietMode"]:
1032 MolName = RDKitUtil.GetMolName(Mol, MolNum)
1033 MiscUtil.PrintWarning(
1034 "Failed to retrieve calculated energy for conformation number %d of molecule %s. Try again after removing any salts or cleaing up the molecule...\n"
1035 % (ConfID, MolName)
1036 )
1037 return (Mol, False)
1038
1039 if ConvergeStatus != 0:
1040 if not OptionsInfo["QuietMode"]:
1041 MiscUtil.PrintWarning(
1042 'Minimization using forcefield failed to converge for molecule %s in %d steps. Try using higher value for "maxIters" in "--confParams" option...\n'
1043 % (MolName, OptionsInfo["ConfGenerationParams"]["MaxIters"])
1044 )
1045
1046 CalcEnergyMap[ConfID] = Energy
1047
1048 SortedConfIDs = sorted(ConfIDs, key=lambda ConfID: CalcEnergyMap[ConfID])
1049 MinEnergyConfID = SortedConfIDs[0]
1050
1051 for ConfID in [Conf.GetId() for Conf in Mol.GetConformers()]:
1052 if ConfID == MinEnergyConfID:
1053 continue
1054 Mol.RemoveConformer(ConfID)
1055
1056 # Set ConfID to 0 for MinEnergyConf...
1057 Mol.GetConformer(MinEnergyConfID).SetId(0)
1058
1059 return (Mol, True)
1060
1061
1062 def ConstrainAndMinimizeMolecule(Mol, TorsionAngle, RefMolCore, RefMolMatches, MolNum=None):
1063 """Constrain and minimize molecule."""
1064
1065 # TorsionMol, CalcStatus, Energy
1066 MolName = RDKitUtil.GetMolName(Mol, MolNum)
1067
1068 # Setup constrained conformers...
1069 MolConfs, MolConfsStatus = ConstrainEmbedAndMinimizeMoleculeUsingRDKit(Mol, RefMolCore, RefMolMatches, MolNum)
1070 if not MolConfsStatus:
1071 if not OptionsInfo["QuietMode"]:
1072 MiscUtil.PrintWarning(
1073 "Conformation generation couldn't be performed for molecule %s: Constrained embedding failed...\n"
1074 % MolName
1075 )
1076 return (Mol, False, None)
1077
1078 # Minimize conformers...
1079 ConfNums = []
1080 CalcEnergyMap = {}
1081 MolConfsMap = {}
1082
1083 for ConfNum, MolConf in enumerate(MolConfs):
1084 if not OptionsInfo["QuietMode"]:
1085 MiscUtil.PrintInfo(
1086 "\nPerforming constrained minimization using Psi4 for molecule, %s, conformer number, %s, at torsion angle %s..."
1087 % (MolName, ConfNum, TorsionAngle)
1088 )
1089
1090 CalcStatus, Energy = ConstrainAndMinimizeMoleculeUsingPsi4(
1091 OptionsInfo["psi4"], MolConf, RefMolCore, RefMolMatches, MolNum
1092 )
1093 if not CalcStatus:
1094 if not OptionsInfo["QuietMode"]:
1095 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s\n" % (MolName))
1096 return (Mol, False, None)
1097
1098 ConfNums.append(ConfNum)
1099 CalcEnergyMap[ConfNum] = Energy
1100 MolConfsMap[ConfNum] = MolConf
1101
1102 SortedConfNums = sorted(ConfNums, key=lambda ConfNum: CalcEnergyMap[ConfNum])
1103 MinEnergyConfNum = SortedConfNums[0]
1104
1105 MinEnergy = CalcEnergyMap[MinEnergyConfNum]
1106 MinEnergyMolConf = MolConfsMap[MinEnergyConfNum]
1107
1108 MinEnergyMolConf.ClearProp("EmbedRMS")
1109
1110 return (MinEnergyMolConf, True, MinEnergy)
1111
1112
1113 def ConstrainAndMinimizeMoleculeUsingPsi4(Psi4Handle, Mol, RefMolCore, RefMolMatches, MolNum, ConfID=-1):
1114 """Minimize molecule using Psi4."""
1115
1116 # Setup a list for constrained atoms...
1117 ConstrainedAtomIndices = SetupConstrainedAtomIndicesForPsi4(Mol, RefMolCore, RefMolMatches)
1118 if len(ConstrainedAtomIndices) == 0:
1119 return (False, None)
1120
1121 # Setup a Psi4Mol...
1122 Psi4Mol = SetupPsi4Mol(Psi4Handle, Mol, MolNum, ConfID)
1123 if Psi4Mol is None:
1124 return (False, None)
1125
1126 # Setup reference wave function...
1127 Reference = SetupReferenceWavefunction(Mol)
1128 Psi4Handle.set_options({"Reference": Reference})
1129
1130 # Setup method name and basis set...
1131 MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
1132
1133 # Setup freeze list for constrained torsion...
1134 FreezeDihedralList = [("%s" % AtomIdex) for AtomIdex in ConstrainedAtomIndices]
1135 FreezeDihedralString = "%s" % " ".join(FreezeDihedralList)
1136 Psi4Handle.set_options({"OPTKING__frozen_dihedral": FreezeDihedralString})
1137
1138 # Optimize geometry...
1139 Status, Energy, WaveFunction = Psi4Util.PerformGeometryOptimization(
1140 Psi4Handle, Psi4Mol, MethodName, BasisSet, ReturnWaveFunction=True, Quiet=OptionsInfo["QuietMode"]
1141 )
1142
1143 if not Status:
1144 PerformPsi4Cleanup(Psi4Handle)
1145 return (False, None)
1146
1147 # Update atom positions...
1148 AtomPositions = Psi4Util.GetAtomPositions(Psi4Handle, WaveFunction, InAngstroms=True)
1149 RDKitUtil.SetAtomPositions(Mol, AtomPositions, ConfID=ConfID)
1150
1151 # Convert energy units...
1152 if OptionsInfo["ApplyEnergyConversionFactor"]:
1153 Energy = Energy * OptionsInfo["EnergyConversionFactor"]
1154
1155 # Clean up
1156 PerformPsi4Cleanup(Psi4Handle)
1157
1158 return (True, Energy)
1159
1160
1161 def ConstrainEmbedAndMinimizeMoleculeUsingRDKit(Mol, RefMolCore, RefMolMatches, MolNum=None):
1162 """Constrain, embed, and minimize molecule."""
1163
1164 # Setup forcefield function to use for constrained minimization...
1165 ForceFieldFunction = None
1166 ForceFieldName = None
1167 if OptionsInfo["ConfGenerationParams"]["UseUFF"]:
1168 ForceFieldFunction = lambda mol, confId=-1: AllChem.UFFGetMoleculeForceField(mol, confId=confId)
1169 ForceFieldName = "UFF"
1170 else:
1171 ForceFieldFunction = lambda mol, confId=-1: AllChem.MMFFGetMoleculeForceField(
1172 mol,
1173 AllChem.MMFFGetMoleculeProperties(
1174 mol, mmffVariant=OptionsInfo["ConfGenerationParams"]["ForceFieldMMFFVariant"]
1175 ),
1176 confId=confId,
1177 )
1178 ForceFieldName = "MMFF"
1179
1180 if ForceFieldFunction is None:
1181 if not OptionsInfo["QuietMode"]:
1182 MiscUtil.PrintWarning(
1183 "Failed to setup forcefield %s for molecule: %s\n" % (ForceFieldName, RDKitUtil.GetMolName(Mol, MolNum))
1184 )
1185 return (None, False)
1186
1187 MaxConfs = OptionsInfo["ConfGenerationParams"]["MaxConfsTorsions"]
1188 EnforceChirality = OptionsInfo["ConfGenerationParams"]["EnforceChirality"]
1189 UseExpTorsionAnglePrefs = OptionsInfo["ConfGenerationParams"]["UseExpTorsionAnglePrefs"]
1190 ETVersion = OptionsInfo["ConfGenerationParams"]["ETVersion"]
1191 UseBasicKnowledge = OptionsInfo["ConfGenerationParams"]["UseBasicKnowledge"]
1192 UseTethers = OptionsInfo["ConfGenerationParams"]["UseTethers"]
1193
1194 MolConfs = []
1195 ConfIDs = [ConfID for ConfID in range(0, MaxConfs)]
1196
1197 for ConfID in ConfIDs:
1198 try:
1199 MolConf = Chem.Mol(Mol)
1200 RDKitUtil.ConstrainAndEmbed(
1201 MolConf,
1202 RefMolCore,
1203 coreMatchesMol=RefMolMatches,
1204 useTethers=UseTethers,
1205 coreConfId=-1,
1206 randomseed=ConfID,
1207 getForceField=ForceFieldFunction,
1208 enforceChirality=EnforceChirality,
1209 useExpTorsionAnglePrefs=UseExpTorsionAnglePrefs,
1210 useBasicKnowledge=UseBasicKnowledge,
1211 ETversion=ETVersion,
1212 )
1213 except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
1214 if not OptionsInfo["QuietMode"]:
1215 MiscUtil.PrintWarning(
1216 "Constrained embedding couldn't be performed for molecule %s:\n%s\n"
1217 % (RDKitUtil.GetMolName(Mol, MolNum), ErrMsg)
1218 )
1219 return (None, False)
1220
1221 MolConfs.append(MolConf)
1222
1223 return FilterConstrainedConformationsByRMSD(Mol, MolConfs, MolNum)
1224
1225
1226 def FilterConstrainedConformationsByRMSD(Mol, MolConfs, MolNum=None):
1227 """Filter contarained conformations by RMSD."""
1228
1229 EmbedRMSDCutoff = OptionsInfo["ConfGenerationParams"]["EmbedRMSDCutoff"]
1230 if EmbedRMSDCutoff is None or EmbedRMSDCutoff <= 0:
1231 if not OptionsInfo["QuietMode"]:
1232 MiscUtil.PrintInfo(
1233 "\nGenerating initial ensemble of constrained conformations by distance geometry and forcefield for %s - EmbedRMSDCutoff: None; Size: %s"
1234 % (RDKitUtil.GetMolName(Mol, MolNum), len(MolConfs))
1235 )
1236 return (MolConfs, True)
1237
1238 FirstMolConf = True
1239 SelectedMolConfs = []
1240 for MolConfIndex, MolConf in enumerate(MolConfs):
1241 if FirstMolConf:
1242 FirstMolConf = False
1243 SelectedMolConfs.append(MolConf)
1244 continue
1245
1246 # Compare RMSD against all selected conformers...
1247 ProbeMolConf = Chem.RemoveHs(Chem.Mol(MolConf))
1248 IgnoreConf = False
1249 for SelectedMolConfIndex, SelectedMolConf in enumerate(SelectedMolConfs):
1250 RefMolConf = Chem.RemoveHs(Chem.Mol(SelectedMolConf))
1251 CalcRMSD = rdMolAlign.AlignMol(ProbeMolConf, RefMolConf)
1252
1253 if CalcRMSD < EmbedRMSDCutoff:
1254 IgnoreConf = True
1255 break
1256
1257 if IgnoreConf:
1258 continue
1259
1260 SelectedMolConfs.append(MolConf)
1261
1262 if not OptionsInfo["QuietMode"]:
1263 MiscUtil.PrintInfo(
1264 "\nGenerating initial ensemble of constrained conformations by distance geometry and forcefield for %s - EmbedRMSDCutoff: %s; Size: %s; Size after RMSD filtering: %s"
1265 % (RDKitUtil.GetMolName(Mol, MolNum), EmbedRMSDCutoff, len(MolConfs), len(SelectedMolConfs))
1266 )
1267
1268 return (SelectedMolConfs, True)
1269
1270
1271 def EmbedMolecule(Mol, MolNum=None):
1272 """Embed conformations."""
1273
1274 ConfIDs = []
1275
1276 MaxConfs = OptionsInfo["ConfGenerationParams"]["MaxConfs"]
1277 RandomSeed = OptionsInfo["ConfGenerationParams"]["RandomSeed"]
1278 EnforceChirality = OptionsInfo["ConfGenerationParams"]["EnforceChirality"]
1279 UseExpTorsionAnglePrefs = OptionsInfo["ConfGenerationParams"]["UseExpTorsionAnglePrefs"]
1280 ETVersion = OptionsInfo["ConfGenerationParams"]["ETVersion"]
1281 UseBasicKnowledge = OptionsInfo["ConfGenerationParams"]["UseBasicKnowledge"]
1282 EmbedRMSDCutoff = OptionsInfo["ConfGenerationParams"]["EmbedRMSDCutoff"]
1283
1284 try:
1285 ConfIDs = AllChem.EmbedMultipleConfs(
1286 Mol,
1287 numConfs=MaxConfs,
1288 randomSeed=RandomSeed,
1289 pruneRmsThresh=EmbedRMSDCutoff,
1290 enforceChirality=EnforceChirality,
1291 useExpTorsionAnglePrefs=UseExpTorsionAnglePrefs,
1292 useBasicKnowledge=UseBasicKnowledge,
1293 ETversion=ETVersion,
1294 )
1295 except ValueError as ErrMsg:
1296 if not OptionsInfo["QuietMode"]:
1297 MolName = RDKitUtil.GetMolName(Mol, MolNum)
1298 MiscUtil.PrintWarning("Embedding failed for molecule %s:\n%s\n" % (MolName, ErrMsg))
1299 ConfIDs = []
1300
1301 return ConfIDs
1302
1303
1304 def MinimizeMoleculeUsingForceField(Mol, MolNum, ConfID=-1):
1305 """Minimize molecule using forcefield available in RDKit."""
1306
1307 try:
1308 if OptionsInfo["ConfGenerationParams"]["UseUFF"]:
1309 ConvergeStatus = AllChem.UFFOptimizeMolecule(
1310 Mol, confId=ConfID, maxIters=OptionsInfo["ConfGenerationParams"]["MaxIters"]
1311 )
1312 elif OptionsInfo["ConfGenerationParams"]["UseMMFF"]:
1313 ConvergeStatus = AllChem.MMFFOptimizeMolecule(
1314 Mol,
1315 confId=ConfID,
1316 maxIters=OptionsInfo["ConfGenerationParams"]["MaxIters"],
1317 mmffVariant=OptionsInfo["ConfGenerationParams"]["ForceFieldMMFFVariant"],
1318 )
1319 else:
1320 MiscUtil.PrintError(
1321 "Minimization couldn't be performed: Specified forcefield, %s, is not supported"
1322 % OptionsInfo["ConfGenerationParams"]["ForceField"]
1323 )
1324 except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
1325 if not OptionsInfo["QuietMode"]:
1326 MolName = RDKitUtil.GetMolName(Mol, MolNum)
1327 MiscUtil.PrintWarning(
1328 "Minimization using forcefield couldn't be performed for molecule %s:\n%s\n" % (MolName, ErrMsg)
1329 )
1330 return (False, None)
1331
1332 return (True, ConvergeStatus)
1333
1334
1335 def CalculateEnergyUsingForceField(Mol, ConfID=None):
1336 """Calculate energy."""
1337
1338 Status = True
1339 Energy = None
1340
1341 if ConfID is None:
1342 ConfID = -1
1343
1344 if OptionsInfo["ConfGenerationParams"]["UseUFF"]:
1345 UFFMoleculeForcefield = AllChem.UFFGetMoleculeForceField(Mol, confId=ConfID)
1346 if UFFMoleculeForcefield is None:
1347 Status = False
1348 else:
1349 Energy = UFFMoleculeForcefield.CalcEnergy()
1350 elif OptionsInfo["ConfGenerationParams"]["UseMMFF"]:
1351 MMFFMoleculeProperties = AllChem.MMFFGetMoleculeProperties(
1352 Mol, mmffVariant=OptionsInfo["ConfGenerationParams"]["ForceFieldMMFFVariant"]
1353 )
1354 MMFFMoleculeForcefield = AllChem.MMFFGetMoleculeForceField(Mol, MMFFMoleculeProperties, confId=ConfID)
1355 if MMFFMoleculeForcefield is None:
1356 Status = False
1357 else:
1358 Energy = MMFFMoleculeForcefield.CalcEnergy()
1359 else:
1360 MiscUtil.PrintError(
1361 "Couldn't retrieve conformer energy: Specified forcefield, %s, is not supported"
1362 % OptionsInfo["ConfGenerationParams"]["ForceField"]
1363 )
1364
1365 return (Status, Energy)
1366
1367
1368 def CalculateEnergyUsingPsi4(Psi4Handle, Psi4Mol, Mol, MolNum=None):
1369 """Calculate single point energy using Psi4."""
1370
1371 Status = False
1372 Energy = None
1373
1374 # Setup reference wave function...
1375 Reference = SetupReferenceWavefunction(Mol)
1376 Psi4Handle.set_options({"Reference": Reference})
1377
1378 # Setup method name and basis set...
1379 MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
1380
1381 Status, Energy = Psi4Util.CalculateSinglePointEnergy(
1382 Psi4Handle, Psi4Mol, MethodName, BasisSet, Quiet=OptionsInfo["QuietMode"]
1383 )
1384
1385 # Convert energy units...
1386 if Status:
1387 if OptionsInfo["ApplyEnergyConversionFactor"]:
1388 Energy = Energy * OptionsInfo["EnergyConversionFactor"]
1389
1390 # Clean up
1391 PerformPsi4Cleanup(Psi4Handle)
1392
1393 return (Status, Energy)
1394
1395
1396 def SetupConstrainedAtomIndicesForPsi4(Mol, RefMolCore, RefMolMatches):
1397 """Setup a list of atom indices to be constrained during Psi4 minimizaiton."""
1398
1399 AtomIndices = []
1400
1401 if RefMolMatches is None:
1402 return AtomIndices
1403 else:
1404 ConstrainAtomIndices = RefMolMatches
1405
1406 # Atom indices start from 1 for Psi4 instead 0 for RDKit...
1407 AtomIndices = [AtomIndex + 1 for AtomIndex in ConstrainAtomIndices]
1408
1409 return AtomIndices
1410
1411
1412 def SetupPsi4Mol(Psi4Handle, Mol, MolNum, ConfID=-1):
1413 """Setup a Psi4 molecule object."""
1414
1415 # Turn off recentering and reorientation to perform optimization in the
1416 # constrained coordinate space...
1417 MolGeometry = RDKitUtil.GetPsi4XYZFormatString(Mol, ConfID=ConfID, NoCom=True, NoReorient=True)
1418
1419 try:
1420 Psi4Mol = Psi4Handle.geometry(MolGeometry)
1421 except Exception as ErrMsg:
1422 Psi4Mol = None
1423 if not OptionsInfo["QuietMode"]:
1424 MiscUtil.PrintWarning("Failed to create Psi4 molecule from geometry string: %s\n" % ErrMsg)
1425 MolName = RDKitUtil.GetMolName(Mol, MolNum)
1426 MiscUtil.PrintWarning("Ignoring molecule: %s" % MolName)
1427
1428 return Psi4Mol
1429
1430
1431 def PerformPsi4Cleanup(Psi4Handle):
1432 """Perform clean up."""
1433
1434 # Clean up after Psi4 run...
1435 Psi4Handle.core.clean()
1436
1437 # Clean up any leftover scratch files...
1438 if OptionsInfo["MPMode"]:
1439 Psi4Util.RemoveScratchFiles(Psi4Handle, OptionsInfo["Psi4RunParams"]["OutputFile"])
1440
1441
1442 def CheckAndValidateMolecule(Mol, MolCount=None):
1443 """Validate molecule for Psi4 calculations."""
1444
1445 if Mol is None:
1446 if not OptionsInfo["QuietMode"]:
1447 MiscUtil.PrintInfo("\nProcessing molecule number %s..." % MolCount)
1448 return False
1449
1450 MolName = RDKitUtil.GetMolName(Mol, MolCount)
1451 if not OptionsInfo["QuietMode"]:
1452 MiscUtil.PrintInfo("\nProcessing molecule %s..." % MolName)
1453
1454 if RDKitUtil.IsMolEmpty(Mol):
1455 if not OptionsInfo["QuietMode"]:
1456 MiscUtil.PrintWarning("Ignoring empty molecule: %s\n" % MolName)
1457 return False
1458
1459 if not RDKitUtil.ValidateElementSymbols(RDKitUtil.GetAtomSymbols(Mol)):
1460 if not OptionsInfo["QuietMode"]:
1461 MiscUtil.PrintWarning("Ignoring molecule containing invalid element symbols: %s\n" % MolName)
1462 return False
1463
1464 if OptionsInfo["Infile3D"]:
1465 if not Mol.GetConformer().Is3D():
1466 if not OptionsInfo["QuietMode"]:
1467 MiscUtil.PrintWarning("3D tag is not set for molecule: %s\n" % MolName)
1468
1469 if OptionsInfo["Infile3D"]:
1470 # Otherwise, Hydrogens are always added...
1471 if RDKitUtil.AreHydrogensMissingInMolecule(Mol):
1472 if not OptionsInfo["QuietMode"]:
1473 MiscUtil.PrintWarning("Missing hydrogens in molecule: %s\n" % MolName)
1474
1475 return True
1476
1477
1478 def SetupMethodNameAndBasisSet(Mol):
1479 """Setup method name and basis set."""
1480
1481 MethodName = OptionsInfo["MethodName"]
1482 if OptionsInfo["MethodNameAuto"]:
1483 MethodName = "B3LYP"
1484
1485 BasisSet = OptionsInfo["BasisSet"]
1486 if OptionsInfo["BasisSetAuto"]:
1487 BasisSet = "6-31+G**" if RDKitUtil.IsAtomSymbolPresentInMol(Mol, "S") else "6-31G**"
1488
1489 return (MethodName, BasisSet)
1490
1491
1492 def SetupReferenceWavefunction(Mol):
1493 """Setup reference wavefunction."""
1494
1495 Reference = OptionsInfo["Reference"]
1496 if OptionsInfo["ReferenceAuto"]:
1497 Reference = "UHF" if (RDKitUtil.GetSpinMultiplicity(Mol) > 1) else "RHF"
1498
1499 return Reference
1500
1501
1502 def SetupEnergyConversionFactor(Psi4Handle):
1503 """Setup converstion factor for energt units. The Psi4 energy units are Hartrees."""
1504
1505 EnergyUnits = OptionsInfo["EnergyUnits"]
1506
1507 ApplyConversionFactor = True
1508 if re.match(r"^kcal\/mol$", EnergyUnits, re.I):
1509 ConversionFactor = Psi4Handle.constants.hartree2kcalmol
1510 elif re.match(r"^kJ\/mol$", EnergyUnits, re.I):
1511 ConversionFactor = Psi4Handle.constants.hartree2kJmol
1512 elif re.match("^eV$", EnergyUnits, re.I):
1513 ConversionFactor = Psi4Handle.constants.hartree2ev
1514 else:
1515 ApplyConversionFactor = False
1516 ConversionFactor = 1.0
1517
1518 OptionsInfo["ApplyEnergyConversionFactor"] = ApplyConversionFactor
1519 OptionsInfo["EnergyConversionFactor"] = ConversionFactor
1520
1521
1522 def GenerateStartingTorsionScanStructureOutfile(Mol, MolNum):
1523 """Write out the structure of molecule used for starting tosion scan."""
1524
1525 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
1526 MolName = GetOutputFileMolName(Mol, MolNum)
1527
1528 Outfile = "%s_%s.%s" % (FileName, MolName, FileExt)
1529
1530 # Set up a molecule writer...
1531 Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
1532 if Writer is None:
1533 MiscUtil.PrintWarning("Failed to setup a writer for output fie %s " % Outfile)
1534 return
1535
1536 Writer.write(Mol)
1537
1538 if Writer is not None:
1539 Writer.close()
1540
1541
1542 def GenerateOutputFiles(Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
1543 """Generate output files."""
1544
1545 StructureOutfile, EnergyTextOutfile, PlotOutfile, ViewerHTMLOutfile = SetupOutputFileNames(
1546 Mol, MolNum, TorsionID, TorsionMatchNum
1547 )
1548
1549 GenerateScannedTorsionsStructureOutfile(
1550 StructureOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles
1551 )
1552 GenerateEnergyTextOutfile(
1553 EnergyTextOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles
1554 )
1555 GeneratePlotOutfile(
1556 PlotOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles
1557 )
1558 GenerateViewerHTMLOutfile(
1559 ViewerHTMLOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles
1560 )
1561
1562
1563 def GenerateScannedTorsionsStructureOutfile(
1564 Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles
1565 ):
1566 """Write out structures generated after torsion scan along with associated data."""
1567
1568 # Set up a molecule writer...
1569 Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
1570 if Writer is None:
1571 MiscUtil.PrintWarning("Failed to setup a writer for output fie %s " % Outfile)
1572 return
1573
1574 MolName = RDKitUtil.GetMolName(Mol, MolNum)
1575
1576 RelativeTorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
1577 for Index, TorsionMol in enumerate(TorsionMols):
1578 TorsionAngle = "%s" % TorsionAngles[Index]
1579 TorsionMol.SetProp("Torsion_Angle", TorsionAngle)
1580
1581 TorsionEnergy = "%.*f" % (OptionsInfo["Precision"], TorsionEnergies[Index])
1582 TorsionMol.SetProp(OptionsInfo["EnergyDataFieldLabel"], TorsionEnergy)
1583
1584 RelativeTorsionEnergy = "%.*f" % (OptionsInfo["Precision"], RelativeTorsionEnergies[Index])
1585 TorsionMol.SetProp(OptionsInfo["EnergyRelativeDataFieldLabel"], RelativeTorsionEnergy)
1586
1587 TorsionMolName = "%s_Deg%s" % (MolName, TorsionAngle)
1588 TorsionMol.SetProp("_Name", TorsionMolName)
1589
1590 Writer.write(TorsionMol)
1591
1592 if Writer is not None:
1593 Writer.close()
1594
1595
1596 def GenerateEnergyTextOutfile(
1597 Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles
1598 ):
1599 """Write out torsion angles and energies."""
1600
1601 # Setup a writer...
1602 Writer = open(Outfile, "w")
1603 if Writer is None:
1604 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
1605
1606 # Write headers...
1607 Writer.write(
1608 "TorsionAngle,%s,%s\n" % (OptionsInfo["EnergyDataFieldLabel"], OptionsInfo["EnergyRelativeDataFieldLabel"])
1609 )
1610
1611 RelativeTorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
1612 for Index, TorsionAngle in enumerate(TorsionAngles):
1613 TorsionEnergy = "%.*f" % (OptionsInfo["Precision"], TorsionEnergies[Index])
1614 RelativeTorsionEnergy = "%.*f" % (OptionsInfo["Precision"], RelativeTorsionEnergies[Index])
1615 Writer.write("%d,%s,%s\n" % (TorsionAngle, TorsionEnergy, RelativeTorsionEnergy))
1616
1617 if Writer is not None:
1618 Writer.close()
1619
1620
1621 def GenerateViewerHTMLOutfile(
1622 Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles
1623 ):
1624 """Write out a HTML file for viewing torsion scan along with associated data."""
1625
1626 # Setup a writer...
1627 Writer = open(Outfile, "w")
1628 if Writer is None:
1629 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
1630
1631 # Setup HTML for torsion scan viewer...
1632 MolName = RDKitUtil.GetMolName(Mol, MolNum)
1633 RelativeTorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
1634
1635 Units = OptionsInfo["EnergyUnits"]
1636 PlotHeight = OptionsInfo["OutPlotTorsionViewerHeight"]
1637
1638 MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
1639 TitleLine = "%s (%s/%s)" % (OptionsInfo["OutPlotParams"]["Title"], MethodName, BasisSet)
1640 if OptionsInfo["OutPlotTitleTorsionSpec"]:
1641 TorsionPattern = OptionsInfo["TorsionsPatternsInfo"]["Pattern"][TorsionID]
1642 TitleLine = "%s: %s" % (TitleLine, TorsionPattern)
1643
1644 ViwerHTMLText = RDKitUtil.SetupHTMLForTorsionScanViewer(
1645 MolName, TorsionMols, TorsionEnergies, RelativeTorsionEnergies, TorsionAngles, Units, PlotHeight, TitleLine
1646 )
1647
1648 # Write out HTML for torsion scan viewer...
1649 Writer.write("%s" % ViwerHTMLText)
1650
1651 Writer.close()
1652
1653
1654 def GeneratePlotOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
1655 """Generate a plot corresponding to torsion angles and energies."""
1656
1657 OutPlotParams = OptionsInfo["OutPlotParams"]
1658
1659 # Initialize seaborn and matplotlib paramaters...
1660 if not OptionsInfo["OutPlotInitialized"]:
1661 OptionsInfo["OutPlotInitialized"] = True
1662 RCParams = {
1663 "figure.figsize": (OutPlotParams["Width"], OutPlotParams["Height"]),
1664 "axes.titleweight": OutPlotParams["TitleWeight"],
1665 "axes.labelweight": OutPlotParams["LabelWeight"],
1666 }
1667 sns.set(
1668 context=OutPlotParams["Context"],
1669 style=OutPlotParams["Style"],
1670 palette=OutPlotParams["Palette"],
1671 font=OutPlotParams["Font"],
1672 font_scale=OutPlotParams["FontScale"],
1673 rc=RCParams,
1674 )
1675
1676 # Create a new figure...
1677 plt.figure()
1678
1679 if OptionsInfo["OutPlotRelativeEnergy"]:
1680 TorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
1681
1682 # Draw plot...
1683 PlotType = OutPlotParams["Type"]
1684 if re.match("linepoint", PlotType, re.I):
1685 Axis = sns.lineplot(x=TorsionAngles, y=TorsionEnergies, marker="o", legend=False)
1686 elif re.match("scatter", PlotType, re.I):
1687 Axis = sns.scatterplot(x=TorsionAngles, y=TorsionEnergies, legend=False)
1688 elif re.match("line", PlotType, re.I):
1689 Axis = sns.lineplot(x=TorsionAngles, y=TorsionEnergies, legend=False)
1690 else:
1691 MiscUtil.PrintError(
1692 'The value, %s, specified for "type" using option "--outPlotParams" is not supported. Valid plot types: linepoint, scatter or line'
1693 % (PlotType)
1694 )
1695
1696 # Setup title and labels...
1697 Title = OutPlotParams["Title"]
1698 if OptionsInfo["OutPlotTitleTorsionSpec"]:
1699 TorsionPattern = OptionsInfo["TorsionsPatternsInfo"]["Pattern"][TorsionID]
1700 Title = "%s: %s" % (OutPlotParams["Title"], TorsionPattern)
1701
1702 # Set labels and title...
1703 Axis.set(xlabel=OutPlotParams["XLabel"], ylabel=OutPlotParams["YLabel"], title=Title)
1704
1705 # Save figure...
1706 plt.savefig(Outfile)
1707
1708 # Close the plot...
1709 plt.close()
1710
1711
1712 def SetupRelativeEnergies(Energies):
1713 """Set up a list of relative energies."""
1714
1715 SortedEnergies = sorted(Energies)
1716 MinEnergy = SortedEnergies[0]
1717 RelativeEnergies = [(Energy - MinEnergy) for Energy in Energies]
1718
1719 return RelativeEnergies
1720
1721
1722 def SetupOutputFileNames(Mol, MolNum, TorsionID, TorsionMatchNum):
1723 """Setup names of output files."""
1724
1725 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
1726 MolName = GetOutputFileMolName(Mol, MolNum)
1727
1728 OutfileRoot = "%s_%s_Torsion%s_Match%s" % (FileName, MolName, TorsionID, TorsionMatchNum)
1729
1730 StructureOutfile = "%s.%s" % (OutfileRoot, FileExt)
1731 EnergyTextOutfile = "%s_Energies.csv" % (OutfileRoot)
1732
1733 PlotExt = OptionsInfo["OutPlotParams"]["OutExt"]
1734 PlotOutfile = "%s_Plot.%s" % (OutfileRoot, PlotExt)
1735
1736 ViewerHTMLOutfile = "%s_Viewer.html" % (OutfileRoot)
1737
1738 return (StructureOutfile, EnergyTextOutfile, PlotOutfile, ViewerHTMLOutfile)
1739
1740
1741 def GetOutputFileMolName(Mol, MolNum):
1742 """Get output file prefix."""
1743
1744 MolName = "Mol%s" % MolNum
1745 if OptionsInfo["OutfileMolName"]:
1746 MolName = re.sub("[^a-zA-Z0-9]", "_", RDKitUtil.GetMolName(Mol, MolNum), flags=re.I)
1747
1748 return MolName
1749
1750
1751 def ProcessTorsionRangeOptions():
1752 """Process tosion range options."""
1753
1754 TosionRangeMode = Options["--torsionRangeMode"]
1755 OptionsInfo["TosionRangeMode"] = TosionRangeMode
1756
1757 if re.match("^Range$", TosionRangeMode, re.I):
1758 ProcessTorsionRangeValues()
1759 elif re.match("^Angles$", TosionRangeMode, re.I):
1760 ProcessTorsionAnglesValues()
1761 else:
1762 MiscUtil.PrintError('The value, %s, option "--torsionRangeMode" is not supported.' % TosionRangeMode)
1763
1764
1765 def ProcessTorsionRangeValues():
1766 """Process tosion range values."""
1767
1768 TorsionRange = Options["--torsionRange"]
1769 if re.match("^auto$", TorsionRange, re.I):
1770 TorsionRange = "0,360,5"
1771 TorsionRangeWords = TorsionRange.split(",")
1772
1773 TorsionStart = int(TorsionRangeWords[0])
1774 TorsionStop = int(TorsionRangeWords[1])
1775 TorsionStep = int(TorsionRangeWords[2])
1776
1777 if TorsionStart >= TorsionStop:
1778 MiscUtil.PrintError(
1779 'The start value, %d, specified for option "--torsionRange" in string "%s" must be less than stop value, %s.'
1780 % (TorsionStart, Options["--torsionRange"], TorsionStop)
1781 )
1782 if TorsionStep == 0:
1783 MiscUtil.PrintError(
1784 'The step value, %d, specified for option "--torsonRange" in string "%s" must be > 0.'
1785 % (TorsionStep, Options["--torsionRange"])
1786 )
1787 if TorsionStep >= (TorsionStop - TorsionStart):
1788 MiscUtil.PrintError(
1789 'The step value, %d, specified for option "--torsonRange" in string "%s" must be less than, %s.'
1790 % (TorsionStep, Options["--torsionRange"], (TorsionStop - TorsionStart))
1791 )
1792
1793 if TorsionStart < 0:
1794 if TorsionStart < -180:
1795 MiscUtil.PrintError(
1796 'The start value, %d, specified for option "--torsionRange" in string "%s" must be >= -180 to use scan range from -180 to 180.'
1797 % (TorsionStart, Options["--torsionRange"])
1798 )
1799 if TorsionStop > 180:
1800 MiscUtil.PrintError(
1801 'The stop value, %d, specified for option "--torsionRange" in string "%s" must be <= 180 to use scan range from -180 to 180.'
1802 % (TorsionStop, Options["--torsionRange"])
1803 )
1804 else:
1805 if TorsionStop > 360:
1806 MiscUtil.PrintError(
1807 'The stop value, %d, specified for option "--torsionRange" in string "%s" must be <= 360 to use scan range from 0 to 360.'
1808 % (TorsionStop, Options["--torsionRange"])
1809 )
1810
1811 TorsionAngles = [Angle for Angle in range(TorsionStart, TorsionStop, TorsionStep)]
1812 TorsionAngles.append(TorsionStop)
1813
1814 OptionsInfo["TorsionRange"] = TorsionRange
1815 OptionsInfo["TorsionStart"] = TorsionStart
1816 OptionsInfo["TorsionStop"] = TorsionStop
1817 OptionsInfo["TorsionStep"] = TorsionStep
1818
1819 OptionsInfo["TorsionAngles"] = TorsionAngles
1820
1821
1822 def ProcessTorsionAnglesValues():
1823 """Process tosion angle values."""
1824
1825 TorsionRange = Options["--torsionRange"]
1826 if re.match("^auto$", TorsionRange, re.I):
1827 MiscUtil.PrintError('The value specified, %s, for option "--torsionRange" is not valid.' % (TorsionRange))
1828
1829 TorsionAngles = []
1830
1831 for TorsionAngle in TorsionRange.split(","):
1832 TorsionAngle = int(TorsionAngle)
1833
1834 if TorsionAngle < -180:
1835 MiscUtil.PrintError(
1836 'The stop value, %d, specified for option "--torsionRange" in string "%s" must be >= -180.'
1837 % (TorsionAngle, TorsionRange)
1838 )
1839
1840 if TorsionAngle > 360:
1841 MiscUtil.PrintError(
1842 'The stop value, %d, specified for option "--torsionRange" in string "%s" must be <= 360.'
1843 % (TorsionAngle, TorsionRange)
1844 )
1845
1846 TorsionAngles.append(TorsionAngle)
1847
1848 OptionsInfo["TorsionRange"] = TorsionRange
1849 OptionsInfo["TorsionStart"] = None
1850 OptionsInfo["TorsionStop"] = None
1851 OptionsInfo["TorsionStep"] = None
1852
1853 OptionsInfo["TorsionAngles"] = sorted(TorsionAngles)
1854
1855
1856 def ProcessOptions():
1857 """Process and validate command line arguments and options."""
1858
1859 MiscUtil.PrintInfo("Processing options...")
1860
1861 # Validate options...
1862 ValidateOptions()
1863
1864 OptionsInfo["ModeMols"] = Options["--modeMols"]
1865 OptionsInfo["FirstMolMode"] = True if re.match("^First$", Options["--modeMols"], re.I) else False
1866
1867 OptionsInfo["ModeTorsions"] = Options["--modeTorsions"]
1868 OptionsInfo["FirstTorsionMode"] = True if re.match("^First$", Options["--modeTorsions"], re.I) else False
1869
1870 OptionsInfo["Infile"] = Options["--infile"]
1871 OptionsInfo["SMILESInfileStatus"] = True if MiscUtil.CheckFileExt(Options["--infile"], "smi csv tsv txt") else False
1872 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
1873 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
1874 "--infileParams",
1875 Options["--infileParams"],
1876 InfileName=Options["--infile"],
1877 ParamsDefaultInfo=ParamsDefaultInfoOverride,
1878 )
1879 OptionsInfo["Infile3D"] = True if re.match("^yes$", Options["--infile3D"], re.I) else False
1880
1881 OptionsInfo["Outfile"] = Options["--outfile"]
1882 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters(
1883 "--outfileParams", Options["--outfileParams"]
1884 )
1885
1886 # Method, basis set, and reference wavefunction...
1887 OptionsInfo["BasisSet"] = Options["--basisSet"]
1888 OptionsInfo["BasisSetAuto"] = True if re.match("^auto$", Options["--basisSet"], re.I) else False
1889
1890 OptionsInfo["MethodName"] = Options["--methodName"]
1891 OptionsInfo["MethodNameAuto"] = True if re.match("^auto$", Options["--methodName"], re.I) else False
1892
1893 OptionsInfo["Reference"] = Options["--reference"]
1894 OptionsInfo["ReferenceAuto"] = True if re.match("^auto$", Options["--reference"], re.I) else False
1895
1896 # Run and options parameters...
1897 OptionsInfo["Psi4OptionsParams"] = Psi4Util.ProcessPsi4OptionsParameters(
1898 "--psi4OptionsParams", Options["--psi4OptionsParams"]
1899 )
1900 OptionsInfo["Psi4RunParams"] = Psi4Util.ProcessPsi4RunParameters(
1901 "--psi4RunParams", Options["--psi4RunParams"], InfileName=OptionsInfo["Infile"]
1902 )
1903
1904 # Conformer generation paramaters...
1905 ParamsDefaultInfoOverride = {"MaxConfs": 250, "MaxConfsTorsions": 50}
1906 OptionsInfo["ConfGenerationParams"] = MiscUtil.ProcessOptionConformerParameters(
1907 "--confParams", Options["--confParams"], ParamsDefaultInfoOverride
1908 )
1909
1910 # Energy units and label...
1911 OptionsInfo["EnergyUnits"] = Options["--energyUnits"]
1912
1913 EnergyDataFieldLabel = Options["--energyDataFieldLabel"]
1914 if re.match("^auto$", EnergyDataFieldLabel, re.I):
1915 EnergyDataFieldLabel = "Psi4_Energy (%s)" % Options["--energyUnits"]
1916 OptionsInfo["EnergyDataFieldLabel"] = EnergyDataFieldLabel
1917
1918 EnergyRelativeDataFieldLabel = Options["--energyRelativeDataFieldLabel"]
1919 if re.match("^auto$", EnergyRelativeDataFieldLabel, re.I):
1920 EnergyRelativeDataFieldLabel = "Psi4_Relative_Energy (%s)" % Options["--energyUnits"]
1921 OptionsInfo["EnergyRelativeDataFieldLabel"] = EnergyRelativeDataFieldLabel
1922
1923 # Plot parameters...
1924 OptionsInfo["OutfileMolName"] = True if re.match("^yes$", Options["--outfileMolName"], re.I) else False
1925 OptionsInfo["OutPlotRelativeEnergy"] = (
1926 True if re.match("^yes$", Options["--outPlotRelativeEnergy"], re.I) else False
1927 )
1928 OptionsInfo["OutPlotTitleTorsionSpec"] = (
1929 True if re.match("^yes$", Options["--outPlotTitleTorsionSpec"], re.I) else False
1930 )
1931 OptionsInfo["OutPlotTorsionViewerHeight"] = int(Options["--outPlotTorsionViewerHeight"])
1932
1933 # The default width and height, 10.0 and 7.5, map to aspect raito of 16/9 (1.778)...
1934 if OptionsInfo["OutPlotRelativeEnergy"]:
1935 EnergyLabel = "Relative Energy (%s)" % OptionsInfo["EnergyUnits"]
1936 else:
1937 EnergyLabel = "Energy (%s)" % OptionsInfo["EnergyUnits"]
1938
1939 DefaultValues = {
1940 "Type": "linepoint",
1941 "Width": 10.0,
1942 "Height": 5.6,
1943 "Title": "Psi4 Torsion Scan",
1944 "XLabel": "Torsion Angle (degrees)",
1945 "YLabel": EnergyLabel,
1946 }
1947 OptionsInfo["OutPlotParams"] = MiscUtil.ProcessOptionSeabornPlotParameters(
1948 "--outPlotParams", Options["--outPlotParams"], DefaultValues
1949 )
1950 if not re.match("^(linepoint|scatter|Line)$", OptionsInfo["OutPlotParams"]["Type"], re.I):
1951 MiscUtil.PrintError(
1952 'The value, %s, specified for "type" using option "--outPlotParams" is not supported. Valid plot types: linepoint, scatter or line'
1953 % (OptionsInfo["OutPlotParams"]["Type"])
1954 )
1955
1956 OptionsInfo["OutPlotInitialized"] = False
1957
1958 OptionsInfo["Overwrite"] = Options["--overwrite"]
1959
1960 OptionsInfo["MaxIters"] = int(Options["--maxIters"])
1961
1962 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
1963 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
1964
1965 # Multiprocessing level...
1966 MPLevelMoleculesMode = False
1967 MPLevelTorsionAnglesMode = False
1968 MPLevel = Options["--mpLevel"]
1969 if re.match("^Molecules$", MPLevel, re.I):
1970 MPLevelMoleculesMode = True
1971 elif re.match("^TorsionAngles$", MPLevel, re.I):
1972 MPLevelTorsionAnglesMode = True
1973 else:
1974 MiscUtil.PrintError('The value, %s, specified for option "--mpLevel" is not valid. ' % MPLevel)
1975 OptionsInfo["MPLevel"] = MPLevel
1976 OptionsInfo["MPLevelMoleculesMode"] = MPLevelMoleculesMode
1977 OptionsInfo["MPLevelTorsionAnglesMode"] = MPLevelTorsionAnglesMode
1978
1979 OptionsInfo["Precision"] = int(Options["--precision"])
1980 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
1981
1982 # Procsss and validate specified SMILES/SMARTS torsion patterns...
1983 TorsionPatterns = Options["--torsions"]
1984 TorsionPatternsList = []
1985 for TorsionPattern in TorsionPatterns.split(","):
1986 TorsionPattern = TorsionPattern.strip()
1987 if not len(TorsionPattern):
1988 MiscUtil.PrintError(
1989 'Empty value specified for SMILES/SMARTS pattern in "-t, --torsions" option: %s' % TorsionPatterns
1990 )
1991
1992 TorsionMol = Chem.MolFromSmarts(TorsionPattern)
1993 if TorsionMol is None:
1994 MiscUtil.PrintError(
1995 'Failed to create torsion pattern molecule. The torsion SMILES/SMARTS pattern, "%s", specified using "-t, --torsions" option, "%s", is not valid.'
1996 % (TorsionPattern, TorsionPatterns)
1997 )
1998 TorsionPatternsList.append(TorsionPattern)
1999
2000 OptionsInfo["TorsionPatterns"] = TorsionPatterns
2001 OptionsInfo["TorsionPatternsList"] = TorsionPatternsList
2002
2003 # Process and validate any specified torsion atom indices for filtering torsion matches...
2004 TorsionsFilterByAtomIndices = Options["--torsionsFilterbyAtomIndices"]
2005 TorsionsFilterByAtomIndicesList = []
2006 if not re.match("^None$", TorsionsFilterByAtomIndices, re.I):
2007 for AtomIndex in TorsionsFilterByAtomIndices.split(","):
2008 AtomIndex = AtomIndex.strip()
2009 if not MiscUtil.IsInteger(AtomIndex):
2010 MiscUtil.PrintError(
2011 'The value specified, %s, for option "--torsionsFilterbyAtomIndices" must be an integer.'
2012 % AtomIndex
2013 )
2014 AtomIndex = int(AtomIndex)
2015 if AtomIndex < 0:
2016 MiscUtil.PrintError(
2017 'The value specified, %s, for option "--torsionsFilterbyAtomIndices" must be >= 0.' % AtomIndex
2018 )
2019 TorsionsFilterByAtomIndicesList.append(AtomIndex)
2020
2021 if len(TorsionsFilterByAtomIndicesList) < 4:
2022 MiscUtil.PrintError(
2023 'The number of values, %s, specified, %s, for option "--torsionsFilterbyAtomIndices" must be >=4.'
2024 % (len(TorsionsFilterByAtomIndicesList), TorsionsFilterByAtomIndices)
2025 )
2026
2027 OptionsInfo["TorsionsFilterByAtomIndices"] = TorsionsFilterByAtomIndices
2028 OptionsInfo["TorsionsFilterByAtomIndicesList"] = TorsionsFilterByAtomIndicesList
2029 OptionsInfo["FilterTorsionsByAtomIndicesMode"] = True if len(TorsionsFilterByAtomIndicesList) > 0 else False
2030
2031 OptionsInfo["TorsionMaxMatches"] = int(Options["--torsionMaxMatches"])
2032 OptionsInfo["TorsionMinimize"] = True if re.match("^yes$", Options["--torsionMinimize"], re.I) else False
2033
2034 ProcessTorsionRangeOptions()
2035
2036 OptionsInfo["UseChirality"] = True if re.match("^yes$", Options["--useChirality"], re.I) else False
2037
2038
2039 def RetrieveOptions():
2040 """Retrieve command line arguments and options."""
2041
2042 # Get options...
2043 global Options
2044 Options = docopt(_docoptUsage_)
2045
2046 # Set current working directory to the specified directory...
2047 WorkingDir = Options["--workingdir"]
2048 if WorkingDir:
2049 os.chdir(WorkingDir)
2050
2051 # Handle examples option...
2052 if "--examples" in Options and Options["--examples"]:
2053 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
2054 sys.exit(0)
2055
2056
2057 def ValidateOptions():
2058 """Validate option values."""
2059
2060 MiscUtil.ValidateOptionTextValue("--energyUnits", Options["--energyUnits"], "Hartrees kcal/mol kJ/mol eV")
2061
2062 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
2063 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi txt csv tsv")
2064 MiscUtil.ValidateOptionTextValue("--infile3D", Options["--infile3D"], "yes no")
2065
2066 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
2067 MiscUtil.ValidateOptionsOutputFileOverwrite(
2068 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"]
2069 )
2070 MiscUtil.ValidateOptionsDistinctFileNames(
2071 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"]
2072 )
2073
2074 if not Options["--overwrite"]:
2075 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
2076 FileNames = glob.glob("%s_*" % FileName)
2077 if len(FileNames):
2078 MiscUtil.PrintError(
2079 'The outfile names, %s_*, generated from file specified, %s, for option "-o, --outfile" already exist. Use option "--overwrite" or "--ov" and try again.\n'
2080 % (FileName, Options["--outfile"])
2081 )
2082
2083 MiscUtil.ValidateOptionTextValue("--outPlotRelativeEnergy", Options["--outPlotRelativeEnergy"], "yes no")
2084 MiscUtil.ValidateOptionTextValue("--outPlotTitleTorsionSpec", Options["--outPlotTitleTorsionSpec"], "yes no")
2085 MiscUtil.ValidateOptionIntegerValue(
2086 "--outPlotTorsionViewerHeight", Options["--outPlotTorsionViewerHeight"], {">": 0}
2087 )
2088
2089 MiscUtil.ValidateOptionTextValue("--outfileMolName ", Options["--outfileMolName"], "yes no")
2090
2091 MiscUtil.ValidateOptionTextValue("--modeMols", Options["--modeMols"], "First All")
2092 MiscUtil.ValidateOptionTextValue("--modeTorsions", Options["--modeTorsions"], "First All")
2093
2094 MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0})
2095
2096 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
2097 MiscUtil.ValidateOptionTextValue("--mpLevel", Options["--mpLevel"], "Molecules TorsionAngles")
2098
2099 MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0})
2100 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
2101
2102 MiscUtil.ValidateOptionIntegerValue("--torsionMaxMatches", Options["--torsionMaxMatches"], {">": 0})
2103 MiscUtil.ValidateOptionTextValue("--torsionMinimize", Options["--torsionMinimize"], "yes no")
2104
2105 MiscUtil.ValidateOptionTextValue("--torsionRangeMode", Options["--torsionRangeMode"], "Range or Angles")
2106 TorsionRange = Options["--torsionRange"]
2107 if re.match("^Range$", Options["--torsionRangeMode"], re.I):
2108 if not re.match("^auto$", TorsionRange, re.I):
2109 MiscUtil.ValidateOptionNumberValues("--torsionRange", TorsionRange, 3, ",", "integer", {})
2110 else:
2111 if re.match("^auto$", TorsionRange, re.I):
2112 MiscUtil.PrintError(
2113 '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.'
2114 % (TorsionRange, Options["--torsionRangeMode"])
2115 )
2116 TorsionAngles = []
2117 for TorsionAngle in TorsionRange.split(","):
2118 TorsionAngle = TorsionAngle.strip()
2119 if not MiscUtil.IsInteger(TorsionAngle):
2120 MiscUtil.PrintError(
2121 'The value specified, %s, for option "--torsionRange" in string "%s" must be an integer.'
2122 % (TorsionAngle, TorsionRange)
2123 )
2124 if TorsionAngle in TorsionAngles:
2125 MiscUtil.PrintError(
2126 'The value specified, %s, for option "--torsionRange" in string "%s" is a duplicate value.'
2127 % (TorsionAngle, TorsionRange)
2128 )
2129 TorsionAngles.append(TorsionAngle)
2130
2131 MiscUtil.ValidateOptionTextValue("--useChirality", Options["--useChirality"], "yes no")
2132
2133
2134 # Setup a usage string for docopt...
2135 _docoptUsage_ = """
2136 Psi4PerformTorsionScan.py - Perform torsion scan
2137
2138 Usage:
2139 Psi4PerformTorsionScan.py [--basisSet <text>] [--confParams <Name,Value,...>] [--energyDataFieldLabel <text>]
2140 [--energyRelativeDataFieldLabel <text>] [--energyUnits <text>] [--infile3D <yes or no>]
2141 [--infileParams <Name,Value,...>] [--maxIters <number>] [--methodName <text>]
2142 [--modeMols <First or All>] [--modeTorsions <First or All>] [--mp <yes or no>]
2143 [--mpLevel <Molecules or TorsionAngles>] [--mpParams <Name,Value,...>]
2144 [--outfileMolName <yes or no>] [--outfileParams <Name,Value,...>] [--outPlotParams <Name,Value,...>]
2145 [--outPlotRelativeEnergy <yes or no>] [--outPlotTitleTorsionSpec <yes or no>] [--outPlotTorsionViewerHeight <number>]
2146 [--overwrite] [--precision <number>] [--psi4OptionsParams <Name,Value,...>] [--psi4RunParams <Name,Value,...>]
2147 [--quiet <yes or no>] [--reference <text>] [--torsionsFilterbyAtomIndices <Index1, Index2, ...>]
2148 [--torsionMaxMatches <number>] [--torsionMinimize <yes or no>] [--torsionRangeMode <Range or Angles>]
2149 [--torsionRange <Start,Stop,Step or Angle1,Angle2,...>] [--useChirality <yes or no>]
2150 [-w <dir>] -t <torsions> -i <infile> -o <outfile>
2151 Psi4PerformTorsionScan.py -h | --help | -e | --examples
2152
2153 Description:
2154 Perform torsion scan for molecules around torsion angles specified using
2155 SMILES/SMARTS patterns. A molecule is optionally minimized before performing
2156 a torsion scan using a forcefield. A set of initial 3D structures are generated for
2157 a molecule by scanning the torsion angle across the specified range and updating
2158 the 3D coordinates of the molecule. A conformation ensemble is optionally generated
2159 for each 3D structure representing a specific torsion angle using a combination of
2160 distance geometry and forcefield followed by constrained geometry optimization
2161 using a quantum chemistry method. The conformation with the lowest energy is
2162 selected to represent the torsion angle. An option is available to skip the generation
2163 of the conformation ensemble and simply calculate the energy for the initial 3D
2164 structure for a specific torsion torsion angle using a quantum chemistry method.
2165
2166 The torsions are specified using SMILES or SMARTS patterns. A substructure match
2167 is performed to select torsion atoms in a molecule. The SMILES pattern match must
2168 correspond to four torsion atoms. The SMARTS patterns containing atom map numbers
2169 may match more than four atoms. The atom map numbers, however, must match
2170 exactly four torsion atoms. For example: [s:1][c:2]([aX2,cH1])!@[CX3:3](O)=[O:4] for
2171 thiophene esters and carboxylates as specified in Torsion Library (TorLib) [Ref 146].
2172
2173 A Psi4 XYZ format geometry string is automatically generated for each molecule
2174 in input file. It contains atom symbols and 3D coordinates for each atom in a
2175 molecule. In addition, the formal charge and spin multiplicity are present in the
2176 the geometry string. These values are either retrieved from molecule properties
2177 named 'FormalCharge' and 'SpinMultiplicty' or dynamically calculated for a
2178 molecule.
2179
2180 A set of five output files is generated for each torsion match in each
2181 molecule. The names of the output files are generated using the root of
2182 the specified output file. They may either contain sequential molecule
2183 numbers or molecule names as shown below:
2184
2185 <OutfileRoot>_Mol<Num>.sdf
2186 <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>.sdf
2187 <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Energies.csv
2188 <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Plot.<ImgExt>
2189 <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Viewer.html
2190
2191 or
2192
2193 <OutfileRoot>_<MolName>.sdf
2194 <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>.sdf
2195 <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Energies.csv
2196 <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Plot.<ImgExt>
2197 <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Viewer.html
2198
2199 The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi,
2200 .csv, .tsv, .txt)
2201
2202 The supported output file formats are: SD (.sdf, .sd)
2203
2204 Options:
2205 -b, --basisSet <text> [default: auto]
2206 Basis set to use for energy calculation or constrained energy minimization.
2207 Default: 6-31+G** for sulfur containing molecules; Otherwise, 6-31G** [ Ref 150 ].
2208 The specified value must be a valid Psi4 basis set. No validation is performed.
2209
2210 The following list shows a representative sample of basis sets available
2211 in Psi4:
2212
2213 STO-3G, 6-31G, 6-31+G, 6-31++G, 6-31G*, 6-31+G*, 6-31++G*,
2214 6-31G**, 6-31+G**, 6-31++G**, 6-311G, 6-311+G, 6-311++G,
2215 6-311G*, 6-311+G*, 6-311++G*, 6-311G**, 6-311+G**, 6-311++G**,
2216 cc-pVDZ, cc-pCVDZ, aug-cc-pVDZ, cc-pVDZ-DK, cc-pCVDZ-DK, def2-SVP,
2217 def2-SVPD, def2-TZVP, def2-TZVPD, def2-TZVPP, def2-TZVPPD
2218
2219 --confParams <Name,Value,...> [default: auto]
2220 A comma delimited list of parameter name and value pairs for generating
2221 initial 3D coordinates for molecules in input file at specific torsion angles. A
2222 conformation ensemble is optionally generated for each 3D structure
2223 representing a specific torsion angle using a combination of distance geometry
2224 and forcefield followed by constrained geometry optimization using a quantum
2225 chemistry method. The conformation with the lowest energy is selected to
2226 represent the torsion angle.
2227
2228 The supported parameter names along with their default values are shown
2229 below:
2230
2231 confMethod,ETKDGv2,
2232 forceField,MMFF, forceFieldMMFFVariant,MMFF94,
2233 enforceChirality,yes,embedRMSDCutoff,0.5,maxConfs,250,
2234 maxConfsTorsions,50,useTethers,yes
2235
2236 confMethod,ETKDGv2 [ Possible values: SDG, KDG, ETDG,
2237 ETKDG , or ETKDGv2]
2238 forceField, MMFF [ Possible values: UFF or MMFF ]
2239 forceFieldMMFFVariant,MMFF94 [ Possible values: MMFF94 or MMFF94s ]
2240 enforceChirality,yes [ Possible values: yes or no ]
2241 useTethers,yes [ Possible values: yes or no ]
2242
2243 confMethod: Conformation generation methodology for generating initial 3D
2244 coordinates. Possible values: Standard Distance Geometry (SDG), Experimental
2245 Torsion-angle preference with Distance Geometry (ETDG), basic Knowledge-terms
2246 with Distance Geometry (KDG) and Experimental Torsion-angle preference
2247 along with basic Knowledge-terms and Distance Geometry (ETKDG or
2248 ETKDGv2) [Ref 129, 167] .
2249
2250 forceField: Forcefield method to use for energy minimization. Possible values:
2251 Universal Force Field (UFF) [ Ref 81 ] or Merck Molecular Mechanics Force
2252 Field [ Ref 83-87 ] .
2253
2254 enforceChirality: Enforce chirality for defined chiral centers during
2255 forcefield minimization.
2256
2257 maxConfs: Maximum number of conformations to generate for each molecule
2258 during the generation of an initial 3D conformation ensemble using a conformation
2259 generation methodology. The conformations are minimized using the specified
2260 forcefield. The lowest energy structure is selected for performing the torsion scan.
2261
2262 maxConfsTorsion: Maximum number of 3D conformations to generate for
2263 conformation ensemble representing a specific torsion. The conformations are
2264 constrained at specific torsions angles and minimized using the specified forcefield
2265 and a quantum chemistry method. The lowest energy conformation is selected to
2266 calculate final torsion energy and written to the output file.
2267
2268 embedRMSDCutoff: RMSD cutoff for retaining initial set of conformers embedded
2269 using distance geometry and forcefield minimization. All embedded conformers
2270 are kept for 'None' value. Otherwise, only those conformers which are different
2271 from each other by the specified RMSD cutoff, 0.5 by default, are kept. The first
2272 embedded conformer is always retained.
2273
2274 useTethers: Use tethers to optimize the final embedded conformation by
2275 applying a series of extra forces to align matching atoms to the positions of
2276 the core atoms. Otherwise, use simple distance constraints during the
2277 optimization.
2278 --energyDataFieldLabel <text> [default: auto]
2279 Energy data field label for writing energy values. Default: Psi4_Energy (<Units>).
2280 --energyRelativeDataFieldLabel <text> [default: auto]
2281 Relative energy data field label for writing energy values. Default:
2282 Psi4_Relative_Energy (<Units>).
2283 --energyUnits <text> [default: kcal/mol]
2284 Energy units. Possible values: Hartrees, kcal/mol, kJ/mol, or eV.
2285 -e, --examples
2286 Print examples.
2287 -h, --help
2288 Print this help message.
2289 -i, --infile <infile>
2290 Input file name.
2291 --infile3D <yes or no> [default: no]
2292 Skip generation and minimization of initial 3D structures for molecules in
2293 input file containing 3D coordinates.
2294 --infileParams <Name,Value,...> [default: auto]
2295 A comma delimited list of parameter name and value pairs for reading
2296 molecules from files. The supported parameter names for different file
2297 formats, along with their default values, are shown below:
2298
2299 SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
2300
2301 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
2302 smilesTitleLine,auto,sanitize,yes
2303
2304 Possible values for smilesDelimiter: space, comma or tab.
2305 --maxIters <number> [default: 50]
2306 Maximum number of iterations to perform for each molecule or conformer
2307 during constrained energy minimization by a quantum chemistry method.
2308 -m, --methodName <text> [default: auto]
2309 Method to use for energy calculation or constrained energy minimization.
2310 Default: B3LYP [ Ref 150 ]. The specified value must be a valid Psi4 method
2311 name. No validation is performed.
2312
2313 The following list shows a representative sample of methods available
2314 in Psi4:
2315
2316 B1LYP, B2PLYP, B2PLYP-D3BJ, B2PLYP-D3MBJ, B3LYP, B3LYP-D3BJ,
2317 B3LYP-D3MBJ, CAM-B3LYP, CAM-B3LYP-D3BJ, HF, HF-D3BJ, HF3c, M05,
2318 M06, M06-2x, M06-HF, M06-L, MN12-L, MN15, MN15-D3BJ,PBE, PBE0,
2319 PBEH3c, PW6B95, PW6B95-D3BJ, WB97, WB97X, WB97X-D, WB97X-D3BJ
2320
2321 --modeMols <First or All> [default: First]
2322 Perform torsion scan for the first molecule or all molecules in input
2323 file.
2324 --modeTorsions <First or All> [default: First]
2325 Perform torsion scan for the first or all specified torsion pattern in
2326 molecules up to a maximum number of matches for each torsion
2327 specification as indicated by '--torsionMaxMatches' option.
2328 --mp <yes or no> [default: no]
2329 Use multiprocessing.
2330
2331 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
2332 function employing lazy RDKit data iterable. This allows processing of
2333 arbitrary large data sets without any additional requirements memory.
2334
2335 All input data may be optionally loaded into memory by mp.Pool.map()
2336 before starting worker processes in a process pool by setting the value
2337 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
2338
2339 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
2340 data mode may adversely impact the performance. The '--mpParams' section
2341 provides additional information to tune the value of 'chunkSize'.
2342 --mpLevel <Molecules or TorsionAngles> [default: Molecules]
2343 Perform multiprocessing at molecules or torsion angles level. Possible values:
2344 Molecules or TorsionAngles. The 'Molecules' value starts a process pool at the
2345 molecules level. All torsion angles of a molecule are processed in a single
2346 process. The 'TorsionAngles' value, however, starts a process pool at the
2347 torsion angles level. Each torsion angle in a torsion match for a molecule is
2348 processed in an individual process in the process pool.
2349 --mpParams <Name,Value,...> [default: auto]
2350 A comma delimited list of parameter name and value pairs to configure
2351 multiprocessing.
2352
2353 The supported parameter names along with their default and possible
2354 values are shown below:
2355
2356 chunkSize, auto
2357 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
2358 numProcesses, auto [ Default: mp.cpu_count() ]
2359
2360 These parameters are used by the following functions to configure and
2361 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
2362 mp.Pool.imap().
2363
2364 The chunkSize determines chunks of input data passed to each worker
2365 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
2366 The default value of chunkSize is dependent on the value of 'inputDataMode'.
2367
2368 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
2369 automatically converts RDKit data iterable into a list, loads all data into
2370 memory, and calculates the default chunkSize using the following method
2371 as shown in its code:
2372
2373 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
2374 if extra: chunkSize += 1
2375
2376 For example, the default chunkSize will be 7 for a pool of 4 worker processes
2377 and 100 data items.
2378
2379 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
2380 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
2381 data into memory. Consequently, the size of input data is not known a priori.
2382 It's not possible to estimate an optimal value for the chunkSize. The default
2383 chunkSize is set to 1.
2384
2385 The default value for the chunkSize during 'Lazy' data mode may adversely
2386 impact the performance due to the overhead associated with exchanging
2387 small chunks of data. It is generally a good idea to explicitly set chunkSize to
2388 a larger value during 'Lazy' input data mode, based on the size of your input
2389 data and number of processes in the process pool.
2390
2391 The mp.Pool.map() function waits for all worker processes to process all
2392 the data and return the results. The mp.Pool.imap() function, however,
2393 returns the the results obtained from worker processes as soon as the
2394 results become available for specified chunks of data.
2395
2396 The order of data in the results returned by both mp.Pool.map() and
2397 mp.Pool.imap() functions always corresponds to the input data.
2398 -o, --outfile <outfile>
2399 Output file name. The output file root is used for generating the names
2400 of the output files corresponding to structures, energies, and plots during
2401 the torsion scan.
2402 --outfileMolName <yes or no> [default: no]
2403 Append molecule name to output file root during the generation of the names
2404 for output files. The default is to use <MolNum>. The non alphabetical
2405 characters in molecule names are replaced by underscores.
2406 --outfileParams <Name,Value,...> [default: auto]
2407 A comma delimited list of parameter name and value pairs for writing
2408 molecules to files. The supported parameter names for different file
2409 formats, along with their default values, are shown below:
2410
2411 SD: kekulize,yes,forceV3000,no
2412
2413 --outPlotParams <Name,Value,...> [default: auto]
2414 A comma delimited list of parameter name and value pairs for generating
2415 plots using Seaborn module. The supported parameter names along with their
2416 default values are shown below:
2417
2418 type,linepoint,outExt,svg,width,10,height,5.6,
2419 title,auto,xlabel,auto,ylabel,auto,titleWeight,bold,labelWeight,bold
2420 style,darkgrid,palette,deep,font,sans-serif,fontScale,1,
2421 context,notebook
2422
2423 Possible values:
2424
2425 type: linepoint, scatter, or line. Both points and lines are drawn
2426 for linepoint plot type.
2427 outExt: Any valid format supported by Python module Matplotlib.
2428 For example: PDF (.pdf), PNG (.png), PS (.ps), SVG (.svg)
2429 titleWeight, labelWeight: Font weight for title and axes labels.
2430 Any valid value.
2431 style: darkgrid, whitegrid, dark, white, ticks
2432 palette: deep, muted, pastel, dark, bright, colorblind
2433 font: Any valid font name
2434 context: paper, notebook, talk, poster, or any valid name
2435
2436 --outPlotRelativeEnergy <yes or no> [default: yes]
2437 Plot relative energies in the torsion plot. The minimum energy value is
2438 subtracted from energy values to calculate relative energies. This option
2439 is not used during the generation of interactive energy plot for torsion
2440 scan viewer, which always plots relative energy.
2441 --outPlotTitleTorsionSpec <yes or no> [default: yes]
2442 Append torsion specification to the title of the torsion plot.
2443 --outPlotTorsionViewerHeight <number> [default: 430]
2444 Plot height in pixels for interactive relative energy plot generated in
2445 torsion scan viewer. This is different from the width and height specified
2446 using '--outPlotParams' for the standalone plots.
2447 --overwrite
2448 Overwrite existing files.
2449 --precision <number> [default: 6]
2450 Floating point precision for writing energy values.
2451 --psi4OptionsParams <Name,Value,...> [default: none]
2452 A comma delimited list of Psi4 option name and value pairs for setting
2453 global and module options. The names are 'option_name' for global options
2454 and 'module_name__option_name' for options local to a module. The
2455 specified option names must be valid Psi4 names. No validation is
2456 performed.
2457
2458 The specified option name and value pairs are processed and passed to
2459 psi4.set_options() as a dictionary. The supported value types are float,
2460 integer, boolean, or string. The float value string is converted into a float.
2461 The valid values for a boolean string are yes, no, true, false, on, or off.
2462 --psi4RunParams <Name,Value,...> [default: auto]
2463 A comma delimited list of parameter name and value pairs for configuring
2464 Psi4 jobs.
2465
2466 The supported parameter names along with their default and possible
2467 values are shown below:
2468
2469 MemoryInGB, 1
2470 NumThreads, 1
2471 OutputFile, auto [ Possible values: stdout, quiet, or FileName ]
2472 ScratchDir, auto [ Possivle values: DirName]
2473 RemoveOutputFile, yes [ Possible values: yes, no, true, or false]
2474
2475 These parameters control the runtime behavior of Psi4.
2476
2477 The default file name for 'OutputFile' is <InFileRoot>_Psi4.out. The PID
2478 is appended to output file name during multiprocessing as shown:
2479 <InFileRoot>_Psi4_<PIDNum>.out. The 'stdout' value for 'OutputType'
2480 sends Psi4 output to stdout. The 'quiet' or 'devnull' value suppresses
2481 all Psi4 output. The 'OutputFile' is set to 'quiet' for 'auto' value during
2482 'Conformers' of '--mpLevel' option.
2483
2484 The default 'Yes' value of 'RemoveOutputFile' option forces the removal
2485 of any existing Psi4 before creating new files to append output from
2486 multiple Psi4 runs.
2487
2488 The option 'ScratchDir' is a directory path to the location of scratch
2489 files. The default value corresponds to Psi4 default. It may be used to
2490 override the deafult path.
2491 -q, --quiet <yes or no> [default: no]
2492 Use quiet mode. The warning and information messages will not be printed.
2493 --reference <text> [default: auto]
2494 Reference wave function to use for energy calculation or constrained energy
2495 minimization. Default: RHF or UHF. The default values are Restricted Hartree-Fock
2496 (RHF) for closed-shell molecules with all electrons paired and Unrestricted
2497 Hartree-Fock (UHF) for open-shell molecules with unpaired electrons.
2498
2499 The specified value must be a valid Psi4 reference wave function. No validation
2500 is performed. For example: ROHF, CUHF, RKS, etc.
2501
2502 The spin multiplicity determines the default value of reference wave function
2503 for input molecules. It is calculated from number of free radical electrons using
2504 Hund's rule of maximum multiplicity defined as 2S + 1 where S is the total
2505 electron spin. The total spin is 1/2 the number of free radical electrons in a
2506 molecule. The value of 'SpinMultiplicity' molecule property takes precedence
2507 over the calculated value of spin multiplicity.
2508 -t, --torsions <SMILES/SMARTS,...,...>
2509 SMILES/SMARTS patterns corresponding to torsion specifications. It's a
2510 comma delimited list of valid SMILES/SMART patterns.
2511
2512 A substructure match is performed to select torsion atoms in a molecule.
2513 The SMILES pattern match must correspond to four torsion atoms. The
2514 SMARTS patterns containing atom map numbers may match more than four
2515 atoms. The atom map numbers, however, must match exactly four torsion
2516 atoms. For example: [s:1][c:2]([aX2,cH1])!@[CX3:3](O)=[O:4] for thiophene
2517 esters and carboxylates as specified in Torsion Library (TorLib) [Ref 146].
2518 --torsionsFilterbyAtomIndices <Index1, Index2, ...> [default: none]
2519 Comma delimited list of atom indices for filtering torsion matches
2520 corresponding to torsion specifications "-t, --torsions". The atom indices
2521 must be valid. No explicit validation is performed. The list must contain at
2522 least 4 atom indices.
2523
2524 The torsion atom indices, matched by "-t, --torsions" specifications, must be
2525 present in the list. Otherwise, the torsion matches are ignored.
2526 --torsionMaxMatches <number> [default: 5]
2527 Maximum number of torsions to match for each torsion specification in a
2528 molecule.
2529 --torsionMinimize <yes or no> [default: no]
2530 Perform constrained energy minimization on a conformation ensemble
2531 for a specific torsion angle and select the lowest energy conformation
2532 representing the torsion angle. A conformation ensemble is generated for
2533 each 3D structure representing a specific torsion angle using a combination
2534 of distance geometry and forcefield followed by constrained geometry
2535 optimization using a quantum chemistry method.
2536 --torsionRangeMode <Range or Angles> [default: Range]
2537 Perform torsion scan using torsion angles corresponding to a torsion range
2538 or an explicit list of torsion angles. Possible values: Range or Angles. You
2539 may use '--torsionRange' option to specify values for torsion angle or
2540 torsion angles.
2541 --torsionRange <Start,Stop,Step or Angle1,Angle2...> [default: auto]
2542 Start, stop, and step size angles or a comma delimited list of angles in
2543 degrees for a torsion scan.
2544
2545 This value is '--torsionRangeMode' specific. It must be a triplet corresponding
2546 to 'start,Stop,Step' for 'Range' value of '--torsionRange' option. Otherwise, it
2547 is comma delimited list of one or more torsion angles for 'Angles' value of
2548 '--torsionRange' option.
2549
2550 The default values, based on '--torsionRangeMode' option, are shown below:
2551
2552 TorsionRangeMode Default value
2553 Range 0,360,5
2554 Angles None
2555
2556 You must explicitly provide a list of torsion angle(s) for 'Angles' of
2557 '--torsionRangeMode' option.
2558 --useChirality <yes or no> [default: no]
2559 Use chirrality during substructure matches for identification of torsions.
2560 -w, --workingdir <dir>
2561 Location of working directory which defaults to the current directory.
2562
2563 Examples:
2564 To perform a torsion scan on the first molecule in a SMILES file using a minimum
2565 energy structure of the molecule selected from an initial ensemble of conformations
2566 generated using distance geometry and forcefield, skip generation of conformation
2567 ensembles for specific torsion angles and constrained energy minimization of the
2568 ensemble, calculating single point at a specific torsion angle energy using B3LYP/6-31G**
2569 and B3LYP/6-31+G** for non-sulfur and sulfur containing molecules, generate output files
2570 corresponding to structure, energy and torsion plot, type:
2571
2572 % Psi4PerformTorsionScan.py -t "CCCC" -i Psi4SampleTorsionScan.smi
2573 -o SampleOut.sdf
2574
2575 To run the previous example for performing a torsion scan using a specific list
2576 of torsion angles, type:
2577
2578 % Psi4PerformTorsionScan.py -t "CCCC" -i Psi4SampleTorsionScan.smi
2579 -o SampleOut.sdf --torsionRangeMode Angles
2580 --torsionRange "0,180,360"
2581
2582 To run the previous example on the first molecule in a SD file containing 3D
2583 coordinates and skip the generations of initial 3D structure, type:
2584
2585 % Psi4PerformTorsionScan.py -t "CCCC" --infile3D yes
2586 -i Psi4SampleTorsionScan3D.sdf -o SampleOut.sdf
2587
2588 To run the first example on all molecules in a SD file, type:
2589
2590 % Psi4PerformTorsionScan.py -t "CCCC" --modeMols All
2591 -i Psi4SampleTorsionScan.sdf -o SampleOut.sdf
2592
2593 To run the first example on all molecules in a SD file containing 3D
2594 coordinates and skip the generation of initial 3D structures, type:
2595
2596 % Psi4PerformTorsionScan.py -t "CCCC" --infile3D yes
2597 --modeMols All -i Psi4SampleTorsionScan3D.sdf -o SampleOut.sdf
2598
2599 To perform a torsion scan on the first molecule in a SMILES file using a minimum
2600 energy structure of the molecule selected from an initial ensemble of conformations
2601 generated using distance geometry and forcefield, generate up to 50 conformations
2602 for specific torsion angles using ETKDGv2 methodology followed by initial MMFF
2603 forcefield minimization and final energy minimization using B3LYP/6-31G** and
2604 B3LYP/6-31+G** for non-sulfur and sulfur containing molecules, generate output files
2605 corresponding to minimum energy structure, energy and torsion plot, type:
2606
2607 % Psi4PerformTorsionScan.py -t "CCCC" --torsionMinimize Yes
2608 -i Psi4SampleTorsionScan.smi -o SampleOut.sdf
2609
2610 To run the previous example on all molecules in a SD file, type:
2611
2612 % Psi4PerformTorsionScan.py -t "CCCC" --modeMols All
2613 --torsionMinimize Yes -i Psi4SampleTorsionScan.sdf -o SampleOut.sdf
2614
2615 To run the previous example on all molecules in a SD file containing 3D
2616 coordinates and skip the generation of initial 3D structures, type:
2617
2618 % Psi4PerformTorsionScan.py -t "CCCC" --modeMols All
2619 --infile3D yes --modeMols All --torsionMinimize Yes
2620 -i Psi4SampleTorsionScan.sdf -o SampleOut.sdf
2621
2622 To run the previous example in multiprocessing mode at molecules level
2623 on all available CPUs without loading all data into memory and write out
2624 a SD file, type:
2625
2626 % Psi4PerformTorsionScan.py -t "CCCC" -i Psi4SampleTorsionScan.smi
2627 -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
2628
2629 To run the previous example in multiprocessing mode at torsion angles level
2630 on all available CPUs without loading all data into memory and write out
2631 a SD file, type:
2632
2633 % Psi4PerformTorsionScan.py -t "CCCC" -i Psi4SampleTorsionScan.smi
2634 -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
2635 --mpLevel TorsionAngles
2636
2637 To run the previous example in multiprocessing mode on all available CPUs
2638 by loading all data into memory and write out a SD file, type:
2639
2640 % Psi4PerformTorsionScan.py -t "CCCC" -i Psi4SampleTorsionScan.smi
2641 -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
2642 --mpParams "inputDataMode,InMemory"
2643
2644 To run the previous example in multiprocessing mode on specific number of
2645 CPUs and chunk size without loading all data into memory and write out a SD file,
2646 type:
2647
2648 % Psi4PerformTorsionScan.py -t "CCCC" -i Psi4SampleTorsionScan.smi
2649 -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
2650 --mpParams "inputDataMode,Lazy,numProcesses,4,chunkSize,8"
2651
2652 Author:
2653 Manish Sud(msud@san.rr.com)
2654
2655 Acknowledgment:
2656 Pat Walters
2657
2658 See also:
2659 Psi4CalculateEnergy.py, Psi4GenerateConformers.py,
2660 Psi4GenerateConstrainedConformers.py, Psi4PerformConstrainedMinimization.py
2661
2662 Copyright:
2663 Copyright (C) 2026 Manish Sud. All rights reserved.
2664
2665 The functionality available in this script is implemented using RDKit, an
2666 open source toolkit for cheminformatics developed by Greg Landrum.
2667
2668 This file is part of MayaChemTools.
2669
2670 MayaChemTools is free software; you can redistribute it and/or modify it under
2671 the terms of the GNU Lesser General Public License as published by the Free
2672 Software Foundation; either version 3 of the License, or (at your option) any
2673 later version.
2674
2675 """
2676
2677 if __name__ == "__main__":
2678 main()
