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