1 #!/bin/env python
2 #
3 # File: Psi4VisualizeDualDescriptors.py
4 # Author: Manish Sud <msud@san.rr.com>
5 #
6 # Copyright (C) 2026 Manish Sud. All rights reserved.
7 #
8 # The functionality available in this script is implemented using Psi4, an
9 # open source quantum chemistry software package, and RDKit, an open
10 # source toolkit for cheminformatics developed by Greg Landrum.
11 #
12 # This file is part of MayaChemTools.
13 #
14 # MayaChemTools is free software; you can redistribute it and/or modify it under
15 # the terms of the GNU Lesser General Public License as published by the Free
16 # Software Foundation; either version 3 of the License, or (at your option) any
17 # later version.
18 #
19 # MayaChemTools is distributed in the hope that it will be useful, but without
20 # any warranty; without even the implied warranty of merchantability of fitness
21 # for a particular purpose. See the GNU Lesser General Public License for more
22 # details.
23 #
24 # You should have received a copy of the GNU Lesser General Public License
25 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
26 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
27 # Boston, MA, 02111-1307, USA.
28 #
29
30 from __future__ import print_function
31
32 import os
33 import sys
34 import time
35 import re
36 import glob
37 import shutil
38 import multiprocessing as mp
39
40 # Psi4 imports...
41 if hasattr(shutil, "which") and shutil.which("psi4") is None:
42 sys.stderr.write("\nWarning: Failed to find 'psi4' in your PATH indicating potential issues with your\n")
43 sys.stderr.write("Psi4 environment. The 'import psi4' directive in the global scope of the script\n")
44 sys.stderr.write("interferes with the multiprocessing functionality. It is imported later in the\n")
45 sys.stderr.write("local scope during the execution of the script and may fail. Check/update your\n")
46 sys.stderr.write("Psi4 environment and try again.\n\n")
47
48 # RDKit imports...
49 try:
50 from rdkit import rdBase
51 from rdkit import Chem
52 except ImportError as ErrMsg:
53 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
54 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
55 sys.exit(1)
56
57 # MayaChemTools imports...
58 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
59 try:
60 from docopt import docopt
61 import MiscUtil
62 import Psi4Util
63 import RDKitUtil
64 except ImportError as ErrMsg:
65 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
66 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
67 sys.exit(1)
68
69 ScriptName = os.path.basename(sys.argv[0])
70 Options = {}
71 OptionsInfo = {}
72
73
74 def main():
75 """Start execution of the script."""
76
77 MiscUtil.PrintInfo(
78 "\n%s (Psi4: Imported later; RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
79 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime())
80 )
81
82 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
83
84 # Retrieve command line arguments and options...
85 RetrieveOptions()
86
87 # Process and validate command line arguments and options...
88 ProcessOptions()
89
90 # Perform actions required by the script...
91 GenerateAndVisualizeDualDescriptors()
92
93 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
94 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
95
96
97 def GenerateAndVisualizeDualDescriptors():
98 """Generate and visualize dual descriptors."""
99
100 if OptionsInfo["GenerateCubeFiles"]:
101 GenerateDualDescriptors()
102
103 if OptionsInfo["VisualizeCubeFiles"]:
104 VisualizeDualDescriptors()
105
106
107 def GenerateDualDescriptors():
108 """Generate cube files for dual descriptors."""
109
110 # Setup a molecule reader...
111 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
112 Mols = RDKitUtil.ReadMolecules(OptionsInfo["InfilePath"], **OptionsInfo["InfileParams"])
113
114 MolCount, ValidMolCount, CubeFilesFailedCount = ProcessMolecules(Mols)
115
116 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
117 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
118 MiscUtil.PrintInfo("Number of molecules failed during generation of cube files: %d" % CubeFilesFailedCount)
119 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + CubeFilesFailedCount))
120
121
122 def ProcessMolecules(Mols):
123 """Process and generate cube files for molecules."""
124
125 if OptionsInfo["MPMode"]:
126 return ProcessMoleculesUsingMultipleProcesses(Mols)
127 else:
128 return ProcessMoleculesUsingSingleProcess(Mols)
129
130
131 def ProcessMoleculesUsingSingleProcess(Mols):
132 """Process and generate dual descriptors cube files for molecules using a single process."""
133
134 # Intialize Psi4...
135 MiscUtil.PrintInfo("\nInitializing Psi4...")
136 Psi4Handle = Psi4Util.InitializePsi4(
137 Psi4RunParams=OptionsInfo["Psi4RunParams"],
138 Psi4OptionsParams=OptionsInfo["Psi4OptionsParams"],
139 PrintVersion=True,
140 PrintHeader=True,
141 )
142 OptionsInfo["psi4"] = Psi4Handle
143
144 MiscUtil.PrintInfo("\nGenerating cube files for dual descriptors...")
145
146 (MolCount, ValidMolCount, CubeFilesFailedCount) = [0] * 3
147 for Mol in Mols:
148 MolCount += 1
149
150 if not CheckAndValidateMolecule(Mol, MolCount):
151 continue
152
153 # Setup a Psi4 molecule...
154 Psi4Mol = SetupPsi4Mol(Psi4Handle, Mol, MolCount)
155 if Psi4Mol is None:
156 continue
157
158 ValidMolCount += 1
159
160 CalcStatus = GenerateMolCubeFiles(Psi4Handle, Psi4Mol, Mol, MolCount)
161
162 if not CalcStatus:
163 if not OptionsInfo["QuietMode"]:
164 MiscUtil.PrintWarning(
165 "Failed to generate cube files for molecule %s" % RDKitUtil.GetMolName(Mol, MolCount)
166 )
167
168 CubeFilesFailedCount += 1
169 continue
170
171 return (MolCount, ValidMolCount, CubeFilesFailedCount)
172
173
174 def ProcessMoleculesUsingMultipleProcesses(Mols):
175 """Process and generate dual descriptors cube files for molecules using multiple processes."""
176
177 MiscUtil.PrintInfo("\nGenerating dual descriptors using multiprocessing...")
178
179 MPParams = OptionsInfo["MPParams"]
180
181 # Setup data for initializing a worker process...
182 InitializeWorkerProcessArgs = (
183 MiscUtil.ObjectToBase64EncodedString(Options),
184 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
185 )
186
187 # Setup a encoded mols data iterable for a worker process...
188 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
189
190 # Setup process pool along with data initialization for each process...
191 if not OptionsInfo["QuietMode"]:
192 MiscUtil.PrintInfo(
193 "\nConfiguring multiprocessing using %s method..."
194 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
195 )
196 MiscUtil.PrintInfo(
197 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
198 % (
199 MPParams["NumProcesses"],
200 MPParams["InputDataMode"],
201 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
202 )
203 )
204
205 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
206
207 # Start processing...
208 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
209 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
210 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
211 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
212 else:
213 MiscUtil.PrintError(
214 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
215 )
216
217 # Print out Psi4 version in the main process...
218 MiscUtil.PrintInfo("\nInitializing Psi4...\n")
219 Psi4Handle = Psi4Util.InitializePsi4(PrintVersion=True, PrintHeader=False)
220 OptionsInfo["psi4"] = Psi4Handle
221
222 (MolCount, ValidMolCount, CubeFilesFailedCount) = [0] * 3
223 for Result in Results:
224 MolCount += 1
225 MolIndex, EncodedMol, CalcStatus = Result
226
227 if EncodedMol is None:
228 continue
229
230 ValidMolCount += 1
231
232 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
233
234 if not CalcStatus:
235 if not OptionsInfo["QuietMode"]:
236 MiscUtil.PrintWarning(
237 "Failed to generate cube files for molecule %s" % RDKitUtil.GetMolName(Mol, MolCount)
238 )
239
240 CubeFilesFailedCount += 1
241 continue
242
243 return (MolCount, ValidMolCount, CubeFilesFailedCount)
244
245
246 def InitializeWorkerProcess(*EncodedArgs):
247 """Initialize data for a worker process."""
248
249 global Options, OptionsInfo
250
251 if not OptionsInfo["QuietMode"]:
252 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
253
254 # Decode Options and OptionInfo...
255 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
256 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
257
258 # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
259 # output files for each process...
260 OptionsInfo["Psi4Initialized"] = False
261
262
263 def InitializePsi4ForWorkerProcess():
264 """Initialize Psi4 for a worker process."""
265
266 if OptionsInfo["Psi4Initialized"]:
267 return
268
269 OptionsInfo["Psi4Initialized"] = True
270
271 # Update output file...
272 OptionsInfo["Psi4RunParams"]["OutputFile"] = Psi4Util.UpdatePsi4OutputFileUsingPID(
273 OptionsInfo["Psi4RunParams"]["OutputFile"], os.getpid()
274 )
275
276 # Intialize Psi4...
277 OptionsInfo["psi4"] = Psi4Util.InitializePsi4(
278 Psi4RunParams=OptionsInfo["Psi4RunParams"],
279 Psi4OptionsParams=OptionsInfo["Psi4OptionsParams"],
280 PrintVersion=False,
281 PrintHeader=True,
282 )
283
284
285 def WorkerProcess(EncodedMolInfo):
286 """Process data for a worker process."""
287
288 if not OptionsInfo["Psi4Initialized"]:
289 InitializePsi4ForWorkerProcess()
290
291 MolIndex, EncodedMol = EncodedMolInfo
292
293 CalcStatus = False
294
295 if EncodedMol is None:
296 return [MolIndex, None, CalcStatus]
297
298 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
299 MolCount = MolIndex + 1
300
301 if not CheckAndValidateMolecule(Mol, MolCount):
302 return [MolIndex, None, CalcStatus]
303
304 # Setup a Psi4 molecule...
305 Psi4Mol = SetupPsi4Mol(OptionsInfo["psi4"], Mol, MolCount)
306 if Psi4Mol is None:
307 return [MolIndex, None, CalcStatus]
308
309 CalcStatus = GenerateMolCubeFiles(OptionsInfo["psi4"], Psi4Mol, Mol, MolCount)
310
311 return [
312 MolIndex,
313 RDKitUtil.MolToBase64EncodedMolString(
314 Mol, PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps
315 ),
316 CalcStatus,
317 ]
318
319
320 def GenerateMolCubeFiles(Psi4Handle, Psi4Mol, Mol, MolNum=None):
321 """Generate cube files for dual descriptors."""
322
323 # Setup reference wave function...
324 Reference = SetupReferenceWavefunction(Mol)
325 Psi4Handle.set_options({"Reference": Reference})
326
327 # Setup method name and basis set...
328 MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
329
330 # Calculate single point energy to setup a wavefunction...
331 Status, Energy, WaveFunction = Psi4Util.CalculateSinglePointEnergy(
332 Psi4Handle, Psi4Mol, MethodName, BasisSet, ReturnWaveFunction=True, Quiet=OptionsInfo["QuietMode"]
333 )
334
335 if not Status:
336 PerformPsi4Cleanup(Psi4Handle)
337 return False
338
339 # Generate cube files...
340 Status = GenerateCubeFilesForDualDescriptors(Psi4Handle, WaveFunction, Mol, MolNum)
341
342 # Clean up
343 PerformPsi4Cleanup(Psi4Handle)
344
345 return True
346
347
348 def GenerateCubeFilesForDualDescriptors(Psi4Handle, WaveFunction, Mol, MolNum):
349 """Generate cube files for dual descriptors."""
350
351 # Setup a temporary local directory to generate cube files...
352 Status, CubeFilesDir, CubeFilesDirPath = SetupMolCubeFilesDir(Mol, MolNum)
353 if not Status:
354 return False
355
356 # Generate cube files using psi4.cubeprop...
357 Status, Psi4CubeFiles = GenerateCubeFilesUsingCubeprop(
358 Psi4Handle, WaveFunction, Mol, MolNum, CubeFilesDir, CubeFilesDirPath
359 )
360 if not Status:
361 return False
362
363 # Copy and rename cube files...
364 if not MoveAndRenameCubeFiles(Mol, MolNum, Psi4CubeFiles):
365 return False
366
367 # Delete temporary local directory...
368 if os.path.isdir(CubeFilesDir):
369 shutil.rmtree(CubeFilesDir)
370
371 if not GenerateMolFile(Mol, MolNum):
372 return False
373
374 return True
375
376
377 def SetupMolCubeFilesDir(Mol, MolNum):
378 """Setup a directory for generating cube files for a molecule."""
379
380 MolPrefix = SetupMolPrefix(Mol, MolNum)
381
382 CubeFilesDir = "%s_CubeFiles" % (MolPrefix)
383 CubeFilesDirPath = os.path.join(os.getcwd(), CubeFilesDir)
384 try:
385 if os.path.isdir(CubeFilesDir):
386 shutil.rmtree(CubeFilesDir)
387 os.mkdir(CubeFilesDir)
388 except Exception as ErrMsg:
389 if not OptionsInfo["QuietMode"]:
390 MiscUtil.PrintWarning("Failed to create cube files: %s\n" % ErrMsg)
391 MiscUtil.PrintWarning("Ignoring molecule: %s" % RDKitUtil.GetMolName(Mol, MolNum))
392 return (False, None, None)
393
394 return (True, CubeFilesDir, CubeFilesDirPath)
395
396
397 def GenerateCubeFilesUsingCubeprop(Psi4Handle, WaveFunction, Mol, MolNum, CubeFilesDir, CubeFilesDirPath):
398 """Generate cube files using cubeprop."""
399
400 Psi4CubeFilesParams = OptionsInfo["Psi4CubeFilesParams"]
401
402 # Generate cube files...
403 GridSpacing = Psi4CubeFilesParams["GridSpacing"]
404 GridOverage = Psi4CubeFilesParams["GridOverage"]
405 IsoContourThreshold = Psi4CubeFilesParams["IsoContourThreshold"]
406 try:
407 Psi4Handle.set_options(
408 {
409 "CUBEPROP_TASKS": ["DUAL_DESCRIPTOR"],
410 "CUBIC_GRID_SPACING": [GridSpacing, GridSpacing, GridSpacing],
411 "CUBIC_GRID_OVERAGE": [GridOverage, GridOverage, GridOverage],
412 "CUBEPROP_FILEPATH": CubeFilesDirPath,
413 "CUBEPROP_ISOCONTOUR_THRESHOLD": IsoContourThreshold,
414 }
415 )
416 Psi4Handle.cubeprop(WaveFunction)
417 except Exception as ErrMsg:
418 shutil.rmtree(CubeFilesDir)
419 if not OptionsInfo["QuietMode"]:
420 MiscUtil.PrintWarning("Failed to create cube files: %s\n" % ErrMsg)
421 MiscUtil.PrintWarning("Ignoring molecule: %s" % RDKitUtil.GetMolName(Mol, MolNum))
422 return (False, None)
423
424 # Collect cube files...
425 Psi4CubeFiles = glob.glob(os.path.join(CubeFilesDir, "DUAL*.cube"))
426 if len(Psi4CubeFiles) == 0:
427 shutil.rmtree(CubeFilesDir)
428 if not OptionsInfo["QuietMode"]:
429 MiscUtil.PrintWarning("Failed to create cube files for dual descriptors...")
430 MiscUtil.PrintWarning("Ignoring molecule: %s" % RDKitUtil.GetMolName(Mol, MolNum))
431 return (False, None)
432
433 return (True, Psi4CubeFiles)
434
435
436 def MoveAndRenameCubeFiles(Mol, MolNum, Psi4CubeFiles):
437 """Move and rename cube files."""
438
439 for Psi4CubeFileName in Psi4CubeFiles:
440 try:
441 NewCubeFileName = SetupMolCubeFileName(Mol, MolNum, Psi4CubeFileName)
442 shutil.move(Psi4CubeFileName, NewCubeFileName)
443 except Exception as ErrMsg:
444 if not OptionsInfo["QuietMode"]:
445 MiscUtil.PrintWarning("Failed to move cube file: %s\n" % ErrMsg)
446 MiscUtil.PrintWarning("Ignoring molecule: %s" % RDKitUtil.GetMolName(Mol, MolNum))
447 return False
448
449 return True
450
451
452 def GenerateMolFile(Mol, MolNum):
453 """Generate a SD file for molecules corresponding to cube files."""
454
455 Outfile = SetupMolFileName(Mol, MolNum)
456
457 OutfileParams = {}
458 Writer = RDKitUtil.MoleculesWriter(Outfile, **OutfileParams)
459 if Writer is None:
460 if not OptionsInfo["QuietMode"]:
461 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
462 return False
463
464 Writer.write(Mol)
465
466 if Writer is not None:
467 Writer.close()
468
469 return True
470
471
472 def VisualizeDualDescriptors():
473 """Visualize cube files for dual descriptors."""
474
475 MiscUtil.PrintInfo("\nVisualizing cube files for dual descriptors...")
476
477 if not OptionsInfo["GenerateCubeFiles"]:
478 MiscUtil.PrintInfo("\nInitializing Psi4...\n")
479 Psi4Handle = Psi4Util.InitializePsi4(PrintVersion=True, PrintHeader=False)
480 OptionsInfo["psi4"] = Psi4Handle
481
482 # Setup a molecule reader...
483 MiscUtil.PrintInfo("Processing file %s..." % OptionsInfo["Infile"])
484 Mols = RDKitUtil.ReadMolecules(OptionsInfo["InfilePath"], **OptionsInfo["InfileParams"])
485
486 # Setup for writing a PyMOL PML file...
487 Outfile = OptionsInfo["Outfile"]
488 OutFH = open(Outfile, "w")
489 if OutFH is None:
490 MiscUtil.PrintError("Failed to open output fie %s " % Outfile)
491 MiscUtil.PrintInfo("\nGenerating file %s..." % Outfile)
492
493 # Setup header...
494 WritePMLHeader(OutFH, ScriptName)
495 WritePyMOLParameters(OutFH)
496
497 # Process cube files for molecules...
498 FirstMol = True
499 (MolCount, ValidMolCount, CubeFilesMissingCount) = [0] * 3
500 for Mol in Mols:
501 MolCount += 1
502
503 if Mol is None:
504 if not OptionsInfo["QuietMode"]:
505 MiscUtil.PrintInfo("\nProcessing molecule number %s..." % MolCount)
506 continue
507 ValidMolCount += 1
508
509 MolName = RDKitUtil.GetMolName(Mol, MolCount)
510 if not OptionsInfo["QuietMode"]:
511 MiscUtil.PrintInfo("\nProcessing molecule %s..." % MolName)
512
513 # Retrieve cube files...
514 CubeFilesStatus, CubeFilesInfo = RetrieveMolCubeFilesInfo(Mol, MolCount)
515 if not CubeFilesStatus:
516 CubeFilesMissingCount += 1
517 if not OptionsInfo["QuietMode"]:
518 MiscUtil.PrintWarning("Ignoring molecule with missing cube file...\n")
519 continue
520
521 # Setup PyMOL object names..
522 PyMOLObjectNames = SetupPyMOLObjectNames(Mol, MolCount, CubeFilesInfo)
523
524 # Setup molecule view...
525 WriteMolView(OutFH, CubeFilesInfo, PyMOLObjectNames)
526
527 # Setup dual descriptor views...
528 WriteDualDescriptorView(OutFH, CubeFilesInfo, PyMOLObjectNames)
529
530 # Setup dual descriptor level group...
531 Enable, Action = [True, "open"]
532 GenerateAndWritePMLForGroup(
533 OutFH, PyMOLObjectNames["DualGroup"], PyMOLObjectNames["DualGroupMembers"], Enable, Action
534 )
535
536 # Setup mol level group...
537 Enable, Action = [False, "close"]
538 if FirstMol:
539 FirstMol = False
540 Enable, Action = [True, "open"]
541 GenerateAndWritePMLForGroup(
542 OutFH, PyMOLObjectNames["MolGroup"], PyMOLObjectNames["MolGroupMembers"], Enable, Action
543 )
544
545 OutFH.close()
546
547 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
548 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
549 MiscUtil.PrintInfo("Number of molecules with missing cube files: %d" % CubeFilesMissingCount)
550 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + CubeFilesMissingCount))
551
552
553 def WritePMLHeader(OutFH, ScriptName):
554 """Write out PML header."""
555
556 HeaderInfo = """\
557 #
558 # This file is automatically generated by the following Psi4 script available in
559 # MayaChemTools: %s
560 #
561 cmd.reinitialize() """ % (ScriptName)
562
563 OutFH.write("%s\n" % HeaderInfo)
564
565
566 def WritePyMOLParameters(OutFH):
567 """Write out PyMOL global parameters."""
568
569 OutFH.write("""\n""\n"Setting up PyMOL gobal parameters..."\n""\n""")
570 PMLCmds = []
571 PMLCmds.append("""cmd.set("mesh_quality", %s)""" % (OptionsInfo["PyMOLViewParams"]["MeshQuality"]))
572 PMLCmds.append("""cmd.set("mesh_width", %.2f)""" % (OptionsInfo["PyMOLViewParams"]["MeshWidth"]))
573
574 PMLCmds.append("""cmd.set("surface_quality", %s)""" % (OptionsInfo["PyMOLViewParams"]["SurfaceQuality"]))
575 PMLCmds.append("""cmd.set("transparency", %.2f)""" % (OptionsInfo["PyMOLViewParams"]["SurfaceTransparency"]))
576 PML = "\n".join(PMLCmds)
577 OutFH.write("%s\n" % PML)
578
579 if OptionsInfo["PyMOLViewParams"]["SetupGlobalVolumeColorRamp"]:
580 PyMOLViewParams = OptionsInfo["PyMOLViewParams"]
581 GenerateAndWritePMLForVolumeColorRamp(
582 OutFH,
583 PyMOLViewParams["GlobalVolumeColorRampName"],
584 PyMOLViewParams["ContourLevel1"],
585 PyMOLViewParams["ContourColor1"],
586 PyMOLViewParams["ContourLevel2"],
587 PyMOLViewParams["ContourColor2"],
588 PyMOLViewParams["VolumeContourWindowFactor"],
589 PyMOLViewParams["VolumeColorRampOpacity"],
590 )
591
592
593 def WriteMolView(OutFH, CubeFilesInfo, PyMOLObjectNames):
594 """Write out PML for viewing molecules."""
595
596 MolFile = CubeFilesInfo["MolFile"]
597
598 OutFH.write("""\n""\n"Loading %s and setting up view for molecule..."\n""\n""" % MolFile)
599
600 PyMOLViewParams = OptionsInfo["PyMOLViewParams"]
601 PMLCmds = []
602
603 # Molecule...
604 Name = PyMOLObjectNames["Mol"]
605 PMLCmds.append("""cmd.load("%s", "%s")""" % (MolFile, Name))
606 PMLCmds.append("""cmd.hide("everything", "%s")""" % (Name))
607 PMLCmds.append("""util.cbag("%s", _self = cmd)""" % (Name))
608 PMLCmds.append("""cmd.show("sticks", "%s")""" % (Name))
609 PMLCmds.append("""cmd.set("stick_radius", %s, "%s")""" % (PyMOLViewParams["DisplayStickRadius"], Name))
610 if PyMOLViewParams["HideHydrogens"]:
611 PMLCmds.append("""cmd.hide("(%s and hydro)")""" % (Name))
612 if re.match("^Sticks$", PyMOLViewParams["DisplayMolecule"]):
613 PMLCmds.append("""cmd.enable("%s")""" % (Name))
614 else:
615 PMLCmds.append("""cmd.disable("%s")""" % (Name))
616
617 # Molecule ball and stick...
618 Name = PyMOLObjectNames["MolBallAndStick"]
619 PMLCmds.append("""cmd.load("%s", "%s")""" % (MolFile, Name))
620 PMLCmds.append("""cmd.hide("everything", "%s")""" % (Name))
621 PMLCmds.append("""util.cbag("%s", _self = cmd)""" % (Name))
622 PMLCmds.append("""cmd.show("sphere", "%s")""" % (Name))
623 PMLCmds.append("""cmd.show("sticks", "%s")""" % (Name))
624 PMLCmds.append("""cmd.set("sphere_scale", %s, "%s")""" % (PyMOLViewParams["DisplaySphereScale"], Name))
625 PMLCmds.append("""cmd.set("stick_radius", %s, "%s")""" % (PyMOLViewParams["DisplayStickRadius"], Name))
626 if PyMOLViewParams["HideHydrogens"]:
627 PMLCmds.append("""cmd.hide("(%s and hydro)")""" % (Name))
628 if re.match("^BallAndStick$", PyMOLViewParams["DisplayMolecule"]):
629 PMLCmds.append("""cmd.enable("%s")""" % (Name))
630 else:
631 PMLCmds.append("""cmd.disable("%s")""" % (Name))
632
633 PML = "\n".join(PMLCmds)
634 OutFH.write("%s\n" % PML)
635
636 # MolAlone group...
637 GenerateAndWritePMLForGroup(
638 OutFH, PyMOLObjectNames["MolAloneGroup"], PyMOLObjectNames["MolAloneGroupMembers"], True, "open"
639 )
640
641
642 def WriteDualDescriptorView(OutFH, CubeFilesInfo, PyMOLObjectNames):
643 """Write out PML for viewing dual descriptors in a molecule."""
644
645 OutFH.write("""\n""\n"Setting up views for dual descriptor..."\n""\n""")
646
647 # Cube...
648 DualCubeFile = CubeFilesInfo["FileName"]
649 DualCubeName = PyMOLObjectNames["DualCube"]
650 PMLCmds = []
651 PMLCmds.append("""cmd.load("%s", "%s")""" % (DualCubeFile, DualCubeName))
652 PMLCmds.append("""cmd.disable("%s")""" % (DualCubeName))
653 PMLCmds.append("")
654 PML = "\n".join(PMLCmds)
655 OutFH.write("%s\n" % PML)
656
657 ContourLevel1 = CubeFilesInfo["ContourLevel1"]
658 ContourLevel2 = CubeFilesInfo["ContourLevel2"]
659
660 ContourColor1 = CubeFilesInfo["ContourColor1"]
661 ContourColor2 = CubeFilesInfo["ContourColor2"]
662
663 ContourWindowFactor = CubeFilesInfo["ContourWindowFactor"]
664 ContourColorOpacity = CubeFilesInfo["ContourColorOpacity"]
665
666 # Volume ramp...
667 if CubeFilesInfo["SetupLocalVolumeColorRamp"]:
668 GenerateAndWritePMLForVolumeColorRamp(
669 OutFH,
670 CubeFilesInfo["VolumeColorRampName"],
671 ContourLevel1,
672 ContourColor1,
673 ContourLevel2,
674 ContourColor2,
675 ContourWindowFactor,
676 ContourColorOpacity,
677 )
678
679 # Volume...
680 DualVolumeName = PyMOLObjectNames["DualVolume"]
681 PMLCmds = []
682 PMLCmds.append(
683 """cmd.volume("%s", "%s", "%s")""" % (DualVolumeName, DualCubeName, CubeFilesInfo["VolumeColorRampName"])
684 )
685 PMLCmds.append("""cmd.disable("%s")""" % (DualVolumeName))
686 PMLCmds.append("")
687 PML = "\n".join(PMLCmds)
688 OutFH.write("%s\n" % PML)
689
690 # Mesh...
691 PMLCmds = []
692 DualMeshName = PyMOLObjectNames["DualMeshPositive"]
693 PMLCmds.append("""cmd.isomesh("%s", "%s", "%s")""" % (DualMeshName, DualCubeName, ContourLevel2))
694 PMLCmds.append("""util.color_deep("%s", "%s")""" % (ContourColor2, DualMeshName))
695 PMLCmds.append("""cmd.enable("%s")""" % (DualMeshName))
696 PMLCmds.append("")
697
698 DualMeshName = PyMOLObjectNames["DualMeshNegative"]
699 PMLCmds.append("""cmd.isomesh("%s", "%s", "%s")""" % (DualMeshName, DualCubeName, ContourLevel1))
700 PMLCmds.append("""util.color_deep("%s", "%s")""" % (ContourColor1, DualMeshName))
701 PMLCmds.append("""cmd.enable("%s")""" % (DualMeshName))
702 PMLCmds.append("")
703 PML = "\n".join(PMLCmds)
704 OutFH.write("%s\n" % PML)
705
706 GenerateAndWritePMLForGroup(
707 OutFH, PyMOLObjectNames["DualMeshGroup"], PyMOLObjectNames["DualMeshGroupMembers"], False, "open"
708 )
709
710 # Surface...
711 PMLCmds = []
712 DualSurfaceName = PyMOLObjectNames["DualSurfacePositive"]
713 PMLCmds.append("""cmd.isosurface("%s", "%s", "%s")""" % (DualSurfaceName, DualCubeName, ContourLevel2))
714 PMLCmds.append("""util.color_deep("%s", "%s")""" % (ContourColor2, DualSurfaceName))
715 PMLCmds.append("""cmd.enable("%s")""" % (DualSurfaceName))
716 PMLCmds.append("")
717
718 DualSurfaceName = PyMOLObjectNames["DualSurfaceNegative"]
719 PMLCmds.append("""cmd.isosurface("%s", "%s", "%s")""" % (DualSurfaceName, DualCubeName, ContourLevel1))
720 PMLCmds.append("""util.color_deep("%s", "%s")""" % (ContourColor1, DualSurfaceName))
721 PMLCmds.append("""cmd.enable("%s")""" % (DualSurfaceName))
722 PMLCmds.append("")
723 PML = "\n".join(PMLCmds)
724 OutFH.write("%s\n" % PML)
725
726 GenerateAndWritePMLForGroup(
727 OutFH, PyMOLObjectNames["DualSurfaceGroup"], PyMOLObjectNames["DualSurfaceGroupMembers"], True, "open"
728 )
729
730
731 def GenerateAndWritePMLForVolumeColorRamp(
732 OutFH,
733 ColorRampName,
734 ContourLevel1,
735 ContourColor1,
736 ContourLevel2,
737 ContourColor2,
738 ContourWindowFactor,
739 ContourColorOpacity,
740 ):
741 """Write out PML for volume color ramp."""
742
743 RampWindowOffset1 = abs(ContourLevel1 * ContourWindowFactor)
744 RampWindowOffset2 = abs(ContourLevel2 * ContourWindowFactor)
745
746 PML = """\
747 cmd.volume_ramp_new("%s", [%.4f, "%s", 0.00, %.4f, "%s", %s, %.4f, "%s", 0.00, %4f, "%s", 0.00, %.4f, "%s", %s, %4f, "%s", 0.00])""" % (
748 ColorRampName,
749 (ContourLevel1 - RampWindowOffset1),
750 ContourColor1,
751 ContourLevel1,
752 ContourColor1,
753 ContourColorOpacity,
754 (ContourLevel1 + RampWindowOffset1),
755 ContourColor1,
756 (ContourLevel2 - RampWindowOffset2),
757 ContourColor2,
758 ContourLevel2,
759 ContourColor2,
760 ContourColorOpacity,
761 (ContourLevel2 + RampWindowOffset2),
762 ContourColor2,
763 )
764
765 OutFH.write("%s\n" % PML)
766
767
768 def GenerateAndWritePMLForGroup(OutFH, GroupName, GroupMembersList, Enable, Action):
769 """Generate and write PML for group."""
770
771 OutFH.write("""\n""\n"Setting up group %s..."\n""\n""" % GroupName)
772
773 PMLCmds = []
774
775 GroupMembers = " ".join(GroupMembersList)
776 PMLCmds.append("""cmd.group("%s", "%s")""" % (GroupName, GroupMembers))
777
778 if Enable is not None:
779 if Enable:
780 PMLCmds.append("""cmd.enable("%s")""" % GroupName)
781 else:
782 PMLCmds.append("""cmd.disable("%s")""" % GroupName)
783
784 if Action is not None:
785 PMLCmds.append("""cmd.group("%s", action="%s")""" % (GroupName, Action))
786
787 PML = "\n".join(PMLCmds)
788
789 OutFH.write("%s\n\n" % PML)
790
791
792 def RetrieveMolCubeFilesInfo(Mol, MolNum):
793 """Retrieve available cube files info for a molecule."""
794
795 # Initialize cube files info...
796 CubeFilesInfo = {}
797 CubeFilesInfo["FileName"] = None
798 CubeFilesInfo["IsocontourRangeMin"] = None
799 CubeFilesInfo["IsocontourRangeMax"] = None
800
801 CubeFilesInfo["ContourLevel1"] = None
802 CubeFilesInfo["ContourLevel2"] = None
803
804 CubeFilesInfo["ContourColor1"] = None
805 CubeFilesInfo["ContourColor2"] = None
806
807 CubeFilesInfo["ContourWindowFactor"] = None
808 CubeFilesInfo["ContourColorOpacity"] = None
809
810 CubeFilesInfo["VolumeColorRampName"] = None
811 CubeFilesInfo["SetupLocalVolumeColorRamp"] = False
812
813 # Setup cube mol info...
814 CubeFilesInfo["MolPrefix"] = SetupMolPrefix(Mol, MolNum)
815 CubeFilesInfo["MolName"] = SetupMolName(Mol, MolNum)
816 CubeFilesInfo["MolFile"] = SetupMolFileName(Mol, MolNum)
817
818 # Retrieve cube files...
819 CubeFiles = glob.glob("%s*DUAL*.cube" % (CubeFilesInfo["MolPrefix"]))
820
821 if len(CubeFiles) == 0:
822 return (False, CubeFilesInfo)
823
824 CubeFileName = CubeFiles[0]
825 if len(CubeFiles) > 1:
826 if not OptionsInfo["QuietMode"]:
827 MiscUtil.PrintWarning(
828 "Multiple cube files available for molecule %s; Using first cube file, %s, to visualize dual descriptor..."
829 % (RDKitUtil.GetMolName(Mol, MolNum), CubeFileName)
830 )
831 CubeFilesInfo["FileName"] = CubeFileName
832
833 IsocontourRangeMin, IsocontourRangeMax, ContourLevel1, ContourLevel2 = SetupIsocontourRangeAndContourLevels(
834 CubeFileName
835 )
836 CubeFilesInfo["IsocontourRangeMin"] = IsocontourRangeMin
837 CubeFilesInfo["IsocontourRangeMax"] = IsocontourRangeMax
838
839 CubeFilesInfo["ContourLevel1"] = ContourLevel1
840 CubeFilesInfo["ContourLevel2"] = ContourLevel2
841 CubeFilesInfo["ContourColor1"] = OptionsInfo["PyMOLViewParams"]["ContourColor1"]
842 CubeFilesInfo["ContourColor2"] = OptionsInfo["PyMOLViewParams"]["ContourColor2"]
843
844 CubeFilesInfo["ContourWindowFactor"] = OptionsInfo["PyMOLViewParams"]["VolumeContourWindowFactor"]
845 CubeFilesInfo["ContourColorOpacity"] = OptionsInfo["PyMOLViewParams"]["VolumeColorRampOpacity"]
846
847 VolumeColorRampName = SetupVolumeColorRampName(CubeFilesInfo)
848 CubeFilesInfo["VolumeColorRampName"] = VolumeColorRampName
849 CubeFilesInfo["SetupLocalVolumeColorRamp"] = (
850 False if OptionsInfo["PyMOLViewParams"]["SetupGlobalVolumeColorRamp"] else True
851 )
852
853 return (True, CubeFilesInfo)
854
855
856 def SetupVolumeColorRampName(CubeFilesInfo):
857 """Setup volume color ramp name and status."""
858
859 PyMOLViewParams = OptionsInfo["PyMOLViewParams"]
860
861 # Setup info for generating local volume color ramps for individual cube files...
862 VolumeColorRampName = PyMOLViewParams["VolumeColorRamp"]
863 if PyMOLViewParams["VolumeColorRampAuto"]:
864 VolumeColorRampName = SetupDefaultVolumeRampName()
865 if PyMOLViewParams["ContourLevel1Auto"] or PyMOLViewParams["ContourLevel2Auto"]:
866 VolumeColorRampName = "%s_%s" % (CubeFilesInfo["MolPrefix"], VolumeColorRampName)
867
868 return VolumeColorRampName
869
870
871 def SetupIsocontourRangeAndContourLevels(CubeFileName):
872 """Setup isocontour range."""
873
874 PyMOLViewParams = OptionsInfo["PyMOLViewParams"]
875
876 # Setup isocontour range and contour levels...
877 IsocontourRangeMin, IsocontourRangeMax = Psi4Util.RetrieveIsocontourRangeFromCubeFile(CubeFileName)
878
879 DefaultIsocontourRange = 0.06
880 if IsocontourRangeMin is None:
881 IsocontourRangeMin = -DefaultIsocontourRange
882 if not OptionsInfo["QuietMode"]:
883 MiscUtil.PrintInfo(
884 "Failed to retrieve isocontour range from the cube file. Setting min isocontour value to %s..."
885 % (IsocontourRangeMin)
886 )
887
888 if IsocontourRangeMax is None:
889 IsocontourRangeMax = DefaultIsocontourRange
890 if not OptionsInfo["QuietMode"]:
891 MiscUtil.PrintInfo(
892 "Failed to retrieve isocontour range from the cube file. Setting max isocontour value to %s..."
893 % (IsocontourRangeMax)
894 )
895
896 # Setup contour levels...
897 ContourLevel = max(abs(IsocontourRangeMin), abs(IsocontourRangeMax)) * PyMOLViewParams["ContourLevelAutoAt"]
898 if ContourLevel >= 0.01:
899 ContourLevel = float("%.2f" % ContourLevel)
900 elif ContourLevel >= 0.001:
901 ContourLevel = float("%.3f" % ContourLevel)
902 elif ContourLevel >= 0.0001:
903 ContourLevel = float("%.4f" % ContourLevel)
904
905 ContourLevel1 = -ContourLevel if PyMOLViewParams["ContourLevel1Auto"] else PyMOLViewParams["ContourLevel1"]
906 ContourLevel2 = ContourLevel if PyMOLViewParams["ContourLevel2Auto"] else PyMOLViewParams["ContourLevel2"]
907
908 if not OptionsInfo["QuietMode"]:
909 if IsocontourRangeMin is not None and IsocontourRangeMax is not None:
910 MiscUtil.PrintInfo(
911 "CubeFileName: %s; Isocontour range for %s percent of the density: %.4f to %.4f; ContourLevel1: %s; ContourLevel2: %s"
912 % (
913 CubeFileName,
914 (100 * OptionsInfo["Psi4CubeFilesParams"]["IsoContourThreshold"]),
915 IsocontourRangeMin,
916 IsocontourRangeMax,
917 ContourLevel1,
918 ContourLevel2,
919 )
920 )
921
922 return (IsocontourRangeMin, IsocontourRangeMax, ContourLevel1, ContourLevel2)
923
924
925 def SetupPyMOLObjectNames(Mol, MolNum, CubeFilesInfo):
926 """Setup PyMOL object names."""
927
928 PyMOLObjectNames = {}
929
930 SetupPyMOLObjectNamesForMol(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames)
931 SetupPyMOLObjectNamesForDualDescriptors(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames)
932
933 return PyMOLObjectNames
934
935
936 def SetupPyMOLObjectNamesForMol(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames):
937 """Setup groups and objects for molecule."""
938
939 MolFileRoot = CubeFilesInfo["MolPrefix"]
940
941 MolGroupName = "%s" % MolFileRoot
942 PyMOLObjectNames["MolGroup"] = MolGroupName
943 PyMOLObjectNames["MolGroupMembers"] = []
944
945 # Molecule alone group...
946 MolAloneGroupName = "%s.Molecule" % (MolGroupName)
947 PyMOLObjectNames["MolAloneGroup"] = MolAloneGroupName
948 PyMOLObjectNames["MolGroupMembers"].append(MolAloneGroupName)
949
950 PyMOLObjectNames["MolAloneGroupMembers"] = []
951
952 # Molecule...
953 MolName = "%s.Molecule" % (MolAloneGroupName)
954 PyMOLObjectNames["Mol"] = MolName
955 PyMOLObjectNames["MolAloneGroupMembers"].append(MolName)
956
957 # BallAndStick...
958 MolBallAndStickName = "%s.BallAndStick" % (MolAloneGroupName)
959 PyMOLObjectNames["MolBallAndStick"] = MolBallAndStickName
960 PyMOLObjectNames["MolAloneGroupMembers"].append(MolBallAndStickName)
961
962
963 def SetupPyMOLObjectNamesForDualDescriptors(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames):
964 """Setup groups and objects for dual descriptors."""
965
966 MolGroupName = PyMOLObjectNames["MolGroup"]
967
968 # Setup dual group..
969 DualGroupName = "%s.Dual" % (MolGroupName)
970 PyMOLObjectNames["DualGroup"] = DualGroupName
971 PyMOLObjectNames["MolGroupMembers"].append(DualGroupName)
972 PyMOLObjectNames["DualGroupMembers"] = []
973
974 # Cube...
975 DualCubeName = "%s.Cube" % DualGroupName
976 PyMOLObjectNames["DualCube"] = DualCubeName
977 PyMOLObjectNames["DualGroupMembers"].append(DualCubeName)
978
979 # Volume...
980 DualVolumeName = "%s.Volume" % DualGroupName
981 PyMOLObjectNames["DualVolume"] = DualVolumeName
982 PyMOLObjectNames["DualGroupMembers"].append(DualVolumeName)
983
984 # Mesh...
985 DualMeshGroupName = "%s.Mesh" % DualGroupName
986 DualMeshPositiveName = "%s.Positive_Electrophilic" % DualMeshGroupName
987 DualMeshNegativeName = "%s.Negative_Nucleophilic" % DualMeshGroupName
988
989 PyMOLObjectNames["DualMeshGroup"] = DualMeshGroupName
990 PyMOLObjectNames["DualGroupMembers"].append(DualMeshGroupName)
991
992 PyMOLObjectNames["DualMeshPositive"] = DualMeshPositiveName
993 PyMOLObjectNames["DualMeshNegative"] = DualMeshNegativeName
994 PyMOLObjectNames["DualMeshGroupMembers"] = []
995 PyMOLObjectNames["DualMeshGroupMembers"].append(DualMeshPositiveName)
996 PyMOLObjectNames["DualMeshGroupMembers"].append(DualMeshNegativeName)
997
998 # Surface....
999 DualSurfaceGroupName = "%s.Surface" % DualGroupName
1000 DualSurfacePositiveName = "%s.Positive_Electrophilic" % DualSurfaceGroupName
1001 DualSurfaceNegativeName = "%s.Negative_Nucleophilic" % DualSurfaceGroupName
1002
1003 PyMOLObjectNames["DualSurfaceGroup"] = DualSurfaceGroupName
1004 PyMOLObjectNames["DualGroupMembers"].append(DualSurfaceGroupName)
1005
1006 PyMOLObjectNames["DualSurfacePositive"] = DualSurfacePositiveName
1007 PyMOLObjectNames["DualSurfaceNegative"] = DualSurfaceNegativeName
1008 PyMOLObjectNames["DualSurfaceGroupMembers"] = []
1009 PyMOLObjectNames["DualSurfaceGroupMembers"].append(DualSurfacePositiveName)
1010 PyMOLObjectNames["DualSurfaceGroupMembers"].append(DualSurfaceNegativeName)
1011
1012
1013 def SetupMolFileName(Mol, MolNum):
1014 """Setup SD file name for a molecule."""
1015
1016 MolPrefix = SetupMolPrefix(Mol, MolNum)
1017 MolFileName = "%s.sdf" % MolPrefix
1018
1019 return MolFileName
1020
1021
1022 def SetupMolCubeFileName(Mol, MolNum, Psi4CubeFileName):
1023 """Setup cube file name for a molecule."""
1024
1025 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Psi4CubeFileName)
1026 CubeFileName = "%s.%s" % (FileName, FileExt)
1027
1028 # Add MolPrefix...
1029 MolPrefix = SetupMolPrefix(Mol, MolNum)
1030 CubeFileName = "%s_%s" % (MolPrefix, CubeFileName)
1031
1032 return CubeFileName
1033
1034
1035 def SetupMolPrefix(Mol, MolNum):
1036 """Get molecule prefix for generating files and directories."""
1037
1038 MolNamePrefix = ""
1039 if Mol.HasProp("_Name"):
1040 MolNamePrefix = re.sub("[^a-zA-Z0-9]", "_", Mol.GetProp("_Name"))
1041
1042 MolNumPrefix = "Mol%s" % MolNum
1043
1044 if OptionsInfo["UseMolNumPrefix"] and OptionsInfo["UseMolNamePrefix"]:
1045 MolPrefix = MolNumPrefix
1046 if len(MolNamePrefix):
1047 MolPrefix = "%s_%s" % (MolNumPrefix, MolNamePrefix)
1048 elif OptionsInfo["UseMolNamePrefix"]:
1049 MolPrefix = MolNamePrefix if len(MolNamePrefix) else MolNumPrefix
1050 elif OptionsInfo["UseMolNumPrefix"]:
1051 MolPrefix = MolNumPrefix
1052
1053 return MolPrefix
1054
1055
1056 def SetupMolName(Mol, MolNum):
1057 """Get molecule name."""
1058
1059 # Test for creating PyMOL object name for molecule...
1060 if Mol.HasProp("_Name"):
1061 MolName = re.sub("[^a-zA-Z0-9]", "_", Mol.GetProp("_Name"))
1062 else:
1063 MolName = "Mol%s" % MolNum
1064
1065 return MolName
1066
1067
1068 def SetupDefaultVolumeRampName():
1069 """Setup a default name for a volume ramp."""
1070
1071 return "psi4_cube_dual"
1072
1073
1074 def SetupPsi4Mol(Psi4Handle, Mol, MolCount=None):
1075 """Setup a Psi4 molecule object."""
1076
1077 MolGeometry = RDKitUtil.GetPsi4XYZFormatString(Mol, NoCom=True, NoReorient=True)
1078
1079 try:
1080 Psi4Mol = Psi4Handle.geometry(MolGeometry)
1081 except Exception as ErrMsg:
1082 Psi4Mol = None
1083 if not OptionsInfo["QuietMode"]:
1084 MiscUtil.PrintWarning("Failed to create Psi4 molecule from geometry string: %s\n" % ErrMsg)
1085 MolName = RDKitUtil.GetMolName(Mol, MolCount)
1086 MiscUtil.PrintWarning("Ignoring molecule: %s" % MolName)
1087
1088 return Psi4Mol
1089
1090
1091 def PerformPsi4Cleanup(Psi4Handle):
1092 """Perform clean up."""
1093
1094 # Clean up after Psi4 run ...
1095 Psi4Handle.core.clean()
1096
1097 # Clean up any leftover scratch files...
1098 if OptionsInfo["MPMode"]:
1099 Psi4Util.RemoveScratchFiles(Psi4Handle, OptionsInfo["Psi4RunParams"]["OutputFile"])
1100
1101
1102 def CheckAndValidateMolecule(Mol, MolCount=None):
1103 """Validate molecule for Psi4 calculations."""
1104
1105 if Mol is None:
1106 if not OptionsInfo["QuietMode"]:
1107 MiscUtil.PrintInfo("\nProcessing molecule number %s..." % MolCount)
1108 return False
1109
1110 MolName = RDKitUtil.GetMolName(Mol, MolCount)
1111 if not OptionsInfo["QuietMode"]:
1112 MiscUtil.PrintInfo("\nProcessing molecule %s..." % MolName)
1113
1114 if RDKitUtil.IsMolEmpty(Mol):
1115 if not OptionsInfo["QuietMode"]:
1116 MiscUtil.PrintWarning("Ignoring empty molecule: %s\n" % MolName)
1117 return False
1118
1119 if not RDKitUtil.ValidateElementSymbols(RDKitUtil.GetAtomSymbols(Mol)):
1120 if not OptionsInfo["QuietMode"]:
1121 MiscUtil.PrintWarning("Ignoring molecule containing invalid element symbols: %s\n" % MolName)
1122 return False
1123
1124 # Check for spin multiplicity to warn about open-shell systems...
1125 if RDKitUtil.GetSpinMultiplicity(Mol) > 1:
1126 if not OptionsInfo["QuietMode"]:
1127 MiscUtil.PrintWarning(
1128 "Generation of dual descriptors not supported for open-shell systems. Ignoring molecule containing unpaired electrons: %s\n"
1129 % MolName
1130 )
1131 return False
1132
1133 # Check for 3D flag...
1134 if not Mol.GetConformer().Is3D():
1135 if not OptionsInfo["QuietMode"]:
1136 MiscUtil.PrintWarning("3D tag is not set for molecule: %s\n" % MolName)
1137
1138 # Check for missing hydrogens...
1139 if RDKitUtil.AreHydrogensMissingInMolecule(Mol):
1140 if not OptionsInfo["QuietMode"]:
1141 MiscUtil.PrintWarning("Missing hydrogens in molecule: %s\n" % MolName)
1142
1143 return True
1144
1145
1146 def SetupMethodNameAndBasisSet(Mol):
1147 """Setup method name and basis set."""
1148
1149 MethodName = OptionsInfo["MethodName"]
1150 if OptionsInfo["MethodNameAuto"]:
1151 MethodName = "B3LYP"
1152
1153 BasisSet = OptionsInfo["BasisSet"]
1154 if OptionsInfo["BasisSetAuto"]:
1155 BasisSet = "6-31+G**" if RDKitUtil.IsAtomSymbolPresentInMol(Mol, "S") else "6-31G**"
1156
1157 return (MethodName, BasisSet)
1158
1159
1160 def SetupReferenceWavefunction(Mol):
1161 """Setup reference wavefunction."""
1162
1163 Reference = OptionsInfo["Reference"]
1164 if OptionsInfo["ReferenceAuto"]:
1165 Reference = "UHF" if (RDKitUtil.GetSpinMultiplicity(Mol) > 1) else "RHF"
1166
1167 return Reference
1168
1169
1170 def CheckAndSetupOutfilesDir():
1171 """Check and setup outfiles directory."""
1172
1173 if OptionsInfo["GenerateCubeFiles"]:
1174 if os.path.isdir(OptionsInfo["OutfilesDir"]):
1175 if not Options["--overwrite"]:
1176 MiscUtil.PrintError(
1177 'The outfiles directory, %s, corresponding to output file specified, %s, for option "-o, --outfile" already exist. Use option "--overwrite" or "--ov" and try again to generate and visualize cube files...\n'
1178 % (OptionsInfo["OutfilesDir"], OptionsInfo["Outfile"])
1179 )
1180
1181 if OptionsInfo["VisualizeCubeFiles"]:
1182 if not Options["--overwrite"]:
1183 if os.path.exists(os.path.join(OptionsInfo["OutfilesDir"], OptionsInfo["Outfile"])):
1184 MiscUtil.PrintError(
1185 'The outfile file specified, %s, for option "-o, --outfile" already exist in outfiles dir, %s. Use option "--overwrite" or "--ov" to overwrite and try again to generate and visualize cube files.\n'
1186 % (OptionsInfo["Outfile"], OptionsInfo["OutfilesDir"])
1187 )
1188
1189 if not OptionsInfo["GenerateCubeFiles"]:
1190 if not os.path.isdir(OptionsInfo["OutfilesDir"]):
1191 MiscUtil.PrintError(
1192 'The outfiles directory, %s, corresponding to output file specified, %s, for option "-o, --outfile" doesn\'t exist. Use value, GenerateCubeFiles or Both, for option "--mode" and try again to generate and visualize cube files...\n'
1193 % (OptionsInfo["OutfilesDir"], OptionsInfo["Outfile"])
1194 )
1195
1196 CubeFiles = glob.glob(os.path.join(OptionsInfo["OutfilesDir"], "*.cube"))
1197 if not len(CubeFiles):
1198 MiscUtil.PrintError(
1199 'The outfiles directory, %s, contains no cube files, "*.cube". Use value, GenerateCubeFiles or Both, for option "--mode" and try again to generate and visualize cube...\n'
1200 % (OptionsInfo["OutfilesDir"])
1201 )
1202
1203 OutfilesDir = OptionsInfo["OutfilesDir"]
1204 if OptionsInfo["GenerateCubeFiles"]:
1205 if not os.path.isdir(OutfilesDir):
1206 MiscUtil.PrintInfo("\nCreating directory %s..." % OutfilesDir)
1207 os.mkdir(OutfilesDir)
1208
1209 MiscUtil.PrintInfo("\nChanging directory to %s..." % OutfilesDir)
1210 os.chdir(OutfilesDir)
1211
1212
1213 def ProcessOptions():
1214 """Process and validate command line arguments and options."""
1215
1216 MiscUtil.PrintInfo("Processing options...")
1217
1218 # Validate options...
1219 ValidateOptions()
1220
1221 OptionsInfo["Infile"] = Options["--infile"]
1222 OptionsInfo["InfilePath"] = os.path.abspath(Options["--infile"])
1223
1224 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
1225 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
1226 "--infileParams",
1227 Options["--infileParams"],
1228 InfileName=Options["--infile"],
1229 ParamsDefaultInfo=ParamsDefaultInfoOverride,
1230 )
1231
1232 Outfile = Options["--outfile"]
1233 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Outfile)
1234
1235 OptionsInfo["Outfile"] = Outfile
1236 OptionsInfo["OutfileRoot"] = FileName
1237 OptionsInfo["OutfileExt"] = FileExt
1238
1239 OutfilesDir = Options["--outfilesDir"]
1240 if re.match("^auto$", OutfilesDir, re.I):
1241 OutfilesDir = "%s_DualDescriptors" % OptionsInfo["OutfileRoot"]
1242 OptionsInfo["OutfilesDir"] = OutfilesDir
1243
1244 OptionsInfo["Overwrite"] = Options["--overwrite"]
1245
1246 # Method, basis set, and reference wavefunction...
1247 OptionsInfo["BasisSet"] = Options["--basisSet"]
1248 OptionsInfo["BasisSetAuto"] = True if re.match("^auto$", Options["--basisSet"], re.I) else False
1249
1250 OptionsInfo["MethodName"] = Options["--methodName"]
1251 OptionsInfo["MethodNameAuto"] = True if re.match("^auto$", Options["--methodName"], re.I) else False
1252
1253 OptionsInfo["Reference"] = Options["--reference"]
1254 OptionsInfo["ReferenceAuto"] = True if re.match("^auto$", Options["--reference"], re.I) else False
1255
1256 # Run, options, and cube files parameters...
1257 OptionsInfo["Psi4OptionsParams"] = Psi4Util.ProcessPsi4OptionsParameters(
1258 "--psi4OptionsParams", Options["--psi4OptionsParams"]
1259 )
1260 OptionsInfo["Psi4RunParams"] = Psi4Util.ProcessPsi4RunParameters(
1261 "--psi4RunParams", Options["--psi4RunParams"], InfileName=OptionsInfo["Infile"]
1262 )
1263
1264 ParamsDefaultInfoOverride = {}
1265 OptionsInfo["Psi4CubeFilesParams"] = Psi4Util.ProcessPsi4CubeFilesParameters(
1266 "--psi4CubeFilesParams", Options["--psi4CubeFilesParams"], ParamsDefaultInfo=ParamsDefaultInfoOverride
1267 )
1268
1269 ParamsDefaultInfoOverride = {
1270 "ContourColor1": "red",
1271 "ContourColor2": "blue",
1272 "ContourLevel1": "auto",
1273 "ContourLevel2": "auto",
1274 "ContourLevelAutoAt": 0.5,
1275 "DisplayMolecule": "BallAndStick",
1276 "DisplaySphereScale": 0.2,
1277 "DisplayStickRadius": 0.1,
1278 "HideHydrogens": True,
1279 "MeshWidth": 0.5,
1280 "MeshQuality": 2,
1281 "SurfaceQuality": 2,
1282 "SurfaceTransparency": 0.25,
1283 "VolumeColorRamp": "auto",
1284 "VolumeColorRampOpacity": 0.2,
1285 "VolumeContourWindowFactor": 0.05,
1286 }
1287 OptionsInfo["PyMOLViewParams"] = MiscUtil.ProcessOptionPyMOLCubeFileViewParameters(
1288 "--pymolViewParams", Options["--pymolViewParams"], ParamsDefaultInfo=ParamsDefaultInfoOverride
1289 )
1290
1291 SetupGlobalVolumeColorRamp = False
1292 GlobalVolumeColorRampName = OptionsInfo["PyMOLViewParams"]["VolumeColorRamp"]
1293 if OptionsInfo["PyMOLViewParams"]["VolumeColorRampAuto"]:
1294 GlobalVolumeColorRampName = SetupDefaultVolumeRampName()
1295 if not (
1296 OptionsInfo["PyMOLViewParams"]["ContourLevel1Auto"] or OptionsInfo["PyMOLViewParams"]["ContourLevel2Auto"]
1297 ):
1298 SetupGlobalVolumeColorRamp = True
1299 OptionsInfo["PyMOLViewParams"]["GlobalVolumeColorRampName"] = GlobalVolumeColorRampName
1300 OptionsInfo["PyMOLViewParams"]["SetupGlobalVolumeColorRamp"] = SetupGlobalVolumeColorRamp
1301
1302 Mode = Options["--mode"]
1303 if re.match("^GenerateCubeFiles$", Mode, re.I):
1304 GenerateCubeFiles = True
1305 VisualizeCubeFiles = False
1306 elif re.match("^VisualizeCubeFiles$", Mode, re.I):
1307 GenerateCubeFiles = False
1308 VisualizeCubeFiles = True
1309 else:
1310 GenerateCubeFiles = True
1311 VisualizeCubeFiles = True
1312 OptionsInfo["Mode"] = Mode
1313 OptionsInfo["GenerateCubeFiles"] = GenerateCubeFiles
1314 OptionsInfo["VisualizeCubeFiles"] = VisualizeCubeFiles
1315
1316 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
1317 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
1318
1319 OutfilesMolPrefix = Options["--outfilesMolPrefix"]
1320 if re.match("^MolName$", OutfilesMolPrefix, re.I):
1321 UseMolNamePrefix = True
1322 UseMolNumPrefix = False
1323 elif re.match("^MolNum$", OutfilesMolPrefix, re.I):
1324 UseMolNamePrefix = False
1325 UseMolNumPrefix = True
1326 else:
1327 UseMolNamePrefix = True
1328 UseMolNumPrefix = True
1329 OptionsInfo["OutfilesMolPrefix"] = OutfilesMolPrefix
1330 OptionsInfo["UseMolNamePrefix"] = UseMolNamePrefix
1331 OptionsInfo["UseMolNumPrefix"] = UseMolNumPrefix
1332
1333 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
1334
1335 CheckAndSetupOutfilesDir()
1336
1337
1338 def RetrieveOptions():
1339 """Retrieve command line arguments and options."""
1340
1341 # Get options...
1342 global Options
1343 Options = docopt(_docoptUsage_)
1344
1345 # Set current working directory to the specified directory...
1346 WorkingDir = Options["--workingdir"]
1347 if WorkingDir:
1348 os.chdir(WorkingDir)
1349
1350 # Handle examples option...
1351 if "--examples" in Options and Options["--examples"]:
1352 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
1353 sys.exit(0)
1354
1355
1356 def ValidateOptions():
1357 """Validate option values."""
1358
1359 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
1360 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol")
1361
1362 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "pml")
1363
1364 MiscUtil.ValidateOptionTextValue("--mode", Options["--mode"], "GenerateCubeFiles VisualizeCubeFiles Both")
1365 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
1366
1367 MiscUtil.ValidateOptionTextValue("--outfilesMolPrefix", Options["--outfilesMolPrefix"], "MolNum MolName Both")
1368 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
1369
1370
1371 # Setup a usage string for docopt...
1372 _docoptUsage_ = """
1373 Psi4VisualizeDualDescriptors.py - Visualize dual descriptors for frontier orbitals
1374
1375 Usage:
1376 Psi4VisualizeDualDescriptors.py [--basisSet <text>] [--infileParams <Name,Value,...>] [--methodName <text>]
1377 [--mode <GenerateCubeFiles, VisualizeCubeFiles, Both>] [--mp <yes or no>] [--mpParams <Name, Value,...>]
1378 [--outfilesDir <text>] [--outfilesMolPrefix <MolNum, MolName, Both> ] [--overwrite]
1379 [--psi4CubeFilesParams <Name,Value,...>] [--psi4OptionsParams <Name,Value,...>]
1380 [--psi4RunParams <Name,Value,...>] [--pymolViewParams <Name,Value,...>] [--quiet <yes or no>]
1381 [--reference <text>] [-w <dir>] -i <infile> -o <outfile>
1382 Psi4VisualizeDualDescriptors.py -h | --help | -e | --examples
1383
1384 Description:
1385 Generate and visualize dual descriptors [ Ref 151 ] corresponding to frontier
1386 molecular orbitals for molecules in input file. A set of cube files, corresponding
1387 to dual descriptors, is generated for molecules. The cube files are used to
1388 create a PyMOL visualization file for viewing volumes, meshes, and surfaces
1389 representing dual descriptors for frontier molecular orbitals. An option is
1390 available to skip the generation of new cube files and use existing cube files
1391 to visualize frontier molecular orbitals.
1392
1393 The dual descriptor, a measure of electrophilicity and nucleophilicity for
1394 a molecule, is calculated by Psi4 from frontier molecular orbitals [ Ref 151]:
1395 f2(r) = rhoLUMO(r) - rhoHOMO(r). The sign of the dual descriptor value
1396 corresponds to electrophilic and nucleophilic sites in molecules as shown
1397 below:
1398
1399 Sign Site Reactivity
1400 Positive Electrophilic Favored for a nucleophilic attack
1401 Negative Nucleophilic Favored for a electrophilic attack
1402
1403 A Psi4 XYZ format geometry string is automatically generated for each molecule
1404 in input file. It contains atom symbols and 3D coordinates for each atom in a
1405 molecule. In addition, the formal charge and spin multiplicity are present in the
1406 the geometry string. These values are either retrieved from molecule properties
1407 named 'FormalCharge' and 'SpinMultiplicty' or dynamically calculated for a
1408 molecule.
1409
1410 A set of cube and SD output files is generated for each closed-shell molecule
1411 in input file as shown below:
1412
1413 Ouput dir: <OutfileRoot>_DualDescriptors or <OutfilesDir>
1414
1415 <MolIDPrefix>.sdf
1416 <MolIDPrefix>*DUAL*.cube
1417
1418 In addition, a <OutfileRoot>.pml is generated containing dual descriptors for
1419 frontier molecular orbitals for all molecules in input file.
1420
1421 The supported input file formats are: Mol (.mol), SD (.sdf, .sd)
1422
1423 The supported output file formats are: PyMOL script file (.pml)
1424
1425 A variety of PyMOL groups and objects are created for visualization of dual
1426 descriptors for closed-shell molecules as shown below:
1427
1428 <MoleculeID>
1429 .Molecule
1430 .Molecule
1431 .BallAndStick
1432 .Dual
1433 .Cube
1434 .Volume
1435 .Mesh
1436 .Positive_Electrophilic
1437 .Negative_Nucleophilic
1438 .Surface
1439 .Positive_Electrophilic
1440 .Negative_Nucleophilic
1441 <MoleculeID>
1442 .Molecule
1443 ... ... ...
1444 .Dual
1445 ... ... ...
1446
1447 Options:
1448 -b, --basisSet <text> [default: auto]
1449 Basis set to use for calculating single point energy before generating
1450 cube files corresponding to dual descriptors for frontier molecular orbitals.
1451 Default: 6-31+G** for sulfur containing molecules; Otherwise, 6-31G**
1452 [ Ref 150 ]. The specified value must be a valid Psi4 basis set. No validation
1453 is performed.
1454
1455 The following list shows a representative sample of basis sets available
1456 in Psi4:
1457
1458 STO-3G, 6-31G, 6-31+G, 6-31++G, 6-31G*, 6-31+G*, 6-31++G*,
1459 6-31G**, 6-31+G**, 6-31++G**, 6-311G, 6-311+G, 6-311++G,
1460 6-311G*, 6-311+G*, 6-311++G*, 6-311G**, 6-311+G**, 6-311++G**,
1461 cc-pVDZ, cc-pCVDZ, aug-cc-pVDZ, cc-pVDZ-DK, cc-pCVDZ-DK, def2-SVP,
1462 def2-SVPD, def2-TZVP, def2-TZVPD, def2-TZVPP, def2-TZVPPD
1463
1464 -e, --examples
1465 Print examples.
1466 -h, --help
1467 Print this help message.
1468 -i, --infile <infile>
1469 Input file name.
1470 --infileParams <Name,Value,...> [default: auto]
1471 A comma delimited list of parameter name and value pairs for reading
1472 molecules from files. The supported parameter names for different file
1473 formats, along with their default values, are shown below:
1474
1475 SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
1476
1477 -m, --methodName <text> [default: auto]
1478 Method to use for calculating single point energy before generating
1479 cube files corresponding to dual descriptors for frontier molecular
1480 orbitals. Default: B3LYP [ Ref 150 ]. The specified value must be a
1481 valid Psi4 method name. No validation is performed.
1482
1483 The following list shows a representative sample of methods available
1484 in Psi4:
1485
1486 B1LYP, B2PLYP, B2PLYP-D3BJ, B2PLYP-D3MBJ, B3LYP, B3LYP-D3BJ,
1487 B3LYP-D3MBJ, CAM-B3LYP, CAM-B3LYP-D3BJ, HF, HF-D3BJ, HF3c, M05,
1488 M06, M06-2x, M06-HF, M06-L, MN12-L, MN15, MN15-D3BJ,PBE, PBE0,
1489 PBEH3c, PW6B95, PW6B95-D3BJ, WB97, WB97X, WB97X-D, WB97X-D3BJ
1490
1491 --mode <GenerateCubeFiles, VisualizeCubeFiles, or Both> [default: Both]
1492 Generate and visualize cube files corresponding to dual descriptors for
1493 frontier molecular orbitals. The 'VisualizeCubes' value skips the generation
1494 of new cube files and uses existing cube files for visualization of dual
1495 descriptors. Multiprocessing is not supported during 'VisualizeCubeFiles'
1496 value of '--mode' option.
1497 --mp <yes or no> [default: no]
1498 Use multiprocessing.
1499
1500 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
1501 function employing lazy RDKit data iterable. This allows processing of
1502 arbitrary large data sets without any additional requirements memory.
1503
1504 All input data may be optionally loaded into memory by mp.Pool.map()
1505 before starting worker processes in a process pool by setting the value
1506 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
1507
1508 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
1509 data mode may adversely impact the performance. The '--mpParams' section
1510 provides additional information to tune the value of 'chunkSize'.
1511 --mpParams <Name,Value,...> [default: auto]
1512 A comma delimited list of parameter name and value pairs to configure
1513 multiprocessing.
1514
1515 The supported parameter names along with their default and possible
1516 values are shown below:
1517
1518 chunkSize, auto
1519 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
1520 numProcesses, auto [ Default: mp.cpu_count() ]
1521
1522 These parameters are used by the following functions to configure and
1523 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
1524 mp.Pool.imap().
1525
1526 The chunkSize determines chunks of input data passed to each worker
1527 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
1528 The default value of chunkSize is dependent on the value of 'inputDataMode'.
1529
1530 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
1531 automatically converts RDKit data iterable into a list, loads all data into
1532 memory, and calculates the default chunkSize using the following method
1533 as shown in its code:
1534
1535 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
1536 if extra: chunkSize += 1
1537
1538 For example, the default chunkSize will be 7 for a pool of 4 worker processes
1539 and 100 data items.
1540
1541 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
1542 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
1543 data into memory. Consequently, the size of input data is not known a priori.
1544 It's not possible to estimate an optimal value for the chunkSize. The default
1545 chunkSize is set to 1.
1546
1547 The default value for the chunkSize during 'Lazy' data mode may adversely
1548 impact the performance due to the overhead associated with exchanging
1549 small chunks of data. It is generally a good idea to explicitly set chunkSize to
1550 a larger value during 'Lazy' input data mode, based on the size of your input
1551 data and number of processes in the process pool.
1552
1553 The mp.Pool.map() function waits for all worker processes to process all
1554 the data and return the results. The mp.Pool.imap() function, however,
1555 returns the the results obtained from worker processes as soon as the
1556 results become available for specified chunks of data.
1557
1558 The order of data in the results returned by both mp.Pool.map() and
1559 mp.Pool.imap() functions always corresponds to the input data.
1560 -o, --outfile <outfile>
1561 Output file name for PyMOL PML file. The PML output file, along with cube
1562 files, is generated in a local directory corresponding to '--outfilesDir'
1563 option.
1564 --outfilesDir <text> [default: auto]
1565 Directory name containing PML and cube files. Default:
1566 <OutfileRoot>_DualDescriptors. This directory must be present during
1567 'VisualizeCubeFiles' value of '--mode' option.
1568 --outfilesMolPrefix <MolNum, MolName, Both> [default: Both]
1569 Molecule prefix to use for the names of cube files. Possible values:
1570 MolNum, MolName, or Both. By default, both molecule number and name
1571 are used. The format of molecule prefix is as follows: MolNum - Mol<Num>;
1572 MolName - <MolName>, Both: Mol<Num>_<MolName>. Empty molecule names
1573 are ignored. Molecule numbers are used for empty molecule names.
1574 --overwrite
1575 Overwrite existing files.
1576 --psi4CubeFilesParams <Name,Value,...> [default: auto]
1577 A comma delimited list of parameter name and value pairs for generating
1578 Psi4 cube files.
1579
1580 The supported parameter names along with their default and possible
1581 values are shown below:
1582
1583 gridSpacing, 0.2, gridOverage, 4.0, isoContourThreshold, 0.85
1584
1585 gridSpacing: Grid spacing for generating cube files. Units: Bohr. A higher
1586 value reduces the size of the cube files on the disk. This option corresponds
1587 to Psi4 option CUBIC_GRID_SPACING.
1588
1589 gridOverage: Grid overage for generating cube files. Units: Bohr.This option
1590 corresponds to Psi4 option CUBIC_GRID_OVERAGE.
1591
1592 isoContourThreshold: IsoContour values for generating cube files that capture
1593 specified percent of the probability density using the least amount of grid
1594 points. Default: 0.85 (85%). This option corresponds to Psi4 option
1595 CUBEPROP_ISOCONTOUR_THRESHOLD.
1596 --psi4OptionsParams <Name,Value,...> [default: none]
1597 A comma delimited list of Psi4 option name and value pairs for setting
1598 global and module options. The names are 'option_name' for global options
1599 and 'module_name__option_name' for options local to a module. The
1600 specified option names must be valid Psi4 names. No validation is
1601 performed.
1602
1603 The specified option name and value pairs are processed and passed to
1604 psi4.set_options() as a dictionary. The supported value types are float,
1605 integer, boolean, or string. The float value string is converted into a float.
1606 The valid values for a boolean string are yes, no, true, false, on, or off.
1607 --psi4RunParams <Name,Value,...> [default: auto]
1608 A comma delimited list of parameter name and value pairs for configuring
1609 Psi4 jobs.
1610
1611 The supported parameter names along with their default and possible
1612 values are shown below:
1613
1614 MemoryInGB, 1
1615 NumThreads, 1
1616 OutputFile, auto [ Possible values: stdout, quiet, or FileName ]
1617 ScratchDir, auto [ Possivle values: DirName]
1618 RemoveOutputFile, yes [ Possible values: yes, no, true, or false]
1619
1620 These parameters control the runtime behavior of Psi4.
1621
1622 The default file name for 'OutputFile' is <InFileRoot>_Psi4.out. The PID
1623 is appended to output file name during multiprocessing as shown:
1624 <InFileRoot>_Psi4_<PIDNum>.out. The 'stdout' value for 'OutputType'
1625 sends Psi4 output to stdout. The 'quiet' or 'devnull' value suppresses
1626 all Psi4 output.
1627
1628 The default 'Yes' value of 'RemoveOutputFile' option forces the removal
1629 of any existing Psi4 before creating new files to append output from
1630 multiple Psi4 runs.
1631
1632 The option 'ScratchDir' is a directory path to the location of scratch
1633 files. The default value corresponds to Psi4 default. It may be used to
1634 override the deafult path.
1635 --pymolViewParams <Name,Value,...> [default: auto]
1636 A comma delimited list of parameter name and value pairs for visualizing
1637 cube files in PyMOL.
1638
1639 contourColor1, red, contourColor2, blue,
1640 contourLevel1, auto, contourLevel2, auto,
1641 contourLevelAutoAt, 0.5,
1642 displayMolecule, BallAndStick, displaySphereScale, 0.2,
1643 displayStickRadius, 0.1, hideHydrogens, yes,
1644 meshWidth, 0.5, meshQuality, 2,
1645 surfaceQuality, 2, surfaceTransparency, 0.25,
1646 volumeColorRamp, auto, volumeColorRampOpacity,0.2
1647 volumeContourWindowFactor,0.05
1648
1649 contourColor1 and contourColor2: Color to use for visualizing volumes,
1650 meshes, and surfaces corresponding to the negative and positive values
1651 in cube files. The specified values must be valid PyMOL color names. No
1652 validation is performed.
1653
1654 contourLevel1 and contourLevel2: Contour levels to use for visualizing
1655 volumes, meshes, and surfaces corresponding to the negative and positive
1656 values in cube files. Default: auto. The specified values for contourLevel1
1657 and contourLevel2 must be negative and positive numbers.
1658
1659 The contour levels are automatically calculated by default. The isocontour
1660 range for specified percent of the density is retrieved from the cube files.
1661 The contour levels are set at 'contourLevelAutoAt' of the absolute maximum
1662 value of the isocontour range. For example: contour levels are set to plus and
1663 minus 0.03 at 'contourLevelAutoAt' of 0.5 for isocontour range of -0.06 to
1664 0.06 covering specified percent of the density.
1665
1666 contourLevelAutoAt: Set contour levels at specified fraction of the absolute
1667 maximum value of the isocontour range retrieved from the cube files. This
1668 option is only used during the automatic calculations of the contour levels.
1669
1670 displayMolecule: Display mode for molecules. Possible values: Sticks or
1671 BallAndStick. Both displays objects are created for molecules.
1672
1673 displaySphereScale: Sphere scale for displaying molecule during
1674 BallAndStick display.
1675
1676 displayStickRadius: Stick radius for displaying molecule during Sticks
1677 and BallAndStick display.
1678
1679 hideHydrogens: Hide hydrogens in molecules. Default: yes. Possible
1680 values: yes or no.
1681
1682 meshQuality: Mesh quality for meshes to visualize cube files. The
1683 higher values represents better quality.
1684
1685 meshWidth: Line width for mesh lines to visualize cube files.
1686
1687 surfaceQuality: Surface quality for surfaces to visualize cube files.
1688 The higher values represents better quality.
1689
1690 surfaceTransparency: Surface transparency for surfaces to visualize cube
1691 files.
1692
1693 volumeColorRamp: Name of a PyMOL volume color ramp to use for visualizing
1694 cube files. Default name(s): <OutfielsMolPrefix>_psi4_cube_dual or psi4_cube_dual
1695 The default volume color ramps are automatically generated using contour
1696 levels and colors during 'auto' value of 'volumeColorRamp'. An explicitly
1697 specified value must be a valid PyMOL volume color ramp. No validation is
1698 preformed.
1699
1700 VolumeColorRampOpacity: Opacity for generating volume color ramps
1701 for visualizing cube files. This value is equivalent to 1 minus Transparency.
1702
1703 volumeContourWindowFactor: Fraction of contour level representing window
1704 widths around contour levels during generation of volume color ramps for
1705 visualizing cube files. For example, the value of 0.05 implies a ramp window
1706 size of 0.0015 at contour level of 0.03.
1707 -q, --quiet <yes or no> [default: no]
1708 Use quiet mode. The warning and information messages will not be printed.
1709 -r, --reference <text> [default: auto]
1710 Reference wave function to use for calculating single point energy before
1711 generating cube files for dual descriptors corresponding to frontier molecular
1712 orbitals. Default: RHF or UHF. The default values are Restricted Hartree-Fock (RHF)
1713 for closed-shell molecules with all electrons paired and Unrestricted artree-Fock (UHF)
1714 for open-shell molecules with unpaired electrons.
1715
1716 The specified value must be a valid Psi4 reference wave function. No validation
1717 is performed. For example: ROHF, CUHF, RKS, etc.
1718
1719 The spin multiplicity determines the default value of reference wave function
1720 for input molecules. It is calculated from number of free radical electrons using
1721 Hund's rule of maximum multiplicity defined as 2S + 1 where S is the total
1722 electron spin. The total spin is 1/2 the number of free radical electrons in a
1723 molecule. The value of 'SpinMultiplicity' molecule property takes precedence
1724 over the calculated value of spin multiplicity.
1725 -w, --workingdir <dir>
1726 Location of working directory which defaults to the current directory.
1727
1728 Examples:
1729 To generate and visualize dual descriptors from frontier orbitals based on a
1730 single point energy calculation using B3LYP/6-31G** and B3LYP/6-31+G** for
1731 non-sulfur and sulfur containing closed-shell molecules in a SD file with 3D
1732 structures, and write a new PML file, type:
1733
1734 % Psi4VisualizeDualDescriptors.py -i Psi4Sample3D.sdf
1735 -o Psi4Sample3DOut.pml
1736
1737 To run the first example to only generate cube files and skip generation of
1738 a PML file to visualize dual descriptors for frontier molecular orbitals, type:
1739
1740 % Psi4VisualizeDualDescriptors.py --mode GenerateCubeFiles
1741 -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1742
1743 To run the first example to skip generation of cube files and use existing cube
1744 files to visualize dual descriptors for frontier molecular orbitals and write out
1745 a PML file, type:
1746
1747 % Psi4VisualizeDualDescriptors.py --mode VisualizeCubeFiles
1748 -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1749
1750 To run the first example in multiprocessing mode on all available CPUs
1751 without loading all data into memory and write out a PML file, type:
1752
1753 % Psi4VisualizeDualDescriptors.py --mp yes -i Psi4Sample3D.sdf
1754 -o Psi4Sample3DOut.pml
1755
1756 To run the first example in multiprocessing mode on all available CPUs
1757 by loading all data into memory and write out a PML file, type:
1758
1759 % Psi4VisualizeFrontierOrbitals.py --mp yes --mpParams "inputDataMode,
1760 InMemory" -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1761
1762 To run the first example in multiprocessing mode on all available CPUs
1763 without loading all data into memory along with multiple threads for each
1764 Psi4 run and write out a SD file, type:
1765
1766 % Psi4VisualizeDualDescriptors.py --mp yes --psi4RunParams
1767 "NumThreads,2" -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1768
1769 To run the first example in using a specific set of parameters to generate and
1770 visualize dual descriptors for frontier molecular orbitals and write out a PML file,
1771 type:
1772
1773 % Psi4VisualizeDualDescriptors.py --mode both -m SCF -b aug-cc-pVDZ
1774 --psi4CubeFilesParams "gridSpacing, 0.2, gridOverage, 4.0"
1775 --psi4RunParams "MemoryInGB, 2" --pymolViewParams "contourColor1,
1776 red, contourColor2, blue,contourLevel1, -0.04, contourLevel2, 0.04,
1777 contourLevelAutoAt, 0.75,volumeColorRamp, auto,
1778 volumeColorRampOpacity,0.25, volumeContourWindowFactor,0.05"
1779 -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1780
1781 Author:
1782 Manish Sud(msud@san.rr.com)
1783
1784 See also:
1785 Psi4PerformMinimization.py, Psi4GenerateConformers.py,
1786 Psi4VisualizeElectrostaticPotential.py , Psi4VisualizeFrontierOrbitals.py
1787
1788 Copyright:
1789 Copyright (C) 2026 Manish Sud. All rights reserved.
1790
1791 The functionality available in this script is implemented using Psi4, an
1792 open source quantum chemistry software package, and RDKit, an open
1793 source toolkit for cheminformatics developed by Greg Landrum.
1794
1795 This file is part of MayaChemTools.
1796
1797 MayaChemTools is free software; you can redistribute it and/or modify it under
1798 the terms of the GNU Lesser General Public License as published by the Free
1799 Software Foundation; either version 3 of the License, or (at your option) any
1800 later version.
1801
1802 """
1803
1804 if __name__ == "__main__":
1805 main()
