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