1 #!/bin/env python
2 #
3 # File: Psi4CalculateEnergy.py
4 # Author: Manish Sud <msud@san.rr.com>
5 #
6 # Copyright (C) 2026 Manish Sud. All rights reserved.
7 #
8 # The functionality available in this script is implemented using Psi4, an
9 # open source quantum chemistry software package, and RDKit, an open
10 # source toolkit for cheminformatics developed by Greg Landrum.
11 #
12 # This file is part of MayaChemTools.
13 #
14 # MayaChemTools is free software; you can redistribute it and/or modify it under
15 # the terms of the GNU Lesser General Public License as published by the Free
16 # Software Foundation; either version 3 of the License, or (at your option) any
17 # later version.
18 #
19 # MayaChemTools is distributed in the hope that it will be useful, but without
20 # any warranty; without even the implied warranty of merchantability of fitness
21 # for a particular purpose. See the GNU Lesser General Public License for more
22 # details.
23 #
24 # You should have received a copy of the GNU Lesser General Public License
25 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
26 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
27 # Boston, MA, 02111-1307, USA.
28 #
29
30 from __future__ import print_function
31
32 import os
33 import sys
34 import time
35 import re
36 import shutil
37 import multiprocessing as mp
38
39 # Psi4 imports...
40 if hasattr(shutil, "which") and shutil.which("psi4") is None:
41 sys.stderr.write("\nWarning: Failed to find 'psi4' in your PATH indicating potential issues with your\n")
42 sys.stderr.write("Psi4 environment. The 'import psi4' directive in the global scope of the script\n")
43 sys.stderr.write("interferes with the multiprocessing functionality. It is imported later in the\n")
44 sys.stderr.write("local scope during the execution of the script and may fail. Check/update your\n")
45 sys.stderr.write("Psi4 environment and try again.\n\n")
46
47 # RDKit imports...
48 try:
49 from rdkit import rdBase
50 from rdkit import Chem
51 except ImportError as ErrMsg:
52 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
53 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
54 sys.exit(1)
55
56 # MayaChemTools imports...
57 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
58 try:
59 from docopt import docopt
60 import MiscUtil
61 import Psi4Util
62 import RDKitUtil
63 except ImportError as ErrMsg:
64 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
65 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
66 sys.exit(1)
67
68 ScriptName = os.path.basename(sys.argv[0])
69 Options = {}
70 OptionsInfo = {}
71
72
73 def main():
74 """Start execution of the script."""
75
76 MiscUtil.PrintInfo(
77 "\n%s (Psi4: Imported later; RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
78 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime())
79 )
80
81 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
82
83 # Retrieve command line arguments and options...
84 RetrieveOptions()
85
86 # Process and validate command line arguments and options...
87 ProcessOptions()
88
89 # Perform actions required by the script...
90 CalculateInteractionEnergy()
91
92 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
93 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
94
95
96 def CalculateInteractionEnergy():
97 """Calculate interaction energy."""
98
99 # Setup a molecule reader...
100 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
101 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
102
103 # Set up a molecule writer...
104 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"])
105 if Writer is None:
106 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"])
107 MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["Outfile"])
108
109 MolCount, ValidMolCount, EnergyFailedCount = ProcessMolecules(Mols, Writer)
110
111 if Writer is not None:
112 Writer.close()
113
114 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
115 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
116 MiscUtil.PrintInfo("Number of molecules failed during energy calculation: %d" % EnergyFailedCount)
117 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + EnergyFailedCount))
118
119
120 def ProcessMolecules(Mols, Writer):
121 """Process and calculate energy of molecules."""
122
123 if OptionsInfo["MPMode"]:
124 return ProcessMoleculesUsingMultipleProcesses(Mols, Writer)
125 else:
126 return ProcessMoleculesUsingSingleProcess(Mols, Writer)
127
128
129 def ProcessMoleculesUsingSingleProcess(Mols, Writer):
130 """Process and calculate energy of molecules using a single process."""
131
132 # Intialize Psi4...
133 MiscUtil.PrintInfo("\nInitializing Psi4...")
134 Psi4Handle = Psi4Util.InitializePsi4(
135 Psi4RunParams=OptionsInfo["Psi4RunParams"],
136 Psi4OptionsParams=OptionsInfo["Psi4OptionsParams"],
137 PrintVersion=True,
138 PrintHeader=True,
139 )
140 OptionsInfo["psi4"] = Psi4Handle
141
142 # Setup conversion factor for energy units...
143 SetupEnergyConversionFactor(Psi4Handle)
144
145 MiscUtil.PrintInfo("\nCalculating interaction energy...")
146
147 (MolCount, ValidMolCount, EnergyFailedCount) = [0] * 3
148 for Mol in Mols:
149 MolCount += 1
150
151 if not CheckAndValidateMolecule(Mol, MolCount):
152 continue
153
154 # Setup a Psi4 molecule...
155 Psi4Mol = SetupPsi4Mol(Psi4Handle, Mol, MolCount)
156 if Psi4Mol is None:
157 continue
158
159 ValidMolCount += 1
160
161 CalcStatus, Energy, EnergyDetails = CalculateMoleculeInteractionEnergy(Psi4Handle, Psi4Mol, Mol, MolCount)
162
163 if not CalcStatus:
164 if not OptionsInfo["QuietMode"]:
165 MiscUtil.PrintWarning(
166 "Failed to calculate interaction energy for molecule %s" % RDKitUtil.GetMolName(Mol, MolCount)
167 )
168
169 EnergyFailedCount += 1
170 continue
171
172 WriteMolecule(Writer, Mol, Energy, EnergyDetails)
173
174 return (MolCount, ValidMolCount, EnergyFailedCount)
175
176
177 def ProcessMoleculesUsingMultipleProcesses(Mols, Writer):
178 """Process and calculate energy of molecules using process."""
179
180 MPParams = OptionsInfo["MPParams"]
181
182 # Setup data for initializing a worker process...
183 InitializeWorkerProcessArgs = (
184 MiscUtil.ObjectToBase64EncodedString(Options),
185 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
186 )
187
188 # Setup a encoded mols data iterable for a worker process...
189 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
190
191 # Setup process pool along with data initialization for each process...
192 if not OptionsInfo["QuietMode"]:
193 MiscUtil.PrintInfo(
194 "\nConfiguring multiprocessing using %s method..."
195 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
196 )
197 MiscUtil.PrintInfo(
198 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
199 % (
200 MPParams["NumProcesses"],
201 MPParams["InputDataMode"],
202 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
203 )
204 )
205
206 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
207
208 # Start processing...
209 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
210 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
211 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
212 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
213 else:
214 MiscUtil.PrintError(
215 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
216 )
217
218 # Print out Psi4 version in the main process...
219 MiscUtil.PrintInfo("\nInitializing Psi4...\n")
220 Psi4Handle = Psi4Util.InitializePsi4(PrintVersion=True, PrintHeader=False)
221 OptionsInfo["psi4"] = Psi4Handle
222
223 (MolCount, ValidMolCount, EnergyFailedCount) = [0] * 3
224 for Result in Results:
225 MolCount += 1
226 MolIndex, EncodedMol, CalcStatus, Energy, EnergyDetails = Result
227
228 if EncodedMol is None:
229 continue
230
231 ValidMolCount += 1
232
233 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
234
235 if not CalcStatus:
236 if not OptionsInfo["QuietMode"]:
237 MolName = RDKitUtil.GetMolName(Mol, MolCount)
238 MiscUtil.PrintWarning("Failed to calculate energy for molecule %s" % MolName)
239
240 EnergyFailedCount += 1
241 continue
242
243 WriteMolecule(Writer, Mol, Energy, EnergyDetails)
244
245 return (MolCount, ValidMolCount, EnergyFailedCount)
246
247
248 def InitializeWorkerProcess(*EncodedArgs):
249 """Initialize data for a worker process."""
250
251 global Options, OptionsInfo
252
253 if not OptionsInfo["QuietMode"]:
254 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
255
256 # Decode Options and OptionInfo...
257 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
258 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
259
260 # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
261 # output files for each process...
262 OptionsInfo["Psi4Initialized"] = False
263
264
265 def InitializePsi4ForWorkerProcess():
266 """Initialize Psi4 for a worker process."""
267
268 if OptionsInfo["Psi4Initialized"]:
269 return
270
271 OptionsInfo["Psi4Initialized"] = True
272
273 # Update output file...
274 OptionsInfo["Psi4RunParams"]["OutputFile"] = Psi4Util.UpdatePsi4OutputFileUsingPID(
275 OptionsInfo["Psi4RunParams"]["OutputFile"], os.getpid()
276 )
277
278 # Intialize Psi4...
279 OptionsInfo["psi4"] = Psi4Util.InitializePsi4(
280 Psi4RunParams=OptionsInfo["Psi4RunParams"],
281 Psi4OptionsParams=OptionsInfo["Psi4OptionsParams"],
282 PrintVersion=False,
283 PrintHeader=True,
284 )
285
286 # Setup conversion factor for energy units...
287 SetupEnergyConversionFactor(OptionsInfo["psi4"])
288
289
290 def WorkerProcess(EncodedMolInfo):
291 """Process data for a worker process."""
292
293 if not OptionsInfo["Psi4Initialized"]:
294 InitializePsi4ForWorkerProcess()
295
296 MolIndex, EncodedMol = EncodedMolInfo
297
298 CalcStatus = False
299 Energy = None
300 EnergyDetails = None
301
302 if EncodedMol is None:
303 return [MolIndex, None, CalcStatus, Energy, EnergyDetails]
304
305 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
306 MolCount = MolIndex + 1
307
308 if not CheckAndValidateMolecule(Mol, MolCount):
309 return [MolIndex, None, CalcStatus, Energy, EnergyDetails]
310
311 # Setup a Psi4 molecule...
312 Psi4Mol = SetupPsi4Mol(OptionsInfo["psi4"], Mol, MolCount)
313 if Psi4Mol is None:
314 return [MolIndex, None, CalcStatus, Energy, EnergyDetails]
315
316 CalcStatus, Energy, EnergyDetails = CalculateMoleculeInteractionEnergy(OptionsInfo["psi4"], Psi4Mol, Mol, MolCount)
317
318 return [
319 MolIndex,
320 RDKitUtil.MolToBase64EncodedMolString(
321 Mol, PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps
322 ),
323 CalcStatus,
324 Energy,
325 EnergyDetails,
326 ]
327
328
329 def CalculateMoleculeInteractionEnergy(Psi4Handle, Psi4Mol, Mol, MolNum=None):
330 """Calculate interaction energy for a molecule containing two components."""
331
332 # Setup reference wave function...
333 Reference = SetupReferenceWavefunction(Mol)
334 Psi4Handle.set_options({"Reference": Reference})
335
336 # Setup method name and basis set...
337 MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
338 MethodAndBasisSet = Psi4Util.JoinMethodNameAndBasisSet(MethodName, BasisSet)
339
340 Status, Energy, EnergyDetails = [False, None, None]
341 try:
342 if OptionsInfo["BSSEType"] is None:
343 Energy = Psi4Handle.energy(MethodAndBasisSet, molecule=Psi4Mol, return_wfn=False)
344 else:
345 Energy = Psi4Handle.energy(
346 MethodAndBasisSet, bsse_type=OptionsInfo["BSSEType"], molecule=Psi4Mol, return_wfn=False
347 )
348 Status = True
349 except Exception as ErrMsg:
350 if not OptionsInfo["QuietMode"]:
351 MiscUtil.PrintWarning("Failed to calculate interaction energy:\n%s\n" % ErrMsg)
352
353 # Retrieve SAPT energy details...
354 if Status and OptionsInfo["SAPTMode"]:
355 EnergyTypeMap = OptionsInfo["SAPTEnergyIDsToTypes"]
356 EnergyDetails = {}
357 for EnergyType in EnergyTypeMap:
358 try:
359 EnergyDetails[EnergyType] = Psi4Handle.variable(EnergyTypeMap[EnergyType])
360 except Exception as ErrMsg:
361 EnergyDetails[EnergyType] = None
362 if not OptionsInfo["QuietMode"]:
363 MiscUtil.PrintWarning(
364 "Failed to retrieve value for interaction energy components %s:\n%s\n" % (EnergyType, ErrMsg)
365 )
366
367 # Convert energy units...
368 if Status and OptionsInfo["ApplyEnergyConversionFactor"]:
369 if Energy is not None:
370 Energy = Energy * OptionsInfo["EnergyConversionFactor"]
371 if EnergyDetails is not None:
372 for EnergyType in EnergyDetails:
373 if EnergyDetails[EnergyType] is not None:
374 EnergyDetails[EnergyType] = EnergyDetails[EnergyType] * OptionsInfo["EnergyConversionFactor"]
375
376 # Clean up
377 PerformPsi4Cleanup(Psi4Handle)
378
379 return (Status, Energy, EnergyDetails)
380
381
382 def SetupPsi4Mol(Psi4Handle, Mol, MolCount=None):
383 """Setup a Psi4 molecule object."""
384
385 MolGeometry = RDKitUtil.GetPsi4XYZFormatString(Mol, NoCom=True, NoReorient=True, CheckFragments=True)
386
387 try:
388 Psi4Mol = Psi4Handle.geometry(MolGeometry)
389 except Exception as ErrMsg:
390 Psi4Mol = None
391 if not OptionsInfo["QuietMode"]:
392 MiscUtil.PrintWarning("Failed to create Psi4 molecule from geometry string: %s\n" % ErrMsg)
393 MolName = RDKitUtil.GetMolName(Mol, MolCount)
394 MiscUtil.PrintWarning("Ignoring molecule: %s" % MolName)
395
396 return Psi4Mol
397
398
399 def PerformPsi4Cleanup(Psi4Handle):
400 """Perform clean up."""
401
402 # Clean up after Psi4 run ...
403 Psi4Handle.core.clean()
404
405 # Clean up any leftover scratch files...
406 if OptionsInfo["MPMode"]:
407 Psi4Util.RemoveScratchFiles(Psi4Handle, OptionsInfo["Psi4RunParams"]["OutputFile"])
408
409
410 def CheckAndValidateMolecule(Mol, MolCount=None):
411 """Validate molecule for Psi4 calculations."""
412
413 if Mol is None:
414 if not OptionsInfo["QuietMode"]:
415 MiscUtil.PrintInfo("\nProcessing molecule number %s..." % MolCount)
416 return False
417
418 MolName = RDKitUtil.GetMolName(Mol, MolCount)
419 if not OptionsInfo["QuietMode"]:
420 MiscUtil.PrintInfo("\nProcessing molecule %s..." % MolName)
421
422 if RDKitUtil.IsMolEmpty(Mol):
423 if not OptionsInfo["QuietMode"]:
424 MiscUtil.PrintWarning("Ignoring empty molecule: %s\n" % MolName)
425 return False
426
427 if not RDKitUtil.ValidateElementSymbols(RDKitUtil.GetAtomSymbols(Mol)):
428 if not OptionsInfo["QuietMode"]:
429 MiscUtil.PrintWarning("Ignoring molecule containing invalid element symbols: %s\n" % MolName)
430 return False
431
432 # Check for number of fragments...
433 NumFragments = RDKitUtil.GetNumFragments(Mol)
434 if OptionsInfo["SAPTMode"] or OptionsInfo["SNSMP2Mode"]:
435 if NumFragments != 2:
436 if not OptionsInfo["QuietMode"]:
437 MiscUtil.PrintWarning(
438 "Ignoring molecule containing invalid number of fragments: %s. It must contain exaclty 2 fragments for calculation of interaction energy using SAPT or SNS-MP2 methodology.\n"
439 % NumFragments
440 )
441 return False
442 else:
443 if NumFragments == 1:
444 if not OptionsInfo["QuietMode"]:
445 MiscUtil.PrintWarning(
446 "Ignoring molecule containing invalid number of fragments: %s. It must contain more than 1 fragment for calculation of interaction energy using counterpoise correctio methodology.\n"
447 % NumFragments
448 )
449 return False
450
451 # Check for 3D flag...
452 if not Mol.GetConformer().Is3D():
453 if not OptionsInfo["QuietMode"]:
454 MiscUtil.PrintWarning("3D tag is not set for molecule: %s\n" % MolName)
455
456 # Check for missing hydrogens...
457 if RDKitUtil.AreHydrogensMissingInMolecule(Mol):
458 if not OptionsInfo["QuietMode"]:
459 MiscUtil.PrintWarning("Missing hydrogens in molecule: %s\n" % MolName)
460
461 return True
462
463
464 def SetupMethodNameAndBasisSet(Mol):
465 """Setup method name and basis set."""
466
467 return (OptionsInfo["MethodName"], OptionsInfo["BasisSet"])
468
469
470 def SetupReferenceWavefunction(Mol):
471 """Setup reference wave function."""
472
473 Reference = OptionsInfo["Reference"]
474 if OptionsInfo["ReferenceAuto"]:
475 Reference = "UHF" if (RDKitUtil.GetSpinMultiplicity(Mol) > 1) else "RHF"
476
477 return Reference
478
479
480 def SetupEnergyConversionFactor(Psi4Handle):
481 """Setup converstion factor for energt units. The Psi4 energy units are Hartrees."""
482
483 EnergyUnits = OptionsInfo["EnergyUnits"]
484
485 ApplyConversionFactor = True
486 if re.match(r"^kcal\/mol$", EnergyUnits, re.I):
487 ConversionFactor = Psi4Handle.constants.hartree2kcalmol
488 elif re.match(r"^kJ\/mol$", EnergyUnits, re.I):
489 ConversionFactor = Psi4Handle.constants.hartree2kJmol
490 elif re.match("^eV$", EnergyUnits, re.I):
491 ConversionFactor = Psi4Handle.constants.hartree2ev
492 else:
493 ApplyConversionFactor = False
494 ConversionFactor = 1.0
495
496 OptionsInfo["ApplyEnergyConversionFactor"] = ApplyConversionFactor
497 OptionsInfo["EnergyConversionFactor"] = ConversionFactor
498
499
500 def WriteMolecule(Writer, Mol, Energy, EnergyDetails):
501 """Write molecule."""
502
503 # Setup energy...
504 EnergyValue = "%.*f" % (OptionsInfo["Precision"], Energy) if Energy is not None else ""
505 Mol.SetProp(OptionsInfo["EnergyDataFieldLabel"], EnergyValue)
506
507 # Setup SAPT energy details...
508 if EnergyDetails is not None and OptionsInfo["SAPTMode"]:
509 for EnergyID in OptionsInfo["SAPTEnergyIDs"]:
510 Energy = EnergyDetails[EnergyID] if EnergyID in EnergyDetails else None
511 EnergyValue = "%.*f" % (OptionsInfo["Precision"], Energy) if Energy is not None else ""
512
513 Mol.SetProp(OptionsInfo["SAPTEnergyIDsToLabels"][EnergyID], EnergyValue)
514
515 Writer.write(Mol)
516
517
518 def ProcessMethodNameAndBasisSetOptions():
519 """Process method name and basis set options."""
520
521 MethodName = Options["--methodName"]
522 MethodNameAuto = True if re.match("^auto$", MethodName, re.I) else False
523 if MethodNameAuto:
524 MethodName = "SAPT0"
525 SAPTMode = True if re.match("^SAPT", MethodName, re.I) else False
526
527 SNSMP2Mode = True if re.match("^SNS-MP2$", MethodName, re.I) else False
528
529 BasisSet = Options["--basisSet"]
530 BasisSetAuto = True if re.match("^auto$", BasisSet, re.I) else False
531 if BasisSetAuto:
532 if SAPTMode:
533 BasisSet = "jun-cc-pVDZ"
534 elif SNSMP2Mode:
535 BasisSet = ""
536 else:
537 MiscUtil.PrintError(
538 'A basis set must be explicitly specified using "-b, --basisSet" option for, "%s", value of "-m, --methodName". '
539 % (MethodName)
540 )
541
542 if SNSMP2Mode and not MiscUtil.IsEmpty(BasisSet):
543 MiscUtil.PrintError(
544 'The basis set value, "%s", specified using "-b, --basisSet" option is not valid. It must be an empty string SNS-MP2 calculations.'
545 % (Options["--basisSet"])
546 )
547
548 OptionsInfo["MethodName"] = MethodName
549 OptionsInfo["MethodNameAuto"] = MethodNameAuto
550 OptionsInfo["SAPTMode"] = SAPTMode
551 OptionsInfo["SNSMP2Mode"] = SNSMP2Mode
552
553 OptionsInfo["BasisSet"] = BasisSet
554 OptionsInfo["BasisSetAuto"] = BasisSetAuto
555
556
557 def ProcessBSSETypeOption():
558 """Process basis set superposition energy correction option."""
559
560 BSSEType = (
561 None
562 if re.match("^None$", Options["--bsseType"], re.I) or MiscUtil.IsEmpty(Options["--bsseType"])
563 else Options["--bsseType"]
564 )
565 BSSETypeAuto = True if re.match("^auto$", Options["--bsseType"], re.I) else False
566
567 if BSSETypeAuto:
568 if OptionsInfo["SAPTMode"] or OptionsInfo["SNSMP2Mode"]:
569 BSSEType = None
570 elif re.match("^HF3c$", OptionsInfo["MethodName"], re.I):
571 BSSEType = "noCP"
572 else:
573 BSSEType = "CP"
574 else:
575 if OptionsInfo["SAPTMode"] or OptionsInfo["SNSMP2Mode"]:
576 if BSSEType is not None:
577 MiscUtil.PrintError(
578 'The BSSE type value, "%s", specified using "--bsseType" option is not valid for SAPT or SNS-MP2 calculations. It must be set to None.'
579 % (Options["--bsseType"])
580 )
581 else:
582 if BSSEType is None:
583 MiscUtil.PrintWarning(
584 'The BSSE type value, "%s", specified using "--bsseType" option is not valid for "%s" calculations. Possible values: CP, noCP or VMFC'
585 % (Options["--bsseType"], OptionsInfo["MethodName"])
586 )
587
588 OptionsInfo["BSSEType"] = BSSEType
589 OptionsInfo["BSSETypeAuto"] = BSSETypeAuto
590
591
592 def ProcessSAPTEnergyDataFieldLabelsOption():
593 """Process data field labels for SAPT interaction energy components."""
594
595 ParamsOptionName = "--energySAPTDataFieldLabels"
596 ParamsOptionValue = Options["--energySAPTDataFieldLabels"]
597
598 ParamsIDs = ["ElectrostaticEnergy", "ExchangeEnergy", "InductionEnergy", "DispersionEnergy"]
599 ParamsIDsToTypes = {
600 "ElectrostaticEnergy": "SAPT ELST ENERGY",
601 "ExchangeEnergy": "SAPT EXCH ENERGY",
602 "InductionEnergy": "SAPT IND ENERGY",
603 "DispersionEnergy": "SAPT DISP ENERGY",
604 }
605
606 UnitsLabel = "(%s)" % Options["--energyUnits"]
607 ParamsIDsToLabels = {
608 "ElectrostaticEnergy": "Psi4_SAPT_Electrostatic_Energy %s" % UnitsLabel,
609 "ExchangeEnergy": "Psi4_SAPT_Exchange_Energy %s" % UnitsLabel,
610 "InductionEnergy": "Psi4_SAPT_Induction_Energy %s" % UnitsLabel,
611 "DispersionEnergy": "Psi4_SAPT_Dispersion_Energy %s" % UnitsLabel,
612 }
613
614 if re.match("^auto$", ParamsOptionValue, re.I):
615 OptionsInfo["SAPTEnergyIDs"] = ParamsIDs
616 OptionsInfo["SAPTEnergyIDsToTypes"] = ParamsIDsToTypes
617 OptionsInfo["SAPTEnergyIDsToLabels"] = ParamsIDsToLabels
618 return
619
620 # Setup a canonical paramater names...
621 ValidParamNames = []
622 CanonicalParamNamesMap = {}
623 for ParamName in sorted(ParamsIDsToLabels):
624 ValidParamNames.append(ParamName)
625 CanonicalParamNamesMap[ParamName.lower()] = ParamName
626
627 ParamsOptionValue = ParamsOptionValue.strip()
628 if not ParamsOptionValue:
629 MiscUtil.PrintError('No valid parameter name and value pairs specified using "%s" option' % ParamsOptionName)
630
631 ParamsOptionValueWords = ParamsOptionValue.split(",")
632 if len(ParamsOptionValueWords) % 2:
633 MiscUtil.PrintError(
634 'The number of comma delimited paramater names and values, %d, specified using "%s" option must be an even number.'
635 % (len(ParamsOptionValueWords), ParamsOptionName)
636 )
637
638 # Validate paramater name and value pairs...
639 for Index in range(0, len(ParamsOptionValueWords), 2):
640 Name = ParamsOptionValueWords[Index].strip()
641 Value = ParamsOptionValueWords[Index + 1].strip()
642
643 CanonicalName = Name.lower()
644 if CanonicalName not in CanonicalParamNamesMap:
645 MiscUtil.PrintError(
646 'The parameter name, %s, specified using "%s" is not a valid name. Supported parameter names: %s'
647 % (Name, ParamsOptionName, " ".join(ValidParamNames))
648 )
649
650 ParamName = CanonicalParamNamesMap[CanonicalName]
651 ParamValue = Value
652
653 # Set value...
654 ParamsIDsToLabels[ParamName] = ParamValue
655
656 OptionsInfo["SAPTEnergyIDs"] = ParamsIDs
657 OptionsInfo["SAPTEnergyIDsToTypes"] = ParamsIDsToTypes
658 OptionsInfo["SAPTEnergyIDsToLabels"] = ParamsIDsToLabels
659
660
661 def ProcessOptions():
662 """Process and validate command line arguments and options."""
663
664 MiscUtil.PrintInfo("Processing options...")
665
666 # Validate options...
667 ValidateOptions()
668
669 OptionsInfo["Infile"] = Options["--infile"]
670 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
671 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
672 "--infileParams",
673 Options["--infileParams"],
674 InfileName=Options["--infile"],
675 ParamsDefaultInfo=ParamsDefaultInfoOverride,
676 )
677
678 OptionsInfo["Outfile"] = Options["--outfile"]
679 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters(
680 "--outfileParams", Options["--outfileParams"]
681 )
682
683 OptionsInfo["Overwrite"] = Options["--overwrite"]
684
685 # Method and basis set...
686 ProcessMethodNameAndBasisSetOptions()
687
688 # BSSEType option...
689 ProcessBSSETypeOption()
690
691 # Reference wavefunction...
692 OptionsInfo["Reference"] = Options["--reference"]
693 OptionsInfo["ReferenceAuto"] = True if re.match("^auto$", Options["--reference"], re.I) else False
694
695 # Run and options parameters...
696 OptionsInfo["Psi4OptionsParams"] = Psi4Util.ProcessPsi4OptionsParameters(
697 "--psi4OptionsParams", Options["--psi4OptionsParams"]
698 )
699 OptionsInfo["Psi4RunParams"] = Psi4Util.ProcessPsi4RunParameters(
700 "--psi4RunParams", Options["--psi4RunParams"], InfileName=OptionsInfo["Infile"]
701 )
702
703 # Interaction energy units...
704 OptionsInfo["EnergyUnits"] = Options["--energyUnits"]
705
706 # Total interaction energy label...
707 EnergyDataFieldLabel = Options["--energyDataFieldLabel"]
708 if re.match("^auto$", EnergyDataFieldLabel, re.I):
709 EnergyLabel = "Psi4_Interaction_Energy"
710 if OptionsInfo["SAPTMode"]:
711 EnergyLabel = "Psi4_SAPT_Interaction_Energy"
712 elif OptionsInfo["SNSMP2Mode"]:
713 EnergyLabel = "Psi4_SNS-MP2_Interaction_Energy"
714 EnergyDataFieldLabel = "%s %s" % (EnergyLabel, Options["--energyUnits"])
715
716 OptionsInfo["EnergyDataFieldLabel"] = EnergyDataFieldLabel
717
718 ProcessSAPTEnergyDataFieldLabelsOption()
719
720 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
721 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
722
723 OptionsInfo["Precision"] = int(Options["--precision"])
724 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
725
726
727 def RetrieveOptions():
728 """Retrieve command line arguments and options."""
729
730 # Get options...
731 global Options
732 Options = docopt(_docoptUsage_)
733
734 # Set current working directory to the specified directory...
735 WorkingDir = Options["--workingdir"]
736 if WorkingDir:
737 os.chdir(WorkingDir)
738
739 # Handle examples option...
740 if "--examples" in Options and Options["--examples"]:
741 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
742 sys.exit(0)
743
744
745 def ValidateOptions():
746 """Validate option values."""
747
748 MiscUtil.ValidateOptionTextValue("--bsseType", Options["--bsseType"], "CP noCP VMFC None auto")
749
750 MiscUtil.ValidateOptionTextValue("--energyUnits", Options["--energyUnits"], "Hartrees kcal/mol kJ/mol eV")
751
752 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
753 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol")
754
755 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
756 MiscUtil.ValidateOptionsOutputFileOverwrite(
757 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"]
758 )
759 MiscUtil.ValidateOptionsDistinctFileNames(
760 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"]
761 )
762
763 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
764 MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0})
765
766 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
767
768
769 # Setup a usage string for docopt...
770 _docoptUsage_ = """
771 Psi4CalculateInteractionEnergy.py - Calculate interaction energy
772
773 Usage:
774 Psi4CalculateInteractionEnergy.py [--basisSet <text>] [--bsseType <CP, noCP, VMFC, or None>]
775 [--energyDataFieldLabel <text>] [--energySAPTDataFieldLabels <Type,Label,...,...>]
776 [--energyUnits <text>] [--infileParams <Name,Value,...>] [--methodName <text>]
777 [--mp <yes or no>] [--mpParams <Name, Value,...>] [ --outfileParams <Name,Value,...> ]
778 [--overwrite] [--precision <number> ] [--psi4OptionsParams <Name,Value,...>]
779 [--psi4RunParams <Name,Value,...>] [--quiet <yes or no>]
780 [--reference <text>] [-w <dir>] -i <infile> -o <outfile>
781 Psi4CalculateInteractionEnergy.py -h | --help | -e | --examples
782
783 Description:
784 Calculate interaction energy for molecules using a specified method name and
785 basis set. The molecules must contain exactly two fragments or disconnected
786 components for Symmetry Adapted Perturbation Theory (SAPT) [ Ref 154-155 ]
787 and Spin-Network-Scaled MP2 (SNS-MP2) [ Ref 156] calculations and more than
788 one fragment for all other calculations. An arbitrary number of fragments may be
789 present in a molecule for Basis Set Superposition Energy (BSSE) correction
790 calculations.
791
792 The interaction energy is calculated at SAPT0 / jun-cc-pVDZ level of theory by
793 default. The SAPT0 calculations are relatively fast but less accurate. You may
794 want to consider calculating interaction energy at WB97X-D3 / aug-cc-pVTZ,
795 B3LYP-D3 / aug-cc-pVTZ, or higher level of theory [ Ref 157 ] to improve the
796 accuracy of your results. The WB97X-D3 and B3LYP-D3 calculations rely on
797 the presence of DFTD3 and gCP Psi4 plugins in your environment.
798
799 The molecules must have 3D coordinates in input file. The molecular geometry
800 is not optimized before the calculation. In addition, hydrogens must be present
801 for all molecules in input file. The 3D coordinates are not modified during the
802 calculation.
803
804 A Psi4 XYZ format geometry string is automatically generated for each molecule
805 in input file. It contains atom symbols and 3D coordinates for each atom in a
806 molecule. In addition, the formal charge and spin multiplicity are present in the
807 the geometry string. These values are either retrieved from molecule properties
808 named 'FormalCharge' and 'SpinMultiplicty' or dynamically calculated for a
809 molecule. A double dash separates each fragment or component in a molecule.
810 The same formal charge and multiplicity values are assigned to each fragment
811 in a molecule.
812
813 The supported input file formats are: Mol (.mol), SD (.sdf, .sd)
814
815 The supported output file formats are: SD (.sdf, .sd)
816
817 Options:
818 -b, --basisSet <text> [default: auto]
819 Basis set to use for interaction energy calculation. Default: jun-cc-pVDZ for
820 SAPT calculations; None for SNS-MP2 calculations to use its default basis
821 set; otherwise, it must be explicitly specified using this option. The specified
822 value must be a valid Psi4 basis set. No validation is performed. You may set
823 an empty string as a value for the basis set.
824
825 The following list shows a representative sample of basis sets available
826 in Psi4:
827
828 STO-3G, 6-31G, 6-31+G, 6-31++G, 6-31G*, 6-31+G*, 6-31++G*,
829 6-31G**, 6-31+G**, 6-31++G**, 6-311G, 6-311+G, 6-311++G,
830 6-311G*, 6-311+G*, 6-311++G*, 6-311G**, 6-311+G**, 6-311++G**,
831 cc-pVDZ, cc-pCVDZ, aug-cc-pVDZ, cc-pVDZ-DK, cc-pCVDZ-DK, def2-SVP,
832 def2-SVPD, def2-TZVP, def2-TZVPD, def2-TZVPP, def2-TZVPPD
833
834 --bsseType <CP, noCP, VMFC, or None> [default: auto]
835 Type of Basis Set Superposition Energy (BSSE) correction to apply during the
836 calculation of interaction energy. Possible values:
837
838 CP: Counterpoise corrected interaction energy
839 noCP: Supramolecular interaction energy without any CP correction
840 VMFC: Valiron-Mayer Function Counterpoise correction
841 None: The Psi4 option 'bsse_type' is not passed to the energy
842 function during the calculation of interaction energy
843
844 Default values:
845
846 None: SAPT and SNS-MP2 calculations. An explicit bsse_type option is not
847 valid for these calculations.
848 HF3c: noCP to use built-in correction
849 CP: All other calculations
850
851 --energyDataFieldLabel <text> [default: auto]
852 Energy data field label for writing interaction energy values. Default:
853 Psi4_SAPT_Interaction_Energy (<Units>) for SAPT calculation;
854 Psi4_SNS-MP2_Interaction_Energy (<Units>) for SNS-MP2 calculation;
855 otherwise, Psi4_Interaction_Energy (<Units>)
856 --energySAPTDataFieldLabels <Type,Label,...,...> [default: auto]
857 A comma delimted interaction energy type and data field label value pairs
858 for writing individual components of total SAPT interaction energy.
859
860 The supported SAPT energy types along with their default data field label
861 values are shown below:
862
863 ElectrostaticEnergy, Psi4_SAPT_Electrostatic_Energy (<Units>),
864 ExchangeEnergy, Psi4_SAPT_Exchange_Energy (<Units>),
865 InductionEnergy, Psi4_SAPT_Induction_Energy (<Units>),
866 DispersionEnergy, Psi4_SAPT_Dispersion_Energy (<Units>)
867
868 --energyUnits <text> [default: kcal/mol]
869 Energy units. Possible values: Hartrees, kcal/mol, kJ/mol, or eV.
870 -e, --examples
871 Print examples.
872 -h, --help
873 Print this help message.
874 -i, --infile <infile>
875 Input file name.
876 --infileParams <Name,Value,...> [default: auto]
877 A comma delimited list of parameter name and value pairs for reading
878 molecules from files. The supported parameter names for different file
879 formats, along with their default values, are shown below:
880
881 SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
882
883 -m, --methodName <text> [default: auto]
884 Method to use for interaction energy calculation. Default: SAPT0. The
885 specified value must be a valid Psi4 method name. No validation is
886 performed.
887
888 The following list shows a representative sample of methods available
889 in Psi4:
890
891 SAPT0, SAPT2, SAPT2+, SAPT2+(CCD), SAPT2+DMP2, SAPT2+(CCD)DMP2
892 SAPT2+(3), SAPT2+(3)(CCD), SAPT2+DMP2, SAPT2+(CCD)DMP2,
893 SAPT2+3, SAPT2+3(CCD), SAPT2+(3)DMP2, SAPT2+3(CCD)DMP2, SNS-MP2,
894 B1LYP, B2PLYP, B2PLYP-D3BJ, B2PLYP-D3MBJ, B3LYP, B3LYP-D3BJ,
895 B3LYP-D3MBJ, CAM-B3LYP, CAM-B3LYP-D3BJ, HF, HF-D3BJ, HF3c, M05,
896 M06, M06-2x, M06-HF, M06-L, MN12-L, MN15, MN15-D3BJ,PBE, PBE0,
897 PBEH3c, PW6B95, PW6B95-D3BJ, WB97, WB97X, WB97X-D, WB97X-D3BJ
898
899 --mp <yes or no> [default: no]
900 Use multiprocessing.
901
902 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
903 function employing lazy RDKit data iterable. This allows processing of
904 arbitrary large data sets without any additional requirements memory.
905
906 All input data may be optionally loaded into memory by mp.Pool.map()
907 before starting worker processes in a process pool by setting the value
908 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
909
910 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
911 data mode may adversely impact the performance. The '--mpParams' section
912 provides additional information to tune the value of 'chunkSize'.
913 --mpParams <Name,Value,...> [default: auto]
914 A comma delimited list of parameter name and value pairs to configure
915 multiprocessing.
916
917 The supported parameter names along with their default and possible
918 values are shown below:
919
920 chunkSize, auto
921 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
922 numProcesses, auto [ Default: mp.cpu_count() ]
923
924 These parameters are used by the following functions to configure and
925 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
926 mp.Pool.imap().
927
928 The chunkSize determines chunks of input data passed to each worker
929 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
930 The default value of chunkSize is dependent on the value of 'inputDataMode'.
931
932 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
933 automatically converts RDKit data iterable into a list, loads all data into
934 memory, and calculates the default chunkSize using the following method
935 as shown in its code:
936
937 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
938 if extra: chunkSize += 1
939
940 For example, the default chunkSize will be 7 for a pool of 4 worker processes
941 and 100 data items.
942
943 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
944 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
945 data into memory. Consequently, the size of input data is not known a priori.
946 It's not possible to estimate an optimal value for the chunkSize. The default
947 chunkSize is set to 1.
948
949 The default value for the chunkSize during 'Lazy' data mode may adversely
950 impact the performance due to the overhead associated with exchanging
951 small chunks of data. It is generally a good idea to explicitly set chunkSize to
952 a larger value during 'Lazy' input data mode, based on the size of your input
953 data and number of processes in the process pool.
954
955 The mp.Pool.map() function waits for all worker processes to process all
956 the data and return the results. The mp.Pool.imap() function, however,
957 returns the the results obtained from worker processes as soon as the
958 results become available for specified chunks of data.
959
960 The order of data in the results returned by both mp.Pool.map() and
961 mp.Pool.imap() functions always corresponds to the input data.
962 -o, --outfile <outfile>
963 Output file name.
964 --outfileParams <Name,Value,...> [default: auto]
965 A comma delimited list of parameter name and value pairs for writing
966 molecules to files. The supported parameter names for different file
967 formats, along with their default values, are shown below:
968
969 SD: kekulize,yes,forceV3000,no
970
971 --overwrite
972 Overwrite existing files.
973 --precision <number> [default: 6]
974 Floating point precision for writing energy values.
975 --psi4OptionsParams <Name,Value,...> [default: none]
976 A comma delimited list of Psi4 option name and value pairs for setting
977 global and module options. The names are 'option_name' for global options
978 and 'module_name__option_name' for options local to a module. The
979 specified option names must be valid Psi4 names. No validation is
980 performed.
981
982 The specified option name and value pairs are processed and passed to
983 psi4.set_options() as a dictionary. The supported value types are float,
984 integer, boolean, or string. The float value string is converted into a float.
985 The valid values for a boolean string are yes, no, true, false, on, or off.
986 --psi4RunParams <Name,Value,...> [default: auto]
987 A comma delimited list of parameter name and value pairs for configuring
988 Psi4 jobs.
989
990 The supported parameter names along with their default and possible
991 values are shown below:
992
993 MemoryInGB, 1
994 NumThreads, 1
995 OutputFile, auto [ Possible values: stdout, quiet, or FileName ]
996 ScratchDir, auto [ Possivle values: DirName]
997 RemoveOutputFile, yes [ Possible values: yes, no, true, or false]
998
999 These parameters control the runtime behavior of Psi4.
1000
1001 The default file name for 'OutputFile' is <InFileRoot>_Psi4.out. The PID
1002 is appended to output file name during multiprocessing as shown:
1003 <InFileRoot>_Psi4_<PIDNum>.out. The 'stdout' value for 'OutputType'
1004 sends Psi4 output to stdout. The 'quiet' or 'devnull' value suppresses
1005 all Psi4 output.
1006
1007 The default 'Yes' value of 'RemoveOutputFile' option forces the removal
1008 of any existing Psi4 before creating new files to append output from
1009 multiple Psi4 runs.
1010
1011 The option 'ScratchDir' is a directory path to the location of scratch
1012 files. The default value corresponds to Psi4 default. It may be used to
1013 override the deafult path.
1014 -q, --quiet <yes or no> [default: no]
1015 Use quiet mode. The warning and information messages will not be printed.
1016 -r, --reference <text> [default: auto]
1017 Reference wave function to use for interaction energy calculation. Default: RHF
1018 or UHF. The default values are Restricted Hartree-Fock (RHF) for closed-shell
1019 molecules with all electrons paired and Unrestricted Hartree-Fock (UHF) for
1020 open-shell molecules with unpaired electrons.
1021
1022 The specified value must be a valid Psi4 reference wave function. No validation
1023 is performed. For example: ROHF, CUHF, RKS, etc.
1024
1025 The spin multiplicity determines the default value of reference wave function
1026 for input molecules. It is calculated from number of free radical electrons using
1027 Hund's rule of maximum multiplicity defined as 2S + 1 where S is the total
1028 electron spin. The total spin is 1/2 the number of free radical electrons in a
1029 molecule. The value of 'SpinMultiplicity' molecule property takes precedence
1030 over the calculated value of spin multiplicity.
1031
1032 The 'SpinMultiplicity' molecule property may contain multiples values
1033 corresponding to the number of fragments in a molecule. This property must
1034 not be set for automatic determination of the reference wave functionn.
1035 -w, --workingdir <dir>
1036 Location of working directory which defaults to the current directory.
1037
1038 Examples:
1039 To calculate interaction energy using SAPT0/aug-cc-pVDZ for molecules in a
1040 SD file, use RHF and UHF for closed-shell and open-shell molecules, and write
1041 a new SD file, type:
1042
1043 % Psi4CalculateInteractionEnergy.py -i Psi4SampleDimers3D.sdf
1044 -o Psi4SampleDimers3DOut.sdf
1045
1046 To run the first example for freezing core electrons and setting SCF type to DF
1047 and write out a new SD file, type:
1048
1049 % Psi4CalculateInteractionEnergy.py --psi4OptionsParams "scf_type, df,
1050 freeze_core, true" -i Psi4SampleDimers3D.sdf -o
1051 Psi4SampleDimers3DOut.sdf
1052
1053 To calculate interaction energy using SNS2-MP methodology for molecules
1054 in a SD containing 3D structures and write a new SD file, type:
1055
1056 % Psi4CalculateInteractionEnergy.py -m "sns-mp2"
1057 -i Psi4SampleDimers3D.sdf -o Psi4SampleDimers3DOut.sdf
1058
1059 To calculate interaction energy at WB97X-D3/aug-cc-pVTZ level of theory,
1060 along with explicit Psi4 run time paramaters, for molecules in a SD containing
1061 3D structures and write a new SD file, type:
1062
1063 % Psi4CalculateInteractionEnergy.py -m WB97X-D3 -b aug-cc-pVTZ
1064 --bsseType CP --psi4RunParams "NumThreads,4,MemoryInGB,6"
1065 -i Psi4SampleDimers3D.sdf -o Psi4SampleDimers3DOut.sdf
1066
1067 To calculate interaction energy at B3LYP-D3/aug-cc-pVTZ level of theory using
1068 default Psi4 run time paramaters for molecules in a SD containing 3D structures
1069 and write a new SD file, type:
1070
1071 % Psi4CalculateInteractionEnergy.py -m B3LYP-D3 -b aug-cc-pVTZ
1072 --bsseType CP -i Psi4SampleDimers3D.sdf -o Psi4SampleDimers3DOut.sdf
1073
1074 To calculate interaction energy at B3LYP-D3/aug-cc-pVTZ level of theory, along
1075 with specifying grid resolution using Psi4 options and explicit Psi4 run time
1076 paramaters, for molecules in a SD containing 3D structures
1077 and write a new SD file, type:
1078
1079 % Psi4CalculateInteractionEnergy.py -m B3LYP-D3 -b aug-cc-pVTZ
1080 --bsseType CP --psi4OptionsParams "dft_spherical_points, 302,
1081 dft_radial_points, 75" --psi4RunParams "NumThreads,4,MemoryInGB,6"
1082 -i Psi4SampleDimers3D.sdf -o Psi4SampleDimers3DOut.sdf
1083
1084 To calculate interaction energy at HF3c level of theory using built-in basis set
1085 for molecules in a SD containing 3D structures and write a new SD file, type:
1086
1087 % Psi4CalculateInteractionEnergy.py -m HF3c -b "" --bsseType noCP
1088 -i Psi4SampleDimers3D.sdf -o Psi4SampleDimers3DOut.sdf
1089
1090 To calculate interaction energy at CCSD(T)/aug-cc-pVDZ level of theory using
1091 default Psi4 run time paramaters for molecules in a SD containing 3D structures
1092 and write a new SD file, type:
1093
1094 % Psi4CalculateInteractionEnergy.py -m "ccsd(t)" -b "aug-cc-pvdz"
1095 -i Psi4SampleDimers3D.sdf -o Psi4SampleDimers3DOut.sdf
1096
1097 To run the first example in multiprocessing mode on all available CPUs
1098 without loading all data into memory and write out a SD file, type:
1099
1100 % Psi4CalculateInteractionEnergy.py --mp yes -i Psi4SampleDimers3D.sdf
1101 -o Psi4SampleDimers3DOut.sdf
1102
1103 To run the first example in multiprocessing mode on all available CPUs
1104 by loading all data into memory and write out a SD file, type:
1105
1106 % Psi4CalculateInteractionEnergy.py --mp yes --mpParams "inputDataMode,
1107 InMemory" -i Psi4SampleDimers3D.sdf -o Psi4SampleDimers3DOut.sdf
1108
1109 To run the first example in multiprocessing mode on all available CPUs
1110 without loading all data into memory along with multiple threads for each
1111 Psi4 run and write out a SD file, type:
1112
1113 % Psi4CalculateInteractionEnergy.py --mp yes --psi4RunParams "NumThreads,2"
1114 -i Psi4SampleDimers3D.sdf -o Psi4SampleDimers3DOut.sdf
1115
1116 Author:
1117 Manish Sud(msud@san.rr.com)
1118
1119 See also:
1120 Psi4CalculateEnergy.py, Psi4CalculatePartialCharges.py, Psi4PerformMinimization.py,
1121 Psi4GenerateConformers.py
1122
1123 Copyright:
1124 Copyright (C) 2026 Manish Sud. All rights reserved.
1125
1126 The functionality available in this script is implemented using Psi4, an
1127 open source quantum chemistry software package, and RDKit, an open
1128 source toolkit for cheminformatics developed by Greg Landrum.
1129
1130 This file is part of MayaChemTools.
1131
1132 MayaChemTools is free software; you can redistribute it and/or modify it under
1133 the terms of the GNU Lesser General Public License as published by the Free
1134 Software Foundation; either version 3 of the License, or (at your option) any
1135 later version.
1136
1137 """
1138
1139 if __name__ == "__main__":
1140 main()
