1 #!/bin/env python
2 #
3 # File: Psi4VisualizeFrontierOrbitals.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 GenerateAndVisualizeFrontierOrbitals()
92
93 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
94 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
95
96
97 def GenerateAndVisualizeFrontierOrbitals():
98 """Generate and visualize frontier molecular orbitals."""
99
100 if OptionsInfo["GenerateCubeFiles"]:
101 GenerateFrontierOrbitals()
102
103 if OptionsInfo["VisualizeCubeFiles"]:
104 VisualizeFrontierOrbitals()
105
106
107 def GenerateFrontierOrbitals():
108 """Generate cube files for frontier orbitals."""
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 frontier orbital 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 frontier orbitals 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 frontier orbitals...")
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 frontier orbitals cube files for molecules using multiple processes."""
176
177 MiscUtil.PrintInfo("\nGenerating frontier orbitals 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 frontier orbitals."""
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 = GenerateCubeFilesForFrontierOrbitals(Psi4Handle, WaveFunction, Mol, MolNum)
341
342 # Clean up
343 PerformPsi4Cleanup(Psi4Handle)
344
345 return True
346
347
348 def GenerateCubeFilesForFrontierOrbitals(Psi4Handle, WaveFunction, Mol, MolNum):
349 """Generate cube files for frontier orbitals."""
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": ["FRONTIER_ORBITALS"],
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, "Psi*.cube"))
426 if len(Psi4CubeFiles) == 0:
427 shutil.rmtree(CubeFilesDir)
428 if not OptionsInfo["QuietMode"]:
429 MiscUtil.PrintWarning("Failed to create cube files for frontier orbitals...")
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 VisualizeFrontierOrbitals():
473 """Visualize cube files for frontier orbitals."""
474
475 MiscUtil.PrintInfo("\nVisualizing cube files for frontier orbitals...")
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 orbital views...
528 FirstOrbital = True
529 for OrbitalID in CubeFilesInfo["OrbitalIDs"]:
530 FirstSpin = True
531 for SpinID in CubeFilesInfo["SpinIDs"][OrbitalID]:
532 WriteOrbitalSpinView(OutFH, CubeFilesInfo, PyMOLObjectNames, OrbitalID, SpinID, FirstSpin)
533 FirstSpin = False
534
535 # Setup orbital level group...
536 Enable, Action = [False, "open"]
537 if FirstOrbital:
538 FirstOrbital = False
539 Enable, Action = [True, "open"]
540 GenerateAndWritePMLForGroup(
541 OutFH,
542 PyMOLObjectNames["Orbitals"][OrbitalID]["OrbitalGroup"],
543 PyMOLObjectNames["Orbitals"][OrbitalID]["OrbitalGroupMembers"],
544 Enable,
545 Action,
546 )
547
548 # Setup mol level group...
549 Enable, Action = [False, "close"]
550 if FirstMol:
551 FirstMol = False
552 Enable, Action = [True, "open"]
553 GenerateAndWritePMLForGroup(
554 OutFH, PyMOLObjectNames["MolGroup"], PyMOLObjectNames["MolGroupMembers"], Enable, Action
555 )
556
557 OutFH.close()
558
559 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
560 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
561 MiscUtil.PrintInfo("Number of molecules with missing cube files: %d" % CubeFilesMissingCount)
562 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + CubeFilesMissingCount))
563
564
565 def WritePMLHeader(OutFH, ScriptName):
566 """Write out PML header."""
567
568 HeaderInfo = """\
569 #
570 # This file is automatically generated by the following Psi4 script available in
571 # MayaChemTools: %s
572 #
573 cmd.reinitialize() """ % (ScriptName)
574
575 OutFH.write("%s\n" % HeaderInfo)
576
577
578 def WritePyMOLParameters(OutFH):
579 """Write out PyMOL global parameters."""
580
581 OutFH.write("""\n""\n"Setting up PyMOL gobal parameters..."\n""\n""")
582 PMLCmds = []
583 PMLCmds.append("""cmd.set("mesh_quality", %s)""" % (OptionsInfo["PyMOLViewParams"]["MeshQuality"]))
584 PMLCmds.append("""cmd.set("mesh_width", %.2f)""" % (OptionsInfo["PyMOLViewParams"]["MeshWidth"]))
585
586 PMLCmds.append("""cmd.set("surface_quality", %s)""" % (OptionsInfo["PyMOLViewParams"]["SurfaceQuality"]))
587 PMLCmds.append("""cmd.set("transparency", %.2f)""" % (OptionsInfo["PyMOLViewParams"]["SurfaceTransparency"]))
588 PML = "\n".join(PMLCmds)
589 OutFH.write("%s\n" % PML)
590
591 if OptionsInfo["PyMOLViewParams"]["SetupGlobalVolumeColorRamp"]:
592 PyMOLViewParams = OptionsInfo["PyMOLViewParams"]
593 GenerateAndWritePMLForVolumeColorRamp(
594 OutFH,
595 PyMOLViewParams["GlobalVolumeColorRampName"],
596 PyMOLViewParams["ContourLevel1"],
597 PyMOLViewParams["ContourColor1"],
598 PyMOLViewParams["ContourLevel2"],
599 PyMOLViewParams["ContourColor2"],
600 PyMOLViewParams["VolumeContourWindowFactor"],
601 PyMOLViewParams["VolumeColorRampOpacity"],
602 )
603
604
605 def WriteMolView(OutFH, CubeFilesInfo, PyMOLObjectNames):
606 """Write out PML for viewing molecules."""
607
608 MolFile = CubeFilesInfo["MolFile"]
609
610 OutFH.write("""\n""\n"Loading %s and setting up view for molecule..."\n""\n""" % MolFile)
611
612 PyMOLViewParams = OptionsInfo["PyMOLViewParams"]
613 PMLCmds = []
614
615 # Molecule...
616 Name = PyMOLObjectNames["Mol"]
617 PMLCmds.append("""cmd.load("%s", "%s")""" % (MolFile, Name))
618 PMLCmds.append("""cmd.hide("everything", "%s")""" % (Name))
619 PMLCmds.append("""util.cbag("%s", _self = cmd)""" % (Name))
620 PMLCmds.append("""cmd.show("sticks", "%s")""" % (Name))
621 PMLCmds.append("""cmd.set("stick_radius", %s, "%s")""" % (PyMOLViewParams["DisplayStickRadius"], Name))
622 if PyMOLViewParams["HideHydrogens"]:
623 PMLCmds.append("""cmd.hide("(%s and hydro)")""" % (Name))
624 if re.match("^Sticks$", PyMOLViewParams["DisplayMolecule"]):
625 PMLCmds.append("""cmd.enable("%s")""" % (Name))
626 else:
627 PMLCmds.append("""cmd.disable("%s")""" % (Name))
628
629 # Molecule ball and stick...
630 Name = PyMOLObjectNames["MolBallAndStick"]
631 PMLCmds.append("""cmd.load("%s", "%s")""" % (MolFile, Name))
632 PMLCmds.append("""cmd.hide("everything", "%s")""" % (Name))
633 PMLCmds.append("""util.cbag("%s", _self = cmd)""" % (Name))
634 PMLCmds.append("""cmd.show("sphere", "%s")""" % (Name))
635 PMLCmds.append("""cmd.show("sticks", "%s")""" % (Name))
636 PMLCmds.append("""cmd.set("sphere_scale", %s, "%s")""" % (PyMOLViewParams["DisplaySphereScale"], Name))
637 PMLCmds.append("""cmd.set("stick_radius", %s, "%s")""" % (PyMOLViewParams["DisplayStickRadius"], Name))
638 if PyMOLViewParams["HideHydrogens"]:
639 PMLCmds.append("""cmd.hide("(%s and hydro)")""" % (Name))
640 if re.match("^BallAndStick$", PyMOLViewParams["DisplayMolecule"]):
641 PMLCmds.append("""cmd.enable("%s")""" % (Name))
642 else:
643 PMLCmds.append("""cmd.disable("%s")""" % (Name))
644
645 PML = "\n".join(PMLCmds)
646 OutFH.write("%s\n" % PML)
647
648 # MolAlone group...
649 GenerateAndWritePMLForGroup(
650 OutFH, PyMOLObjectNames["MolAloneGroup"], PyMOLObjectNames["MolAloneGroupMembers"], True, "open"
651 )
652
653
654 def WriteOrbitalSpinView(OutFH, CubeFilesInfo, PyMOLObjectNames, OrbitalID, SpinID, FirstSpin):
655 """Write out PML for viewing spins in an orbital."""
656
657 OutFH.write("""\n""\n"Setting up views for spin %s in orbital %s..."\n""\n""" % (SpinID, OrbitalID))
658
659 # Cube...
660 SpinCubeFile = CubeFilesInfo["FileName"][OrbitalID][SpinID]
661 SpinCubeName = PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinCube"]
662 PMLCmds = []
663 PMLCmds.append("""cmd.load("%s", "%s")""" % (SpinCubeFile, SpinCubeName))
664 PMLCmds.append("""cmd.disable("%s")""" % (SpinCubeName))
665 PMLCmds.append("")
666 PML = "\n".join(PMLCmds)
667 OutFH.write("%s\n" % PML)
668
669 ContourLevel1 = CubeFilesInfo["ContourLevel1"][OrbitalID][SpinID]
670 ContourLevel2 = CubeFilesInfo["ContourLevel2"][OrbitalID][SpinID]
671
672 ContourColor1 = CubeFilesInfo["ContourColor1"][OrbitalID][SpinID]
673 ContourColor2 = CubeFilesInfo["ContourColor2"][OrbitalID][SpinID]
674
675 ContourWindowFactor = CubeFilesInfo["ContourWindowFactor"][OrbitalID][SpinID]
676 ContourColorOpacity = CubeFilesInfo["ContourColorOpacity"][OrbitalID][SpinID]
677
678 SetupLocalVolumeColorRamp = CubeFilesInfo["SetupLocalVolumeColorRamp"][OrbitalID][SpinID]
679 VolumeColorRampName = CubeFilesInfo["VolumeColorRampName"][OrbitalID][SpinID]
680
681 # Volume ramp...
682 if SetupLocalVolumeColorRamp:
683 GenerateAndWritePMLForVolumeColorRamp(
684 OutFH,
685 VolumeColorRampName,
686 ContourLevel1,
687 ContourColor1,
688 ContourLevel2,
689 ContourColor2,
690 ContourWindowFactor,
691 ContourColorOpacity,
692 )
693
694 # Volume...
695 SpinVolumeName = PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinVolume"]
696 PMLCmds = []
697 PMLCmds.append("""cmd.volume("%s", "%s", "%s")""" % (SpinVolumeName, SpinCubeName, VolumeColorRampName))
698 PMLCmds.append("""cmd.disable("%s")""" % (SpinVolumeName))
699 PMLCmds.append("")
700 PML = "\n".join(PMLCmds)
701 OutFH.write("%s\n" % PML)
702
703 # Mesh...
704 SpinMeshName = PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshNegative"]
705 PMLCmds = []
706 PMLCmds.append("""cmd.isomesh("%s", "%s", "%s")""" % (SpinMeshName, SpinCubeName, ContourLevel1))
707 PMLCmds.append("""util.color_deep("%s", "%s")""" % (ContourColor1, SpinMeshName))
708 PMLCmds.append("""cmd.enable("%s")""" % (SpinMeshName))
709 PMLCmds.append("")
710
711 SpinMeshName = PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshPositive"]
712 PMLCmds.append("""cmd.isomesh("%s", "%s", "%s")""" % (SpinMeshName, SpinCubeName, ContourLevel2))
713 PMLCmds.append("""util.color_deep("%s", "%s")""" % (ContourColor2, SpinMeshName))
714 PMLCmds.append("""cmd.enable("%s")""" % (SpinMeshName))
715 PMLCmds.append("")
716 PML = "\n".join(PMLCmds)
717 OutFH.write("%s\n" % PML)
718
719 GenerateAndWritePMLForGroup(
720 OutFH,
721 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshGroup"],
722 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshGroupMembers"],
723 False,
724 "close",
725 )
726
727 # Surface...
728 SpinSurfaceName = PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfaceNegative"]
729 PMLCmds = []
730 PMLCmds.append("""cmd.isosurface("%s", "%s", "%s")""" % (SpinSurfaceName, SpinCubeName, ContourLevel1))
731 PMLCmds.append("""util.color_deep("%s", "%s")""" % (ContourColor1, SpinSurfaceName))
732 PMLCmds.append("""cmd.enable("%s")""" % (SpinSurfaceName))
733 PMLCmds.append("")
734
735 SpinSurfaceName = PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfacePositive"]
736 PMLCmds.append("""cmd.isosurface("%s", "%s", "%s")""" % (SpinSurfaceName, SpinCubeName, ContourLevel2))
737 PMLCmds.append("""util.color_deep("%s", "%s")""" % (ContourColor2, SpinSurfaceName))
738 PMLCmds.append("""cmd.enable("%s")""" % (SpinSurfaceName))
739 PMLCmds.append("")
740 PML = "\n".join(PMLCmds)
741 OutFH.write("%s\n" % PML)
742
743 GenerateAndWritePMLForGroup(
744 OutFH,
745 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfaceGroup"],
746 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfaceGroupMembers"],
747 True,
748 "close",
749 )
750
751 Enable = True if FirstSpin else False
752 Action = "open" if FirstSpin else "close"
753 GenerateAndWritePMLForGroup(
754 OutFH,
755 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinGroup"],
756 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinGroupMembers"],
757 Enable,
758 Action,
759 )
760
761
762 def GenerateAndWritePMLForVolumeColorRamp(
763 OutFH,
764 ColorRampName,
765 ContourLevel1,
766 ContourColor1,
767 ContourLevel2,
768 ContourColor2,
769 ContourWindowFactor,
770 ContourColorOpacity,
771 ):
772 """Write out PML for volume color ramp."""
773
774 RampWindowOffset1 = abs(ContourLevel1 * ContourWindowFactor)
775 RampWindowOffset2 = abs(ContourLevel2 * ContourWindowFactor)
776
777 PML = """\
778 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])""" % (
779 ColorRampName,
780 (ContourLevel1 - RampWindowOffset1),
781 ContourColor1,
782 ContourLevel1,
783 ContourColor1,
784 ContourColorOpacity,
785 (ContourLevel1 + RampWindowOffset1),
786 ContourColor1,
787 (ContourLevel2 - RampWindowOffset2),
788 ContourColor2,
789 ContourLevel2,
790 ContourColor2,
791 ContourColorOpacity,
792 (ContourLevel2 + RampWindowOffset2),
793 ContourColor2,
794 )
795
796 OutFH.write("%s\n" % PML)
797
798
799 def GenerateAndWritePMLForGroup(OutFH, GroupName, GroupMembersList, Enable, Action):
800 """Generate and write PML for group."""
801
802 OutFH.write("""\n""\n"Setting up group %s..."\n""\n""" % GroupName)
803
804 PMLCmds = []
805
806 GroupMembers = " ".join(GroupMembersList)
807 PMLCmds.append("""cmd.group("%s", "%s")""" % (GroupName, GroupMembers))
808
809 if Enable is not None:
810 if Enable:
811 PMLCmds.append("""cmd.enable("%s")""" % GroupName)
812 else:
813 PMLCmds.append("""cmd.disable("%s")""" % GroupName)
814
815 if Action is not None:
816 PMLCmds.append("""cmd.group("%s", action="%s")""" % (GroupName, Action))
817
818 PML = "\n".join(PMLCmds)
819
820 OutFH.write("%s\n" % PML)
821
822
823 def RetrieveMolCubeFilesInfo(Mol, MolNum):
824 """Retrieve available cube files info for a molecule."""
825
826 CubeFilesInfo = {}
827 CubeFilesInfo["OrbitalIDs"] = []
828 CubeFilesInfo["SpinIDs"] = {}
829 CubeFilesInfo["FileName"] = {}
830
831 CubeFilesInfo["IsocontourRangeMin"] = {}
832 CubeFilesInfo["IsocontourRangeMax"] = {}
833
834 CubeFilesInfo["ContourLevel1"] = {}
835 CubeFilesInfo["ContourLevel2"] = {}
836 CubeFilesInfo["ContourColor1"] = {}
837 CubeFilesInfo["ContourColor2"] = {}
838
839 CubeFilesInfo["ContourWindowFactor"] = {}
840 CubeFilesInfo["ContourColorOpacity"] = {}
841
842 CubeFilesInfo["VolumeColorRampName"] = {}
843 CubeFilesInfo["SetupLocalVolumeColorRamp"] = {}
844
845 # Setup cube mol info...
846 CubeFilesInfo["MolPrefix"] = SetupMolPrefix(Mol, MolNum)
847 CubeFilesInfo["MolName"] = SetupMolName(Mol, MolNum)
848 CubeFilesInfo["MolFile"] = SetupMolFileName(Mol, MolNum)
849
850 MolPrefix = CubeFilesInfo["MolPrefix"]
851 CubeFilesCount = 0
852 for OrbitalID in ["HOMO", "LUMO", "DOMO", "LVMO", "SOMO"]:
853 CubeFiles = glob.glob("%s*%s.cube" % (MolPrefix, OrbitalID))
854 if not len(CubeFiles):
855 continue
856
857 CubeFilesInfo["OrbitalIDs"].append(OrbitalID)
858 CubeFilesInfo["SpinIDs"][OrbitalID] = []
859 CubeFilesInfo["FileName"][OrbitalID] = {}
860
861 CubeFilesInfo["IsocontourRangeMin"][OrbitalID] = {}
862 CubeFilesInfo["IsocontourRangeMax"][OrbitalID] = {}
863
864 CubeFilesInfo["ContourLevel1"][OrbitalID] = {}
865 CubeFilesInfo["ContourLevel2"][OrbitalID] = {}
866 CubeFilesInfo["ContourColor1"][OrbitalID] = {}
867 CubeFilesInfo["ContourColor2"][OrbitalID] = {}
868
869 CubeFilesInfo["ContourWindowFactor"][OrbitalID] = {}
870 CubeFilesInfo["ContourColorOpacity"][OrbitalID] = {}
871
872 CubeFilesInfo["VolumeColorRampName"][OrbitalID] = {}
873 CubeFilesInfo["SetupLocalVolumeColorRamp"][OrbitalID] = {}
874
875 AlphaSpinCount, BetaSpinCount, SpinCount, CurrentSpinCount = [0] * 4
876 for CubeFile in CubeFiles:
877 if re.match("%s_a_" % MolPrefix, CubeFile):
878 SpinID = "Alpha"
879 AlphaSpinCount += 1
880 CurrentSpinCount = AlphaSpinCount
881 elif re.match("%s_b_" % MolPrefix, CubeFile):
882 SpinID = "Beta"
883 BetaSpinCount += 1
884 CurrentSpinCount = BetaSpinCount
885 else:
886 SpinID = "Spin"
887 SpinCount += 1
888 CurrentSpinCount = SpinCount
889
890 # Set up a unique spin ID...
891 if SpinID in CubeFilesInfo["FileName"][OrbitalID]:
892 SpinID = "%s_%s" % (SpinID, CurrentSpinCount)
893
894 CubeFilesCount += 1
895 CubeFilesInfo["SpinIDs"][OrbitalID].append(SpinID)
896 CubeFilesInfo["FileName"][OrbitalID][SpinID] = CubeFile
897
898 IsocontourRangeMin, IsocontourRangeMax, ContourLevel1, ContourLevel2 = SetupIsocontourRangeAndContourLevels(
899 CubeFile
900 )
901
902 CubeFilesInfo["IsocontourRangeMin"][OrbitalID][SpinID] = IsocontourRangeMin
903 CubeFilesInfo["IsocontourRangeMax"][OrbitalID][SpinID] = IsocontourRangeMax
904
905 CubeFilesInfo["ContourLevel1"][OrbitalID][SpinID] = ContourLevel1
906 CubeFilesInfo["ContourLevel2"][OrbitalID][SpinID] = ContourLevel2
907 CubeFilesInfo["ContourColor1"][OrbitalID][SpinID] = OptionsInfo["PyMOLViewParams"]["ContourColor1"]
908 CubeFilesInfo["ContourColor2"][OrbitalID][SpinID] = OptionsInfo["PyMOLViewParams"]["ContourColor2"]
909
910 CubeFilesInfo["ContourWindowFactor"][OrbitalID][SpinID] = OptionsInfo["PyMOLViewParams"][
911 "VolumeContourWindowFactor"
912 ]
913 CubeFilesInfo["ContourColorOpacity"][OrbitalID][SpinID] = OptionsInfo["PyMOLViewParams"][
914 "VolumeColorRampOpacity"
915 ]
916
917 VolumeColorRampName = SetupVolumeColorRampName(CubeFilesInfo)
918 CubeFilesInfo["VolumeColorRampName"][OrbitalID][SpinID] = VolumeColorRampName
919 CubeFilesInfo["SetupLocalVolumeColorRamp"][OrbitalID][SpinID] = (
920 False if OptionsInfo["PyMOLViewParams"]["SetupGlobalVolumeColorRamp"] else True
921 )
922
923 Status = True if CubeFilesCount > 0 else False
924
925 return (Status, CubeFilesInfo)
926
927
928 def SetupVolumeColorRampName(CubeFilesInfo):
929 """Setup volume color ramp name and status."""
930
931 PyMOLViewParams = OptionsInfo["PyMOLViewParams"]
932
933 # Setup info for generating local volume color ramps for individual cube files...
934 VolumeColorRampName = PyMOLViewParams["VolumeColorRamp"]
935 if PyMOLViewParams["VolumeColorRampAuto"]:
936 VolumeColorRampName = SetupDefaultVolumeRampName()
937 if PyMOLViewParams["ContourLevel1Auto"] or PyMOLViewParams["ContourLevel2Auto"]:
938 VolumeColorRampName = "%s_%s" % (CubeFilesInfo["MolPrefix"], VolumeColorRampName)
939
940 return VolumeColorRampName
941
942
943 def SetupIsocontourRangeAndContourLevels(CubeFileName):
944 """Setup isocontour range."""
945
946 PyMOLViewParams = OptionsInfo["PyMOLViewParams"]
947
948 # Setup isocontour range and contour levels...
949 IsocontourRangeMin, IsocontourRangeMax = Psi4Util.RetrieveIsocontourRangeFromCubeFile(CubeFileName)
950
951 DefaultIsocontourRange = 0.06
952 if IsocontourRangeMin is None:
953 IsocontourRangeMin = -DefaultIsocontourRange
954 if not OptionsInfo["QuietMode"]:
955 MiscUtil.PrintInfo(
956 "Failed to retrieve isocontour range from the cube file. Setting min isocontour value to %s..."
957 % (IsocontourRangeMin)
958 )
959
960 if IsocontourRangeMax is None:
961 IsocontourRangeMax = DefaultIsocontourRange
962 if not OptionsInfo["QuietMode"]:
963 MiscUtil.PrintInfo(
964 "Failed to retrieve isocontour range from the cube file. Setting max isocontour value to %s..."
965 % (IsocontourRangeMax)
966 )
967
968 # Setup contour levels...
969 ContourLevel = max(abs(IsocontourRangeMin), abs(IsocontourRangeMax)) * PyMOLViewParams["ContourLevelAutoAt"]
970 if ContourLevel >= 0.01:
971 ContourLevel = float("%.2f" % ContourLevel)
972 elif ContourLevel >= 0.001:
973 ContourLevel = float("%.3f" % ContourLevel)
974 elif ContourLevel >= 0.0001:
975 ContourLevel = float("%.4f" % ContourLevel)
976
977 ContourLevel1 = -ContourLevel if PyMOLViewParams["ContourLevel1Auto"] else PyMOLViewParams["ContourLevel1"]
978 ContourLevel2 = ContourLevel if PyMOLViewParams["ContourLevel2Auto"] else PyMOLViewParams["ContourLevel2"]
979
980 if not OptionsInfo["QuietMode"]:
981 if IsocontourRangeMin is not None and IsocontourRangeMax is not None:
982 MiscUtil.PrintInfo(
983 "CubeFileName: %s; Isocontour range for %s percent of the density: %.4f to %.4f; ContourLevel1: %s; ContourLevel2: %s"
984 % (
985 CubeFileName,
986 (100 * OptionsInfo["Psi4CubeFilesParams"]["IsoContourThreshold"]),
987 IsocontourRangeMin,
988 IsocontourRangeMax,
989 ContourLevel1,
990 ContourLevel2,
991 )
992 )
993
994 return (IsocontourRangeMin, IsocontourRangeMax, ContourLevel1, ContourLevel2)
995
996
997 def SetupPyMOLObjectNames(Mol, MolNum, CubeFilesInfo):
998 """Setup PyMOL object names."""
999
1000 PyMOLObjectNames = {}
1001 PyMOLObjectNames["Orbitals"] = {}
1002 PyMOLObjectNames["Spins"] = {}
1003
1004 SetupPyMOLObjectNamesForMol(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames)
1005
1006 for OrbitalID in CubeFilesInfo["OrbitalIDs"]:
1007 SetupPyMOLObjectNamesForOrbital(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames, OrbitalID)
1008 for SpinID in CubeFilesInfo["SpinIDs"][OrbitalID]:
1009 SetupPyMOLObjectNamesForSpin(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames, OrbitalID, SpinID)
1010
1011 return PyMOLObjectNames
1012
1013
1014 def SetupPyMOLObjectNamesForMol(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames):
1015 """Setup groups and objects for molecule."""
1016
1017 MolFileRoot = CubeFilesInfo["MolPrefix"]
1018
1019 MolGroupName = "%s" % MolFileRoot
1020 PyMOLObjectNames["MolGroup"] = MolGroupName
1021 PyMOLObjectNames["MolGroupMembers"] = []
1022
1023 # Molecule alone group...
1024 MolAloneGroupName = "%s.Molecule" % (MolGroupName)
1025 PyMOLObjectNames["MolAloneGroup"] = MolAloneGroupName
1026 PyMOLObjectNames["MolGroupMembers"].append(MolAloneGroupName)
1027
1028 PyMOLObjectNames["MolAloneGroupMembers"] = []
1029
1030 # Molecule...
1031 MolName = "%s.Molecule" % (MolAloneGroupName)
1032 PyMOLObjectNames["Mol"] = MolName
1033 PyMOLObjectNames["MolAloneGroupMembers"].append(MolName)
1034
1035 # BallAndStick...
1036 MolBallAndStickName = "%s.BallAndStick" % (MolAloneGroupName)
1037 PyMOLObjectNames["MolBallAndStick"] = MolBallAndStickName
1038 PyMOLObjectNames["MolAloneGroupMembers"].append(MolBallAndStickName)
1039
1040
1041 def SetupPyMOLObjectNamesForOrbital(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames, OrbitalID):
1042 """Setup groups and objects for orbital."""
1043
1044 PyMOLObjectNames["Orbitals"][OrbitalID] = {}
1045 PyMOLObjectNames["Spins"][OrbitalID] = {}
1046
1047 MolGroupName = PyMOLObjectNames["MolGroup"]
1048
1049 # Setup orbital group...
1050 OrbitalGroupName = "%s.%s" % (MolGroupName, OrbitalID)
1051 PyMOLObjectNames["Orbitals"][OrbitalID]["OrbitalGroup"] = OrbitalGroupName
1052 PyMOLObjectNames["MolGroupMembers"].append(OrbitalGroupName)
1053 PyMOLObjectNames["Orbitals"][OrbitalID]["OrbitalGroupMembers"] = []
1054
1055
1056 def SetupPyMOLObjectNamesForSpin(Mol, MolNum, CubeFilesInfo, PyMOLObjectNames, OrbitalID, SpinID):
1057 """Setup groups and objects for spin."""
1058
1059 PyMOLObjectNames["Spins"][OrbitalID][SpinID] = {}
1060
1061 OrbitalGroupName = PyMOLObjectNames["Orbitals"][OrbitalID]["OrbitalGroup"]
1062
1063 # Setup spin group at orbital level...
1064 SpinGroupName = "%s.%s" % (OrbitalGroupName, SpinID)
1065 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinGroup"] = SpinGroupName
1066 PyMOLObjectNames["Orbitals"][OrbitalID]["OrbitalGroupMembers"].append(SpinGroupName)
1067 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinGroupMembers"] = []
1068
1069 # Cube...
1070 SpinCubeName = "%s.Cube" % SpinGroupName
1071 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinCube"] = SpinCubeName
1072 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinGroupMembers"].append(SpinCubeName)
1073
1074 # Volume...
1075 SpinVolumeName = "%s.Volume" % SpinGroupName
1076 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinVolume"] = SpinVolumeName
1077 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinGroupMembers"].append(SpinVolumeName)
1078
1079 # Mesh...
1080 SpinMeshGroupName = "%s.Mesh" % SpinGroupName
1081 SpinMeshNegativeName = "%s.Negative_Phase" % SpinMeshGroupName
1082 SpinMeshPositiveName = "%s.Positive_Phase" % SpinMeshGroupName
1083
1084 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshGroup"] = SpinMeshGroupName
1085 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinGroupMembers"].append(SpinMeshGroupName)
1086
1087 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshNegative"] = SpinMeshNegativeName
1088 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshPositive"] = SpinMeshPositiveName
1089 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshGroupMembers"] = []
1090 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshGroupMembers"].append(SpinMeshNegativeName)
1091 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinMeshGroupMembers"].append(SpinMeshPositiveName)
1092
1093 # Surface...
1094 SpinSurfaceGroupName = "%s.Surface" % SpinGroupName
1095 SpinSurfaceNegativeName = "%s.Negative_Phase" % SpinSurfaceGroupName
1096 SpinSurfacePositiveName = "%s.Positive_Phase" % SpinSurfaceGroupName
1097
1098 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfaceGroup"] = SpinSurfaceGroupName
1099 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinGroupMembers"].append(SpinSurfaceGroupName)
1100
1101 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfaceNegative"] = SpinSurfaceNegativeName
1102 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfacePositive"] = SpinSurfacePositiveName
1103 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfaceGroupMembers"] = []
1104 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfaceGroupMembers"].append(SpinSurfaceNegativeName)
1105 PyMOLObjectNames["Spins"][OrbitalID][SpinID]["SpinSurfaceGroupMembers"].append(SpinSurfacePositiveName)
1106
1107
1108 def SetupMolFileName(Mol, MolNum):
1109 """Setup SD file name for a molecule."""
1110
1111 MolPrefix = SetupMolPrefix(Mol, MolNum)
1112 MolFileName = "%s.sdf" % MolPrefix
1113
1114 return MolFileName
1115
1116
1117 def SetupMolCubeFileName(Mol, MolNum, Psi4CubeFileName):
1118 """Setup cube file name for a molecule."""
1119
1120 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Psi4CubeFileName)
1121 CubeFileName = "%s.%s" % (FileName, FileExt)
1122
1123 # Replace Psi by MolPrefix...
1124 MolPrefix = SetupMolPrefix(Mol, MolNum)
1125 CubeFileName = re.sub("^Psi", MolPrefix, CubeFileName)
1126
1127 # Replace prime (') and dprime ('') in the cube file names..
1128 CubeFileName = re.sub('"', "dblprime", CubeFileName)
1129 CubeFileName = re.sub("'", "prime", CubeFileName)
1130
1131 return CubeFileName
1132
1133
1134 def SetupMolPrefix(Mol, MolNum):
1135 """Get molecule prefix for generating files and directories."""
1136
1137 MolNamePrefix = ""
1138 if Mol.HasProp("_Name"):
1139 MolNamePrefix = re.sub("[^a-zA-Z0-9]", "_", Mol.GetProp("_Name"))
1140
1141 MolNumPrefix = "Mol%s" % MolNum
1142
1143 if OptionsInfo["UseMolNumPrefix"] and OptionsInfo["UseMolNamePrefix"]:
1144 MolPrefix = MolNumPrefix
1145 if len(MolNamePrefix):
1146 MolPrefix = "%s_%s" % (MolNumPrefix, MolNamePrefix)
1147 elif OptionsInfo["UseMolNamePrefix"]:
1148 MolPrefix = MolNamePrefix if len(MolNamePrefix) else MolNumPrefix
1149 elif OptionsInfo["UseMolNumPrefix"]:
1150 MolPrefix = MolNumPrefix
1151
1152 return MolPrefix
1153
1154
1155 def SetupMolName(Mol, MolNum):
1156 """Get molecule name."""
1157
1158 # Test for creating PyMOL object name for molecule...
1159 if Mol.HasProp("_Name"):
1160 MolName = re.sub("[^a-zA-Z0-9]", "_", Mol.GetProp("_Name"))
1161 else:
1162 MolName = "Mol%s" % MolNum
1163
1164 return MolName
1165
1166
1167 def SetupDefaultVolumeRampName():
1168 """Setup a default name for a volume ramp."""
1169
1170 return "psi4_cube_mo"
1171
1172
1173 def SetupPsi4Mol(Psi4Handle, Mol, MolCount=None):
1174 """Setup a Psi4 molecule object."""
1175
1176 MolGeometry = RDKitUtil.GetPsi4XYZFormatString(Mol, NoCom=True, NoReorient=True)
1177
1178 try:
1179 Psi4Mol = Psi4Handle.geometry(MolGeometry)
1180 except Exception as ErrMsg:
1181 Psi4Mol = None
1182 if not OptionsInfo["QuietMode"]:
1183 MiscUtil.PrintWarning("Failed to create Psi4 molecule from geometry string: %s\n" % ErrMsg)
1184 MolName = RDKitUtil.GetMolName(Mol, MolCount)
1185 MiscUtil.PrintWarning("Ignoring molecule: %s" % MolName)
1186
1187 return Psi4Mol
1188
1189
1190 def PerformPsi4Cleanup(Psi4Handle):
1191 """Perform clean up."""
1192
1193 # Clean up after Psi4 run...
1194 Psi4Handle.core.clean()
1195
1196 # Clean up any leftover scratch files...
1197 if OptionsInfo["MPMode"]:
1198 Psi4Util.RemoveScratchFiles(Psi4Handle, OptionsInfo["Psi4RunParams"]["OutputFile"])
1199
1200
1201 def CheckAndValidateMolecule(Mol, MolCount=None):
1202 """Validate molecule for Psi4 calculations."""
1203
1204 if Mol is None:
1205 if not OptionsInfo["QuietMode"]:
1206 MiscUtil.PrintInfo("\nProcessing molecule number %s..." % MolCount)
1207 return False
1208
1209 MolName = RDKitUtil.GetMolName(Mol, MolCount)
1210 if not OptionsInfo["QuietMode"]:
1211 MiscUtil.PrintInfo("\nProcessing molecule %s..." % MolName)
1212
1213 if RDKitUtil.IsMolEmpty(Mol):
1214 if not OptionsInfo["QuietMode"]:
1215 MiscUtil.PrintWarning("Ignoring empty molecule: %s\n" % MolName)
1216 return False
1217
1218 if not RDKitUtil.ValidateElementSymbols(RDKitUtil.GetAtomSymbols(Mol)):
1219 if not OptionsInfo["QuietMode"]:
1220 MiscUtil.PrintWarning("Ignoring molecule containing invalid element symbols: %s\n" % MolName)
1221 return False
1222
1223 # Check for 3D flag...
1224 if not Mol.GetConformer().Is3D():
1225 if not OptionsInfo["QuietMode"]:
1226 MiscUtil.PrintWarning("3D tag is not set for molecule: %s\n" % MolName)
1227
1228 # Check for missing hydrogens...
1229 if RDKitUtil.AreHydrogensMissingInMolecule(Mol):
1230 if not OptionsInfo["QuietMode"]:
1231 MiscUtil.PrintWarning("Missing hydrogens in molecule: %s\n" % MolName)
1232
1233 return True
1234
1235
1236 def SetupMethodNameAndBasisSet(Mol):
1237 """Setup method name and basis set."""
1238
1239 MethodName = OptionsInfo["MethodName"]
1240 if OptionsInfo["MethodNameAuto"]:
1241 MethodName = "B3LYP"
1242
1243 BasisSet = OptionsInfo["BasisSet"]
1244 if OptionsInfo["BasisSetAuto"]:
1245 BasisSet = "6-31+G**" if RDKitUtil.IsAtomSymbolPresentInMol(Mol, "S") else "6-31G**"
1246
1247 return (MethodName, BasisSet)
1248
1249
1250 def SetupReferenceWavefunction(Mol):
1251 """Setup reference wavefunction."""
1252
1253 Reference = OptionsInfo["Reference"]
1254 if OptionsInfo["ReferenceAuto"]:
1255 Reference = "UHF" if (RDKitUtil.GetSpinMultiplicity(Mol) > 1) else "RHF"
1256
1257 return Reference
1258
1259
1260 def CheckAndSetupOutfilesDir():
1261 """Check and setup outfiles directory."""
1262
1263 if OptionsInfo["GenerateCubeFiles"]:
1264 if os.path.isdir(OptionsInfo["OutfilesDir"]):
1265 if not Options["--overwrite"]:
1266 MiscUtil.PrintError(
1267 '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'
1268 % (OptionsInfo["OutfilesDir"], OptionsInfo["Outfile"])
1269 )
1270
1271 if OptionsInfo["VisualizeCubeFiles"]:
1272 if not Options["--overwrite"]:
1273 if os.path.exists(os.path.join(OptionsInfo["OutfilesDir"], OptionsInfo["Outfile"])):
1274 MiscUtil.PrintError(
1275 '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'
1276 % (OptionsInfo["Outfile"], OptionsInfo["OutfilesDir"])
1277 )
1278
1279 if not OptionsInfo["GenerateCubeFiles"]:
1280 if not os.path.isdir(OptionsInfo["OutfilesDir"]):
1281 MiscUtil.PrintError(
1282 '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'
1283 % (OptionsInfo["OutfilesDir"], OptionsInfo["Outfile"])
1284 )
1285
1286 CubeFiles = glob.glob(os.path.join(OptionsInfo["OutfilesDir"], "*.cube"))
1287 if not len(CubeFiles):
1288 MiscUtil.PrintError(
1289 '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'
1290 % (OptionsInfo["OutfilesDir"])
1291 )
1292
1293 OutfilesDir = OptionsInfo["OutfilesDir"]
1294 if OptionsInfo["GenerateCubeFiles"]:
1295 if not os.path.isdir(OutfilesDir):
1296 MiscUtil.PrintInfo("\nCreating directory %s..." % OutfilesDir)
1297 os.mkdir(OutfilesDir)
1298
1299 MiscUtil.PrintInfo("\nChanging directory to %s..." % OutfilesDir)
1300 os.chdir(OutfilesDir)
1301
1302
1303 def ProcessOptions():
1304 """Process and validate command line arguments and options."""
1305
1306 MiscUtil.PrintInfo("Processing options...")
1307
1308 # Validate options...
1309 ValidateOptions()
1310
1311 OptionsInfo["Infile"] = Options["--infile"]
1312 OptionsInfo["InfilePath"] = os.path.abspath(Options["--infile"])
1313
1314 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
1315 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
1316 "--infileParams",
1317 Options["--infileParams"],
1318 InfileName=Options["--infile"],
1319 ParamsDefaultInfo=ParamsDefaultInfoOverride,
1320 )
1321
1322 Outfile = Options["--outfile"]
1323 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Outfile)
1324
1325 OptionsInfo["Outfile"] = Outfile
1326 OptionsInfo["OutfileRoot"] = FileName
1327 OptionsInfo["OutfileExt"] = FileExt
1328
1329 OutfilesDir = Options["--outfilesDir"]
1330 if re.match("^auto$", OutfilesDir, re.I):
1331 OutfilesDir = "%s_FrontierOrbitals" % OptionsInfo["OutfileRoot"]
1332 OptionsInfo["OutfilesDir"] = OutfilesDir
1333
1334 OptionsInfo["Overwrite"] = Options["--overwrite"]
1335
1336 # Method, basis set, and reference wavefunction...
1337 OptionsInfo["BasisSet"] = Options["--basisSet"]
1338 OptionsInfo["BasisSetAuto"] = True if re.match("^auto$", Options["--basisSet"], re.I) else False
1339
1340 OptionsInfo["MethodName"] = Options["--methodName"]
1341 OptionsInfo["MethodNameAuto"] = True if re.match("^auto$", Options["--methodName"], re.I) else False
1342
1343 OptionsInfo["Reference"] = Options["--reference"]
1344 OptionsInfo["ReferenceAuto"] = True if re.match("^auto$", Options["--reference"], re.I) else False
1345
1346 # Run, options, and cube files parameters...
1347 OptionsInfo["Psi4OptionsParams"] = Psi4Util.ProcessPsi4OptionsParameters(
1348 "--psi4OptionsParams", Options["--psi4OptionsParams"]
1349 )
1350 OptionsInfo["Psi4RunParams"] = Psi4Util.ProcessPsi4RunParameters(
1351 "--psi4RunParams", Options["--psi4RunParams"], InfileName=OptionsInfo["Infile"]
1352 )
1353
1354 ParamsDefaultInfoOverride = {}
1355 OptionsInfo["Psi4CubeFilesParams"] = Psi4Util.ProcessPsi4CubeFilesParameters(
1356 "--psi4CubeFilesParams", Options["--psi4CubeFilesParams"], ParamsDefaultInfo=ParamsDefaultInfoOverride
1357 )
1358
1359 ParamsDefaultInfoOverride = {
1360 "ContourColor1": "red",
1361 "ContourColor2": "blue",
1362 "ContourLevel1": "auto",
1363 "ContourLevel2": "auto",
1364 "ContourLevelAutoAt": 0.75,
1365 "DisplayMolecule": "BallAndStick",
1366 "DisplaySphereScale": 0.2,
1367 "DisplayStickRadius": 0.1,
1368 "HideHydrogens": True,
1369 "MeshWidth": 0.5,
1370 "MeshQuality": 2,
1371 "SurfaceQuality": 2,
1372 "SurfaceTransparency": 0.25,
1373 "VolumeColorRamp": "auto",
1374 "VolumeColorRampOpacity": 0.2,
1375 "VolumeContourWindowFactor": 0.05,
1376 }
1377 OptionsInfo["PyMOLViewParams"] = MiscUtil.ProcessOptionPyMOLCubeFileViewParameters(
1378 "--pymolViewParams", Options["--pymolViewParams"], ParamsDefaultInfo=ParamsDefaultInfoOverride
1379 )
1380
1381 SetupGlobalVolumeColorRamp = False
1382 GlobalVolumeColorRampName = OptionsInfo["PyMOLViewParams"]["VolumeColorRamp"]
1383 if OptionsInfo["PyMOLViewParams"]["VolumeColorRampAuto"]:
1384 GlobalVolumeColorRampName = SetupDefaultVolumeRampName()
1385 if not (
1386 OptionsInfo["PyMOLViewParams"]["ContourLevel1Auto"] or OptionsInfo["PyMOLViewParams"]["ContourLevel2Auto"]
1387 ):
1388 SetupGlobalVolumeColorRamp = True
1389 OptionsInfo["PyMOLViewParams"]["GlobalVolumeColorRampName"] = GlobalVolumeColorRampName
1390 OptionsInfo["PyMOLViewParams"]["SetupGlobalVolumeColorRamp"] = SetupGlobalVolumeColorRamp
1391
1392 Mode = Options["--mode"]
1393 if re.match("^GenerateCubeFiles$", Mode, re.I):
1394 GenerateCubeFiles = True
1395 VisualizeCubeFiles = False
1396 elif re.match("^VisualizeCubeFiles$", Mode, re.I):
1397 GenerateCubeFiles = False
1398 VisualizeCubeFiles = True
1399 else:
1400 GenerateCubeFiles = True
1401 VisualizeCubeFiles = True
1402 OptionsInfo["Mode"] = Mode
1403 OptionsInfo["GenerateCubeFiles"] = GenerateCubeFiles
1404 OptionsInfo["VisualizeCubeFiles"] = VisualizeCubeFiles
1405
1406 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
1407 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
1408
1409 OutfilesMolPrefix = Options["--outfilesMolPrefix"]
1410 if re.match("^MolName$", OutfilesMolPrefix, re.I):
1411 UseMolNamePrefix = True
1412 UseMolNumPrefix = False
1413 elif re.match("^MolNum$", OutfilesMolPrefix, re.I):
1414 UseMolNamePrefix = False
1415 UseMolNumPrefix = True
1416 else:
1417 UseMolNamePrefix = True
1418 UseMolNumPrefix = True
1419 OptionsInfo["OutfilesMolPrefix"] = OutfilesMolPrefix
1420 OptionsInfo["UseMolNamePrefix"] = UseMolNamePrefix
1421 OptionsInfo["UseMolNumPrefix"] = UseMolNumPrefix
1422
1423 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
1424
1425 CheckAndSetupOutfilesDir()
1426
1427
1428 def RetrieveOptions():
1429 """Retrieve command line arguments and options."""
1430
1431 # Get options...
1432 global Options
1433 Options = docopt(_docoptUsage_)
1434
1435 # Set current working directory to the specified directory...
1436 WorkingDir = Options["--workingdir"]
1437 if WorkingDir:
1438 os.chdir(WorkingDir)
1439
1440 # Handle examples option...
1441 if "--examples" in Options and Options["--examples"]:
1442 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
1443 sys.exit(0)
1444
1445
1446 def ValidateOptions():
1447 """Validate option values."""
1448
1449 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
1450 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol")
1451
1452 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "pml")
1453
1454 MiscUtil.ValidateOptionTextValue("--mode", Options["--mode"], "GenerateCubeFiles VisualizeCubeFiles Both")
1455 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
1456
1457 MiscUtil.ValidateOptionTextValue("--outfilesMolPrefix", Options["--outfilesMolPrefix"], "MolNum MolName Both")
1458 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
1459
1460
1461 # Setup a usage string for docopt...
1462 _docoptUsage_ = """
1463 Psi4VisualizeFrontierOrbitals.py - Visualize frontier molecular orbitals
1464
1465 Usage:
1466 Psi4VisualizeFrontierOrbitals.py [--basisSet <text>] [--infileParams <Name,Value,...>] [--methodName <text>]
1467 [--mode <GenerateCubeFiles, VisualizeCubeFiles, Both>] [--mp <yes or no>] [--mpParams <Name, Value,...>]
1468 [--outfilesDir <text>] [--outfilesMolPrefix <MolNum, MolName, Both> ] [--overwrite]
1469 [--psi4CubeFilesParams <Name,Value,...>] [--psi4OptionsParams <Name,Value,...>]
1470 [--psi4RunParams <Name,Value,...>] [--pymolViewParams <Name,Value,...>] [--quiet <yes or no>]
1471 [--reference <text>] [-w <dir>] -i <infile> -o <outfile>
1472 Psi4VisualizeFrontierOrbitals.py -h | --help | -e | --examples
1473
1474 Description:
1475 Generate and visualize frontier molecular orbitals for molecules in input
1476 file. A set of cube files, corresponding to frontier orbitals, is generated for
1477 molecules. The cube files are used to create a PyMOL visualization file for
1478 viewing volumes, meshes, and surfaces representing frontier molecular
1479 orbitals. An option is available to skip the generation of new cube files
1480 and use existing cube files to visualize frontier molecular orbitals.
1481
1482 A Psi4 XYZ format geometry string is automatically generated for each molecule
1483 in input file. It contains atom symbols and 3D coordinates for each atom in a
1484 molecule. In addition, the formal charge and spin multiplicity are present in the
1485 the geometry string. These values are either retrieved from molecule properties
1486 named 'FormalCharge' and 'SpinMultiplicty' or dynamically calculated for a
1487 molecule.
1488
1489 A set of cube and SD output files is generated for each molecule in input file
1490 as shown below:
1491
1492 Ouput dir: <OutfileRoot>_FrontierOrbitals or <OutfilesDir>
1493
1494 Closed-shell systems:
1495
1496 <MolIDPrefix>.sdf
1497 <MolIDPrefix>*HOMO.cube
1498 <MolIDPrefix>*LUMO.cube
1499
1500 Open-shell systems:
1501
1502 <MolIDPrefix>.sdf
1503 <MolIDPrefix>*DOMO.cube
1504 <MolIDPrefix>*LVMO.cube
1505 <MolIDPrefix>*SOMO.cube
1506
1507 <MolIDPrefix>: Mol<Num>_<MolName>, Mol<Num>, or <MolName>
1508
1509 Abbreviations:
1510
1511 HOMO: Highest Occupied Molecular Orbital
1512 LUMO: Lowest Unoccupied Molecular Orbital
1513
1514 DOMO: Double Occupied Molecular Orbital
1515 SOMO: Singly Occupied Molecular Orbital
1516 LVMO: Lowest Virtual (Unoccupied) Molecular Orbital
1517
1518 In addition, a <OutfileRoot>.pml is generated containing frontier molecular
1519 orbitals for all molecules in input file.
1520
1521 The supported input file formats are: Mol (.mol), SD (.sdf, .sd)
1522
1523 The supported output file formats are: PyMOL script file (.pml)
1524
1525 A variety of PyMOL groups and objects are created for visualization of frontier
1526 molecular orbitals as shown below:
1527
1528 Closed-shell systems:
1529
1530 <MoleculeID>
1531 .Molecule
1532 .Molecule
1533 .BallAndStick
1534 .HOMO
1535 .Alpha
1536 .Cube
1537 .Volume
1538 .Mesh
1539 .Negative_Phase
1540 .Positive_Phase
1541 .Surface
1542 .Negative_Phase
1543 .Positive_Phase
1544 .LUMO
1545 .Alpha
1546 .Cube
1547 .Volume
1548 .Mesh
1549 .Negative_Phase
1550 .Positive_Phase
1551 .Surface
1552 .Negative_Phase
1553 .Positive_Phase
1554 <MoleculeID>
1555 .Molecule
1556 ... ... ...
1557 .HOMO
1558 ... ... ...
1559 .LUMO
1560 ... ... ...
1561
1562 Open-shell systems:
1563
1564 <MoleculeID>
1565 .Molecule
1566 .Molecule
1567 .BallAndStick
1568 .DOMO
1569 .Alpha
1570 .Cube
1571 .Volume
1572 .Mesh
1573 .Negative_Phase
1574 .Positive_Phase
1575 .Surface
1576 .Negative_Phase
1577 .Positive_Phase
1578 .Beta
1579 .Cube
1580 .Volume
1581 .Mesh
1582 .Negative_Phase
1583 .Positive_Phase
1584 .Surface
1585 .Negative_Phase
1586 .Positive_Phase
1587 .LVMO
1588 .Alpha
1589 .Cube
1590 .Volume
1591 .Mesh
1592 .Negative_Phase
1593 .Positive_Phase
1594 .Surface
1595 .Negative_Phase
1596 .Positive_Phase
1597 .Beta
1598 .Cube
1599 .Volume
1600 .Mesh
1601 .Negative_Phase
1602 .Positive_Phase
1603 .Surface
1604 .Negative_Phase
1605 .Positive_Phase
1606 .SOMO
1607 .Alpha
1608 .Cube
1609 .Volume
1610 .Mesh
1611 .Negative_Phase
1612 .Positive_Phase
1613 .Surface
1614 .Negative_Phase
1615 .Positive_Phase
1616 .Alpha_<Num>
1617 ... ... ...
1618 <MoleculeID>
1619 .Molecule
1620 ... ... ...
1621 .DOMO
1622 ... ... ...
1623 .LVMO
1624 ... ... ...
1625 .SOMO
1626 ... ... ...
1627
1628 Options:
1629 -b, --basisSet <text> [default: auto]
1630 Basis set to use for calculating single point energy before generating
1631 cube files for frontier molecular orbitals. Default: 6-31+G** for sulfur
1632 containing molecules; Otherwise, 6-31G** [ Ref 150 ]. The specified
1633 value must be a valid Psi4 basis set. No validation is performed.
1634
1635 The following list shows a representative sample of basis sets available
1636 in Psi4:
1637
1638 STO-3G, 6-31G, 6-31+G, 6-31++G, 6-31G*, 6-31+G*, 6-31++G*,
1639 6-31G**, 6-31+G**, 6-31++G**, 6-311G, 6-311+G, 6-311++G,
1640 6-311G*, 6-311+G*, 6-311++G*, 6-311G**, 6-311+G**, 6-311++G**,
1641 cc-pVDZ, cc-pCVDZ, aug-cc-pVDZ, cc-pVDZ-DK, cc-pCVDZ-DK, def2-SVP,
1642 def2-SVPD, def2-TZVP, def2-TZVPD, def2-TZVPP, def2-TZVPPD
1643
1644 -e, --examples
1645 Print examples.
1646 -h, --help
1647 Print this help message.
1648 -i, --infile <infile>
1649 Input file name.
1650 --infileParams <Name,Value,...> [default: auto]
1651 A comma delimited list of parameter name and value pairs for reading
1652 molecules from files. The supported parameter names for different file
1653 formats, along with their default values, are shown below:
1654
1655 SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
1656
1657 -m, --methodName <text> [default: auto]
1658 Method to use for calculating single point energy before generating
1659 cube files for frontier molecular orbitals. Default: B3LYP [ Ref 150 ].
1660 The specified value must be a valid Psi4 method name. No validation is
1661 performed.
1662
1663 The following list shows a representative sample of methods available
1664 in Psi4:
1665
1666 B1LYP, B2PLYP, B2PLYP-D3BJ, B2PLYP-D3MBJ, B3LYP, B3LYP-D3BJ,
1667 B3LYP-D3MBJ, CAM-B3LYP, CAM-B3LYP-D3BJ, HF, HF-D3BJ, HF3c, M05,
1668 M06, M06-2x, M06-HF, M06-L, MN12-L, MN15, MN15-D3BJ,PBE, PBE0,
1669 PBEH3c, PW6B95, PW6B95-D3BJ, WB97, WB97X, WB97X-D, WB97X-D3BJ
1670
1671 --mode <GenerateCubeFiles, VisualizeCubeFiles, or Both> [default: Both]
1672 Generate and visualize cube files for frontier molecular orbitals. The
1673 'VisualizeCubes' value skips the generation of new cube files and uses
1674 existing cube files for visualization of molecular orbitals. Multiprocessing
1675 is not supported during 'VisualizeCubeFiles' value of '--mode' option.
1676 --mp <yes or no> [default: no]
1677 Use multiprocessing.
1678
1679 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
1680 function employing lazy RDKit data iterable. This allows processing of
1681 arbitrary large data sets without any additional requirements memory.
1682
1683 All input data may be optionally loaded into memory by mp.Pool.map()
1684 before starting worker processes in a process pool by setting the value
1685 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
1686
1687 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
1688 data mode may adversely impact the performance. The '--mpParams' section
1689 provides additional information to tune the value of 'chunkSize'.
1690 --mpParams <Name,Value,...> [default: auto]
1691 A comma delimited list of parameter name and value pairs to configure
1692 multiprocessing.
1693
1694 The supported parameter names along with their default and possible
1695 values are shown below:
1696
1697 chunkSize, auto
1698 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
1699 numProcesses, auto [ Default: mp.cpu_count() ]
1700
1701 These parameters are used by the following functions to configure and
1702 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
1703 mp.Pool.imap().
1704
1705 The chunkSize determines chunks of input data passed to each worker
1706 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
1707 The default value of chunkSize is dependent on the value of 'inputDataMode'.
1708
1709 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
1710 automatically converts RDKit data iterable into a list, loads all data into
1711 memory, and calculates the default chunkSize using the following method
1712 as shown in its code:
1713
1714 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
1715 if extra: chunkSize += 1
1716
1717 For example, the default chunkSize will be 7 for a pool of 4 worker processes
1718 and 100 data items.
1719
1720 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
1721 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
1722 data into memory. Consequently, the size of input data is not known a priori.
1723 It's not possible to estimate an optimal value for the chunkSize. The default
1724 chunkSize is set to 1.
1725
1726 The default value for the chunkSize during 'Lazy' data mode may adversely
1727 impact the performance due to the overhead associated with exchanging
1728 small chunks of data. It is generally a good idea to explicitly set chunkSize to
1729 a larger value during 'Lazy' input data mode, based on the size of your input
1730 data and number of processes in the process pool.
1731
1732 The mp.Pool.map() function waits for all worker processes to process all
1733 the data and return the results. The mp.Pool.imap() function, however,
1734 returns the the results obtained from worker processes as soon as the
1735 results become available for specified chunks of data.
1736
1737 The order of data in the results returned by both mp.Pool.map() and
1738 mp.Pool.imap() functions always corresponds to the input data.
1739 -o, --outfile <outfile>
1740 Output file name for PyMOL PML file. The PML output file, along with cube
1741 files, is generated in a local directory corresponding to '--outfilesDir'
1742 option.
1743 --outfilesDir <text> [default: auto]
1744 Directory name containing PML and cube files. Default:
1745 <OutfileRoot>_FrontierOrbitals. This directory must be present during
1746 'VisualizeCubeFiles' value of '--mode' option.
1747 --outfilesMolPrefix <MolNum, MolName, Both> [default: Both]
1748 Molecule prefix to use for the names of cube files. Possible values:
1749 MolNum, MolName, or Both. By default, both molecule number and name
1750 are used. The format of molecule prefix is as follows: MolNum - Mol<Num>;
1751 MolName - <MolName>, Both: Mol<Num>_<MolName>. Empty molecule names
1752 are ignored. Molecule numbers are used for empty molecule names.
1753 --overwrite
1754 Overwrite existing files.
1755 --psi4CubeFilesParams <Name,Value,...> [default: auto]
1756 A comma delimited list of parameter name and value pairs for generating
1757 Psi4 cube files.
1758
1759 The supported parameter names along with their default and possible
1760 values are shown below:
1761
1762 gridSpacing, 0.2, gridOverage, 4.0, isoContourThreshold, 0.85
1763
1764 gridSpacing: Grid spacing for generating cube files. Units: Bohr. A higher
1765 value reduces the size of the cube files on the disk. This option corresponds
1766 to Psi4 option CUBIC_GRID_SPACING.
1767
1768 gridOverage: Grid overage for generating cube files. Units: Bohr.This option
1769 corresponds to Psi4 option CUBIC_GRID_OVERAGE.
1770
1771 isoContourThreshold: IsoContour values for generating cube files that capture
1772 specified percent of the probability density using the least amount of grid
1773 points. Default: 0.85 (85%). This option corresponds to Psi4 option
1774 CUBEPROP_ISOCONTOUR_THRESHOLD.
1775 --psi4OptionsParams <Name,Value,...> [default: none]
1776 A comma delimited list of Psi4 option name and value pairs for setting
1777 global and module options. The names are 'option_name' for global options
1778 and 'module_name__option_name' for options local to a module. The
1779 specified option names must be valid Psi4 names. No validation is
1780 performed.
1781
1782 The specified option name and value pairs are processed and passed to
1783 psi4.set_options() as a dictionary. The supported value types are float,
1784 integer, boolean, or string. The float value string is converted into a float.
1785 The valid values for a boolean string are yes, no, true, false, on, or off.
1786 --psi4RunParams <Name,Value,...> [default: auto]
1787 A comma delimited list of parameter name and value pairs for configuring
1788 Psi4 jobs.
1789
1790 The supported parameter names along with their default and possible
1791 values are shown below:
1792
1793 MemoryInGB, 1
1794 NumThreads, 1
1795 OutputFile, auto [ Possible values: stdout, quiet, or FileName ]
1796 ScratchDir, auto [ Possivle values: DirName]
1797 RemoveOutputFile, yes [ Possible values: yes, no, true, or false]
1798
1799 These parameters control the runtime behavior of Psi4.
1800
1801 The default file name for 'OutputFile' is <InFileRoot>_Psi4.out. The PID
1802 is appended to output file name during multiprocessing as shown:
1803 <InFileRoot>_Psi4_<PIDNum>.out. The 'stdout' value for 'OutputType'
1804 sends Psi4 output to stdout. The 'quiet' or 'devnull' value suppresses
1805 all Psi4 output.
1806
1807 The default 'Yes' value of 'RemoveOutputFile' option forces the removal
1808 of any existing Psi4 before creating new files to append output from
1809 multiple Psi4 runs.
1810
1811 The option 'ScratchDir' is a directory path to the location of scratch
1812 files. The default value corresponds to Psi4 default. It may be used to
1813 override the deafult path.
1814 --pymolViewParams <Name,Value,...> [default: auto]
1815 A comma delimited list of parameter name and value pairs for visualizing
1816 cube files in PyMOL.
1817
1818 contourColor1, red, contourColor2, blue,
1819 contourLevel1, auto, contourLevel2, auto,
1820 contourLevelAutoAt, 0.75,
1821 displayMolecule, BallAndStick, displaySphereScale, 0.2,
1822 displayStickRadius, 0.1, hideHydrogens, yes,
1823 meshWidth, 0.5, meshQuality, 2,
1824 surfaceQuality, 2, surfaceTransparency, 0.25,
1825 volumeColorRamp, auto, volumeColorRampOpacity,0.2
1826 volumeContourWindowFactor,0.05
1827
1828 contourColor1 and contourColor2: Color to use for visualizing volumes,
1829 meshes, and surfaces corresponding to the negative and positive values
1830 in cube files. The specified values must be valid PyMOL color names. No
1831 validation is performed.
1832
1833 contourLevel1 and contourLevel2: Contour levels to use for visualizing
1834 volumes, meshes, and surfaces corresponding to the negative and positive
1835 values in cube files. Default: auto. The specified values for contourLevel1
1836 and contourLevel2 must be negative and positive numbers.
1837
1838 The contour levels are automatically calculated by default. The isocontour
1839 range for specified percent of the density is retrieved from the cube files.
1840 The contour levels are set at 'contourLevelAutoAt' of the absolute maximum
1841 value of the isocontour range. For example: contour levels are set to plus and
1842 minus 0.03 at 'contourLevelAutoAt' of 0.5 for isocontour range of -0.06 to
1843 0.06 covering specified percent of the density.
1844
1845 contourLevelAutoAt: Set contour levels at specified fraction of the absolute
1846 maximum value of the isocontour range retrieved from the cube files. This
1847 option is only used during the automatic calculations of the contour levels.
1848
1849 hideHydrogens: Hide hydrogens in molecules. Default: yes. Possible
1850 values: yes or no.
1851
1852 displayMolecule: Display mode for molecules. Possible values: Sticks or
1853 BallAndStick. Both displays objects are created for molecules.
1854
1855 displaySphereScale: Sphere scale for displaying molecule during
1856 BallAndStick display.
1857
1858 displayStickRadius: Stick radius for displaying molecule during Sticks
1859 and BallAndStick display.
1860
1861 hideHydrogens: Hide hydrogens in molecules. Default: yes. Possible
1862 values: yes or no.
1863 meshWidth: Line width for mesh lines to visualize cube files.
1864
1865 meshQuality: Mesh quality for meshes to visualize cube files. The
1866 higher values represents better quality.
1867
1868 meshWidth: Line width for mesh lines to visualize cube files.
1869
1870 surfaceQuality: Surface quality for surfaces to visualize cube files.
1871 The higher values represents better quality.
1872
1873 surfaceTransparency: Surface transparency for surfaces to visualize cube
1874 files.
1875
1876 volumeColorRamp: Name of a PyMOL volume color ramp to use for visualizing
1877 cube files. Default name(s): <OutfielsMolPrefix>_psi4_cube_mo or psi4_cube_mo
1878 The default volume color ramps are automatically generated using contour
1879 levels and colors during 'auto' value of 'volumeColorRamp'. An explicitly
1880 specified value must be a valid PyMOL volume color ramp. No validation is
1881 preformed.
1882
1883 VolumeColorRampOpacity: Opacity for generating volume color ramps
1884 for visualizing cube files. This value is equivalent to 1 minus Transparency.
1885
1886 volumeContourWindowFactor: Fraction of contour level representing window
1887 widths around contour levels during generation of volume color ramps for
1888 visualizing cube files. For example, the value of 0.05 implies a ramp window
1889 size of 0.0015 at contour level of 0.03.
1890 -q, --quiet <yes or no> [default: no]
1891 Use quiet mode. The warning and information messages will not be printed.
1892 -r, --reference <text> [default: auto]
1893 Reference wave function to use for calculating single point energy before
1894 generating cube files for frontier molecular orbitals. Default: RHF or UHF.
1895 The default values are Restricted Hartree-Fock (RHF) for closed-shell molecules
1896 with all electrons paired and Unrestricted artree-Fock (UHF) for open-shell
1897 molecules with unpaired electrons.
1898
1899 The specified value must be a valid Psi4 reference wave function. No validation
1900 is performed. For example: ROHF, CUHF, RKS, etc.
1901
1902 The spin multiplicity determines the default value of reference wave function
1903 for input molecules. It is calculated from number of free radical electrons using
1904 Hund's rule of maximum multiplicity defined as 2S + 1 where S is the total
1905 electron spin. The total spin is 1/2 the number of free radical electrons in a
1906 molecule. The value of 'SpinMultiplicity' molecule property takes precedence
1907 over the calculated value of spin multiplicity.
1908 -w, --workingdir <dir>
1909 Location of working directory which defaults to the current directory.
1910
1911 Examples:
1912 To generate and visualize frontier molecular orbitals based on a single point
1913 energy calculation using B3LYP/6-31G** and B3LYP/6-31+G** for non-sulfur
1914 and sulfur containing molecules in a SD file with 3D structures, use RHF and
1915 UHF for closed-shell and open-shell molecules, and write a new PML file, type:
1916
1917 % Psi4VisualizeFrontierOrbitals.py -i Psi4Sample3D.sdf
1918 -o Psi4Sample3DOut.pml
1919
1920 To run the first example to only generate cube files and skip generation of
1921 a PML file to visualize frontier molecular orbitals, type:
1922
1923 % Psi4VisualizeFrontierOrbitals.py --mode GenerateCubeFiles
1924 -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1925
1926 To run the first example to skip generation of cube files and use existing cube
1927 files to visualize frontier molecular orbitals and write out a PML file, type:
1928
1929 % Psi4VisualizeFrontierOrbitals.py --mode VisualizeCubeFiles
1930 -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1931
1932 To run the first example in multiprocessing mode on all available CPUs
1933 without loading all data into memory and write out a PML file, type:
1934
1935 % Psi4VisualizeFrontierOrbitals.py --mp yes -i Psi4Sample3D.sdf
1936 -o Psi4Sample3DOut.pml
1937
1938 To run the first example in multiprocessing mode on all available CPUs
1939 by loading all data into memory and write out a PML file, type:
1940
1941 % Psi4VisualizeFrontierOrbitals.py --mp yes --mpParams "inputDataMode,
1942 InMemory" -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1943
1944 To run the first example in multiprocessing mode on all available CPUs
1945 without loading all data into memory along with multiple threads for each
1946 Psi4 run and write out a SD file, type:
1947
1948 % Psi4VisualizeFrontierOrbitals.py --mp yes --psi4RunParams
1949 "NumThreads,2" -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1950
1951 To run the first example in using a specific set of parameters to generate and
1952 visualize frontier molecular orbitals and write out a PML file, type:
1953
1954 % Psi4VisualizeFrontierOrbitals.py --mode both -m SCF -b aug-cc-pVDZ
1955 --psi4CubeFilesParams "gridSpacing, 0.2, gridOverage, 4.0"
1956 --psi4RunParams "MemoryInGB, 2" --pymolViewParams "contourColor1,
1957 red, contourColor2, blue,contourLevel1, -0.04, contourLevel2, 0.04,
1958 contourLevelAutoAt, 0.75,volumeColorRamp, auto,
1959 volumeColorRampOpacity,0.25, volumeContourWindowFactor,0.05"
1960 -i Psi4Sample3D.sdf -o Psi4Sample3DOut.pml
1961
1962 Author:
1963 Manish Sud(msud@san.rr.com)
1964
1965 See also:
1966 Psi4PerformMinimization.py, Psi4GenerateConformers.py,
1967 Psi4VisualizeDualDescriptors.py, Psi4VisualizeElectrostaticPotential.py
1968
1969 Copyright:
1970 Copyright (C) 2026 Manish Sud. All rights reserved.
1971
1972 The functionality available in this script is implemented using Psi4, an
1973 open source quantum chemistry software package, and RDKit, an open
1974 source toolkit for cheminformatics developed by Greg Landrum.
1975
1976 This file is part of MayaChemTools.
1977
1978 MayaChemTools is free software; you can redistribute it and/or modify it under
1979 the terms of the GNU Lesser General Public License as published by the Free
1980 Software Foundation; either version 3 of the License, or (at your option) any
1981 later version.
1982
1983 """
1984
1985 if __name__ == "__main__":
1986 main()
