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