MayaChemTools

    1 #!/bin/env python
    2 #
    3 # File: Psi4PerformTorsionScan.py
    4 # Author: Manish Sud <msud@san.rr.com>
    5 #
    6 # Copyright (C) 2021 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 
   31 from __future__ import print_function
   32 
   33 # Add local python path to the global path and import standard library modules...
   34 import os
   35 import sys;  sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
   36 import time
   37 import re
   38 import glob
   39 import shutil
   40 import multiprocessing as mp
   41 
   42 import matplotlib.pyplot as plt
   43 import seaborn as sns
   44 
   45 # Psi4 imports...
   46 if (hasattr(shutil, 'which') and shutil.which("psi4") is None):
   47     sys.stderr.write("\nWarning: Failed to find 'psi4' in your PATH indicating potential issues with your\n")
   48     sys.stderr.write("Psi4 environment. The 'import psi4' directive in the global scope of the script\n")
   49     sys.stderr.write("interferes with the multiprocessing functionality. It is imported later in the\n")
   50     sys.stderr.write("local scope during the execution of the script and may fail. Check/update your\n")
   51     sys.stderr.write("Psi4 environment and try again.\n\n")
   52 
   53 # RDKit imports...
   54 try:
   55     from rdkit import rdBase
   56     from rdkit import Chem
   57     from rdkit.Chem import AllChem
   58     from rdkit.Chem import rdMolAlign
   59     from rdkit.Chem import rdMolTransforms
   60 except ImportError as ErrMsg:
   61     sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
   62     sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
   63     sys.exit(1)
   64 
   65 # MayaChemTools imports...
   66 try:
   67     from docopt import docopt
   68     import MiscUtil
   69     import Psi4Util
   70     import RDKitUtil
   71 except ImportError as ErrMsg:
   72     sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
   73     sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
   74     sys.exit(1)
   75 
   76 ScriptName = os.path.basename(sys.argv[0])
   77 Options = {}
   78 OptionsInfo = {}
   79 
   80 def main():
   81     """Start execution of the script"""
   82     
   83     MiscUtil.PrintInfo("\n%s (Psi4: Imported later; RDK v%s; %s): Starting...\n" % (ScriptName, rdBase.rdkitVersion, time.asctime()))
   84     
   85     (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
   86     
   87     # Retrieve command line arguments and options...
   88     RetrieveOptions()
   89     
   90     # Process and validate command line arguments and options...
   91     ProcessOptions()
   92     
   93     # Perform actions required by the script...
   94     PerformTorsionScan()
   95 
   96     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
   97     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
   98 
   99 def PerformTorsionScan():
  100     """Perform torsion scan."""
  101     
  102     # Setup a molecule reader for input file...
  103     MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
  104     OptionsInfo["InfileParams"]["AllowEmptyMols"] = True
  105     Mols  = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
  106     
  107     PlotExt = OptionsInfo["OutPlotParams"]["OutExt"]
  108     FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
  109     MiscUtil.PrintInfo("Generating output files %s_*.sdf, %s_*Torsion*Match*.sdf, %s_*Torsion*Match*Energies.csv, %s_*Torsion*Match*Plot.%s..." % (FileName, FileName, FileName, FileName, PlotExt))
  110 
  111     MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount = ProcessMolecules(Mols)
  112     
  113     MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
  114     MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
  115     MiscUtil.PrintInfo("Number of molecules failed during initial minimization: %d" % MinimizationFailedCount)
  116     MiscUtil.PrintInfo("Number of molecules without any matched torsions: %d" % TorsionsMissingCount)
  117     MiscUtil.PrintInfo("Number of molecules failed during torsion scan: %d" % TorsionsScanFailedCount)
  118     MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + TorsionsMissingCount + MinimizationFailedCount + TorsionsScanFailedCount))
  119 
  120 def ProcessMolecules(Mols):
  121     """Process molecules to perform torsion scan. """
  122 
  123     if OptionsInfo["MPMode"]:
  124         return ProcessMoleculesUsingMultipleProcesses(Mols)
  125     else:
  126         return ProcessMoleculesUsingSingleProcess(Mols)
  127 
  128 def ProcessMoleculesUsingSingleProcess(Mols):
  129     """Process molecules to perform torsion scan using a single process."""
  130 
  131     # Intialize Psi4...
  132     MiscUtil.PrintInfo("\nInitializing Psi4...")
  133     Psi4Handle = Psi4Util.InitializePsi4(Psi4RunParams = OptionsInfo["Psi4RunParams"], Psi4OptionsParams = OptionsInfo["Psi4OptionsParams"], PrintVersion = True, PrintHeader = True)
  134     OptionsInfo["psi4"] = Psi4Handle
  135 
  136     # Setup max iterations global variable...
  137     Psi4Util.UpdatePsi4OptionsParameters(Psi4Handle, {'GEOM_MAXITER': OptionsInfo["MaxIters"]})
  138     
  139     # Setup conversion factor for energy units...
  140     SetupEnergyConversionFactor(Psi4Handle)
  141     
  142     MolInfoText = "first molecule"
  143     if not OptionsInfo["FirstMolMode"]:
  144         MolInfoText = "all molecules"
  145 
  146     if OptionsInfo["TorsionMinimize"]:
  147         MiscUtil.PrintInfo("\nPeforming torsion scan on %s by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  148     else:
  149         MiscUtil.PrintInfo("\nPeforming torsion scan on %s by skipping generation of conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  150     
  151     SetupTorsionsPatternsInfo()
  152     
  153     (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
  154 
  155     for Mol in Mols:
  156         MolCount += 1
  157 
  158         if OptionsInfo["FirstMolMode"] and MolCount > 1:
  159             MolCount -= 1
  160             break
  161         
  162         if not CheckAndValidateMolecule(Mol, MolCount):
  163             continue
  164 
  165         # Setup 2D coordinates for SMILES input file...
  166         if OptionsInfo["SMILESInfileStatus"]:
  167             AllChem.Compute2DCoords(Mol)
  168         
  169         ValidMolCount += 1
  170 
  171         Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(Mol, MolCount)
  172         
  173         if not MinimizationCalcStatus:
  174             MinimizationFailedCount += 1
  175             continue
  176         
  177         if not TorsionsMatchStatus:
  178             TorsionsMissingCount += 1
  179             continue
  180 
  181         if not TorsionsScanCalcStatus:
  182             TorsionsScanFailedCount += 1
  183             continue
  184 
  185     return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
  186 
  187 def ProcessMoleculesUsingMultipleProcesses(Mols):
  188     """Process and minimize molecules using multiprocessing."""
  189     
  190     if OptionsInfo["MPLevelTorsionAnglesMode"]:
  191         return ProcessMoleculesUsingMultipleProcessesAtTorsionAnglesLevel(Mols)
  192     elif OptionsInfo["MPLevelMoleculesMode"]:
  193         return ProcessMoleculesUsingMultipleProcessesAtMoleculesLevel(Mols)
  194     else:
  195         MiscUtil.PrintError("The value, %s,  option \"--mpLevel\" is not supported." % (OptionsInfo["MPLevel"]))
  196         
  197 def ProcessMoleculesUsingMultipleProcessesAtMoleculesLevel(Mols):
  198     """Process molecules to perform torsion scan using multiprocessing at molecules level."""
  199 
  200     MolInfoText = "first molecule"
  201     if not OptionsInfo["FirstMolMode"]:
  202         MolInfoText = "all molecules"
  203 
  204     if OptionsInfo["TorsionMinimize"]:
  205         MiscUtil.PrintInfo("\nPeforming torsion scan on %s using multiprocessing at molecules level by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  206     else:
  207         MiscUtil.PrintInfo("\nPeforming torsion scan %s using multiprocessing at molecules level by skipping generation of conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  208         
  209     MPParams = OptionsInfo["MPParams"]
  210     
  211     # Setup data for initializing a worker process...
  212     InitializeWorkerProcessArgs = (MiscUtil.ObjectToBase64EncodedString(Options), MiscUtil.ObjectToBase64EncodedString(OptionsInfo))
  213 
  214     if OptionsInfo["FirstMolMode"]:
  215         Mol = Mols[0]
  216         Mols = [Mol]
  217 
  218     # Setup a encoded mols data iterable for a worker process...
  219     WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
  220 
  221     # Setup process pool along with data initialization for each process...
  222     MiscUtil.PrintInfo("\nConfiguring multiprocessing using %s method..." % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()"))
  223     MiscUtil.PrintInfo("NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" % (MPParams["NumProcesses"], MPParams["InputDataMode"], ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"])))
  224     
  225     ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
  226     
  227     # Start processing...
  228     if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
  229         Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
  230     elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
  231         Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
  232     else:
  233         MiscUtil.PrintError("The value, %s, specified for \"--inputDataMode\" is not supported." % (MPParams["InputDataMode"]))
  234 
  235     (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
  236     
  237     for Result in Results:
  238         MolCount += 1
  239         
  240         MolIndex, EncodedMol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = Result
  241         
  242         if EncodedMol is None:
  243             continue
  244         ValidMolCount += 1
  245     
  246         if not MinimizationCalcStatus:
  247             MinimizationFailedCount += 1
  248             continue
  249         
  250         if not TorsionsMatchStatus:
  251             TorsionsMissingCount += 1
  252             continue
  253 
  254         if not TorsionsScanCalcStatus:
  255             TorsionsScanFailedCount += 1
  256             continue
  257         
  258         Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
  259         
  260     return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
  261     
  262 def InitializeWorkerProcess(*EncodedArgs):
  263     """Initialize data for a worker process."""
  264     
  265     global Options, OptionsInfo
  266 
  267     if not OptionsInfo["QuietMode"]:
  268         MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
  269 
  270     # Decode Options and OptionInfo...
  271     Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
  272     OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
  273 
  274     # Initialize torsion patterns info...
  275     SetupTorsionsPatternsInfo()
  276     
  277     # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
  278     # output files for each process...
  279     OptionsInfo["Psi4Initialized"]  = False
  280 
  281 def WorkerProcess(EncodedMolInfo):
  282     """Process data for a worker process."""
  283 
  284     if not OptionsInfo["Psi4Initialized"]:
  285         InitializePsi4ForWorkerProcess()
  286     
  287     MolIndex, EncodedMol = EncodedMolInfo
  288 
  289     MolNum = MolIndex + 1
  290     (MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus) = [False] * 3
  291     
  292     if EncodedMol is None:
  293         return [MolIndex, None, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
  294     
  295     Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
  296     if not CheckAndValidateMolecule(Mol, MolNum):
  297         return [MolIndex, None, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
  298     
  299     # Setup 2D coordinates for SMILES input file...
  300     if OptionsInfo["SMILESInfileStatus"]:
  301         AllChem.Compute2DCoords(Mol)
  302     
  303     Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(Mol, MolNum)
  304     
  305     return [MolIndex, RDKitUtil.MolToBase64EncodedMolString(Mol, PropertyPickleFlags = Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps), MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
  306 
  307 def ProcessMoleculesUsingMultipleProcessesAtTorsionAnglesLevel(Mols):
  308     """Process molecules to perform torsion scan using multiprocessing at torsion angles level."""
  309 
  310     MolInfoText = "first molecule"
  311     if not OptionsInfo["FirstMolMode"]:
  312         MolInfoText = "all molecules"
  313 
  314     if OptionsInfo["TorsionMinimize"]:
  315         MiscUtil.PrintInfo("\nPeforming torsion scan on %s using multiprocessing at torsion angles level by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  316     else:
  317         MiscUtil.PrintInfo("\nPeforming torsion scan %s using multiprocessing at torsion angles level by skipping generation of conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  318     
  319     SetupTorsionsPatternsInfo()
  320     
  321     (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
  322 
  323     for Mol in Mols:
  324         MolCount += 1
  325 
  326         if OptionsInfo["FirstMolMode"] and MolCount > 1:
  327             MolCount -= 1
  328             break
  329         
  330         if not CheckAndValidateMolecule(Mol, MolCount):
  331             continue
  332 
  333         # Setup 2D coordinates for SMILES input file...
  334         if OptionsInfo["SMILESInfileStatus"]:
  335             AllChem.Compute2DCoords(Mol)
  336         
  337         ValidMolCount += 1
  338 
  339         Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(Mol, MolCount, UseMultiProcessingAtTorsionAnglesLevel = True)
  340         
  341         if not MinimizationCalcStatus:
  342             MinimizationFailedCount += 1
  343             continue
  344         
  345         if not TorsionsMatchStatus:
  346             TorsionsMissingCount += 1
  347             continue
  348 
  349         if not TorsionsScanCalcStatus:
  350             TorsionsScanFailedCount += 1
  351             continue
  352 
  353     return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
  354 
  355 def ScanSingleTorsionInMolUsingMultipleProcessesAtTorsionAnglesLevel(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum):
  356     """Perform torsion scan for a molecule using multiple processses at torsion angles
  357     level along with constrained energy minimization.
  358     """
  359     
  360     if OptionsInfo["MPLevelMoleculesMode"]:
  361         MiscUtil.PrintError("Single torison scanning for a molecule is not allowed in multiprocessing mode at molecules level.\n")
  362 
  363     Mols, Angles = SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum)
  364 
  365     MPParams = OptionsInfo["MPParams"]
  366     
  367     # Setup data for initializing a worker process...
  368     
  369     # Track and avoid encoding TorsionsPatternsInfo as it contains RDKit molecule object...
  370     TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
  371     OptionsInfo["TorsionsPatternsInfo"] = None
  372     
  373     InitializeWorkerProcessArgs = (MiscUtil.ObjectToBase64EncodedString(Options), MiscUtil.ObjectToBase64EncodedString(OptionsInfo))
  374     
  375     # Restore TorsionsPatternsInfo...
  376     OptionsInfo["TorsionsPatternsInfo"] = TorsionsPatternsInfo
  377 
  378     # Setup a encoded mols data iterable for a worker process...
  379     WorkerProcessDataIterable = GenerateBase64EncodedMolStringsWithTorsionScanInfo(Mol, (MolNum -1), Mols, Angles, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches)
  380     
  381     # Setup process pool along with data initialization for each process...
  382     MiscUtil.PrintInfo("\nConfiguring multiprocessing using %s method..." % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()"))
  383     MiscUtil.PrintInfo("NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" % (MPParams["NumProcesses"], MPParams["InputDataMode"], ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"])))
  384     
  385     ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeTorsionAngleWorkerProcess, InitializeWorkerProcessArgs)
  386     
  387     # Start processing...
  388     if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
  389         Results = ProcessPool.imap(TorsionAngleWorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
  390     elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
  391         Results = ProcessPool.map(TorsionAngleWorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
  392     else:
  393         MiscUtil.PrintError("The value, %s, specified for \"--inputDataMode\" is not supported." % (MPParams["InputDataMode"]))
  394 
  395     TorsionMols = []
  396     TorsionEnergies = []
  397     TorsionAngles = []
  398     
  399     for Result in Results:
  400         EncodedTorsionMol, CalcStatus, Angle, Energy = Result
  401         
  402         if not CalcStatus:
  403             return (Mol, False, None, None, None)
  404 
  405         if EncodedTorsionMol is None:
  406             return (Mol, False, None, None, None)
  407         TorsionMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionMol)
  408         
  409         TorsionMols.append(TorsionMol)
  410         TorsionEnergies.append(Energy)
  411         TorsionAngles.append(Angle)
  412     
  413     return (Mol, True, TorsionMols, TorsionEnergies, TorsionAngles)
  414 
  415 def InitializeTorsionAngleWorkerProcess(*EncodedArgs):
  416     """Initialize data for a worker process."""
  417     
  418     global Options, OptionsInfo
  419 
  420     if not OptionsInfo["QuietMode"]:
  421         MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
  422 
  423     # Decode Options and OptionInfo...
  424     Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
  425     OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
  426 
  427     # Initialize torsion patterns info...
  428     SetupTorsionsPatternsInfo()
  429     
  430     # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
  431     # output files for each process...
  432     OptionsInfo["Psi4Initialized"]  = False
  433 
  434 def TorsionAngleWorkerProcess(EncodedMolInfo):
  435     """Process data for a worker process."""
  436 
  437     if not OptionsInfo["Psi4Initialized"]:
  438         InitializePsi4ForWorkerProcess()
  439     
  440     MolIndex, EncodedMol, EncodedTorsionMol, TorsionAngle, TorsionID, TorsionPattern, EncodedTorsionPatternMol, TorsionMatches = EncodedMolInfo
  441 
  442     if EncodedMol is None or EncodedTorsionMol is None or EncodedTorsionPatternMol is None:
  443         return (None, False, None, None)
  444     
  445     Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
  446     TorsionMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionMol)
  447     TorsionPatternMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedTorsionPatternMol)
  448     
  449     TorsionMol, CalcStatus, Energy = MinimizeCalculateEnergyForTorsionMol(Mol, TorsionMol, TorsionAngle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, (MolIndex + 1))
  450 
  451     return (RDKitUtil.MolToBase64EncodedMolString(TorsionMol, PropertyPickleFlags = Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps), CalcStatus, TorsionAngle, Energy)
  452 
  453 def GenerateBase64EncodedMolStringsWithTorsionScanInfo(Mol, MolIndex, TorsionMols, TorsionAngles, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, PropertyPickleFlags = Chem.PropertyPickleOptions.AllProps):
  454     """Set up an iterator for generating base64 encoded molecule string for
  455     a torsion in a molecule along with appropriate trosion scan information.
  456     """
  457     
  458     for Index, TorsionMol in enumerate(TorsionMols):
  459         yield [MolIndex, None, None, TorsionAngles[Index], TorsionID, TorsionPattern, None, TorsionMatches] if (Mol is None or TorsionMol is None) else [MolIndex, RDKitUtil.MolToBase64EncodedMolString(Mol, PropertyPickleFlags), RDKitUtil.MolToBase64EncodedMolString(TorsionMol, PropertyPickleFlags), TorsionAngles[Index], TorsionID, TorsionPattern, RDKitUtil.MolToBase64EncodedMolString(TorsionPatternMol, PropertyPickleFlags), TorsionMatches]
  460 
  461 def InitializePsi4ForWorkerProcess():
  462     """Initialize Psi4 for a worker process."""
  463     
  464     if OptionsInfo["Psi4Initialized"]:
  465         return
  466 
  467     OptionsInfo["Psi4Initialized"] = True
  468 
  469     if OptionsInfo["MPLevelTorsionAnglesMode"] and re.match("auto", OptionsInfo["Psi4RunParams"]["OutputFileSpecified"], re.I):
  470         # Run Psi4 in quiet mode during multiprocessing at Torsions level for 'auto' OutputFile...
  471         OptionsInfo["Psi4RunParams"]["OutputFile"] = "quiet"
  472     else:
  473         # Update output file...
  474         OptionsInfo["Psi4RunParams"]["OutputFile"] = Psi4Util.UpdatePsi4OutputFileUsingPID(OptionsInfo["Psi4RunParams"]["OutputFile"], os.getpid())
  475             
  476     # Intialize Psi4...
  477     OptionsInfo["psi4"] = Psi4Util.InitializePsi4(Psi4RunParams = OptionsInfo["Psi4RunParams"], Psi4OptionsParams = OptionsInfo["Psi4OptionsParams"], PrintVersion = False, PrintHeader = True)
  478     
  479     # Setup max iterations global variable...
  480     Psi4Util.UpdatePsi4OptionsParameters(OptionsInfo["psi4"], {'GEOM_MAXITER': OptionsInfo["MaxIters"]})
  481     
  482     # Setup conversion factor for energy units...
  483     SetupEnergyConversionFactor(OptionsInfo["psi4"])
  484 
  485 def PerformMinimizationAndTorsionScan(Mol, MolNum, UseMultiProcessingAtTorsionAnglesLevel = False):
  486     """Peform minimization and torsions scan."""
  487 
  488     if not OptionsInfo["Infile3D"]:
  489         # Add hydrogens...
  490         Mol = Chem.AddHs(Mol, addCoords = True)
  491 
  492         Mol, MinimizationCalcStatus = MinimizeMolecule(Mol, MolNum)
  493         if not MinimizationCalcStatus:
  494             return (Mol, False, False, False)
  495         
  496     TorsionsMolInfo = SetupTorsionsMolInfo(Mol, MolNum)
  497     if TorsionsMolInfo["NumOfMatches"] == 0:
  498         return (Mol, True, False, False)
  499     
  500     Mol, ScanCalcStatus = ScanAllTorsionsInMol(Mol, TorsionsMolInfo, MolNum, UseMultiProcessingAtTorsionAnglesLevel)
  501     if not ScanCalcStatus:
  502         return (Mol, True, True, False)
  503     
  504     return (Mol, True, True, True)
  505 
  506 def ScanAllTorsionsInMol(Mol, TorsionsMolInfo,  MolNum, UseMultiProcessingAtTorsionAnglesLevel = False):
  507     """Peform scans on all torsions in a molecule."""
  508 
  509     if TorsionsMolInfo["NumOfMatches"] == 0:
  510         return Mol, True
  511 
  512     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  513     
  514     FirstTorsionMode = OptionsInfo["FirstTorsionMode"]
  515     TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
  516 
  517     TorsionPatternCount, TorsionScanCount, TorsionMatchCount = [0] * 3
  518     TorsionMaxMatches = OptionsInfo["TorsionMaxMatches"]
  519     
  520     for TorsionID in TorsionsPatternsInfo["IDs"]:
  521         TorsionPatternCount +=  1
  522         TorsionPattern = TorsionsPatternsInfo["Pattern"][TorsionID]
  523         TorsionPatternMol = TorsionsPatternsInfo["Mol"][TorsionID]
  524         
  525         TorsionsMatches = TorsionsMolInfo["Matches"][TorsionID]
  526 
  527         if TorsionsMatches is None:
  528             continue
  529         
  530         if FirstTorsionMode and TorsionPatternCount > 1:
  531             if not OptionsInfo["QuietMode"]:
  532                 MiscUtil.PrintWarning("Already scaned first torsion pattern, \"%s\" for molecule %s during \"%s\" value of \"--modeTorsions\" option . Abandoning torsion scan...\n" % (TorsionPattern, MolName, OptionsInfo["ModeTorsions"]))
  533             break
  534 
  535         for Index, TorsionMatches in enumerate(TorsionsMatches):
  536             TorsionMatchNum = Index + 1
  537             TorsionMatchCount +=  1
  538 
  539             if TorsionMatchCount > TorsionMaxMatches:
  540                 if not OptionsInfo["QuietMode"]:
  541                     MiscUtil.PrintWarning("Already scaned a maximum of %s torsion matches for molecule %s specified by \"--torsionMaxMatches\" option. Abandoning torsion scan...\n" % (TorsionMaxMatches, MolName))
  542                 break
  543 
  544             TmpMol, TorsionScanStatus,  TorsionMols, TorsionEnergies, TorsionAngles = ScanSingleTorsionInMol(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches,  TorsionMatchNum, MolNum, UseMultiProcessingAtTorsionAnglesLevel)
  545             if not TorsionScanStatus:
  546                 continue
  547             
  548             TorsionScanCount +=  1
  549             GenerateOutputFiles(Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
  550         
  551         if TorsionMatchCount > TorsionMaxMatches:
  552             break
  553     
  554     if TorsionScanCount:
  555         GenerateStartingTorsionScanStructureOutfile(Mol, MolNum)
  556 
  557     Status = True if TorsionScanCount else False
  558     
  559     return (Mol, Status)
  560 
  561 def ScanSingleTorsionInMol(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, TorsionMatchNum, MolNum, UseMultiProcessingAtTorsionAnglesLevel):
  562     """Perform torsion scan for a molecule along with constrained energy minimization."""
  563     
  564     if not OptionsInfo["QuietMode"]:
  565         MiscUtil.PrintInfo("\nProcessing torsion pattern, %s, match number, %s, in molecule %s..." % (TorsionPattern, TorsionMatchNum, RDKitUtil.GetMolName(Mol, MolNum)))
  566         
  567         if OptionsInfo["TorsionMinimize"]:
  568             MiscUtil.PrintInfo("Generating initial ensemble of constrained conformations by distance geometry and forcefield followed by Psi4 constraned minimization to select the lowest energy structure at specific torsion angles for molecule %s..." % (RDKitUtil.GetMolName(Mol, MolNum)))
  569         else:
  570             MiscUtil.PrintInfo("Calculating single point energy using Psi4 for molecule, %s, at specific torsion angles..." % (RDKitUtil.GetMolName(Mol, MolNum)))
  571     
  572     if UseMultiProcessingAtTorsionAnglesLevel:
  573         return ScanSingleTorsionInMolUsingMultipleProcessesAtTorsionAnglesLevel(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum)
  574     else:
  575         return ScanSingleTorsionInMolUsingSingleProcess(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum)
  576 
  577 def ScanSingleTorsionInMolUsingSingleProcess(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum):
  578     """Perform torsion scan for a molecule using single processs along with constrained
  579     energy minimization."""
  580 
  581     TorsionMols = []
  582     TorsionEnergies = []
  583     TorsionAngles = []
  584 
  585     Mols, Angles = SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum)
  586     
  587     for Index, Angle in enumerate(Angles):
  588         TorsionMol = Mols[Index]
  589         TorsionMol, CalcStatus, Energy = MinimizeCalculateEnergyForTorsionMol(Mol, TorsionMol, Angle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum)
  590         
  591         if not CalcStatus:
  592             return (Mol, False, None, None, None)
  593         
  594         TorsionMols.append(TorsionMol)
  595         TorsionEnergies.append(Energy)
  596         TorsionAngles.append(Angle)
  597     
  598     return (Mol, True, TorsionMols, TorsionEnergies, TorsionAngles)
  599 
  600 def SetupMolsForSingleTorsionScanInMol(Mol, TorsionMatches, MolNum = None):
  601     """Setup molecules corresponding to all torsion angles in a molecule."""
  602 
  603     StartAngle = OptionsInfo["TorsionStart"]
  604     StopAngle = OptionsInfo["TorsionStop"]
  605     StepSize = OptionsInfo["TorsionStep"]
  606     
  607     AtomIndex1, AtomIndex2, AtomIndex3, AtomIndex4 = TorsionMatches
  608 
  609     TorsionMols = []
  610     
  611     TorsionAngles = [Angle for Angle in range(StartAngle, StopAngle, StepSize)]
  612     TorsionAngles.append(StopAngle)
  613 
  614     for Angle in TorsionAngles:
  615         TorsionMol = Chem.Mol(Mol)
  616         TorsionMolConf = TorsionMol.GetConformer(0)
  617         
  618         rdMolTransforms.SetDihedralDeg(TorsionMolConf, AtomIndex1, AtomIndex2, AtomIndex3, AtomIndex4, Angle)
  619         TorsionMols.append(TorsionMol)
  620 
  621     return (TorsionMols, TorsionAngles)
  622 
  623 def MinimizeCalculateEnergyForTorsionMol(Mol, TorsionMol, TorsionAngle, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum):
  624     """"Calculate energy of a torsion molecule by performing an optional constrained
  625     energy minimzation.
  626     """
  627 
  628     if OptionsInfo["TorsionMinimize"]:
  629         if not OptionsInfo["QuietMode"]:
  630             MolName = RDKitUtil.GetMolName(Mol, MolNum)
  631             MiscUtil.PrintInfo("\nProcessing torsion angle %s for molecule %s..." % (TorsionAngle, MolName))
  632         
  633         # Perform constrained minimization...
  634         TorsionMatchesMol = RDKitUtil.MolFromSubstructureMatch(TorsionMol, TorsionPatternMol, TorsionMatches)
  635         TorsionMol, CalcStatus, Energy = ConstrainAndMinimizeMolecule(TorsionMol, TorsionAngle, TorsionMatchesMol, TorsionMatches, MolNum)
  636         
  637         if not CalcStatus:
  638             if not OptionsInfo["QuietMode"]:
  639                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  640                 MiscUtil.PrintWarning("Failed to perform constrained minimization for molecule %s with torsion angle set to %s during torsion scan for torsion pattern %s. Abandoning torsion scan..." % (MolName, TorsionAngle, TorsionPattern))
  641             return (TorsionMol, False, None)
  642     else:
  643         # Setup a Psi4 molecule...
  644         Psi4Mol = SetupPsi4Mol(OptionsInfo["psi4"], TorsionMol, MolNum)
  645         if Psi4Mol is None:
  646             return (TorsionMol, False, None)
  647             
  648         # Calculate single point Psi4 energy...
  649         CalcStatus, Energy = CalculateEnergyUsingPsi4(OptionsInfo["psi4"], Psi4Mol, TorsionMol, MolNum)
  650         if not CalcStatus:
  651             if not OptionsInfo["QuietMode"]:
  652                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  653                 MiscUtil.PrintWarning("Failed to calculate Psi4 energy for molecule %s with torsion angle set to %s during torsion scan for torsion pattern %s. Abandoning torsion scan..." % (MolName, TorsionAngle, TorsionPattern))
  654             return (TorsionMol, False, None)
  655         
  656     return (TorsionMol, CalcStatus, Energy)
  657     
  658 def SetupTorsionsMolInfo(Mol, MolNum = None):
  659     """Setup torsions info for a molecule."""
  660 
  661     TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
  662     
  663     # Initialize...
  664     TorsionsMolInfo = {}
  665     TorsionsMolInfo["IDs"] = []
  666     TorsionsMolInfo["NumOfMatches"] = 0
  667     TorsionsMolInfo["Matches"] = {}
  668     for TorsionID in TorsionsPatternsInfo["IDs"]:
  669         TorsionsMolInfo["IDs"].append(TorsionID)
  670         TorsionsMolInfo["Matches"][TorsionID] = None
  671     
  672     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  673     UseChirality = OptionsInfo["UseChirality"]
  674     
  675     for TorsionID in TorsionsPatternsInfo["IDs"]:
  676         # Match torsions..
  677         TorsionPattern = TorsionsPatternsInfo["Pattern"][TorsionID]
  678         TorsionPatternMol = TorsionsPatternsInfo["Mol"][TorsionID]
  679         TorsionsMatches = RDKitUtil.FilterSubstructureMatchesByAtomMapNumbers(Mol, TorsionPatternMol, Mol.GetSubstructMatches(TorsionPatternMol, useChirality = UseChirality))
  680         
  681         # Validate tosion matches...
  682         ValidTorsionsMatches = []
  683         for Index, TorsionMatch in enumerate(TorsionsMatches):
  684             if len(TorsionMatch) != 4:
  685                 if not OptionsInfo["QuietMode"]:
  686                     MiscUtil.PrintWarning("Ignoring invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: It must match exactly 4 atoms." % (TorsionMatch, TorsionPattern, MolName))
  687                 continue
  688 
  689             if not RDKitUtil.AreAtomIndicesSequentiallyConnected(Mol, TorsionMatch):
  690                 if not OptionsInfo["QuietMode"]:
  691                     MiscUtil.PrintWarning("Ignoring invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices must be sequentially connected." % (TorsionMatch, TorsionPattern, MolName))
  692                 continue
  693 
  694             Bond = Mol.GetBondBetweenAtoms(TorsionMatch[1], TorsionMatch[2])
  695             if Bond.IsInRing():
  696                 if not OptionsInfo["QuietMode"]:
  697                     MiscUtil.PrintWarning("Ignoring invalid torsion match to atom indices, %s, for torsion pattern, %s, in molecule %s: Matched atom indices, %s and %s, are not allowed to be in a ring." % (TorsionMatch, TorsionPattern, MolName, TorsionMatch[1], TorsionMatch[2]))
  698                 continue
  699             
  700             ValidTorsionsMatches.append(TorsionMatch)
  701         
  702         # Track valid matches...
  703         if len(ValidTorsionsMatches):
  704             TorsionsMolInfo["NumOfMatches"] += len(ValidTorsionsMatches)
  705             TorsionsMolInfo["Matches"][TorsionID] = ValidTorsionsMatches
  706         
  707     if TorsionsMolInfo["NumOfMatches"] == 0:
  708         if not OptionsInfo["QuietMode"]:
  709             MiscUtil.PrintWarning("Failed to match any torsions  in molecule %s" % (MolName))
  710 
  711     return TorsionsMolInfo
  712 
  713 def SetupTorsionsPatternsInfo():
  714     """Setup torsions patterns info."""
  715 
  716     TorsionsPatternsInfo = {}
  717     TorsionsPatternsInfo["IDs"] = []
  718     TorsionsPatternsInfo["Pattern"] = {}
  719     TorsionsPatternsInfo["Mol"] = {}
  720 
  721     TorsionID = 0
  722     for TorsionPattern in OptionsInfo["TorsionPatternsList"]:
  723         TorsionID += 1
  724         
  725         TorsionMol = Chem.MolFromSmarts(TorsionPattern)
  726         if TorsionMol is None:
  727             MiscUtil.PrintError("Failed to create torsion pattern molecule. The torsion SMILES/SMARTS pattern, \"%s\", specified using \"-t, --torsions\" option is not valid." % (TorsionPattern))
  728         
  729         TorsionsPatternsInfo["IDs"].append(TorsionID)
  730         TorsionsPatternsInfo["Pattern"][TorsionID] = TorsionPattern
  731         TorsionsPatternsInfo["Mol"][TorsionID] = TorsionMol
  732 
  733     OptionsInfo["TorsionsPatternsInfo"] = TorsionsPatternsInfo
  734     
  735 def MinimizeMolecule(Mol, MolNum = None):
  736     """Minimize molecule."""
  737 
  738     return GenerateAndMinimizeConformersUsingForceField(Mol, MolNum)
  739 
  740 def GenerateAndMinimizeConformersUsingForceField(Mol, MolNum = None):
  741     """Generate and minimize conformers for a molecule to get the lowest energy conformer
  742     as the minimized structure."""
  743 
  744     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  745     
  746     # Setup conformers...
  747     ConfIDs = EmbedMolecule(Mol, MolNum)
  748     if not len(ConfIDs):
  749         if not OptionsInfo["QuietMode"]:
  750             MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s: Embedding failed...\n" % MolName)
  751         return (Mol, False)
  752     
  753     if not OptionsInfo["QuietMode"]:
  754         MiscUtil.PrintInfo("Performing initial minimization of molecule, %s, using forcefield by generating a conformation ensemble and selecting the lowest energy conformer - EmbedRMSDCutoff: %s; Size: %s; Size after RMSD filtering: %s" % (MolName, OptionsInfo["ConfGenerationParams"]["EmbedRMSDCutoff"], OptionsInfo["ConfGenerationParams"]["MaxConfs"], len(ConfIDs)))
  755     
  756     # Minimize conformers...
  757     CalcEnergyMap = {}
  758     for ConfID in ConfIDs:
  759         # Perform forcefield minimization...
  760         Status, ConvergeStatus = MinimizeMoleculeUsingForceField(Mol, MolNum, ConfID)
  761         if not Status:
  762             return (Mol, False)
  763         
  764         EnergyStatus, Energy = CalculateEnergyUsingForceField(Mol, ConfID)
  765         if not EnergyStatus:
  766             if not OptionsInfo["QuietMode"]:
  767                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  768                 MiscUtil.PrintWarning("Failed to retrieve calculated energy for conformation number %d of molecule %s. Try again after removing any salts or cleaing up the molecule...\n" % (ConfID, MolName))
  769             return (Mol, False)
  770         
  771         if ConvergeStatus != 0:
  772             if not OptionsInfo["QuietMode"]:
  773                 MiscUtil.PrintWarning("Minimization using forcefield failed to converge for molecule %s in %d steps. Try using higher value for \"maxIters\" in \"--confParams\" option...\n" % (MolName, OptionsInfo["ConfGenerationParams"]["MaxIters"]))
  774         
  775         CalcEnergyMap[ConfID] = Energy
  776     
  777     SortedConfIDs = sorted(ConfIDs, key = lambda ConfID: CalcEnergyMap[ConfID])
  778     MinEnergyConfID = SortedConfIDs[0]
  779         
  780     for ConfID in [Conf.GetId() for Conf in Mol.GetConformers()]:
  781         if ConfID == MinEnergyConfID:
  782             continue
  783         Mol.RemoveConformer(ConfID)
  784     
  785     # Set ConfID to 0 for MinEnergyConf...
  786     Mol.GetConformer(MinEnergyConfID).SetId(0)
  787 
  788     return (Mol, True)
  789 
  790 def ConstrainAndMinimizeMolecule(Mol, TorsionAngle, RefMolCore, RefMolMatches, MolNum = None):
  791     "Constrain and Minimize molecule."
  792 
  793     # TorsionMol, CalcStatus, Energy
  794     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  795 
  796     # Setup constrained conformers...
  797     MolConfs, MolConfsStatus = ConstrainEmbedAndMinimizeMoleculeUsingRDKit(Mol, RefMolCore, RefMolMatches, MolNum)
  798     if not MolConfsStatus:
  799         if not OptionsInfo["QuietMode"]:
  800             MiscUtil.PrintWarning("Conformation generation couldn't be performed for molecule %s: Constrained embedding failed...\n" % MolName)
  801         return (Mol, False, None)
  802 
  803     # Minimize conformers...
  804     ConfNums = []
  805     CalcEnergyMap = {}
  806     MolConfsMap = {}
  807     
  808     for ConfNum, MolConf in enumerate(MolConfs):
  809         if not OptionsInfo["QuietMode"]:
  810             MiscUtil.PrintInfo("\nPerforming constrained minimization using Psi4 for molecule, %s, conformer number, %s, at torsion angle %s..." % (MolName, ConfNum, TorsionAngle))
  811 
  812         CalcStatus, Energy = ConstrainAndMinimizeMoleculeUsingPsi4(OptionsInfo["psi4"], MolConf, RefMolCore, RefMolMatches, MolNum)
  813         if not CalcStatus:
  814             if not OptionsInfo["QuietMode"]:
  815                 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s\n" % (MolName))
  816             return (Mol, False, None)
  817         
  818         ConfNums.append(ConfNum)
  819         CalcEnergyMap[ConfNum] = Energy
  820         MolConfsMap[ConfNum] = MolConf
  821     
  822     SortedConfNums = sorted(ConfNums, key = lambda ConfNum: CalcEnergyMap[ConfNum])
  823     MinEnergyConfNum = SortedConfNums[0]
  824     
  825     MinEnergy = CalcEnergyMap[MinEnergyConfNum]
  826     MinEnergyMolConf = MolConfsMap[MinEnergyConfNum]
  827     
  828     MinEnergyMolConf.ClearProp('EmbedRMS')
  829     
  830     return (MinEnergyMolConf, True, MinEnergy)
  831 
  832 def ConstrainAndMinimizeMoleculeUsingPsi4(Psi4Handle, Mol, RefMolCore, RefMolMatches, MolNum, ConfID = -1):
  833     "Minimize molecule using Psi4."
  834 
  835     # Setup a list for constrained atoms...
  836     ConstrainedAtomIndices = SetupConstrainedAtomIndicesForPsi4(Mol, RefMolCore, RefMolMatches)
  837     if len(ConstrainedAtomIndices) == 0:
  838         return (False, None)
  839     
  840     # Setup a Psi4Mol...
  841     Psi4Mol = SetupPsi4Mol(Psi4Handle, Mol, MolNum, ConfID)
  842     if Psi4Mol is None:
  843         return (False, None)
  844         
  845     #  Setup reference wave function...
  846     Reference = SetupReferenceWavefunction(Mol)
  847     Psi4Handle.set_options({'Reference': Reference})
  848     
  849     # Setup method name and basis set...
  850     MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
  851 
  852     # Setup freeze list for constrained atoms...
  853     FreezeXYZList = [("%s xyz" % AtomIdex) for AtomIdex in ConstrainedAtomIndices]
  854     FreezeXYZString = " %s " % " ".join(FreezeXYZList)
  855     Psi4Handle.set_module_options("OPTKING", {"frozen_cartesian": FreezeXYZString})
  856     
  857     # Optimize geometry...
  858     Status, Energy, WaveFunction = Psi4Util.PerformGeometryOptimization(Psi4Handle, Psi4Mol, MethodName, BasisSet, ReturnWaveFunction = True, Quiet = OptionsInfo["QuietMode"])
  859     
  860     if not Status:
  861         PerformPsi4Cleanup(Psi4Handle)
  862         return (False, None)
  863 
  864     # Update atom positions...
  865     AtomPositions = Psi4Util.GetAtomPositions(Psi4Handle, WaveFunction, InAngstroms = True)
  866     RDKitUtil.SetAtomPositions(Mol, AtomPositions, ConfID = ConfID)
  867 
  868     # Convert energy units...
  869     if OptionsInfo["ApplyEnergyConversionFactor"]:
  870         Energy = Energy * OptionsInfo["EnergyConversionFactor"]
  871     
  872     # Clean up
  873     PerformPsi4Cleanup(Psi4Handle)
  874     
  875     return (True, Energy)
  876     
  877 def ConstrainEmbedAndMinimizeMoleculeUsingRDKit(Mol, RefMolCore, RefMolMatches, MolNum = None):
  878     "Constrain, embed, and minimize molecule."
  879 
  880     # Setup forcefield function to use for constrained minimization...
  881     ForceFieldFunction = None
  882     ForceFieldName = None
  883     if OptionsInfo["ConfGenerationParams"]["UseUFF"]:
  884         ForceFieldFunction = lambda mol, confId = -1 : AllChem.UFFGetMoleculeForceField(mol, confId = confId)
  885         ForeceFieldName = "UFF"
  886     else:
  887         ForceFieldFunction = lambda mol, confId = -1 : AllChem.MMFFGetMoleculeForceField(mol, AllChem.MMFFGetMoleculeProperties(mol, mmffVariant = OptionsInfo["ConfGenerationParams"]["MMFFVariant"]) , confId = confId)
  888         ForeceFieldName = "MMFF"
  889 
  890     if ForceFieldFunction is None:
  891         if not OptionsInfo["QuietMode"]:
  892             MiscUtil.PrintWarning("Failed to setup forcefield %s for molecule: %s\n" % (ForceFieldName, RDKitUtil.GetMolName(Mol, MolNum)))
  893         return (None, False)
  894     
  895     MaxConfs = OptionsInfo["ConfGenerationParams"]["MaxConfsTorsions"]
  896     EnforceChirality = OptionsInfo["ConfGenerationParams"]["EnforceChirality"]
  897     UseExpTorsionAnglePrefs = OptionsInfo["ConfGenerationParams"]["UseExpTorsionAnglePrefs"]
  898     UseBasicKnowledge = OptionsInfo["ConfGenerationParams"]["UseBasicKnowledge"]
  899     UseTethers = OptionsInfo["ConfGenerationParams"]["UseTethers"]
  900 
  901     MolConfs = []
  902     ConfIDs = [ConfID for ConfID in range(0, MaxConfs)]
  903     
  904     for ConfID in ConfIDs:
  905         try:
  906             MolConf = Chem.Mol(Mol)
  907             RDKitUtil.ConstrainAndEmbed(MolConf, RefMolCore, coreMatchesMol = RefMolMatches, useTethers = UseTethers, coreConfId = -1, randomseed = ConfID, getForceField = ForceFieldFunction, enforceChirality = EnforceChirality, useExpTorsionAnglePrefs = UseExpTorsionAnglePrefs, useBasicKnowledge = UseBasicKnowledge)
  908         except (ValueError, RuntimeError, Chem.rdchem.KekulizeException)  as ErrMsg:
  909             if not OptionsInfo["QuietMode"]:
  910                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  911                 MiscUtil.PrintWarning("Constrained embedding couldn't  be performed for molecule %s:\n%s\n" % (RDKitUtil.GetMolName(Mol, MolNum), ErrMsg))
  912             return (None, False)
  913         
  914         MolConfs.append(MolConf)
  915 
  916     return FilterConstrainedConformationsByRMSD(Mol, MolConfs, MolNum)
  917 
  918 def FilterConstrainedConformationsByRMSD(Mol, MolConfs, MolNum = None):
  919     """Filter contarained conformations by RMSD."""
  920     
  921     EmbedRMSDCutoff = OptionsInfo["ConfGenerationParams"]["EmbedRMSDCutoff"]
  922     if EmbedRMSDCutoff is None or EmbedRMSDCutoff <= 0:
  923         if not OptionsInfo["QuietMode"]:
  924             MiscUtil.PrintInfo("\nGenerating initial ensemble of  constrained conformations by distance geometry  and forcefield for %s - EmbedRMSDCutoff: None; Size: %s" % (RDKitUtil.GetMolName(Mol, MolNum), len(MolConfs)))
  925         return (MolConfs, True)
  926 
  927     FirstMolConf = True
  928     SelectedMolConfs = []
  929     for MolConfIndex, MolConf in enumerate(MolConfs):
  930         if FirstMolConf:
  931             FirstMolConf = False
  932             SelectedMolConfs.append(MolConf)
  933             continue
  934         
  935         # Compare RMSD against all selected conformers...
  936         ProbeMolConf = Chem.RemoveHs(Chem.Mol(MolConf))
  937         IgnoreConf = False
  938         for SelectedMolConfIndex, SelectedMolConf in enumerate(SelectedMolConfs):
  939             RefMolConf = Chem.RemoveHs(Chem.Mol(SelectedMolConf))
  940             CalcRMSD = rdMolAlign.AlignMol(ProbeMolConf, RefMolConf)
  941             
  942             if CalcRMSD < EmbedRMSDCutoff:
  943                 IgnoreConf = True
  944                 break
  945         
  946         if IgnoreConf:
  947             continue
  948         
  949         SelectedMolConfs.append(MolConf)
  950         
  951     if not OptionsInfo["QuietMode"]:
  952         MiscUtil.PrintInfo("\nGenerating initial ensemble of constrained conformations by distance geometry and forcefield for %s - EmbedRMSDCutoff: %s; Size: %s; Size after RMSD filtering: %s" % (RDKitUtil.GetMolName(Mol, MolNum), EmbedRMSDCutoff, len(MolConfs), len(SelectedMolConfs)))
  953     
  954     return (SelectedMolConfs, True)
  955 
  956 def EmbedMolecule(Mol, MolNum = None):
  957     "Embed conformations"
  958 
  959     ConfIDs = []
  960     
  961     MaxConfs = OptionsInfo["ConfGenerationParams"]["MaxConfs"]
  962     RandomSeed = OptionsInfo["ConfGenerationParams"]["RandomSeed"]
  963     EnforceChirality = OptionsInfo["ConfGenerationParams"]["EnforceChirality"]
  964     UseExpTorsionAnglePrefs = OptionsInfo["ConfGenerationParams"]["UseExpTorsionAnglePrefs"]
  965     UseBasicKnowledge = OptionsInfo["ConfGenerationParams"]["UseBasicKnowledge"]
  966     EmbedRMSDCutoff = OptionsInfo["ConfGenerationParams"]["EmbedRMSDCutoff"]
  967 
  968     try:
  969         ConfIDs = AllChem.EmbedMultipleConfs(Mol, numConfs = MaxConfs, randomSeed = RandomSeed, pruneRmsThresh = EmbedRMSDCutoff, enforceChirality = EnforceChirality, useExpTorsionAnglePrefs = UseExpTorsionAnglePrefs, useBasicKnowledge = UseBasicKnowledge)
  970     except ValueError as ErrMsg:
  971         if not OptionsInfo["QuietMode"]:
  972             MolName = RDKitUtil.GetMolName(Mol, MolNum)
  973             MiscUtil.PrintWarning("Embedding failed  for molecule %s:\n%s\n" % (MolName, ErrMsg))
  974         ConfIDs = []
  975 
  976     return ConfIDs
  977 
  978 def MinimizeMoleculeUsingForceField(Mol, MolNum, ConfID = -1):
  979     "Minimize molecule using forcefield available in RDKit."
  980     
  981     try:
  982         if OptionsInfo["ConfGenerationParams"]["UseUFF"]:
  983             ConvergeStatus = AllChem.UFFOptimizeMolecule(Mol, confId = ConfID, maxIters = OptionsInfo["ConfGenerationParams"]["MaxIters"])
  984         elif OptionsInfo["ConfGenerationParams"]["UseMMFF"]:
  985             ConvergeStatus = AllChem.MMFFOptimizeMolecule(Mol, confId = ConfID, maxIters = OptionsInfo["ConfGenerationParams"]["MaxIters"], mmffVariant = OptionsInfo["ConfGenerationParams"]["MMFFVariant"])
  986         else:
  987             MiscUtil.PrintError("Minimization couldn't be performed: Specified forcefield, %s, is not supported" % OptionsInfo["ConfGenerationParams"]["ForceField"])
  988     except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
  989         if not OptionsInfo["QuietMode"]:
  990             MolName = RDKitUtil.GetMolName(Mol, MolNum)
  991             MiscUtil.PrintWarning("Minimization using forcefield couldn't be performed for molecule %s:\n%s\n" % (MolName, ErrMsg))
  992         return (False, None)
  993     
  994     return (True, ConvergeStatus)
  995 
  996 def CalculateEnergyUsingForceField(Mol, ConfID = None):
  997     "Calculate energy."
  998 
  999     Status = True
 1000     Energy = None
 1001 
 1002     if ConfID is None:
 1003         ConfID = -1
 1004     
 1005     if OptionsInfo["ConfGenerationParams"]["UseUFF"]:
 1006         UFFMoleculeForcefield = AllChem.UFFGetMoleculeForceField(Mol, confId = ConfID)
 1007         if UFFMoleculeForcefield is None:
 1008             Status = False
 1009         else:
 1010             Energy = UFFMoleculeForcefield.CalcEnergy()
 1011     elif OptionsInfo["ConfGenerationParams"]["UseMMFF"]:
 1012         MMFFMoleculeProperties = AllChem.MMFFGetMoleculeProperties(Mol, mmffVariant = OptionsInfo["ConfGenerationParams"]["MMFFVariant"])
 1013         MMFFMoleculeForcefield = AllChem.MMFFGetMoleculeForceField(Mol, MMFFMoleculeProperties, confId = ConfID)
 1014         if MMFFMoleculeForcefield is None:
 1015             Status = False
 1016         else:
 1017             Energy = MMFFMoleculeForcefield.CalcEnergy()
 1018     else:
 1019         MiscUtil.PrintError("Couldn't retrieve conformer energy: Specified forcefield, %s, is not supported" % OptionsInfo["ConfGenerationParams"]["ForceField"])
 1020     
 1021     return (Status, Energy)
 1022 
 1023 def CalculateEnergyUsingPsi4(Psi4Handle, Psi4Mol, Mol, MolNum = None):
 1024     "Calculate single point energy using Psi4."
 1025     
 1026     Status = False
 1027     Energy = None
 1028     
 1029     #  Setup reference wave function...
 1030     Reference = SetupReferenceWavefunction(Mol)
 1031     Psi4Handle.set_options({'Reference': Reference})
 1032     
 1033     # Setup method name and basis set...
 1034     MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
 1035     
 1036     Status, Energy = Psi4Util.CalculateSinglePointEnergy(Psi4Handle, Psi4Mol, MethodName, BasisSet, Quiet = OptionsInfo["QuietMode"])
 1037     
 1038     # Convert energy units...
 1039     if Status:
 1040         if OptionsInfo["ApplyEnergyConversionFactor"]:
 1041             Energy = Energy * OptionsInfo["EnergyConversionFactor"]
 1042 
 1043     # Clean up
 1044     PerformPsi4Cleanup(Psi4Handle)
 1045     
 1046     return (Status, Energy)
 1047 
 1048 def SetupConstrainedAtomIndicesForPsi4(Mol, RefMolCore, RefMolMatches, ConstrainHydrogens = False):
 1049     """Setup a list of atom indices to be constrained during Psi4 minimizaiton."""
 1050 
 1051     AtomIndices = []
 1052 
 1053     if RefMolMatches is None:
 1054         ConstrainAtomIndices = Mol.GetSubstructMatch(RefMolCore)
 1055     else:
 1056         ConstrainAtomIndices = RefMolMatches
 1057         
 1058     # Collect matched heavy atoms along with attached hydrogens...
 1059     for AtomIndex in ConstrainAtomIndices:
 1060         Atom = Mol.GetAtomWithIdx(AtomIndex)
 1061         if Atom.GetAtomicNum() == 1:
 1062             continue
 1063         
 1064         AtomIndices.append(AtomIndex)
 1065         for AtomNbr in Atom.GetNeighbors():
 1066             if AtomNbr.GetAtomicNum() == 1:
 1067                 if ConstrainHydrogens:
 1068                     AtomNbrIndex = AtomNbr.GetIdx()
 1069                     AtomIndices.append(AtomNbrIndex)
 1070     
 1071     AtomIndices = sorted(AtomIndices)
 1072 
 1073     # Atom indices start from 1 for Psi4 instead 0 for RDKit...
 1074     AtomIndices = [ AtomIndex + 1 for AtomIndex in AtomIndices]
 1075     
 1076     return AtomIndices
 1077     
 1078 def SetupPsi4Mol(Psi4Handle, Mol, MolNum, ConfID = -1):
 1079     """Setup a Psi4 molecule object."""
 1080 
 1081     # Turn off recentering and reorientation to perform optimization in the
 1082     # constrained coordinate space...
 1083     MolGeometry = RDKitUtil.GetPsi4XYZFormatString(Mol, ConfID = ConfID, NoCom = True, NoReorient = True)
 1084 
 1085     try:
 1086         Psi4Mol = Psi4Handle.geometry(MolGeometry)
 1087     except Exception as ErrMsg:
 1088         Psi4Mol = None
 1089         if not OptionsInfo["QuietMode"]:
 1090             MiscUtil.PrintWarning("Failed to create Psi4 molecule from geometry string: %s\n" % ErrMsg)
 1091             MolName = RDKitUtil.GetMolName(Mol, MolNum)
 1092             MiscUtil.PrintWarning("Ignoring molecule: %s" % MolName)
 1093 
 1094     return Psi4Mol
 1095 
 1096 def PerformPsi4Cleanup(Psi4Handle):
 1097     """Perform clean up. """
 1098 
 1099     # No need to perform any explicit clean by calling
 1100     # Psi4Handle.core.opt_clean(). It's already done by Psi4 after
 1101     # optimization...
 1102     
 1103     # Clean up any leftover scratch files...
 1104     if OptionsInfo["MPMode"]:
 1105         Psi4Util.RemoveScratchFiles(Psi4Handle, OptionsInfo["Psi4RunParams"]["OutputFile"])
 1106 
 1107 def CheckAndValidateMolecule(Mol, MolCount = None):
 1108     """Validate molecule for Psi4 calculations."""
 1109     
 1110     if Mol is None:
 1111         if not OptionsInfo["QuietMode"]:
 1112             MiscUtil.PrintInfo("\nProcessing molecule number %s..." % MolCount)
 1113         return False
 1114     
 1115     MolName = RDKitUtil.GetMolName(Mol, MolCount)
 1116     if not OptionsInfo["QuietMode"]:
 1117         MiscUtil.PrintInfo("\nProcessing molecule %s..." % MolName)
 1118     
 1119     if RDKitUtil.IsMolEmpty(Mol):
 1120         if not OptionsInfo["QuietMode"]:
 1121             MiscUtil.PrintWarning("Ignoring empty molecule: %s\n" % MolName)
 1122         return False
 1123     
 1124     if not RDKitUtil.ValidateElementSymbols(RDKitUtil.GetAtomSymbols(Mol)):
 1125         if not OptionsInfo["QuietMode"]:
 1126             MiscUtil.PrintWarning("Ignoring molecule containing invalid element symbols: %s\n" % MolName)
 1127         return False
 1128 
 1129     if OptionsInfo["Infile3D"]:
 1130         if not Mol.GetConformer().Is3D():
 1131             if not OptionsInfo["QuietMode"]:
 1132                 MiscUtil.PrintWarning("3D tag is not set for molecule: %s\n" % MolName)
 1133     
 1134     if OptionsInfo["Infile3D"]:
 1135         # Otherwise, Hydrogens are always added...
 1136         if RDKitUtil.AreHydrogensMissingInMolecule(Mol):
 1137             if not OptionsInfo["QuietMode"]:
 1138                 MiscUtil.PrintWarning("Missing hydrogens in molecule: %s\n" % MolName)
 1139         
 1140     return True
 1141 
 1142 def SetupMethodNameAndBasisSet(Mol):
 1143     """Setup method name and basis set. """
 1144     
 1145     MethodName = OptionsInfo["MethodName"]
 1146     if OptionsInfo["MethodNameAuto"]:
 1147         MethodName = "B3LYP"
 1148     
 1149     BasisSet = OptionsInfo["BasisSet"]
 1150     if OptionsInfo["BasisSetAuto"]:
 1151         BasisSet = "6-31+G**" if RDKitUtil.IsAtomSymbolPresentInMol(Mol, "S") else "6-31G**"
 1152     
 1153     return (MethodName, BasisSet)
 1154 
 1155 def SetupReferenceWavefunction(Mol):
 1156     """Setup method name and basis set. """
 1157     
 1158     Reference = OptionsInfo["Reference"]
 1159     if OptionsInfo["ReferenceAuto"]:
 1160         Reference = 'UHF' if (RDKitUtil.GetSpinMultiplicity(Mol) > 1) else 'RHF'
 1161     
 1162     return Reference
 1163 
 1164 def SetupEnergyConversionFactor(Psi4Handle):
 1165     """Setup converstion factor for energt units. The Psi4 energy units are Hartrees."""
 1166     
 1167     EnergyUnits = OptionsInfo["EnergyUnits"]
 1168     
 1169     ApplyConversionFactor = True
 1170     if re.match("^kcal\/mol$", EnergyUnits, re.I):
 1171         ConversionFactor = Psi4Handle.constants.hartree2kcalmol
 1172     elif re.match("^kJ\/mol$", EnergyUnits, re.I):
 1173         ConversionFactor = Psi4Handle.constants.hartree2kJmol
 1174     elif re.match("^eV$", EnergyUnits, re.I):
 1175         ConversionFactor = Psi4Handle.constants.hartree2ev
 1176     else:
 1177         ApplyConversionFactor = False
 1178         ConversionFactor = 1.0
 1179     
 1180     OptionsInfo["ApplyEnergyConversionFactor"] = ApplyConversionFactor
 1181     OptionsInfo["EnergyConversionFactor"] = ConversionFactor
 1182 
 1183 def GenerateStartingTorsionScanStructureOutfile(Mol, MolNum):
 1184     """Write out the structure of molecule used for starting tosion scan. """
 1185     
 1186     FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
 1187     MolName = GetOutputFileMolName(Mol, MolNum)
 1188     
 1189     Outfile  = "%s_%s.%s" % (FileName, MolName, FileExt)
 1190     
 1191     # Set up a molecule writer...
 1192     Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
 1193     if Writer is None:
 1194         MiscUtil.PrintWarning("Failed to setup a writer for output fie %s " % Outfile)
 1195         return
 1196     
 1197     Writer.write(Mol)
 1198     
 1199     if Writer is not None:
 1200         Writer.close()
 1201 
 1202 def GenerateOutputFiles(Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
 1203     """Generate output files. """
 1204     
 1205     StructureOutfile, EnergyTextOutfile, PlotOutfile = SetupOutputFileNames(Mol, MolNum, TorsionID, TorsionMatchNum)
 1206     
 1207     GenerateScannedTorsionsStructureOutfile(StructureOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
 1208     GenerateEnergyTextOutfile(EnergyTextOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
 1209     GeneratePlotOutfile(PlotOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
 1210 
 1211 def GenerateScannedTorsionsStructureOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
 1212     """Write out structures generated after torsion scan along with associated data."""
 1213 
 1214     # Set up a molecule writer...
 1215     Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
 1216     if Writer is None:
 1217         MiscUtil.PrintWarning("Failed to setup a writer for output fie %s " % Outfile)
 1218         return
 1219     
 1220     MolName = RDKitUtil.GetMolName(Mol, MolNum)
 1221     
 1222     RelativeTorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
 1223     for Index, TorsionMol in enumerate(TorsionMols):
 1224         TorsionAngle = "%s" % TorsionAngles[Index]
 1225         TorsionMol.SetProp("Torsion_Angle", TorsionAngle)
 1226         
 1227         TorsionEnergy = "%.*f" % (OptionsInfo["Precision"], TorsionEnergies[Index])
 1228         TorsionMol.SetProp(OptionsInfo["EnergyDataFieldLabel"], TorsionEnergy)
 1229 
 1230         RelativeTorsionEnergy = "%.*f" % (OptionsInfo["Precision"], RelativeTorsionEnergies[Index])
 1231         TorsionMol.SetProp(OptionsInfo["EnergyRelativeDataFieldLabel"], RelativeTorsionEnergy)
 1232 
 1233         TorsionMolName = "%s_Deg%s" % (MolName, TorsionAngle)
 1234         TorsionMol.SetProp("_Name", TorsionMolName)
 1235         
 1236         Writer.write(TorsionMol)
 1237         
 1238     if Writer is not None:
 1239         Writer.close()
 1240     
 1241 def GenerateEnergyTextOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
 1242     """Write out torsion angles and energies."""
 1243 
 1244     # Setup a writer...
 1245     Writer = open(Outfile, "w")
 1246     if Writer is None:
 1247         MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
 1248 
 1249     # Write headers...
 1250     Writer.write("TorsionAngle,%s,%s\n" % (OptionsInfo["EnergyDataFieldLabel"], OptionsInfo["EnergyRelativeDataFieldLabel"]))
 1251 
 1252     RelativeTorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
 1253     for Index, TorsionAngle in enumerate(TorsionAngles):
 1254         TorsionEnergy = "%.*f" % (OptionsInfo["Precision"], TorsionEnergies[Index])
 1255         RelativeTorsionEnergy = "%.*f" % (OptionsInfo["Precision"], RelativeTorsionEnergies[Index])
 1256         Writer.write("%d,%s,%s\n" % (TorsionAngle, TorsionEnergy, RelativeTorsionEnergy))
 1257 
 1258     if Writer is not None:
 1259         Writer.close()
 1260     
 1261 def GeneratePlotOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
 1262     """Generate a plot corresponding to torsion angles and energies."""
 1263 
 1264     OutPlotParams = OptionsInfo["OutPlotParams"]
 1265 
 1266     # Initialize seaborn and matplotlib paramaters...
 1267     if not OptionsInfo["OutPlotInitialized"]:
 1268         OptionsInfo["OutPlotInitialized"] = True
 1269         RCParams = {"figure.figsize":(OutPlotParams["Width"], OutPlotParams["Height"]),
 1270                     "axes.titleweight": OutPlotParams["TitleWeight"],
 1271                     "axes.labelweight": OutPlotParams["LabelWeight"]}
 1272         sns.set(context = OutPlotParams["Context"], style = OutPlotParams["Style"], palette = OutPlotParams["Palette"], font = OutPlotParams["Font"], font_scale = OutPlotParams["FontScale"], rc = RCParams)
 1273 
 1274     # Create a new figure...
 1275     plt.figure()
 1276 
 1277     if OptionsInfo["OutPlotRelativeEnergy"]: 
 1278         TorsionEnergies = SetupRelativeEnergies(TorsionEnergies)
 1279     
 1280     # Draw plot...
 1281     PlotType = OutPlotParams["Type"]
 1282     if re.match("linepoint", PlotType, re.I):
 1283         Axis = sns.lineplot(x = TorsionAngles, y = TorsionEnergies, marker = "o",  legend = False)
 1284     elif re.match("scatter", PlotType, re.I):
 1285         Axis = sns.scatterplot(x = TorsionAngles, y = TorsionEnergies, legend = False)
 1286     elif re.match("line", PlotType, re.I):
 1287         Axis = sns.lineplot(x = TorsionAngles, y = TorsionEnergies, legend = False)
 1288     else:
 1289         MiscUtil.PrintError("The value, %s, specified for \"type\" using option \"--outPlotParams\" is not supported. Valid plot types: linepoint, scatter or line" % (PlotType))
 1290 
 1291     # Setup title and labels...
 1292     Title = OutPlotParams["Title"]
 1293     if OptionsInfo["OutPlotTitleTorsionSpec"]:
 1294         TorsionPattern = OptionsInfo["TorsionsPatternsInfo"]["Pattern"][TorsionID]
 1295         Title = "%s: %s" % (OutPlotParams["Title"], TorsionPattern)
 1296 
 1297     # Set labels and title...
 1298     Axis.set(xlabel = OutPlotParams["XLabel"], ylabel = OutPlotParams["YLabel"], title = Title)
 1299     
 1300     # Save figure...
 1301     plt.savefig(Outfile)
 1302 
 1303     # Close the plot...
 1304     plt.close()
 1305 
 1306 def SetupRelativeEnergies(Energies):
 1307     """Set up a list of relative energies. """
 1308     
 1309     SortedEnergies = sorted(Energies)
 1310     MinEnergy = SortedEnergies[0]
 1311     RelativeEnergies = [(Energy - MinEnergy) for Energy in Energies]
 1312 
 1313     return RelativeEnergies
 1314     
 1315 def SetupOutputFileNames(Mol, MolNum, TorsionID, TorsionMatchNum):
 1316     """Setup names of output files. """
 1317     
 1318     FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
 1319     MolName = GetOutputFileMolName(Mol, MolNum)
 1320     
 1321     OutfileRoot  = "%s_%s_Torsion%s_Match%s" % (FileName, MolName, TorsionID, TorsionMatchNum)
 1322     
 1323     StructureOutfile = "%s.%s" % (OutfileRoot, FileExt)
 1324     EnergyTextOutfile = "%s_Energies.csv" % (OutfileRoot)
 1325     PlotExt = OptionsInfo["OutPlotParams"]["OutExt"]
 1326     PlotOutfile = "%s_Plot.%s" % (OutfileRoot, PlotExt)
 1327 
 1328     return (StructureOutfile, EnergyTextOutfile, PlotOutfile)
 1329 
 1330 def GetOutputFileMolName(Mol, MolNum):
 1331     """Get output file prefix. """
 1332     
 1333     MolName = "Mol%s" % MolNum
 1334     if OptionsInfo["OutfileMolName"]:
 1335         MolName = re.sub("[^a-zA-Z0-9]", "_", RDKitUtil.GetMolName(Mol, MolNum), re.I)
 1336 
 1337     return MolName
 1338 
 1339 def ProcessOptions():
 1340     """Process and validate command line arguments and options."""
 1341     
 1342     MiscUtil.PrintInfo("Processing options...")
 1343 
 1344     # Validate options...
 1345     ValidateOptions()
 1346 
 1347     OptionsInfo["ModeMols"] = Options["--modeMols"]
 1348     OptionsInfo["FirstMolMode"] = True if re.match("^First$", Options["--modeMols"], re.I) else False
 1349     
 1350     OptionsInfo["ModeTorsions"] = Options["--modeTorsions"]
 1351     OptionsInfo["FirstTorsionMode"] = True if re.match("^First$", Options["--modeTorsions"], re.I) else False
 1352     
 1353     OptionsInfo["Infile"] = Options["--infile"]
 1354     OptionsInfo["SMILESInfileStatus"] = True if  MiscUtil.CheckFileExt(Options["--infile"], "smi csv tsv txt") else False
 1355     ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
 1356     OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"], InfileName = Options["--infile"], ParamsDefaultInfo = ParamsDefaultInfoOverride)
 1357     OptionsInfo["Infile3D"] = True if re.match("^yes$", Options["--infile3D"], re.I) else False
 1358     
 1359     OptionsInfo["Outfile"] = Options["--outfile"]
 1360     OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters("--outfileParams", Options["--outfileParams"])
 1361     
 1362     # Method, basis set, and reference wavefunction...
 1363     OptionsInfo["BasisSet"] = Options["--basisSet"]
 1364     OptionsInfo["BasisSetAuto"] = True if re.match("^auto$", Options["--basisSet"], re.I) else False
 1365     
 1366     OptionsInfo["MethodName"] = Options["--methodName"]
 1367     OptionsInfo["MethodNameAuto"] = True if re.match("^auto$", Options["--methodName"], re.I) else False
 1368     
 1369     OptionsInfo["Reference"] = Options["--reference"]
 1370     OptionsInfo["ReferenceAuto"] = True if re.match("^auto$", Options["--reference"], re.I) else False
 1371     
 1372     # Run and options parameters...
 1373     OptionsInfo["Psi4OptionsParams"] = Psi4Util.ProcessPsi4OptionsParameters("--psi4OptionsParams", Options["--psi4OptionsParams"])
 1374     OptionsInfo["Psi4RunParams"] = Psi4Util.ProcessPsi4RunParameters("--psi4RunParams", Options["--psi4RunParams"], InfileName = OptionsInfo["Infile"])
 1375 
 1376     # Conformer generation paramaters...
 1377     ParamsDefaultInfoOverride = {"MaxConfs": 250, "MaxConfsTorsions": 50}
 1378     OptionsInfo["ConfGenerationParams"] = MiscUtil.ProcessOptionConformerParameters("--confParams", Options["--confParams"], ParamsDefaultInfoOverride)
 1379     
 1380     # Energy units and label...
 1381     OptionsInfo["EnergyUnits"] = Options["--energyUnits"]
 1382     
 1383     EnergyDataFieldLabel = Options["--energyDataFieldLabel"]
 1384     if re.match("^auto$", EnergyDataFieldLabel, re.I):
 1385         EnergyDataFieldLabel = "Psi4_Energy (%s)" % Options["--energyUnits"]
 1386     OptionsInfo["EnergyDataFieldLabel"] = EnergyDataFieldLabel
 1387     
 1388     EnergyRelativeDataFieldLabel = Options["--energyRelativeDataFieldLabel"]
 1389     if re.match("^auto$", EnergyRelativeDataFieldLabel, re.I):
 1390         EnergyRelativeDataFieldLabel = "Psi4_Relative_Energy (%s)" % Options["--energyUnits"]
 1391     OptionsInfo["EnergyRelativeDataFieldLabel"] = EnergyRelativeDataFieldLabel
 1392 
 1393     # Plot parameters...
 1394     OptionsInfo["OutfileMolName"] = True if re.match("^yes$", Options["--outfileMolName"], re.I) else False
 1395     OptionsInfo["OutPlotRelativeEnergy"] = True if re.match("^yes$", Options["--outPlotRelativeEnergy"], re.I) else False
 1396     OptionsInfo["OutPlotTitleTorsionSpec"] = True if re.match("^yes$", Options["--outPlotTitleTorsionSpec"], re.I) else False
 1397     
 1398     # The default width and height, 10.0 and 7.5, map to aspect raito of 16/9 (1.778)...
 1399     if OptionsInfo["OutPlotRelativeEnergy"]:
 1400         EnergyLabel = "Relative Energy (%s)" % OptionsInfo["EnergyUnits"]
 1401     else:
 1402         EnergyLabel = "Energy (%s)" % OptionsInfo["EnergyUnits"]
 1403         
 1404     DefaultValues = {'Type': 'linepoint', 'Width': 10.0, 'Height': 5.6, 'Title': 'Psi4 Torsion Scan', 'XLabel': 'Torsion Angle (degrees)', 'YLabel': EnergyLabel}
 1405     OptionsInfo["OutPlotParams"] = MiscUtil.ProcessOptionSeabornPlotParameters("--outPlotParams", Options["--outPlotParams"], DefaultValues)
 1406     if not re.match("^(linepoint|scatter|Line)$", OptionsInfo["OutPlotParams"]["Type"], re.I):
 1407         MiscUtil.PrintError("The value, %s, specified for \"type\" using option \"--outPlotParams\" is not supported. Valid plot types: linepoint, scatter or line" % (OptionsInfo["OutPlotParams"]["Type"]))
 1408     
 1409     OptionsInfo["OutPlotInitialized"] = False
 1410     
 1411     OptionsInfo["Overwrite"] = Options["--overwrite"]
 1412 
 1413     OptionsInfo["MaxIters"] = int(Options["--maxIters"])
 1414     
 1415     OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
 1416     OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
 1417     
 1418     # Multiprocessing level...
 1419     MPLevelMoleculesMode = False
 1420     MPLevelTorsionAnglesMode = False
 1421     MPLevel = Options["--mpLevel"]
 1422     if re.match("^Molecules$", MPLevel, re.I):
 1423         MPLevelMoleculesMode = True
 1424     elif re.match("^TorsionAngles$", MPLevel, re.I):
 1425         MPLevelTorsionAnglesMode = True
 1426     else:
 1427         MiscUtil.PrintError("The value, %s, specified for option \"--mpLevel\" is not valid. " % MPLevel)
 1428     OptionsInfo["MPLevel"] = MPLevel
 1429     OptionsInfo["MPLevelMoleculesMode"] = MPLevelMoleculesMode
 1430     OptionsInfo["MPLevelTorsionAnglesMode"] = MPLevelTorsionAnglesMode
 1431     
 1432     OptionsInfo["Precision"] = int(Options["--precision"])
 1433     OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
 1434     
 1435     # Procsss and validate specified SMILES/SMARTS torsion patterns...
 1436     TorsionPatterns = Options["--torsions"]
 1437     TorsionPatternsList = []
 1438     for TorsionPattern in TorsionPatterns.split(","):
 1439         TorsionPattern = TorsionPattern.strip()
 1440         if not len(TorsionPattern):
 1441             MiscUtil.PrintError("Empty value specified for SMILES/SMARTS pattern in  \"-t, --torsions\" option: %s" % TorsionPatterns)
 1442         
 1443         TorsionMol = Chem.MolFromSmarts(TorsionPattern)
 1444         if TorsionMol is None:
 1445             MiscUtil.PrintError("Failed to create torsion pattern molecule. The torsion SMILES/SMARTS pattern, \"%s\", specified using \"-t, --torsions\" option, \"%s\",  is not valid." % (TorsionPattern, TorsionPatterns))
 1446         TorsionPatternsList.append(TorsionPattern)
 1447     
 1448     OptionsInfo["TorsionPatterns"] = TorsionPatterns
 1449     OptionsInfo["TorsionPatternsList"] = TorsionPatternsList
 1450     
 1451     OptionsInfo["TorsionMaxMatches"] = int(Options["--torsionMaxMatches"])
 1452     OptionsInfo["TorsionMinimize"] = True if re.match("^yes$", Options["--torsionMinimize"], re.I) else False
 1453 
 1454     TorsionRange = Options["--torsionRange"]
 1455     TorsionRangeWords = TorsionRange.split(",")
 1456     
 1457     TorsionStart = int(TorsionRangeWords[0])
 1458     TorsionStop = int(TorsionRangeWords[1])
 1459     TorsionStep = int(TorsionRangeWords[2])
 1460     
 1461     if TorsionStart >= TorsionStop:
 1462         MiscUtil.PrintError("The start value, %d, specified for option \"--torsionRange\" in string \"%s\" must be less than stop value, %s." % (TorsionStart, Options["--torsionRange"], TorsionStop))
 1463     if TorsionStep == 0:
 1464         MiscUtil.PrintError("The step value, %d, specified for option \"--torsonRange\" in string \"%s\" must be > 0." % (TorsionStep, Options["--torsionRange"]))
 1465     if TorsionStep >= (TorsionStop - TorsionStart):
 1466         MiscUtil.PrintError("The step value, %d, specified for option \"--torsonRange\" in string \"%s\" must be less than, %s." % (TorsionStep, Options["--torsionRange"], (TorsionStop - TorsionStart)))
 1467     
 1468     if TorsionStart < 0:
 1469         if TorsionStart < -180:
 1470             MiscUtil.PrintError("The start value, %d, specified for option \"--torsionRange\" in string \"%s\" must be  >= -180 to use scan range from -180 to 180." % (TorsionStart, Options["--torsionRange"]))
 1471         if TorsionStop > 180:
 1472             MiscUtil.PrintError("The stop value, %d, specified for option \"--torsionRange\" in string \"%s\" must be <= 180 to use scan range from -180 to 180." % (TorsionStop, Options["--torsionRange"]))
 1473     else:
 1474         if TorsionStop > 360:
 1475             MiscUtil.PrintError("The stop value, %d, specified for option \"--torsionRange\" in string \"%s\" must be  <= 360 to use scan range from 0 to 360." % (TorsionStop, Options["--torsionRange"]))
 1476     
 1477     OptionsInfo["TorsionRange"] = TorsionRange
 1478     OptionsInfo["TorsionStart"] = TorsionStart
 1479     OptionsInfo["TorsionStop"] = TorsionStop
 1480     OptionsInfo["TorsionStep"] = TorsionStep
 1481     
 1482     OptionsInfo["UseChirality"] = True if re.match("^yes$", Options["--useChirality"], re.I) else False
 1483     
 1484 def RetrieveOptions():
 1485     """Retrieve command line arguments and options"""
 1486     
 1487     # Get options...
 1488     global Options
 1489     Options = docopt(_docoptUsage_)
 1490     
 1491     # Set current working directory to the specified directory...
 1492     WorkingDir = Options["--workingdir"]
 1493     if WorkingDir:
 1494         os.chdir(WorkingDir)
 1495     
 1496     # Handle examples option...
 1497     if "--examples" in Options and Options["--examples"]:
 1498         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
 1499         sys.exit(0)
 1500 
 1501 def ValidateOptions():
 1502     """Validate option values"""
 1503 
 1504     MiscUtil.ValidateOptionTextValue("--energyUnits", Options["--energyUnits"], "Hartrees kcal/mol kJ/mol eV")
 1505     
 1506     MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
 1507     MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi txt csv tsv")
 1508     MiscUtil.ValidateOptionTextValue("--infile3D", Options["--infile3D"], "yes no")
 1509 
 1510     MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
 1511     MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
 1512     MiscUtil.ValidateOptionsDistinctFileNames("-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"])
 1513     
 1514     if not Options["--overwrite"]:
 1515         FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
 1516         FileNames = glob.glob("%s_*" % FileName)
 1517         if len(FileNames):
 1518             MiscUtil.PrintError("The outfile names, %s_*, generated from file specified, %s, for option \"-o, --outfile\" already exist. Use option \"--overwrite\" or \"--ov\"  and try again.\n" % (FileName, Options["--outfile"]))
 1519     
 1520     MiscUtil.ValidateOptionTextValue("--outPlotRelativeEnergy", Options["--outPlotRelativeEnergy"], "yes no")
 1521     MiscUtil.ValidateOptionTextValue("--outPlotTitleTorsionSpec", Options["--outPlotTitleTorsionSpec"], "yes no")
 1522     
 1523     MiscUtil.ValidateOptionTextValue("--outfileMolName ", Options["--outfileMolName"], "yes no")
 1524     
 1525     MiscUtil.ValidateOptionTextValue("--modeMols", Options["--modeMols"], "First All")
 1526     MiscUtil.ValidateOptionTextValue("--modeTorsions", Options["--modeTorsions"], "First All")
 1527     
 1528     MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0})
 1529     
 1530     MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
 1531     MiscUtil.ValidateOptionTextValue("--mpLevel", Options["--mpLevel"], "Molecules TorsionAngles")
 1532     
 1533     MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0})
 1534     MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
 1535     
 1536     MiscUtil.ValidateOptionIntegerValue("--torsionMaxMatches", Options["--torsionMaxMatches"], {">": 0})
 1537     MiscUtil.ValidateOptionTextValue("--torsionMinimize", Options["--torsionMinimize"], "yes no")
 1538     MiscUtil.ValidateOptionNumberValues("--torsionRange", Options["--torsionRange"], 3, ",", "integer", {})
 1539     
 1540     MiscUtil.ValidateOptionTextValue("--useChirality", Options["--useChirality"], "yes no")
 1541     
 1542 # Setup a usage string for docopt...
 1543 _docoptUsage_ = """
 1544 Psi4PerformTorsionScan.py - Perform torsion scan
 1545 
 1546 Usage:
 1547     Psi4PerformTorsionScan.py [--basisSet <text>] [--confParams <Name,Value,...>] [--energyDataFieldLabel <text>]
 1548                               [--energyRelativeDataFieldLabel <text>] [--energyUnits <text>] [--infile3D <yes or no>]
 1549                               [--infileParams <Name,Value,...>] [--maxIters <number>] [--methodName <text>]
 1550                               [--modeMols <First or All>] [--modeTorsions <First or All>] [--mp <yes or no>]
 1551                               [--mpLevel <Molecules or TorsionAngles>] [--mpParams <Name,Value,...>]
 1552                               [--outfileMolName <yes or no>] [--outfileParams <Name,Value,...>] [--outPlotParams <Name,Value,...>]
 1553                               [--outPlotRelativeEnergy <yes or no>] [--outPlotTitleTorsionSpec <yes or no>] [--overwrite]
 1554                               [--precision <number>] [--psi4OptionsParams <Name,Value,...>] [--psi4RunParams <Name,Value,...>]
 1555                               [--quiet <yes or no>] [--reference <text>] [--torsionMaxMatches <number>] [--torsionMinimize <yes or no>]
 1556                               [--torsionRange <Start,Stop,Step>] [--useChirality <yes or no>] [-w <dir>] -t <torsions> -i <infile>  -o <outfile> 
 1557     Psi4PerformTorsionScan.py -h | --help | -e | --examples
 1558 
 1559 Description:
 1560     Perform torsion scan for molecules around torsion angles specified using
 1561     SMILES/SMARTS patterns. A molecule is optionally minimized before performing
 1562     a torsion scan using a forcefield. A set of initial 3D structures are generated for
 1563     a molecule by scanning the torsion angle across the specified range and updating
 1564     the 3D coordinates of the molecule. A conformation ensemble is optionally generated
 1565     for each 3D structure representing a specific torsion angle using a combination of
 1566     distance geometry and forcefield followed by constrained geometry optimization
 1567     using a quantum chemistry method. The conformation with the lowest energy is
 1568     selected to represent the torsion angle. An option is available to skip the generation
 1569     of the conformation ensemble and simply calculate the energy for the initial 3D
 1570     structure for a specific torsion torsion angle using a quantum chemistry method.
 1571     
 1572     The torsions are specified using SMILES or SMARTS patterns. A substructure match
 1573     is performed to select torsion atoms in a molecule. The SMILES pattern match must
 1574     correspond to four torsion atoms. The SMARTS patterns containing atom indices may
 1575     match  more than four atoms. The atoms indices, however, must match exactly four
 1576     torsion atoms. For example: [s:1][c:2]([aX2,cH1])!@[CX3:3](O)=[O:4] for thiophene
 1577     esters and carboxylates as specified in Torsion Library (TorLib) [Ref 146].
 1578 
 1579     A Psi4 XYZ format geometry string is automatically generated for each molecule
 1580     in input file. It contains atom symbols and 3D coordinates for each atom in a
 1581     molecule. In addition, the formal charge and spin multiplicity are present in the
 1582     the geometry string. These values are either retrieved from molecule properties
 1583     named 'FormalCharge' and 'SpinMultiplicty' or dynamically calculated for a
 1584     molecule.
 1585     
 1586     A set of four output files is generated for each torsion match in each
 1587     molecule. The names of the output files are generated using the root of
 1588     the specified output file. They may either contain sequential molecule
 1589     numbers or molecule names as shown below:
 1590         
 1591         <OutfileRoot>_Mol<Num>.sdf
 1592         <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>.sdf
 1593         <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Energies.csv
 1594         <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Plot.<ImgExt>
 1595         
 1596         or
 1597         
 1598         <OutfileRoot>_<MolName>.sdf
 1599         <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>.sdf
 1600         <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Energies.csv
 1601         <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Plot.<ImgExt>
 1602         
 1603     The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi,
 1604     .csv, .tsv, .txt)
 1605 
 1606     The supported output file formats are: SD (.sdf, .sd)
 1607 
 1608 Options:
 1609     -b, --basisSet <text>  [default: auto]
 1610         Basis set to use for energy calculation or constrained energy minimization.
 1611         Default: 6-31+G** for sulfur containing molecules; Otherwise, 6-31G** [ Ref 150 ].
 1612         The specified value must be a valid Psi4 basis set. No validation is performed.
 1613         
 1614         The following list shows a representative sample of basis sets available
 1615         in Psi4:
 1616             
 1617             STO-3G, 6-31G, 6-31+G, 6-31++G, 6-31G*, 6-31+G*,  6-31++G*, 
 1618             6-31G**, 6-31+G**, 6-31++G**, 6-311G, 6-311+G, 6-311++G,
 1619             6-311G*, 6-311+G*, 6-311++G*, 6-311G**, 6-311+G**, 6-311++G**,
 1620             cc-pVDZ, cc-pCVDZ, aug-cc-pVDZ, cc-pVDZ-DK, cc-pCVDZ-DK, def2-SVP,
 1621             def2-SVPD, def2-TZVP, def2-TZVPD, def2-TZVPP, def2-TZVPPD
 1622             
 1623     --confParams <Name,Value,...>  [default: auto]
 1624         A comma delimited list of parameter name and value pairs for generating
 1625         initial 3D coordinates for molecules in input file at specific torsion angles. A
 1626         conformation ensemble is optionally generated for each 3D structure
 1627         representing a specific torsion angle using a combination of distance geometry
 1628         and forcefield followed by constrained geometry optimization using a quantum
 1629         chemistry method. The conformation with the lowest energy is selected to
 1630         represent the torsion angle.
 1631         
 1632         The supported parameter names along with their default values are shown
 1633         below:
 1634             
 1635             confMethod,ETKDG,
 1636             forceField,MMFF, forceFieldMMFFVariant,MMFF94,
 1637             enforceChirality,yes,embedRMSDCutoff,0.5,maxConfs,250,
 1638             maxConfsTorsions,50,useTethers,yes
 1639             
 1640             confMethod,ETKDG   [ Possible values: SDG, ETDG, KDG, ETKDG ]
 1641             forceField, MMFF   [ Possible values: UFF or MMFF ]
 1642             forceFieldMMFFVariant,MMFF94   [ Possible values: MMFF94 or MMFF94s ]
 1643             enforceChirality,yes   [ Possible values: yes or no ]
 1644             useTethers,yes   [ Possible values: yes or no ]
 1645             
 1646         confMethod: Conformation generation methodology for generating initial 3D
 1647         coordinates. Possible values: Standard Distance Geometry (SDG), Experimental
 1648         Torsion-angle preference with Distance Geometry (ETDG), basic Knowledge-terms
 1649         with Distance Geometry (KDG) and Experimental Torsion-angle preference
 1650         along with basic Knowledge-terms with Distance Geometry (ETKDG) [Ref 129] .
 1651         
 1652         forceField: Forcefield method to use for energy minimization. Possible values:
 1653         Universal Force Field (UFF) [ Ref 81 ] or Merck Molecular Mechanics Force
 1654         Field [ Ref 83-87 ] .
 1655         
 1656         enforceChirality: Enforce chirality for defined chiral centers during
 1657         forcefield minimization.
 1658         
 1659         maxConfs: Maximum number of conformations to generate for each molecule
 1660         during the generation of an initial 3D conformation ensemble using a conformation
 1661         generation methodology. The conformations are minimized using the specified
 1662         forcefield. The lowest energy structure is selected for performing the torsion scan.
 1663         
 1664         maxConfsTorsion: Maximum number of 3D conformations to generate for
 1665         conformation ensemble representing a specific torsion. The conformations are
 1666         constrained at specific torsions angles and minimized using the specified forcefield
 1667         and a quantum chemistry method. The lowest energy conformation is selected to
 1668         calculate final torsion energy and written to the output file.
 1669         
 1670         embedRMSDCutoff: RMSD cutoff for retaining initial set of conformers embedded
 1671         using distance geometry and forcefield minimization. All embedded conformers
 1672         are kept for 'None' value. Otherwise, only those conformers which are different
 1673         from each other by the specified RMSD cutoff, 0.5 by default, are kept. The first
 1674         embedded conformer is always retained.
 1675         
 1676         useTethers: Use tethers to optimize the final embedded conformation by
 1677         applying a series of extra forces to align matching atoms to the positions of
 1678         the core atoms. Otherwise, use simple distance constraints during the
 1679         optimization.
 1680     --energyDataFieldLabel <text>  [default: auto]
 1681         Energy data field label for writing energy values. Default: Psi4_Energy (<Units>). 
 1682     --energyRelativeDataFieldLabel <text>  [default: auto]
 1683         Relative energy data field label for writing energy values. Default:
 1684         Psi4_Relative_Energy (<Units>). 
 1685     --energyUnits <text>  [default: kcal/mol]
 1686         Energy units. Possible values: Hartrees, kcal/mol, kJ/mol, or eV.
 1687     -e, --examples
 1688         Print examples.
 1689     -h, --help
 1690         Print this help message.
 1691     -i, --infile <infile>
 1692         Input file name.
 1693     --infile3D <yes or no>  [default: no]
 1694         Skip generation and minimization of initial 3D structures for molecules in
 1695         input file containing 3D coordinates.
 1696     --infileParams <Name,Value,...>  [default: auto]
 1697         A comma delimited list of parameter name and value pairs for reading
 1698         molecules from files. The supported parameter names for different file
 1699         formats, along with their default values, are shown below:
 1700             
 1701             SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
 1702             
 1703             SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
 1704                 smilesTitleLine,auto,sanitize,yes
 1705             
 1706         Possible values for smilesDelimiter: space, comma or tab.
 1707     --maxIters <number>  [default: 50]
 1708         Maximum number of iterations to perform for each molecule or conformer
 1709         during constrained energy minimization by a quantum chemistry method.
 1710     -m, --methodName <text>  [default: auto]
 1711         Method to use for energy calculation or constrained energy minimization.
 1712         Default: B3LYP [ Ref 150 ]. The specified value must be a valid Psi4 method
 1713         name. No validation is performed.
 1714         
 1715         The following list shows a representative sample of methods available
 1716         in Psi4:
 1717             
 1718             B1LYP, B2PLYP, B2PLYP-D3BJ, B2PLYP-D3MBJ, B3LYP, B3LYP-D3BJ,
 1719             B3LYP-D3MBJ, CAM-B3LYP, CAM-B3LYP-D3BJ, HF, HF-D3BJ,  HF3c, M05,
 1720             M06, M06-2x, M06-HF, M06-L, MN12-L, MN15, MN15-D3BJ,PBE, PBE0,
 1721             PBEH3c, PW6B95, PW6B95-D3BJ, WB97, WB97X, WB97X-D, WB97X-D3BJ
 1722             
 1723     --modeMols <First or All>  [default: First]
 1724         Perform torsion scan for the first molecule or all molecules in input
 1725         file.
 1726     --modeTorsions <First or All>  [default: First]
 1727         Perform torsion scan for the first or all specified torsion pattern in
 1728         molecules up to a maximum number of matches for each torsion
 1729         specification as indicated by '--torsionMaxMatches' option. 
 1730     --mp <yes or no>  [default: no]
 1731         Use multiprocessing.
 1732          
 1733         By default, input data is retrieved in a lazy manner via mp.Pool.imap()
 1734         function employing lazy RDKit data iterable. This allows processing of
 1735         arbitrary large data sets without any additional requirements memory.
 1736         
 1737         All input data may be optionally loaded into memory by mp.Pool.map()
 1738         before starting worker processes in a process pool by setting the value
 1739         of 'inputDataMode' to 'InMemory' in '--mpParams' option.
 1740         
 1741         A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
 1742         data mode may adversely impact the performance. The '--mpParams' section
 1743         provides additional information to tune the value of 'chunkSize'.
 1744     --mpLevel <Molecules or TorsionAngles>  [default: Molecules]
 1745         Perform multiprocessing at molecules or torsion angles level. Possible values:
 1746         Molecules or TorsionAngles. The 'Molecules' value starts a process pool at the
 1747         molecules level. All torsion angles of a molecule are processed in a single
 1748         process. The 'TorsionAngles' value, however, starts a process pool at the 
 1749         torsion angles level. Each torsion angle in a torsion match for a molecule is
 1750         processed in an individual process in the process pool.
 1751     --mpParams <Name,Value,...>  [default: auto]
 1752         A comma delimited list of parameter name and value pairs to configure
 1753         multiprocessing.
 1754         
 1755         The supported parameter names along with their default and possible
 1756         values are shown below:
 1757         
 1758             chunkSize, auto
 1759             inputDataMode, Lazy   [ Possible values: InMemory or Lazy ]
 1760             numProcesses, auto   [ Default: mp.cpu_count() ]
 1761         
 1762         These parameters are used by the following functions to configure and
 1763         control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
 1764         mp.Pool.imap().
 1765         
 1766         The chunkSize determines chunks of input data passed to each worker
 1767         process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
 1768         The default value of chunkSize is dependent on the value of 'inputDataMode'.
 1769         
 1770         The mp.Pool.map() function, invoked during 'InMemory' input data mode,
 1771         automatically converts RDKit data iterable into a list, loads all data into
 1772         memory, and calculates the default chunkSize using the following method
 1773         as shown in its code:
 1774         
 1775             chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
 1776             if extra: chunkSize += 1
 1777         
 1778         For example, the default chunkSize will be 7 for a pool of 4 worker processes
 1779         and 100 data items.
 1780         
 1781         The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
 1782         'lazy' RDKit data iterable to retrieve data as needed, without loading all the
 1783         data into memory. Consequently, the size of input data is not known a priori.
 1784         It's not possible to estimate an optimal value for the chunkSize. The default 
 1785         chunkSize is set to 1.
 1786         
 1787         The default value for the chunkSize during 'Lazy' data mode may adversely
 1788         impact the performance due to the overhead associated with exchanging
 1789         small chunks of data. It is generally a good idea to explicitly set chunkSize to
 1790         a larger value during 'Lazy' input data mode, based on the size of your input
 1791         data and number of processes in the process pool.
 1792         
 1793         The mp.Pool.map() function waits for all worker processes to process all
 1794         the data and return the results. The mp.Pool.imap() function, however,
 1795         returns the the results obtained from worker processes as soon as the
 1796         results become available for specified chunks of data.
 1797         
 1798         The order of data in the results returned by both mp.Pool.map() and 
 1799         mp.Pool.imap() functions always corresponds to the input data.
 1800     -o, --outfile <outfile>
 1801         Output file name. The output file root is used for generating the names
 1802         of the output files corresponding to structures, energies, and plots during
 1803         the torsion scan.
 1804     --outfileMolName <yes or no>  [default: no]
 1805         Append molecule name to output file root during the generation of the names
 1806         for output files. The default is to use <MolNum>. The non alphabetical
 1807         characters in molecule names are replaced by underscores.
 1808     --outfileParams <Name,Value,...>  [default: auto]
 1809         A comma delimited list of parameter name and value pairs for writing
 1810         molecules to files. The supported parameter names for different file
 1811         formats, along with their default values, are shown below:
 1812             
 1813             SD: kekulize,yes
 1814             
 1815     --outPlotParams <Name,Value,...>  [default: auto]
 1816         A comma delimited list of parameter name and value pairs for generating
 1817         plots using Seaborn module. The supported parameter names along with their
 1818         default values are shown below:
 1819             
 1820             type,linepoint,outExt,svg,width,10,height,5.6,
 1821             title,auto,xlabel,auto,ylabel,auto,titleWeight,bold,labelWeight,bold
 1822             style,darkgrid,palette,deep,font,sans-serif,fontScale,1,
 1823             context,notebook
 1824             
 1825         Possible values:
 1826             
 1827             type: linepoint, scatter, or line. Both points and lines are drawn
 1828                 for linepoint plot type.
 1829             outExt: Any valid format supported by Python module Matplotlib.
 1830                 For example: PDF (.pdf), PNG (.png), PS (.ps), SVG (.svg)
 1831             titleWeight, labelWeight: Font weight for title and axes labels.
 1832                 Any valid value.
 1833             style: darkgrid, whitegrid, dark, white, ticks
 1834             palette: deep, muted, pastel, dark, bright, colorblind
 1835             font: Any valid font name
 1836             
 1837     --outPlotRelativeEnergy <yes or no>  [default: yes]
 1838         Plot relative energies in the torsion plot. The minimum energy value is
 1839         subtracted from energy values to calculate relative energies.
 1840     --outPlotTitleTorsionSpec <yes or no>  [default: yes]
 1841         Append torsion specification to the title of the torsion plot.
 1842     --overwrite
 1843         Overwrite existing files.
 1844     --precision <number>  [default: 6]
 1845         Floating point precision for writing energy values.
 1846     --psi4OptionsParams <Name,Value,...>  [default: none]
 1847         A comma delimited list of Psi4 option name and value pairs for setting
 1848         global and module options. The names are 'option_name' for global options
 1849         and 'module_name__option_name' for options local to a module. The
 1850         specified option names must be valid Psi4 names. No validation is
 1851         performed.
 1852         
 1853         The specified option name and  value pairs are processed and passed to
 1854         psi4.set_options() as a dictionary. The supported value types are float,
 1855         integer, boolean, or string. The float value string is converted into a float.
 1856         The valid values for a boolean string are yes, no, true, false, on, or off. 
 1857     --psi4RunParams <Name,Value,...>  [default: auto]
 1858         A comma delimited list of parameter name and value pairs for configuring
 1859         Psi4 jobs.
 1860         
 1861         The supported parameter names along with their default and possible
 1862         values are shown below:
 1863              
 1864             MemoryInGB, 1
 1865             NumThreads, 1
 1866             OutputFile, auto   [ Possible  values: stdout, quiet, or FileName ]
 1867             ScratchDir, auto   [ Possivle values: DirName]
 1868             RemoveOutputFile, yes   [ Possible values: yes, no, true, or false]
 1869             
 1870         These parameters control the runtime behavior of Psi4.
 1871         
 1872         The default file name for 'OutputFile' is <InFileRoot>_Psi4.out. The PID
 1873         is appended to output file name during multiprocessing as shown:
 1874         <InFileRoot>_Psi4_<PIDNum>.out. The 'stdout' value for 'OutputType'
 1875         sends Psi4 output to stdout. The 'quiet' or 'devnull' value suppresses
 1876         all Psi4 output. The 'OutputFile' is set to 'quiet' for 'auto' value during 
 1877         'Conformers' of '--mpLevel' option.
 1878         
 1879         The default 'Yes' value of 'RemoveOutputFile' option forces the removal
 1880         of any existing Psi4 before creating new files to append output from
 1881         multiple Psi4 runs.
 1882         
 1883         The option 'ScratchDir' is a directory path to the location of scratch
 1884         files. The default value corresponds to Psi4 default. It may be used to
 1885         override the deafult path.
 1886     -q, --quiet <yes or no>  [default: no]
 1887         Use quiet mode. The warning and information messages will not be printed.
 1888     --reference <text>  [default: auto]
 1889         Reference wave function to use for energy calculation or constrained energy
 1890         minimization. Default: RHF or UHF. The default values are Restricted Hartree-Fock
 1891         (RHF) for closed-shell molecules with all electrons paired and Unrestricted
 1892         Hartree-Fock (UHF) for open-shell molecules with unpaired electrons.
 1893         
 1894         The specified value must be a valid Psi4 reference wave function. No validation
 1895         is performed. For example: ROHF, CUHF, RKS, etc.
 1896         
 1897         The spin multiplicity determines the default value of reference wave function
 1898         for input molecules. It is calculated from number of free radical electrons using
 1899         Hund's rule of maximum multiplicity defined as 2S + 1 where S is the total
 1900         electron spin. The total spin is 1/2 the number of free radical electrons in a 
 1901         molecule. The value of 'SpinMultiplicity' molecule property takes precedence
 1902         over the calculated value of spin multiplicity.
 1903     -t, --torsions <SMILES/SMARTS,...,...>
 1904         SMILES/SMARTS patterns corresponding to torsion specifications. It's a 
 1905         comma delimited list of valid SMILES/SMART patterns.
 1906         
 1907         A substructure match is performed to select torsion atoms in a molecule.
 1908         The SMILES pattern match must correspond to four torsion atoms. The
 1909         SMARTS patterns contain atom indices may match  more than four atoms.
 1910         The atoms indices, however, must match exactly four torsion atoms. For example:
 1911         [s:1][c:2]([aX2,cH1])!@[CX3:3](O)=[O:4] for thiophene esters and carboxylates
 1912         as specified in Torsion Library (TorLib) [Ref 146].
 1913     --torsionMaxMatches <number>  [default: 5]
 1914         Maximum number of torsions to match for each torsion specification in a
 1915         molecule.
 1916     --torsionMinimize <yes or no>  [default: no]
 1917         Perform constrained energy minimization on a conformation ensemble
 1918         for  a specific torsion angle and select the lowest energy conformation
 1919         representing the torsion angle. A conformation ensemble is generated for
 1920         each 3D structure representing a specific torsion angle using a combination
 1921         of distance geometry and forcefield followed by constrained geometry
 1922         optimization using a quantum chemistry method.
 1923     --torsionRange <Start,Stop,Step>  [default: 0,360,5]
 1924         Start, stop, and step size angles in degrees for a torsion scan. In addition,
 1925         you may specify values using start and stop angles from -180 to 180.
 1926     --useChirality <yes or no>  [default: no]
 1927         Use chirrality during substructure matches for identification of torsions.
 1928     -w, --workingdir <dir>
 1929         Location of working directory which defaults to the current directory.
 1930 
 1931 Examples:
 1932     To perform a torsion scan on the first molecule in a SMILES file using a minimum
 1933     energy structure of the molecule selected from an initial ensemble of conformations
 1934     generated using distance geometry and forcefield, skip generation of conformation
 1935     ensembles for specific torsion angles and constrained energy minimization of the
 1936     ensemble, calculating single point at a specific torsion angle energy using B3LYP/6-31G**
 1937     and B3LYP/6-31+G** for non-sulfur and sulfur containing molecules, generate output files
 1938     corresponding to structure, energy and torsion plot, type:
 1939     
 1940         % Psi4PerformTorsionScan.py  -t "CCCC" -i Psi4SampleTorsionScan.smi 
 1941           -o SampleOut.sdf
 1942 
 1943     To run the previous example on the first molecule in a SD file containing 3D
 1944     coordinates and skip the generations of initial 3D structure, type: 
 1945     
 1946         % Psi4PerformTorsionScan.py  -t "CCCC"  --infile3D yes
 1947           -i Psi4SampleTorsionScan3D.sdf  -o SampleOut.sdf
 1948 
 1949     To run the first example on all molecules in a SD file, type:
 1950     
 1951         % Psi4PerformTorsionScan.py  -t "CCCC" --modeMols All
 1952           -i Psi4SampleTorsionScan.sdf -o SampleOut.sdf
 1953 
 1954     To run the first example on all molecules in a SD file containing 3D
 1955     coordinates and skip the generation of initial 3D structures, type: 
 1956     
 1957         % Psi4PerformTorsionScan.py  -t "CCCC"  --infile3D yes
 1958           --modeMols All -i Psi4SampleTorsionScan3D.sdf  -o SampleOut.sdf
 1959 
 1960     To perform a torsion scan on the first molecule in a SMILES file using a minimum
 1961     energy structure of the molecule selected from an initial ensemble of conformations
 1962     generated using distance geometry and forcefield,  generate up to 50 conformations
 1963     for specific torsion angles using ETKDG methodology followed by initial MMFF
 1964     forcefield minimization and final energy minimization using B3LYP/6-31G** and
 1965     B3LYP/6-31+G** for non-sulfur and sulfur containing molecules, generate output files
 1966     corresponding to minimum energy structure, energy and torsion plot, type:
 1967 
 1968         % Psi4PerformTorsionScan.py  -t "CCCC" --torsionMinimize Yes
 1969            -i Psi4SampleTorsionScan.smi -o SampleOut.sdf
 1970 
 1971     To run the previous example on all molecules in a SD file, type:
 1972     
 1973         % Psi4PerformTorsionScan.py  -t "CCCC" --modeMols All
 1974            --torsionMinimize Yes -i Psi4SampleTorsionScan.sdf -o SampleOut.sdf
 1975 
 1976     To run the previous example on all molecules in a SD file containing 3D
 1977     coordinates and skip the generation of initial 3D structures, type:
 1978     
 1979         % Psi4PerformTorsionScan.py  -t "CCCC" --modeMols All
 1980            --infile3D yes --modeMols All  --torsionMinimize Yes
 1981            -i Psi4SampleTorsionScan.sdf -o SampleOut.sdf
 1982 
 1983     To run the previous example in multiprocessing mode at molecules level
 1984     on all available CPUs without loading all data into memory and write out
 1985     a SD file, type:
 1986 
 1987         % Psi4PerformTorsionScan.py  -t "CCCC" -i Psi4SampleTorsionScan.smi 
 1988           -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
 1989     
 1990     To run the previous example in multiprocessing mode at torsion angles level
 1991     on all available CPUs without loading all data into memory and write out
 1992     a SD file, type:
 1993 
 1994         % Psi4PerformTorsionScan.py  -t "CCCC" -i Psi4SampleTorsionScan.smi 
 1995           -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
 1996           --mpLevel TorsionAngles
 1997     
 1998     To run the previous example in multiprocessing mode on all available CPUs
 1999     by loading all data into memory and write out a SD file, type:
 2000 
 2001         % Psi4PerformTorsionScan.py  -t "CCCC" -i Psi4SampleTorsionScan.smi 
 2002           -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
 2003           --mpParams "inputDataMode,InMemory"
 2004     
 2005     To run the previous example in multiprocessing mode on specific number of
 2006     CPUs and chunk size without loading all data into memory and write out a SD file,
 2007     type:
 2008 
 2009         % Psi4PerformTorsionScan.py  -t "CCCC" -i Psi4SampleTorsionScan.smi 
 2010           -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
 2011           --mpParams "inputDataMode,Lazy,numProcesses,4,chunkSize,8"
 2012 
 2013 Author:
 2014     Manish Sud(msud@san.rr.com)
 2015 
 2016 See also:
 2017     Psi4CalculateEnergy.py, Psi4GenerateConformers.py,
 2018     Psi4GenerateConstrainedConformers.py, Psi4PerformConstrainedMinimization.py
 2019 
 2020 Copyright:
 2021     Copyright (C) 2021 Manish Sud. All rights reserved.
 2022 
 2023     The functionality available in this script is implemented using RDKit, an
 2024     open source toolkit for cheminformatics developed by Greg Landrum.
 2025 
 2026     This file is part of MayaChemTools.
 2027 
 2028     MayaChemTools is free software; you can redistribute it and/or modify it under
 2029     the terms of the GNU Lesser General Public License as published by the Free
 2030     Software Foundation; either version 3 of the License, or (at your option) any
 2031     later version.
 2032 
 2033 """
 2034 
 2035 if __name__ == "__main__":
 2036     main()