MayaChemTools

    1 #!/bin/env python
    2 #
    3 # File: Psi4CalculatePartialCharges.py
    4 # Author: Manish Sud <msud@san.rr.com>
    5 #
    6 # Copyright (C) 2022 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 shutil
   38 import multiprocessing as mp
   39 
   40 # Psi4 imports...
   41 if (hasattr(shutil, 'which') and shutil.which("psi4") is None):
   42     sys.stderr.write("\nWarning: Failed to find 'psi4' in your PATH indicating potential issues with your\n")
   43     sys.stderr.write("Psi4 environment. The 'import psi4' directive in the global scope of the script\n")
   44     sys.stderr.write("interferes with the multiprocessing functionality. It is imported later in the\n")
   45     sys.stderr.write("local scope during the execution of the script and may fail. Check/update your\n")
   46     sys.stderr.write("Psi4 environment and try again.\n\n")
   47 
   48 # RDKit imports...
   49 try:
   50     from rdkit import rdBase
   51     from rdkit import Chem
   52     from rdkit.Chem import AllChem
   53 except ImportError as ErrMsg:
   54     sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
   55     sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
   56     sys.exit(1)
   57 
   58 # MayaChemTools imports...
   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 def main():
   74     """Start execution of the script"""
   75     
   76     MiscUtil.PrintInfo("\n%s (Psi4: Imported later; RDKit v%s; MayaChemTools v%s; %s): Starting...\n" % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime()))
   77     
   78     (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
   79     
   80     # Retrieve command line arguments and options...
   81     RetrieveOptions()
   82     
   83     # Process and validate command line arguments and options...
   84     ProcessOptions()
   85     
   86     # Perform actions required by the script...
   87     CalculatePartialCharges()
   88     
   89     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
   90     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
   91 
   92 def CalculatePartialCharges():
   93     """Calculate partial atomic charges."""
   94     
   95     CheckSupportForRESPCalculations()
   96     
   97     # Setup a molecule reader...
   98     MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
   99     Mols  = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
  100     
  101     # Set up a molecule writer...
  102     Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"])
  103     if Writer is None:
  104         MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"])
  105     MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["Outfile"])
  106 
  107     MolCount, ValidMolCount, CalcFailedCount = ProcessMolecules(Mols, Writer)
  108 
  109     if Writer is not None:
  110         Writer.close()
  111     
  112     MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
  113     MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
  114     MiscUtil.PrintInfo("Number of molecules failed during calculation of partial charges: %d" % CalcFailedCount)
  115     MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + CalcFailedCount))
  116 
  117 def ProcessMolecules(Mols, Writer):
  118     """Process molecules and calculate partial charges."""
  119     
  120     if OptionsInfo["MPMode"]:
  121         return ProcessMoleculesUsingMultipleProcesses(Mols, Writer)
  122     else:
  123         return ProcessMoleculesUsingSingleProcess(Mols, Writer)
  124 
  125 def ProcessMoleculesUsingSingleProcess(Mols, Writer):
  126     """Process molecules and calculate partial charges using a single process. """
  127     
  128     # Intialize Psi4...
  129     MiscUtil.PrintInfo("\nInitializing Psi4...")
  130     Psi4Handle = Psi4Util.InitializePsi4(Psi4RunParams = OptionsInfo["Psi4RunParams"], Psi4OptionsParams = OptionsInfo["Psi4OptionsParams"], PrintVersion = True, PrintHeader = True)
  131     OptionsInfo["psi4"] = Psi4Handle
  132 
  133     # Initialize RESP...
  134     OptionsInfo["resp"] = SetupRESP()
  135     
  136     MiscUtil.PrintInfo("\nCalculating partial atomic charges...")
  137     
  138     (MolCount, ValidMolCount, CalcFailedCount) = [0] * 3
  139     for Mol in Mols:
  140         MolCount += 1
  141         
  142         if not CheckAndValidateMolecule(Mol, MolCount):
  143             continue
  144         
  145         # Setup a Psi4 molecule...
  146         Psi4Mol = SetupPsi4Mol(Psi4Handle, Mol, MolCount)
  147         if Psi4Mol is None:
  148             continue
  149         
  150         ValidMolCount += 1
  151         
  152         # Retrieve charges...
  153         CalcStatus, PartialCharges = CalculateMolPartialCharges(Psi4Handle, Psi4Mol, Mol, MolCount)
  154 
  155         if not CalcStatus:
  156             if not OptionsInfo["QuietMode"]:
  157                 MolName = RDKitUtil.GetMolName(Mol, MolCount)
  158                 MiscUtil.PrintWarning("Failed to calculate partial atomic charges for molecule %s" % MolName)
  159             
  160             CalcFailedCount += 1
  161             continue
  162         
  163         WriteMolPartialCharges(Writer, Mol, PartialCharges)
  164     
  165     return (MolCount, ValidMolCount, CalcFailedCount)
  166 
  167 def ProcessMoleculesUsingMultipleProcesses(Mols, Writer):
  168     """Process molecules and calculate partial charges using a multiprocessing. """
  169 
  170     MiscUtil.PrintInfo("\nCalculating partial atomic charges using multiprocessing...")
  171     
  172     MPParams = OptionsInfo["MPParams"]
  173     
  174     # Setup data for initializing a worker process...
  175     InitializeWorkerProcessArgs = (MiscUtil.ObjectToBase64EncodedString(Options), MiscUtil.ObjectToBase64EncodedString(OptionsInfo))
  176     
  177     # Setup a encoded mols data iterable for a worker process...
  178     WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
  179     
  180     # Setup process pool along with data initialization for each process...
  181     if not OptionsInfo["QuietMode"]:
  182         MiscUtil.PrintInfo("\nConfiguring multiprocessing using %s method..." % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()"))
  183         MiscUtil.PrintInfo("NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" % (MPParams["NumProcesses"], MPParams["InputDataMode"], ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"])))
  184     
  185     ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
  186     
  187     # Start processing...
  188     if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
  189         Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
  190     elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
  191         Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
  192     else:
  193         MiscUtil.PrintError("The value, %s, specified for \"--inputDataMode\" is not supported." % (MPParams["InputDataMode"]))
  194 
  195     # Print out Psi4 version in the main process...
  196     MiscUtil.PrintInfo("\nInitializing Psi4...\n")
  197     Psi4Handle  = Psi4Util.InitializePsi4(PrintVersion = True, PrintHeader = False)
  198     OptionsInfo["psi4"] = Psi4Handle
  199     
  200     (MolCount, ValidMolCount, CalcFailedCount) = [0] * 3
  201     for Result in Results:
  202         MolCount += 1
  203         MolIndex, EncodedMol, CalcStatus, PartialCharges = Result
  204         
  205         if EncodedMol is None:
  206             continue
  207         
  208         ValidMolCount += 1
  209 
  210         Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
  211         
  212         if not CalcStatus:
  213             if not OptionsInfo["QuietMode"]:
  214                 MolName = RDKitUtil.GetMolName(Mol, MolCount)
  215                 MiscUtil.PrintWarning("Failed to calculate partial atomic charges for molecule %s" % MolName)
  216             
  217             CalcFailedCount += 1
  218             continue
  219         
  220         WriteMolPartialCharges(Writer, Mol, PartialCharges)
  221     
  222     return (MolCount, ValidMolCount, CalcFailedCount)
  223 
  224 def InitializeWorkerProcess(*EncodedArgs):
  225     """Initialize data for a worker process."""
  226     
  227     global Options, OptionsInfo
  228     
  229     if not OptionsInfo["QuietMode"]:
  230         MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
  231     
  232     # Decode Options and OptionInfo...
  233     Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
  234     OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
  235 
  236     # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
  237     # output files for each process...
  238     OptionsInfo["Psi4Initialized"]  = False
  239     
  240 def InitializePsi4ForWorkerProcess():
  241     """Initialize Psi4 for a worker process."""
  242     
  243     if OptionsInfo["Psi4Initialized"]:
  244         return
  245 
  246     OptionsInfo["Psi4Initialized"] = True
  247     
  248     # Update output file...
  249     OptionsInfo["Psi4RunParams"]["OutputFile"] = Psi4Util.UpdatePsi4OutputFileUsingPID(OptionsInfo["Psi4RunParams"]["OutputFile"], os.getpid())
  250     
  251     # Intialize Psi4...
  252     OptionsInfo["psi4"] = Psi4Util.InitializePsi4(Psi4RunParams = OptionsInfo["Psi4RunParams"], Psi4OptionsParams = OptionsInfo["Psi4OptionsParams"], PrintVersion = False, PrintHeader = True)
  253     
  254 def WorkerProcess(EncodedMolInfo):
  255     """Process data for a worker process."""
  256     
  257     if not OptionsInfo["Psi4Initialized"]:
  258         InitializePsi4ForWorkerProcess()
  259     
  260     MolIndex, EncodedMol = EncodedMolInfo
  261     
  262     CalcStatus = False
  263     PartialCharges = None
  264     
  265     if EncodedMol is None:
  266         return [MolIndex, None, CalcStatus, PartialCharges]
  267     
  268     Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
  269     MolCount = MolIndex + 1
  270     
  271     if not CheckAndValidateMolecule(Mol, MolCount):
  272         return [MolIndex, None, CalcStatus, PartialCharges]
  273     
  274     # Setup a Psi4 molecule...
  275     Psi4Mol = SetupPsi4Mol(OptionsInfo["psi4"], Mol, MolCount)
  276     if Psi4Mol is None:
  277         return [MolIndex, None, CalcStatus, PartialCharges]
  278     
  279     CalcStatus, PartialCharges = CalculateMolPartialCharges(OptionsInfo["psi4"], Psi4Mol, Mol, MolCount)
  280     
  281     return [MolIndex, RDKitUtil.MolToBase64EncodedMolString(Mol, PropertyPickleFlags = Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps), CalcStatus, PartialCharges]
  282 
  283 def CalculateMolPartialCharges(Psi4Handle, Psi4Mol, Mol, MolNum):
  284     """Calculate partial atomic charges for a molecule."""
  285 
  286     if OptionsInfo["RESPChargesTypeMode"]:
  287         return CalculateRespMolPartialCharges(Psi4Handle, Psi4Mol, Mol, MolNum)
  288     else:
  289         return CalculateNonRespMolPartialCharges(Psi4Handle, Psi4Mol, Mol, MolNum)
  290     
  291 def CalculateNonRespMolPartialCharges(Psi4Handle, Psi4Mol, Mol, MolNum):
  292     """Calculate non-RESP partial atomic charges for a molecule."""
  293     
  294     Status = False
  295     PartialCharges = []
  296     
  297     #  Setup reference wave function...
  298     Reference = SetupReferenceWavefunction(Mol)
  299     Psi4Handle.set_options({'Reference': Reference})
  300     
  301     # Setup method name and basis set...
  302     MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
  303 
  304     # Calculate single point energy to setup a wavefunction...
  305     Status, Energy, WaveFunction = Psi4Util.CalculateSinglePointEnergy(Psi4Handle, Psi4Mol, MethodName, BasisSet, ReturnWaveFunction = True, Quiet = OptionsInfo["QuietMode"])
  306 
  307     if not Status:
  308         PerformPsi4Cleanup(Psi4Handle)
  309         return (False, PartialCharges)
  310     
  311     # Calculate atomic point charges...
  312     try:
  313         Psi4Handle.oeprop(WaveFunction, OptionsInfo["ChargesPropType"])
  314     except Exception as ErrMsg:
  315         if not OptionsInfo["QuietMode"]:
  316             MiscUtil.PrintWarning("Failed to calculate partial atomic charges:\n%s" % ErrMsg)
  317         PerformPsi4Cleanup(Psi4Handle)
  318         return (False, PartialCharges)
  319     
  320     AtomicPointCharges = WaveFunction.atomic_point_charges().np.tolist()
  321     
  322     # Clean up...
  323     PerformPsi4Cleanup(Psi4Handle)
  324     
  325     # Format charges...
  326     PartialCharges = ["%.*f" % (OptionsInfo["Precision"], float(Value)) for Value in AtomicPointCharges]
  327 
  328     return (Status, PartialCharges)
  329 
  330 def CalculateRespMolPartialCharges(Psi4Handle, Psi4Mol, Mol, MolNum):
  331     """Calculate RESP partial atomic charges for a molecule."""
  332 
  333     PartialCharges = []
  334 
  335     RESPHandle = OptionsInfo["resp"]
  336     RESPOptions = OptionsInfo["ChargesRespOtions"]
  337     
  338     # Setup method name and basis set...
  339     MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
  340     RESPOptions["METHOD_ESP"] = MethodName
  341     RESPOptions["BASIS_ESP"] = BasisSet
  342 
  343     # Calculate RESP charges...
  344     try:
  345         RESPCalcResults =  RESPHandle.resp([Psi4Mol], RESPOptions)
  346     except Exception as ErrMsg:
  347         if not OptionsInfo["QuietMode"]:
  348             MiscUtil.PrintWarning("Failed to calculate RESP partial atomic charges:\n%s" % ErrMsg)
  349         RemoveRESPGridFiles(MolNum)
  350         return (False, PartialCharges)
  351     
  352     ESPCharges = RESPCalcResults[0]
  353     RESPCharges = RESPCalcResults[1]
  354     
  355     PerformRESPCleanup(MolNum)
  356     
  357     # Format charges...
  358     PartialCharges = ["%.*f" % (OptionsInfo["Precision"], float(Value)) for Value in RESPCharges]
  359     
  360     return (True, PartialCharges)
  361 
  362 def PerformRESPCleanup(MolNum):
  363     """Peform RESP cleanup."""
  364 
  365     if OptionsInfo["ChargesRespParams"]["RemoveGridFiles"]:
  366         RemoveRESPGridFiles(MolNum)
  367     else:
  368         RenameRESPGridFiles(MolNum)
  369 
  370 def RemoveRESPGridFiles(MolNum):
  371     """Remove RESP grid files."""
  372 
  373     try:
  374         for GridFile in ["1_default_grid.dat", "1_default_grid_esp.dat", "results.out"]:
  375             if os.path.isfile(GridFile):
  376                 os.remove(GridFile)
  377     except Exception as ErrMsg:
  378         if not OptionsInfo["QuietMode"]:
  379             MiscUtil.PrintWarning("Failed to remove RESP results/grid files: %s\n" % ErrMsg)
  380 
  381 def RenameRESPGridFiles(MolNum):
  382     """Rename RESP grid files. """
  383 
  384     try:
  385         MolPrefix = "Mol%s" % MolNum
  386         GridFiles = ["1_default_grid.dat", "1_default_grid_esp.dat", "results.out"]
  387         NewGridFiles = ["%s_grid.dat" % MolPrefix, "%s_grid_esp.dat" % MolPrefix, "%s_resp_results.out" % MolPrefix]
  388         for GridFile, NewGridFile in zip(GridFiles, NewGridFiles):
  389             if os.path.isfile(GridFile):
  390                 shutil.move(GridFile, NewGridFile)
  391     except Exception as ErrMsg:
  392         if not OptionsInfo["QuietMode"]:
  393             MiscUtil.PrintWarning("Failed to move RESP results/grid files: %s\n" % ErrMsg)
  394 
  395 def WriteMolPartialCharges(Writer, Mol, PartialCharges):
  396     """Write out partial atomic charges for a molecule."""
  397 
  398     if PartialCharges is None:
  399         return
  400     
  401     if OptionsInfo["AtomAliasesFormatMode"]:
  402         for Atom, PartialCharge in zip(Mol.GetAtoms(), PartialCharges):
  403             Atom.SetProp('molFileAlias', PartialCharge)
  404     else:
  405         ChargesValues = "\n".join(PartialCharges)
  406         Mol.SetProp(OptionsInfo["DataFieldLabel"], ChargesValues)
  407     
  408     Writer.write(Mol)
  409     
  410 def SetupPsi4Mol(Psi4Handle, Mol, MolCount = None):
  411     """Setup a Psi4 molecule object."""
  412     
  413     MolGeometry = RDKitUtil.GetPsi4XYZFormatString(Mol, NoCom = True, NoReorient = True)
  414     
  415     try:
  416         Psi4Mol = Psi4Handle.geometry(MolGeometry)
  417     except Exception as ErrMsg:
  418         Psi4Mol = None
  419         if not OptionsInfo["QuietMode"]:
  420             MiscUtil.PrintWarning("Failed to create Psi4 molecule from geometry string: %s\n" % ErrMsg)
  421             MolName = RDKitUtil.GetMolName(Mol, MolCount)
  422             MiscUtil.PrintWarning("Ignoring molecule: %s" % MolName)
  423     
  424     return Psi4Mol
  425 
  426 def PerformPsi4Cleanup(Psi4Handle):
  427     """Perform clean up. """
  428     
  429     # Clean up after Psi4 run ...
  430     Psi4Handle.core.clean()
  431     
  432     # Clean up any leftover scratch files...
  433     if OptionsInfo["MPMode"]:
  434         Psi4Util.RemoveScratchFiles(Psi4Handle, OptionsInfo["Psi4RunParams"]["OutputFile"])
  435 
  436 def SetupRESP():
  437     """Load resp and return its handle."""
  438     
  439     if not OptionsInfo["RESPChargesTypeMode"]:
  440         return None
  441         
  442     try:
  443         import resp
  444     except ImportError as ErrMsg:
  445         sys.stderr.write("\nFailed to import Psi4 module/package resp: %s\n" % ErrMsg)
  446         sys.stderr.write("Check/update your Psi4 environment and try again.\n\n")
  447         sys.exit(1)
  448 
  449     return resp
  450     
  451 def CheckAndValidateMolecule(Mol, MolCount = None):
  452     """Validate molecule for Psi4 calculations."""
  453     
  454     if Mol is None:
  455         if not OptionsInfo["QuietMode"]:
  456             MiscUtil.PrintInfo("\nProcessing molecule number %s..." % MolCount)
  457         return False
  458     
  459     MolName = RDKitUtil.GetMolName(Mol, MolCount)
  460     if not OptionsInfo["QuietMode"]:
  461         MiscUtil.PrintInfo("\nProcessing molecule %s..." % MolName)
  462     
  463     if RDKitUtil.IsMolEmpty(Mol):
  464         if not OptionsInfo["QuietMode"]:
  465             MiscUtil.PrintWarning("Ignoring empty molecule: %s\n" % MolName)
  466         return False
  467     
  468     if not RDKitUtil.ValidateElementSymbols(RDKitUtil.GetAtomSymbols(Mol)):
  469         if not OptionsInfo["QuietMode"]:
  470             MiscUtil.PrintWarning("Ignoring molecule containing invalid element symbols: %s\n" % MolName)
  471         return False
  472     
  473     # Check for 3D flag...
  474     if not Mol.GetConformer().Is3D():
  475         if not OptionsInfo["QuietMode"]:
  476             MiscUtil.PrintWarning("3D tag is not set for molecule: %s\n" % MolName)
  477     
  478     # Check for missing hydrogens...
  479     if RDKitUtil.AreHydrogensMissingInMolecule(Mol):
  480         if not OptionsInfo["QuietMode"]:
  481             MiscUtil.PrintWarning("Missing hydrogens in molecule: %s\n" % MolName)
  482     
  483     return True
  484 
  485 def SetupMethodNameAndBasisSet(Mol):
  486     """Setup method name and basis set. """
  487     
  488     MethodName = OptionsInfo["MethodName"]
  489     if OptionsInfo["MethodNameAuto"]:
  490         MethodName = "B3LYP"
  491     
  492     BasisSet = OptionsInfo["BasisSet"]
  493     if OptionsInfo["BasisSetAuto"]:
  494         BasisSet = "6-31+G**" if RDKitUtil.IsAtomSymbolPresentInMol(Mol, "S") else "6-31G**"
  495     
  496     return (MethodName, BasisSet)
  497 
  498 def SetupReferenceWavefunction(Mol):
  499     """Setup method name and basis set. """
  500     
  501     Reference = OptionsInfo["Reference"]
  502     if OptionsInfo["ReferenceAuto"]:
  503         Reference = 'UHF' if (RDKitUtil.GetSpinMultiplicity(Mol) > 1) else 'RHF'
  504     
  505     return Reference
  506 
  507 def ProcessOptionChargesRespParameters():
  508     """Process charges RSEP parameters option."""
  509 
  510     ParamsOptionName = "--chargesRespParams"
  511     ParamsOptionValue = Options["--chargesRespParams"]
  512 
  513     VDWScaleFactors = [1.4, 1.6, 1.8, 2.0]
  514     VDWRadii = {'H': 1.20, 'HE': 1.20, 'LI': 1.37, 'BE': 1.45, 'B': 1.45, 'C': 1.50, 'N': 1.50, 'O': 1.40, 'F': 1.35, 'NE': 1.30, 'NA': 1.57, 'MG': 1.36, 'AL': 1.24, 'SI': 1.17, 'P': 1.80, 'S': 1.75, 'CL': 1.70}
  515     
  516     ParamsInfo = {"MaxIter": 25, "RestrainHydrogens": False, "RemoveGridFiles": True, "RespA": 0.0005, "RespB": 0.1, "Tolerance": 1e-5, "VDWRadii": VDWRadii, "VDWScaleFactors": VDWScaleFactors, "VDWPointDensity": 1.0}
  517 
  518     if re.match("^auto$", ParamsOptionValue, re.I):
  519         SetupChargesRespOptions(ParamsInfo)
  520         return
  521     
  522     # Setup a canonical paramater names...
  523     ValidParamNames = []
  524     CanonicalParamNamesMap = {}
  525     for ParamName in sorted(ParamsInfo):
  526         ValidParamNames.append(ParamName)
  527         CanonicalParamNamesMap[ParamName.lower()] = ParamName
  528     
  529     ParamsOptionValue = ParamsOptionValue.strip()
  530     if not ParamsOptionValue:
  531         MiscUtil.PrintError("No valid parameter name and value pairs specified using \"%s\" option" % ParamsOptionName)
  532     
  533     ParamsOptionValueWords = ParamsOptionValue.split(",")
  534     if len(ParamsOptionValueWords) % 2:
  535         MiscUtil.PrintError("The number of comma delimited paramater names and values, %d, specified using \"%s\" option must be an even number." % (len(ParamsOptionValueWords), ParamsOptionName))
  536     
  537     # Validate paramater name and value pairs...
  538     for Index in range(0, len(ParamsOptionValueWords), 2):
  539         Name = ParamsOptionValueWords[Index].strip()
  540         Value = ParamsOptionValueWords[Index + 1].strip()
  541 
  542         CanonicalName = Name.lower()
  543         if  not CanonicalName in CanonicalParamNamesMap:
  544             MiscUtil.PrintError("The parameter name, %s, specified using \"%s\" is not a valid name. Supported parameter names: %s" % (Name, ParamsOptionName, " ".join(ValidParamNames)))
  545 
  546         ParamName = CanonicalParamNamesMap[CanonicalName]
  547         ParamValue = Value
  548         
  549         if re.match("^MaxIter$", ParamName, re.I):
  550             if not MiscUtil.IsInteger(Value):
  551                 MiscUtil.PrintError("The parameter value, %s, specified for parameter name, %s, using \"%s\" option must be an integer." % (Value, Name, ParamsOptionName))
  552             Value = int(Value)
  553             if Value <= 0:
  554                 MiscUtil.PrintError("The parameter value, %s, specified for parameter name, %s, using \"%s\" option is not a valid value. Supported values: > 0" % (Value, Name, ParamsOptionName))
  555             ParamValue = Value
  556         elif re.match("^(RestrainHydrogens|RemoveGridFiles)$", ParamName, re.I):
  557             if not re.match("^(yes|no)$", Value, re.I):
  558                 MiscUtil.PrintError("The parameter value, %s, specified for parameter name, %s, using \"%s\" option is not a valid value. Supported values: yes or no" % (Value, Name, ParamsOptionName))
  559             ParamValue = True if re.match("^yes$", Value, re.I) else False
  560         elif re.match("^(RespA|RespB|VDWPointDensity)$", ParamName, re.I):
  561             if not MiscUtil.IsNumber(Value):
  562                 MiscUtil.PrintError("The parameter value, %s, specified for parameter name, %s, using \"%s\" option must be a float." % (Value, Name, ParamsOptionName))
  563             Value = float(Value)
  564             if Value <= 0:
  565                 MiscUtil.PrintError("The parameter value, %s, specified for parameter name, %s, using \"%s\" option is not a valid value. Supported values: > 0" % (Value, Name, ParamsOptionName))
  566             ParamValue = Value
  567         elif re.match("^Tolerance$", ParamName, re.I):
  568             if not MiscUtil.IsNumber(Value):
  569                 MiscUtil.PrintError("The parameter value, %s, specified for parameter name, %s, using \"%s\" option must be a float." % (Value, Name, ParamsOptionName))
  570             Value = float(Value)
  571             if Value < 0:
  572                 MiscUtil.PrintError("The parameter value, %s, specified for parameter name, %s, using \"%s\" option is not a valid value. Supported values: >= 0" % (Value, Name, ParamsOptionName))
  573             ParamValue = Value
  574         elif re.match("^VDWScaleFactors$", ParamName, re.I):
  575             ScaleFactorsValue = Value.strip()
  576             if not ScaleFactorsValue:
  577                 MiscUtil.PrintError("No parameter value specified for parameter name, %s, using \"%s\" option." % (Name, ParamsOptionName))
  578             ScaleFactorsWords = ScaleFactorsValue.split()
  579             
  580             ScaleFactors = []
  581             LastScaleFactor = 0.0
  582             for ScaleFactor in ScaleFactorsWords:
  583                 if not MiscUtil.IsNumber(ScaleFactor):
  584                     MiscUtil.PrintError("The value, %s, in parameter value, %s, specified for parameter name, %s, using \"%s\" option must be a float." % (ScaleFactor, Value, Name, ParamsOptionName))
  585                 ScaleFactor = float(ScaleFactor)
  586                 if ScaleFactor <= 0:
  587                     MiscUtil.PrintError("The value, %s, in parameter value, %s, specified for parameter name, %s, using \"%s\" option is not a valid value. Supported values: > 0" % (ScaleFactor, Value, Name, ParamsOptionName))
  588                 if len(ScaleFactors):
  589                     if ScaleFactor <= LastScaleFactor:
  590                         MiscUtil.PrintError("The value, %s, in parameter value, %s, specified for parameter name, %s, using \"%s\" option is not a valid value. It must be greater than the previous value, %s, specified for the scale factor." % (ScaleFactor, Value, Name, ParamsOptionName, LastScaleFactor))
  591                         
  592                 LastScaleFactor = ScaleFactor
  593                 ScaleFactors.append(ScaleFactor)
  594             
  595             ParamValue = ScaleFactors
  596         elif re.match("^VDWRadii$", ParamName, re.I):
  597             RadiiValue = Value.strip()
  598             if not RadiiValue:
  599                 MiscUtil.PrintError("No parameter value specified for parameter name, %s, using \"%s\" option." % (Name, ParamsOptionName))
  600             RadiiWords = RadiiValue.split()
  601             if len(RadiiWords) % 2:
  602                 MiscUtil.PrintError("The number of space delimited values, %s, in parameter value, %s, specified for parameter name, %s, using \"%s\" option is not valid. It must be an even number." % (len(RadiiWords), Value, Name, ParamsOptionName))
  603             
  604             for RadiiWordsIndex in range(0, len(RadiiWords), 2):
  605                 ElementSymbol = RadiiWords[RadiiWordsIndex].upper()
  606                 VDWRadius = RadiiWords[RadiiWordsIndex + 1]
  607 
  608                 if not MiscUtil.IsNumber(VDWRadius):
  609                     MiscUtil.PrintError("The vdw radius value, %s, in parameter value, %s, specified for parameter name, %s, using \"%s\" option is not a valid." % (VDWRadius, Value, Name, ParamsOptionName))
  610                     
  611                 if not RDKitUtil.IsValidElementSymbol(ElementSymbol.capitalize()):
  612                     MiscUtil.PrintWarning("The element symbol, %s, in parameter value, %s, specified for parameter name, %s, using \"%s\" option is not a valid." % (ElementSymbol, Value, Name, ParamsOptionName))
  613                 VDWRadii[ElementSymbol] = float(VDWRadius)
  614             
  615             ParamValue = VDWRadii
  616         else:
  617             ParamValue = Value
  618             
  619         # Set value...
  620         ParamsInfo[ParamName] = ParamValue
  621     
  622     SetupChargesRespOptions(ParamsInfo)
  623 
  624 def SetupChargesRespOptions(ParamsInfo):
  625     """Setup options for calculating RESP charges. """
  626 
  627     # Initialize ESP options...
  628     ChargesRespOtions = {}
  629     
  630     ChargesRespOtions["METHOD_ESP"] = None
  631     ChargesRespOtions["BASIS_ESP"] = None
  632 
  633     # Setup ESP options...
  634     ParamNameToESPOptionID = {"MaxIter": "MAX_IT", "RespA": "RESP_A", "RespB": "RESP_B", "Tolerance": "TOLER", "VDWRadii": "VDW_RADII", "VDWScaleFactors": "VDW_SCALE_FACTORS", "VDWPointDensity": "VDW_POINT_DENSITY"}
  635     
  636     for ParamName in ParamNameToESPOptionID:
  637         ESPOptionID = ParamNameToESPOptionID[ParamName]
  638         ChargesRespOtions[ESPOptionID] = ParamsInfo[ParamName]
  639 
  640     # Setup IHFREE option...
  641     ChargesRespOtions["IHFREE"] = False if ParamsInfo["RestrainHydrogens"] else True
  642 
  643     OptionsInfo["ChargesRespParams"] = ParamsInfo
  644     OptionsInfo["ChargesRespOtions"] = ChargesRespOtions
  645 
  646 def CheckSupportForRESPCalculations():
  647     """Check support for RESP calculations. """
  648 
  649     if not OptionsInfo["MPMode"]:
  650         return
  651     
  652     if not OptionsInfo["RESPChargesTypeMode"]:
  653         return
  654     
  655     MiscUtil.PrintInfo("")
  656     MiscUtil.PrintError("Multiprocessing is not supported during the calculation of RSEP partial atomic\ncharges. The RESP module is not conducive for multiprocessing. The names of the results\nand grid output files are not unique for molecules during the RESP calculations spread\nacross multiple processes.")
  657 
  658 def ProcessOptions():
  659     """Process and validate command line arguments and options"""
  660     
  661     MiscUtil.PrintInfo("Processing options...")
  662     
  663     # Validate options...
  664     ValidateOptions()
  665     
  666     OptionsInfo["Infile"] = Options["--infile"]
  667     ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
  668     OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"], InfileName = Options["--infile"], ParamsDefaultInfo = ParamsDefaultInfoOverride)
  669     
  670     OptionsInfo["Outfile"] = Options["--outfile"]
  671     OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters("--outfileParams", Options["--outfileParams"])
  672     
  673     OptionsInfo["Overwrite"] = Options["--overwrite"]
  674     
  675     # Method, basis set, and reference wavefunction...
  676     OptionsInfo["BasisSet"] = Options["--basisSet"]
  677     OptionsInfo["BasisSetAuto"] = True if re.match("^auto$", Options["--basisSet"], re.I) else False
  678     
  679     OptionsInfo["MethodName"] = Options["--methodName"]
  680     OptionsInfo["MethodNameAuto"] = True if re.match("^auto$", Options["--methodName"], re.I) else False
  681     
  682     OptionsInfo["Reference"] = Options["--reference"]
  683     OptionsInfo["ReferenceAuto"] = True if re.match("^auto$", Options["--reference"], re.I) else False
  684     
  685     # Run and options parameters...
  686     OptionsInfo["Psi4OptionsParams"] = Psi4Util.ProcessPsi4OptionsParameters("--psi4OptionsParams", Options["--psi4OptionsParams"])
  687     OptionsInfo["Psi4RunParams"] = Psi4Util.ProcessPsi4RunParameters("--psi4RunParams", Options["--psi4RunParams"], InfileName = OptionsInfo["Infile"])
  688 
  689     ChargesType = Options["--chargesType"]
  690     ChargesPropType = None
  691     RESPChargesTypeMode = False
  692     if re.match("^Mulliken$", ChargesType, re.I):
  693         ChargesType = 'Mulliken'
  694         ChargesPropType = 'MULLIKEN_CHARGES'
  695     elif re.match("^Lowdin$", ChargesType, re.I):
  696         ChargesType = 'Lowdin'
  697         ChargesPropType = 'LOWDIN_CHARGES'
  698     elif re.match("^RESP$", ChargesType, re.I):
  699         ChargesType = 'RESP'
  700         ChargesPropType = None
  701         RESPChargesTypeMode = True
  702     else:
  703         MiscUtil.PrintError("The value, %s, specified for charge mode is not supported. " % Options["--chargesType"])
  704         
  705     OptionsInfo["ChargesType"] = ChargesType
  706     OptionsInfo["ChargesPropType"] = ChargesPropType
  707     OptionsInfo["RESPChargesTypeMode"] = RESPChargesTypeMode
  708 
  709     ProcessOptionChargesRespParameters()
  710     
  711     AtomAliasesFormatMode = True
  712     if re.match("^DataField", Options["--chargesSDFormat"], re.I):
  713         AtomAliasesFormatMode = False
  714     OptionsInfo["AtomAliasesFormatMode"] = AtomAliasesFormatMode
  715 
  716     DataFieldLabel = Options["--dataFieldLabel"]
  717     if re.match("^auto$", DataFieldLabel, re.I):
  718         DataFieldLabel = "Psi4_%s_Charges (a.u.)" % Options["--chargesType"]
  719     OptionsInfo["DataFieldLabel"] = DataFieldLabel
  720     
  721     OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
  722     OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
  723     
  724     OptionsInfo["Precision"] = int(Options["--precision"])
  725     
  726     OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
  727 
  728 def RetrieveOptions():
  729     """Retrieve command line arguments and options"""
  730     
  731     # Get options...
  732     global Options
  733     Options = docopt(_docoptUsage_)
  734     
  735     # Set current working directory to the specified directory...
  736     WorkingDir = Options["--workingdir"]
  737     if WorkingDir:
  738         os.chdir(WorkingDir)
  739     
  740     # Handle examples option...
  741     if "--examples" in Options and Options["--examples"]:
  742         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
  743         sys.exit(0)
  744 
  745 def ValidateOptions():
  746     """Validate option values"""
  747 
  748     MiscUtil.ValidateOptionTextValue("-c, --chargesType", Options["--chargesType"], "Mulliken Lowdin RESP")
  749     MiscUtil.ValidateOptionTextValue("--chargesSDFormat", Options["--chargesSDFormat"], "AtomAliases DataField")
  750     
  751     MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
  752     MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol")
  753     
  754     MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
  755     MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
  756     MiscUtil.ValidateOptionsDistinctFileNames("-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"])
  757     
  758     MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
  759     MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
  760     
  761     MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0})
  762     
  763 # Setup a usage string for docopt...
  764 _docoptUsage_ = """
  765 Psi4CalculatePartialCharges.py - Calculate partial atomic charges
  766 
  767 Usage:
  768     Psi4CalculatePartialCharges.py [--basisSet <text>] [--chargesType <Mulliken or Lowdin>] [--chargesRespParams <Name,Value,...>]
  769                                    [--chargesSDFormat <AtomAliases or DataField>] [--dataFieldLabel <text>] [--infileParams <Name,Value,...>]
  770                                    [--methodName <text>] [--mp <yes or no>] [--mpParams <Name, Value,...>] [ --outfileParams <Name,Value,...> ]
  771                                    [--overwrite] [--precision <number>] [--psi4OptionsParams <Name,Value,...>] [--psi4RunParams <Name,Value,...>] 
  772                                    [--quiet <yes or no>] [--reference <text>] [-w <dir>] -i <infile> -o <outfile> 
  773     Psi4CalculatePartialCharges.py -h | --help | -e | --examples
  774 
  775 Description:
  776     Calculate partial atomic charges for molecules using a specified method name
  777     and basis set. The molecules must have 3D coordinates in input file. The molecular
  778     geometry is not optimized before the calculation. In addition, hydrogens must
  779     be present for all molecules in input file. A single point energy calculation is 
  780     performed before calculating the partial atomic charges. The 3D coordinates
  781     are not modified during the calculation.
  782     
  783     A Psi4 XYZ format geometry string is automatically generated for each molecule
  784     in input file. It contains atom symbols and 3D coordinates for each atom in a
  785     molecule. In addition, the formal charge and spin multiplicity are present in the
  786     the geometry string. These values are either retrieved from molecule properties
  787     named 'FormalCharge' and 'SpinMultiplicty' or dynamically calculated for a
  788     molecule.
  789 
  790     The supported input file formats are: Mol (.mol), SD (.sdf, .sd)
  791 
  792     The supported output file formats are: SD (.sdf, .sd)
  793 
  794 Options:
  795     -b, --basisSet <text>  [default: auto]
  796         Basis set to use for calculating single point energy before partial atomic
  797         charges. Default: 6-31+G** for sulfur containing molecules; Otherwise,
  798         6-31G** [ Ref 150 ]. The specified value must be a valid Psi4 basis set.
  799         No validation is performed.
  800         
  801         The following list shows a representative sample of basis sets available
  802         in Psi4:
  803             
  804             STO-3G, 6-31G, 6-31+G, 6-31++G, 6-31G*, 6-31+G*,  6-31++G*, 
  805             6-31G**, 6-31+G**, 6-31++G**, 6-311G, 6-311+G, 6-311++G,
  806             6-311G*, 6-311+G*, 6-311++G*, 6-311G**, 6-311+G**, 6-311++G**,
  807             cc-pVDZ, cc-pCVDZ, aug-cc-pVDZ, cc-pVDZ-DK, cc-pCVDZ-DK, def2-SVP,
  808             def2-SVPD, def2-TZVP, def2-TZVPD, def2-TZVPP, def2-TZVPPD
  809             
  810     -c, --chargesType <Mulliken, Lowdin, or RESP>  [default: Mulliken]
  811         Type of partial atomic charges to calculate. Possible values: Mulliken, Lowdin,
  812         or RESP [ Ref 158 ]. Multiprocessing is not supported during the calculation
  813         of RSEP charges. In addition, the RSEP calculation relies on the presence of
  814         the RESP Psi4 Plugin in your environment.
  815     --chargesRespParams <Name,Value,...>  [default: auto]
  816         A comma delimited list of parameter name and value  pairs for calculating
  817         RESP [ Ref 158 ] charges. A space is used as a delimiter for multiple values in
  818         a name and value pair. The supported parameter names, along with
  819         their default values, are shown below:
  820             
  821             maxIter, 25
  822             restrainHydrogens, no
  823             removeGridFiles, yes
  824             respA, 0.0005
  825             respB, 0.1
  826             tolerance, 1e-5
  827             vdwRadii, auto
  828             vdwScaleFactors, 1.4 1.6 1.8 2.0
  829             vdwPointDensity, 1.0
  830             
  831         maxIter: Maximum number of iterations to perform during charge fitting.
  832          
  833         restrainHydrogens: Restrain hydrogens during charge fitting.
  834         
  835         removeGridFiles: Keep or remove the following ESP grid and output files:
  836         1_default_grid.dat, 1_default_grid_esp.dat, results.out. The output
  837         files are removed by default. You may optionally keep the output files. The
  838         output files are automatically renamed to the following file for 'No' value of
  839         'removeGridFiles': Mol<MolNum>_grid.dat, Mol<MolNum>_grid_esp.dat,
  840         Mol<MolNum>_resp_results.out.
  841         
  842         respA: Scale factor defining the  asymptotic limits of the strength of the
  843         restraint.
  844         
  845         respB: The 'tightness' of the hyperbola around its minimum for the
  846         restraint.
  847         
  848         tolerance: Tolerance for charges during charge fitting to the ESP.
  849         
  850         vdwRadii: vdw radii for elements in angstroms. It's a space delimited list of
  851         element symbol and radius value pairs. The default list is shown below:
  852             
  853             H 1.20 He 1.20 Li 1.37 Be 1.45 B 1.45 C 1.50 N 1.50 O 1.40 F 1.35
  854             Ne 1.30 Na 1.57 Mg 1.36 Al 1.24 Si 1.17P 1.80 S 1.75 Cl 1.7
  855             
  856         You may specify all or a subset of element symbol and vdw radius pairs to
  857         update the default values.
  858         
  859         vdwScaleFactors: The vdw radii are scaled by the scale factors to set the
  860         grid points at the shells for calculating the ESP using quantum methodology.
  861         The default number of shells is 4 and corresponds to the number of vdw
  862         scale factors.The 'shell' points are written to a grid file for calculating the ESP.
  863         
  864         vdwPointDensity: Approximate number of points to generate per square
  865         angstrom surface area.
  866     --chargesSDFormat <AtomAliases or DataField>  [default: AtomAliases]
  867         Format for writing out partial atomic charges to SD file. Possible values:
  868         AtomAliases or DataField.
  869         
  870         The charges are stored as atom property named 'molFileAlias' for
  871         'AtomAliases' format and may be retrieved using the RDKit function
  872         'GetProp' for atoms: Aotm.GetProp('molFileAliases').
  873         
  874         The charges are stored under a data field label specified using
  875         '-d, --dataFieldLabel' for 'DataField' format and may be retrieved using the
  876         RDKit function 'GetProp' for molecules.
  877     -d, --dataFieldLabel <text>  [default: auto]
  878         Data field label to use for storing charged in SD file during 'DataField' value
  879         of '-c, --chargesSDFormat'. Default: Psi4_<ChargesType>_Charges (a.u.)
  880     -e, --examples
  881         Print examples.
  882     -h, --help
  883         Print this help message.
  884     -i, --infile <infile>
  885         Input file name.
  886     --infileParams <Name,Value,...>  [default: auto]
  887         A comma delimited list of parameter name and value pairs for reading
  888         molecules from files. The supported parameter names for different file
  889         formats, along with their default values, are shown below:
  890             
  891             SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
  892             
  893     -m, --methodName <text>  [default: auto]
  894         Method to use for calculating single point energy before partial atomic
  895         charges. 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
  977             
  978     --overwrite
  979         Overwrite existing files.
  980     --precision <number>  [default: 4]
  981         Floating point precision for writing energy values.
  982     --psi4OptionsParams <Name,Value,...>  [default: none]
  983         A comma delimited list of Psi4 option name and value pairs for setting
  984         global and module options. The names are 'option_name' for global options
  985         and 'module_name__option_name' for options local to a module. The
  986         specified option names must be valid Psi4 names. No validation is
  987         performed.
  988         
  989         The specified option name and  value pairs are processed and passed to
  990         psi4.set_options() as a dictionary. The supported value types are float,
  991         integer, boolean, or string. The float value string is converted into a float.
  992         The valid values for a boolean string are yes, no, true, false, on, or off. 
  993     --psi4RunParams <Name,Value,...>  [default: auto]
  994         A comma delimited list of parameter name and value pairs for configuring
  995         Psi4 jobs.
  996         
  997         The supported parameter names along with their default and possible
  998         values are shown below:
  999              
 1000             MemoryInGB, 1
 1001             NumThreads, 1
 1002             OutputFile, auto   [ Possible  values: stdout, quiet, or FileName ]
 1003             ScratchDir, auto   [ Possivle values: DirName]
 1004             RemoveOutputFile, yes   [ Possible values: yes, no, true, or false]
 1005             
 1006         These parameters control the runtime behavior of Psi4.
 1007         
 1008         The default file name for 'OutputFile' is <InFileRoot>_Psi4.out. The PID
 1009         is appended to output file name during multiprocessing as shown:
 1010         <InFileRoot>_Psi4_<PIDNum>.out. The 'stdout' value for 'OutputType'
 1011         sends Psi4 output to stdout. The 'quiet' or 'devnull' value suppresses
 1012         all Psi4 output.
 1013         
 1014         The default 'Yes' value of 'RemoveOutputFile' option forces the removal
 1015         of any existing Psi4 before creating new files to append output from
 1016         multiple Psi4 runs.
 1017         
 1018         The option 'ScratchDir' is a directory path to the location of scratch
 1019         files. The default value corresponds to Psi4 default. It may be used to
 1020         override the deafult path.
 1021     -q, --quiet <yes or no>  [default: no]
 1022         Use quiet mode. The warning and information messages will not be printed.
 1023     -r, --reference <text>  [default: auto]
 1024         Reference wave function to use for calculating single point energy before
 1025         partial atomic charges. Default: RHF or UHF. The default values are Restricted
 1026         Hartree-Fock (RHF) for closed-shell molecules with all electrons paired and
 1027         Unrestricted Hartree-Fock (UHF) for open-shell molecules with unpaired electrons.
 1028         
 1029         The specified value must be a valid Psi4 reference wave function. No validation
 1030         is performed. For example: ROHF, CUHF, RKS, etc.
 1031         
 1032         The spin multiplicity determines the default value of reference wave function
 1033         for input molecules. It is calculated from number of free radical electrons using
 1034         Hund's rule of maximum multiplicity defined as 2S + 1 where S is the total
 1035         electron spin. The total spin is 1/2 the number of free radical electrons in a 
 1036         molecule. The value of 'SpinMultiplicity' molecule property takes precedence
 1037         over the calculated value of spin multiplicity.
 1038     -w, --workingdir <dir>
 1039         Location of working directory which defaults to the current directory.
 1040 
 1041 Examples:
 1042     To calculate Mulliken partial atomic charges using  B3LYP/6-31G** and
 1043     B3LYP/6-31+G** for non-sulfur and sulfur containing molecules in a SD
 1044     file with 3D structures, use RHF and UHF for closed-shell and open-shell
 1045     molecules, and write a new SD file, type:
 1046 
 1047         % Psi4CalculatePartialCharges.py  -i Psi4Sample3D.sdf 
 1048           -o Psi4Sample3DOut.sdf
 1049 
 1050     To run the first example for calculating RESP charges using a default set of
 1051     parameters for the RESP calculation and write out a SD file, type:
 1052 
 1053         % Psi4CalculatePartialCharges.py  --chargesType RESP
 1054            -i Psi4Sample3D.sdf -o Psi4Sample3DOut.sdf
 1055 
 1056     To run the first example for calculating RESP charges using an explicit set
 1057     of specific parameters for the RESP calculation and write out a SD file, type:
 1058 
 1059         % Psi4CalculatePartialCharges.py  --chargesType RESP
 1060            --chargesRespParams "respA, 0.0005, respB, 0.1, vdwScaleFactors,
 1061            1.4 1.6 1.8 2.0" -i Psi4Sample3D.sdf -o Psi4Sample3DOut.sdf
 1062 
 1063     To run the first example in multiprocessing mode on all available CPUs
 1064     without loading all data into memory and write out a SD file, type:
 1065 
 1066         % Psi4CalculatePartialCharges.py --mp yes -i Psi4Sample3D.sdf
 1067           -o Psi4Sample3DOut.sdf
 1068 
 1069     To run the first example in multiprocessing mode on all available CPUs
 1070     by loading all data into memory and write out a SD file, type:
 1071 
 1072         % Psi4CalculatePartialCharges.py  --mp yes --mpParams "inputDataMode,
 1073           InMemory" -i Psi4Sample3D.sdf  -o Psi4Sample3DOut.sdf
 1074 
 1075     To run the first example in multiprocessing mode on all available CPUs
 1076     without loading all data into memory along with multiple threads for each
 1077     Psi4 run and write out a SD file, type:
 1078 
 1079         % Psi4CalculatePartialCharges.py --mp yes --psi4RunParams "NumThreads,2"
 1080            -i Psi4Sample3D.sdf -o Psi4Sample3DOut.sdf
 1081 
 1082     To run the first example for writing out charges to a new SD file under a
 1083     datafield instead of storing them as atom property, type:
 1084 
 1085         % Psi4CalculatePartialCharges.py  --chargesSDFormat DataField
 1086           -i Psi4Sample3D.sdf  -o Psi4Sample3DOut.sdf
 1087 
 1088     To calculate specific partial atomic charges using a specific method  and basis
 1089     set for molecules in a SD ontaining 3D structures and write them out to a specific
 1090     datafield in a new SD file, type:
 1091 
 1092         % Psi4CalculatePartialCharges.py  -c Lowdin -m SCF -b aug-cc-pVDZ
 1093           --chargesSDFormat DataField --dataFieldLabel "Lowdin_Charges"
 1094           -i Psi4Sample3D.sdf  -o Psi4Sample3DOut.sdf
 1095 
 1096 Author:
 1097     Manish Sud(msud@san.rr.com)
 1098 
 1099 See also:
 1100     Psi4CalculateEnergy.py, Psi4PerformMinimization.py, Psi4GenerateConformers.py
 1101 
 1102 Copyright:
 1103     Copyright (C) 2022 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()