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