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