MayaChemTools

    1 #!/bin/env python
    2 #
    3 # File: RDKitPerformTorsionScan.py
    4 # Author: Manish Sud <msud@san.rr.com>
    5 #
    6 # Copyright (C) 2020 Manish Sud. All rights reserved.
    7 #
    8 # The functionality available in this script is implemented using RDKit, an
    9 # open source toolkit for cheminformatics developed by Greg Landrum.
   10 #
   11 # This file is part of MayaChemTools.
   12 #
   13 # MayaChemTools is free software; you can redistribute it and/or modify it under
   14 # the terms of the GNU Lesser General Public License as published by the Free
   15 # Software Foundation; either version 3 of the License, or (at your option) any
   16 # later version.
   17 #
   18 # MayaChemTools is distributed in the hope that it will be useful, but without
   19 # any warranty; without even the implied warranty of merchantability of fitness
   20 # for a particular purpose.  See the GNU Lesser General Public License for more
   21 # details.
   22 #
   23 # You should have received a copy of the GNU Lesser General Public License
   24 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
   25 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
   26 # Boston, MA, 02111-1307, USA.
   27 #
   28 
   29 from __future__ import print_function
   30 
   31 # Add local python path to the global path and import standard library modules...
   32 import os
   33 import sys;  sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
   34 import time
   35 import re
   36 import glob
   37 import multiprocessing as mp
   38 
   39 import matplotlib.pyplot as plt
   40 import seaborn as sns
   41 
   42 # RDKit imports...
   43 try:
   44     from rdkit import rdBase
   45     from rdkit import Chem
   46     from rdkit.Chem import AllChem
   47     from rdkit.Chem import rdFMCS
   48     from rdkit.Chem import rdMolTransforms
   49 except ImportError as ErrMsg:
   50     sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
   51     sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
   52     sys.exit(1)
   53 
   54 # MayaChemTools imports...
   55 try:
   56     from docopt import docopt
   57     import MiscUtil
   58     import RDKitUtil
   59 except ImportError as ErrMsg:
   60     sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
   61     sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
   62     sys.exit(1)
   63 
   64 ScriptName = os.path.basename(sys.argv[0])
   65 Options = {}
   66 OptionsInfo = {}
   67 
   68 def main():
   69     """Start execution of the script"""
   70     
   71     MiscUtil.PrintInfo("\n%s (RDK v%s; %s): Starting...\n" % (ScriptName, rdBase.rdkitVersion, time.asctime()))
   72     
   73     (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
   74     
   75     # Retrieve command line arguments and options...
   76     RetrieveOptions()
   77     
   78     # Process and validate command line arguments and options...
   79     ProcessOptions()
   80     
   81     # Perform actions required by the script...
   82     PerformTorsionScan()
   83 
   84     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
   85     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
   86 
   87 def PerformTorsionScan():
   88     """Perform torsion scan."""
   89     
   90     # Setup a molecule reader for input file...
   91     MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
   92     OptionsInfo["InfileParams"]["AllowEmptyMols"] = True
   93     Mols  = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
   94     
   95     PlotExt = OptionsInfo["OutPlotParams"]["OutExt"]
   96     FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
   97     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))
   98 
   99     MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount = ProcessMolecules(Mols)
  100     
  101     MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
  102     MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
  103     MiscUtil.PrintInfo("Number of molecules failed during initial minimization: %d" % MinimizationFailedCount)
  104     MiscUtil.PrintInfo("Number of molecules without any matched torsions: %d" % TorsionsMissingCount)
  105     MiscUtil.PrintInfo("Number of molecules failed during torsion scan: %d" % TorsionsScanFailedCount)
  106     MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + TorsionsMissingCount + MinimizationFailedCount + TorsionsScanFailedCount))
  107 
  108 def ProcessMolecules(Mols):
  109     """Process molecules to perform torsion scan. """
  110 
  111     if OptionsInfo["MPMode"]:
  112         return ProcessMoleculesUsingMultipleProcesses(Mols)
  113     else:
  114         return ProcessMoleculesUsingSingleProcess( Mols)
  115 
  116 def ProcessMoleculesUsingSingleProcess(Mols):
  117     """Process molecules to perform torsion scan using a single process."""
  118 
  119     MolInfoText = "first molecule"
  120     if not OptionsInfo["FirstMolMode"]:
  121         MolInfoText = "all molecules"
  122 
  123     if OptionsInfo["TorsionMinimize"]:
  124         MiscUtil.PrintInfo("\nPeforming torsion scan on %s by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  125     else:
  126         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))
  127     
  128     SetupTorsionsPatternsInfo()
  129     
  130     (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
  131 
  132     for Mol in Mols:
  133         MolCount += 1
  134 
  135         if OptionsInfo["FirstMolMode"] and MolCount > 1:
  136             MolCount -= 1
  137             break
  138         
  139         if Mol is None:
  140             continue
  141         
  142         if RDKitUtil.IsMolEmpty(Mol):
  143             if not OptionsInfo["QuietMode"]:
  144                 MolName = RDKitUtil.GetMolName(Mol, MolCount)
  145                 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
  146             continue
  147         ValidMolCount += 1
  148 
  149         Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(Mol, MolCount)
  150         
  151         if not MinimizationCalcStatus:
  152             MinimizationFailedCount += 1
  153             continue
  154         
  155         if not TorsionsMatchStatus:
  156             TorsionsMissingCount += 1
  157             continue
  158 
  159         if not TorsionsScanCalcStatus:
  160             TorsionsScanFailedCount += 1
  161             continue
  162 
  163     return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
  164 
  165 def ProcessMoleculesUsingMultipleProcesses(Mols):
  166     """Process molecules to perform torsion scan using multiprocessing."""
  167 
  168     MolInfoText = "first molecule"
  169     if not OptionsInfo["FirstMolMode"]:
  170         MolInfoText = "all molecules"
  171 
  172     if OptionsInfo["TorsionMinimize"]:
  173         MiscUtil.PrintInfo("\nPeforming torsion scan on %s using multiprocessing by generating conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  174     else:
  175         MiscUtil.PrintInfo("\nPeforming torsion scan %s using multiprocessing by skipping generation of conformation ensembles for specific torsion angles and constrained energy minimization of the ensembles..." % (MolInfoText))
  176         
  177     
  178     MPParams = OptionsInfo["MPParams"]
  179     
  180     # Setup data for initializing a worker process...
  181     MiscUtil.PrintInfo("Encoding options info...")
  182     
  183     InitializeWorkerProcessArgs = (MiscUtil.ObjectToBase64EncodedString(Options), MiscUtil.ObjectToBase64EncodedString(OptionsInfo))
  184 
  185     if OptionsInfo["FirstMolMode"]:
  186         Mol = Mols[0]
  187         Mols = [Mol]
  188 
  189     # Setup a encoded mols data iterable for a worker process...
  190     WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
  191 
  192     # Setup process pool along with data initialization for each process...
  193     MiscUtil.PrintInfo("\nConfiguring multiprocessing using %s method..." % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()"))
  194     MiscUtil.PrintInfo("NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" % (MPParams["NumProcesses"], MPParams["InputDataMode"], ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"])))
  195     
  196     ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
  197     
  198     # Start processing...
  199     if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
  200         Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
  201     elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
  202         Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
  203     else:
  204         MiscUtil.PrintError("The value, %s, specified for \"--inputDataMode\" is not supported." % (MPParams["InputDataMode"]))
  205 
  206     (MolCount, ValidMolCount, TorsionsMissingCount, MinimizationFailedCount, TorsionsScanFailedCount) = [0] * 5
  207     
  208     for Result in Results:
  209         MolCount += 1
  210         
  211         MolIndex, EncodedMol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = Result
  212         
  213         if EncodedMol is None:
  214             continue
  215         ValidMolCount += 1
  216     
  217         if not MinimizationCalcStatus:
  218             MinimizationFailedCount += 1
  219             continue
  220         
  221         if not TorsionsMatchStatus:
  222             TorsionsMissingCount += 1
  223             continue
  224 
  225         if not TorsionsScanCalcStatus:
  226             TorsionsScanFailedCount += 1
  227             continue
  228         
  229         Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
  230         
  231     return (MolCount, ValidMolCount, MinimizationFailedCount, TorsionsMissingCount, TorsionsScanFailedCount)
  232     
  233 def InitializeWorkerProcess(*EncodedArgs):
  234     """Initialize data for a worker process."""
  235     
  236     global Options, OptionsInfo
  237 
  238     MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
  239 
  240     # Decode Options and OptionInfo...
  241     Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
  242     OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
  243 
  244     # Initialize torsion patterns info...
  245     SetupTorsionsPatternsInfo()
  246 
  247 def WorkerProcess(EncodedMolInfo):
  248     """Process data for a worker process."""
  249 
  250     MolIndex, EncodedMol = EncodedMolInfo
  251 
  252     (MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus) = [False] * 3
  253     
  254     if EncodedMol is None:
  255         return [MolIndex, None, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
  256     
  257     Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
  258     if RDKitUtil.IsMolEmpty(Mol):
  259         if not OptionsInfo["QuietMode"]:
  260             MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1))
  261             MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
  262         return [MolIndex, None, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
  263     
  264     Mol, MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus = PerformMinimizationAndTorsionScan(Mol, (MolIndex + 1))
  265     
  266     return [MolIndex, RDKitUtil.MolToBase64EncodedMolString(Mol, PropertyPickleFlags = Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps), MinimizationCalcStatus, TorsionsMatchStatus, TorsionsScanCalcStatus]
  267 
  268 def PerformMinimizationAndTorsionScan(Mol, MolNum = None):
  269     """Peform minimization and torsions scan."""
  270     
  271     if not OptionsInfo["Infile3D"]:
  272         Mol, MinimizationCalcStatus = MinimizeMolecule(Mol, MolNum)
  273         if not MinimizationCalcStatus:
  274             return (Mol, False, False, False)
  275         
  276     TorsionsMolInfo = SetupTorsionsMolInfo(Mol, MolNum)
  277     if TorsionsMolInfo["NumOfMatches"] == 0:
  278         return (Mol, True, False, False)
  279     
  280     Mol, ScanCalcStatus = ScanAllTorsionsInMol(Mol, TorsionsMolInfo, MolNum)
  281     if not ScanCalcStatus:
  282         return (Mol, True, True, False)
  283     
  284     return (Mol, True, True, True)
  285 
  286 def ScanAllTorsionsInMol(Mol, TorsionsMolInfo,  MolNum = None):
  287     """Peform scans on all torsions in a molecule."""
  288     
  289     if TorsionsMolInfo["NumOfMatches"] == 0:
  290         return Mol, True
  291 
  292     Mol = AddHydrogens(Mol)
  293     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  294     
  295     FirstTorsionMode = OptionsInfo["FirstTorsionMode"]
  296     TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
  297 
  298     TorsionPatternCount, TorsionScanCount, TorsionMatchCount = [0] * 3
  299     TorsionMaxMatches = OptionsInfo["TorsionMaxMatches"]
  300     
  301     for TorsionID in TorsionsPatternsInfo["IDs"]:
  302         TorsionPatternCount +=  1
  303         TorsionPattern = TorsionsPatternsInfo["Pattern"][TorsionID]
  304         TorsionPatternMol = TorsionsPatternsInfo["Mol"][TorsionID]
  305         
  306         TorsionsMatches = TorsionsMolInfo["Matches"][TorsionID]
  307 
  308         if TorsionsMatches is None:
  309             continue
  310         
  311         if FirstTorsionMode and TorsionPatternCount > 1:
  312             if not OptionsInfo["QuietMode"]:
  313                 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"]))
  314             break
  315 
  316         for Index, TorsionMatches in enumerate(TorsionsMatches):
  317             TorsionMatchNum = Index + 1
  318             TorsionMatchCount +=  1
  319 
  320             if TorsionMatchCount > TorsionMaxMatches:
  321                 if not OptionsInfo["QuietMode"]:
  322                     MiscUtil.PrintWarning("Already scaned a maximum of %s torsion matches for molecule %s specified by \"--torsionMaxMatches\" option. Abandoning torsion scan...\n" % (TorsionMaxMatches, MolName))
  323                 break
  324             
  325             Mol, TorsionScanStatus,  TorsionMols, TorsionEnergies, TorsionAngles = ScanSingleTorsionInMol(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches,  MolNum)
  326             if not TorsionScanStatus:
  327                 continue
  328             
  329             TorsionScanCount +=  1
  330             GenerateOutputFiles(Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
  331         
  332         if TorsionMatchCount > TorsionMaxMatches:
  333             break
  334     
  335     if  OptionsInfo["RemoveHydrogens"]:
  336         Mol = Chem.RemoveHs(Mol)
  337 
  338     if TorsionScanCount:
  339         GenerateStartingTorsionScanStructureOutfile(Mol, MolNum)
  340 
  341     Status = True if TorsionScanCount else False
  342     
  343     return (Mol, Status)
  344     
  345 def ScanSingleTorsionInMol(Mol, TorsionID, TorsionPattern, TorsionPatternMol, TorsionMatches, MolNum = None):
  346     """Perform torsion scan for a molecule along with constrained energy minimization."""
  347 
  348     StartAngle = OptionsInfo["TorsionStart"]
  349     StopAngle = OptionsInfo["TorsionStop"]
  350     StepSize = OptionsInfo["TorsionStep"]
  351     TorsionMinimize = OptionsInfo["TorsionMinimize"]
  352     
  353     AtomIndex1, AtomIndex2, AtomIndex3, AtomIndex4 = TorsionMatches
  354 
  355     TorsionMols = []
  356     TorsionEnergies = []
  357     TorsionAngles = []
  358     
  359     Angles = [Angle for Angle in range(StartAngle, StopAngle, StepSize)]
  360     Angles.append(StopAngle)
  361 
  362     for Angle in Angles:
  363         TorsionMol = Chem.Mol(Mol)
  364         TorsionMolConf = TorsionMol.GetConformer(0)
  365 
  366         rdMolTransforms.SetDihedralDeg(TorsionMolConf, AtomIndex1, AtomIndex2, AtomIndex3, AtomIndex4, Angle)
  367 
  368         if TorsionMinimize:
  369             # Perform constrained minimization...
  370             TorsionMatchesMol = RDKitUtil.MolFromSubstructureMatch(TorsionMol, TorsionPatternMol, TorsionMatches)
  371             TorsionMol, CalcStatus, Energy = ConstrainAndMinimizeMolecule(TorsionMol, TorsionMatchesMol, MolNum)
  372             if not CalcStatus:
  373                 if not OptionsInfo["QuietMode"]:
  374                     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  375                     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, Angle, TorsionPattern))
  376                 return (Mol, False, None, None, None)
  377         else:
  378             # Calculate energy...
  379             CalcStatus, Energy = GetEnergy(TorsionMol)
  380             if not CalcStatus:
  381                 if not OptionsInfo["QuietMode"]:
  382                     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  383                     MiscUtil.PrintWarning("Failed to retrieve calculated energy for molecule %s with torsion angle set to %s during torsion scan for torsion pattern %s. Abandoning torsion scan..." % (MolName, Angle, TorsionPattern))
  384                 return (Mol, False, None, None, None)
  385         
  386         if  OptionsInfo["RemoveHydrogens"]:
  387             TorsionMol = Chem.RemoveHs(TorsionMol)
  388         
  389         TorsionMols.append(TorsionMol)
  390         TorsionEnergies.append(Energy)
  391         TorsionAngles.append(Angle)
  392     
  393     return (Mol, True, TorsionMols, TorsionEnergies, TorsionAngles)
  394 
  395 def SetupTorsionsMolInfo(Mol, MolNum = None):
  396     """Setup torsions info for a molecule."""
  397 
  398     TorsionsPatternsInfo = OptionsInfo["TorsionsPatternsInfo"]
  399     
  400     # Initialize...
  401     TorsionsMolInfo = {}
  402     TorsionsMolInfo["IDs"] = []
  403     TorsionsMolInfo["NumOfMatches"] = 0
  404     TorsionsMolInfo["Matches"] = {}
  405     for TorsionID in TorsionsPatternsInfo["IDs"]:
  406         TorsionsMolInfo["IDs"].append(TorsionID)
  407         TorsionsMolInfo["Matches"][TorsionID] = None
  408     
  409     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  410     UseChirality = OptionsInfo["UseChirality"]
  411     
  412     for TorsionID in TorsionsPatternsInfo["IDs"]:
  413         # Match torsions..
  414         TorsionPattern = TorsionsPatternsInfo["Pattern"][TorsionID]
  415         TorsionPatternMol = TorsionsPatternsInfo["Mol"][TorsionID]
  416         TorsionsMatches = RDKitUtil.FilterSubstructureMatchesByAtomMapNumbers(Mol, TorsionPatternMol, Mol.GetSubstructMatches(TorsionPatternMol, useChirality = UseChirality))
  417         
  418         # Validate tosion matches...
  419         ValidTorsionsMatches = []
  420         for Index, TorsionMatch in enumerate(TorsionsMatches):
  421             if len(TorsionMatch) != 4:
  422                 if not OptionsInfo["QuietMode"]:
  423                     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))
  424                 continue
  425 
  426             if not RDKitUtil.AreAtomIndicesSequentiallyConnected(Mol, TorsionMatch):
  427                 if not OptionsInfo["QuietMode"]:
  428                     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))
  429                 continue
  430 
  431             Bond = Mol.GetBondBetweenAtoms(TorsionMatch[1], TorsionMatch[2])
  432             if Bond.IsInRing():
  433                 if not OptionsInfo["QuietMode"]:
  434                     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]))
  435                 continue
  436             
  437             ValidTorsionsMatches.append(TorsionMatch)
  438         
  439         # Track valid matches...
  440         if len(ValidTorsionsMatches):
  441             TorsionsMolInfo["NumOfMatches"] += len(ValidTorsionsMatches)
  442             TorsionsMolInfo["Matches"][TorsionID] = ValidTorsionsMatches
  443         
  444     if TorsionsMolInfo["NumOfMatches"] == 0:
  445         if not OptionsInfo["QuietMode"]:
  446             MiscUtil.PrintWarning("Failed to match any torsions  in molecule %s" % (MolName))
  447 
  448     return TorsionsMolInfo
  449 
  450 def SetupTorsionsPatternsInfo():
  451     """Setup torsions patterns info."""
  452 
  453     TorsionsPatternsInfo = {}
  454     TorsionsPatternsInfo["IDs"] = []
  455     TorsionsPatternsInfo["Pattern"] = {}
  456     TorsionsPatternsInfo["Mol"] = {}
  457 
  458     TorsionID = 0
  459     for TorsionPattern in OptionsInfo["TorsionPatternsList"]:
  460         TorsionID += 1
  461         
  462         TorsionMol = Chem.MolFromSmarts(TorsionPattern)
  463         if TorsionMol is None:
  464             MiscUtil.PrintError("Failed to create torsion pattern molecule. The torsion SMILES/SMARTS pattern, \"%s\", specified using \"-t, --torsions\" option is not valid." % (TorsionPattern))
  465         
  466         TorsionsPatternsInfo["IDs"].append(TorsionID)
  467         TorsionsPatternsInfo["Pattern"][TorsionID] = TorsionPattern
  468         TorsionsPatternsInfo["Mol"][TorsionID] = TorsionMol
  469 
  470     OptionsInfo["TorsionsPatternsInfo"] = TorsionsPatternsInfo
  471     
  472 def MinimizeMolecule(Mol, MolNum = None):
  473     "Generate and minimize conformers for a molecule to get the lowest energy conformer."
  474 
  475     # Add hydrogens before  minimization. No need to remove them after minimization
  476     # as they will be used during torsion match and constrained minimization...
  477     #
  478     Mol = AddHydrogens(Mol)
  479     
  480     ConfIDs = EmbedMolecule(Mol, MolNum)
  481     if not len(ConfIDs):
  482         if not OptionsInfo["QuietMode"]:
  483             MolName = RDKitUtil.GetMolName(Mol, MolNum)
  484             MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s: Embedding failed...\n" % MolName)
  485         return (Mol, False)
  486 
  487     CalcEnergyMap = {}
  488     for ConfID in ConfIDs:
  489         try:
  490             if OptionsInfo["UseUFF"]:
  491                 Status = AllChem.UFFOptimizeMolecule(Mol, confId = ConfID, maxIters = OptionsInfo["MaxIters"])
  492             elif OptionsInfo["UseMMFF"]:
  493                 Status = AllChem.MMFFOptimizeMolecule(Mol, confId = ConfID, maxIters = OptionsInfo["MaxIters"], mmffVariant = OptionsInfo["MMFFVariant"])
  494             else:
  495                 MiscUtil.PrintError("Minimization couldn't be performed: Specified forcefield, %s, is not supported" % OptionsInfo["ForceField"])
  496         except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
  497             if not OptionsInfo["QuietMode"]:
  498                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  499                 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s:\n%s\n" % (MolName, ErrMsg))
  500             return (Mol, False)
  501         
  502         EnergyStatus, Energy = GetEnergy(Mol, ConfID)
  503         if not EnergyStatus:
  504             if not OptionsInfo["QuietMode"]:
  505                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  506                 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))
  507             return (Mol, False)
  508         
  509         if Status != 0:
  510             if not OptionsInfo["QuietMode"]:
  511                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  512                 MiscUtil.PrintWarning("Minimization failed to converge for conformation number %d of molecule %s in %d steps. Try using higher value for \"--maxIters\" option...\n" % (ConfID, MolName, OptionsInfo["MaxIters"]))
  513             
  514         CalcEnergyMap[ConfID] = Energy
  515     
  516     SortedConfIDs = sorted(ConfIDs, key = lambda ConfID: CalcEnergyMap[ConfID])
  517     MinEnergyConfID = SortedConfIDs[0]
  518         
  519     for ConfID in [Conf.GetId() for Conf in Mol.GetConformers()]:
  520         if ConfID == MinEnergyConfID:
  521             continue
  522         Mol.RemoveConformer(ConfID)
  523     
  524     # Set ConfID to 0 for MinEnergyConf...
  525     Mol.GetConformer(MinEnergyConfID).SetId(0)
  526 
  527     return (Mol, True)
  528 
  529 def ConstrainAndMinimizeMolecule(Mol, RefMolCore, MolNum = None):
  530     "Constrain and Minimize molecule."
  531 
  532     # Setup forcefield function to use for constrained minimization...
  533     ForceFieldFunction = None
  534     ForceFieldName = None
  535     if OptionsInfo["UseUFF"]:
  536         ForceFieldFunction = lambda mol, confId = -1 : AllChem.UFFGetMoleculeForceField(mol, confId = confId)
  537         ForeceFieldName = "UFF"
  538     else:
  539         ForceFieldFunction = lambda mol, confId = -1 : AllChem.MMFFGetMoleculeForceField(mol, AllChem.MMFFGetMoleculeProperties(mol, mmffVariant = OptionsInfo["MMFFVariant"]) , confId = confId)
  540         ForeceFieldName = "MMFF"
  541 
  542     if ForceFieldFunction is None:
  543         if not OptionsInfo["QuietMode"]:
  544             MiscUtil.PrintWarning("Failed to setup forcefield %s for molecule: %s\n" % (ForceFieldName, RDKitUtil.GetMolName(Mol, MolNum)))
  545         return (None, False, None)
  546         
  547     MaxConfs = OptionsInfo["MaxConfsTorsion"]
  548     EnforceChirality = OptionsInfo["EnforceChirality"]
  549     UseExpTorsionAnglePrefs = OptionsInfo["UseExpTorsionAnglePrefs"]
  550     UseBasicKnowledge = OptionsInfo["UseBasicKnowledge"]
  551     UseTethers = OptionsInfo["UseTethers"]
  552 
  553     CalcEnergyMap = {}
  554     MolConfsMap = {}
  555     ConfIDs = [ConfID for ConfID in range(0, MaxConfs)]
  556 
  557     for ConfID in ConfIDs:
  558         try:
  559             MolConf = Chem.Mol(Mol)
  560             AllChem.ConstrainedEmbed(MolConf, RefMolCore, useTethers = UseTethers, coreConfId = -1, randomseed = ConfID, getForceField = ForceFieldFunction, enforceChirality = EnforceChirality, useExpTorsionAnglePrefs = UseExpTorsionAnglePrefs, useBasicKnowledge = UseBasicKnowledge)
  561         except (ValueError, RuntimeError, Chem.rdchem.KekulizeException)  as ErrMsg:
  562             if not OptionsInfo["QuietMode"]:
  563                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  564                 MiscUtil.PrintWarning("Constrained embedding coupldn't  be performed for molecule %s:\n%s\n" % (RDKitUtil.GetMolName(Mol, MolNum), ErrMsg))
  565             return (None, False, None)
  566         
  567         EnergyStatus, Energy = GetEnergy(MolConf)
  568         
  569         if not EnergyStatus:
  570             if not OptionsInfo["QuietMode"]:
  571                 MolName = RDKitUtil.GetMolName(Mol, MolNum)
  572                 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))
  573             return (None, False, None)
  574         
  575         CalcEnergyMap[ConfID] = Energy
  576         MolConfsMap[ConfID] = MolConf
  577 
  578     SortedConfIDs = sorted(ConfIDs, key = lambda ConfID: CalcEnergyMap[ConfID])
  579     MinEnergyConfID = SortedConfIDs[0]
  580     
  581     MinEnergy = CalcEnergyMap[MinEnergyConfID]
  582     MinEnergyMolConf = MolConfsMap[MinEnergyConfID]
  583     
  584     MinEnergyMolConf.ClearProp('EmbedRMS')
  585     
  586     return (MinEnergyMolConf, True, MinEnergy)
  587 
  588 def GetEnergy(Mol, ConfID = None):
  589     "Calculate energy."
  590 
  591     Status = True
  592     Energy = None
  593 
  594     if ConfID is None:
  595         ConfID = -1
  596     
  597     if OptionsInfo["UseUFF"]:
  598         UFFMoleculeForcefield = AllChem.UFFGetMoleculeForceField(Mol, confId = ConfID)
  599         if UFFMoleculeForcefield is None:
  600             Status = False
  601         else:
  602             Energy = UFFMoleculeForcefield.CalcEnergy()
  603     elif OptionsInfo["UseMMFF"]:
  604         MMFFMoleculeProperties = AllChem.MMFFGetMoleculeProperties(Mol, mmffVariant = OptionsInfo["MMFFVariant"])
  605         MMFFMoleculeForcefield = AllChem.MMFFGetMoleculeForceField(Mol, MMFFMoleculeProperties, confId = ConfID)
  606         if MMFFMoleculeForcefield is None:
  607             Status = False
  608         else:
  609             Energy = MMFFMoleculeForcefield.CalcEnergy()
  610     else:
  611         MiscUtil.PrintError("Couldn't retrieve conformer energy: Specified forcefield, %s, is not supported" % OptionsInfo["ForceField"])
  612     
  613     return (Status, Energy)
  614 
  615 def EmbedMolecule(Mol, MolNum = None):
  616     "Embed conformations"
  617 
  618     ConfIDs = []
  619     
  620     MaxConfs = OptionsInfo["MaxConfs"]
  621     RandomSeed = OptionsInfo["RandomSeed"]
  622     EnforceChirality = OptionsInfo["EnforceChirality"]
  623     UseExpTorsionAnglePrefs = OptionsInfo["UseExpTorsionAnglePrefs"]
  624     UseBasicKnowledge = OptionsInfo["UseBasicKnowledge"]
  625     
  626     try:
  627         ConfIDs = AllChem.EmbedMultipleConfs(Mol, numConfs = MaxConfs, randomSeed = RandomSeed, enforceChirality = EnforceChirality, useExpTorsionAnglePrefs = UseExpTorsionAnglePrefs, useBasicKnowledge = UseBasicKnowledge)
  628     except ValueError as ErrMsg:
  629         if not OptionsInfo["QuietMode"]:
  630             MolName = RDKitUtil.GetMolName(Mol, MolNum)
  631             MiscUtil.PrintWarning("Embedding failed  for molecule %s:\n%s\n" % (MolName, ErrMsg))
  632         ConfIDs = []
  633     
  634     return ConfIDs
  635 
  636 def GenerateStartingTorsionScanStructureOutfile(Mol, MolNum):
  637     """Write out the structure of molecule used for starting tosion scan. """
  638     
  639     FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
  640     MolName = GetOutputFileMolName(Mol, MolNum)
  641     
  642     Outfile  = "%s_%s.%s" % (FileName, MolName, FileExt)
  643     
  644     # Set up a molecule writer...
  645     Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
  646     if Writer is None:
  647         MiscUtil.PrintWarning("Failed to setup a writer for output fie %s " % Outfile)
  648         return
  649     
  650     Writer.write(Mol)
  651     
  652     if Writer is not None:
  653         Writer.close()
  654 
  655 def GenerateOutputFiles(Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
  656     """Generate output files. """
  657     
  658     StructureOutfile, EnergyTextOutfile, PlotOutfile = SetupOutputFileNames(Mol, MolNum, TorsionID, TorsionMatchNum)
  659     
  660     GenerateScannedTorsionsStructureOutfile(StructureOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
  661     GenerateEnergyTextOutfile(EnergyTextOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
  662     GeneratePlotOutfile(PlotOutfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles)
  663 
  664 def GenerateScannedTorsionsStructureOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
  665     """Write out structures generated after torsion scan along with associated data."""
  666 
  667     # Set up a molecule writer...
  668     Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
  669     if Writer is None:
  670         MiscUtil.PrintWarning("Failed to setup a writer for output fie %s " % Outfile)
  671         return
  672     
  673     MolName = RDKitUtil.GetMolName(Mol, MolNum)
  674     
  675     for Index, TorsionMol in enumerate(TorsionMols):
  676         TorsionAngle = "%s" % TorsionAngles[Index]
  677         TorsionMol.SetProp("Torsion_Angle", TorsionAngle)
  678         
  679         TorsionEnergy = "%.2f" % TorsionEnergies[Index]
  680         TorsionMol.SetProp(OptionsInfo["EnergyLabel"], TorsionEnergy)
  681 
  682         TorsionMolName = "%s_Deg%s" % (MolName, TorsionAngle)
  683         TorsionMol.SetProp("_Name", TorsionMolName)
  684         
  685         Writer.write(TorsionMol)
  686         
  687     if Writer is not None:
  688         Writer.close()
  689     
  690 def GenerateEnergyTextOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
  691     """Write out torsion angles and energies."""
  692 
  693     # Setup a writer...
  694     Writer = open(Outfile, "w")
  695     if Writer is None:
  696         MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
  697     
  698     # Write headers...
  699     Writer.write("TorsionAngle,%s\n" % (OptionsInfo["EnergyLabel"]))
  700 
  701     for Index, TorsionAngle in enumerate(TorsionAngles):
  702         Writer.write("%d,%.2f\n" % (TorsionAngle, TorsionEnergies[Index]))
  703 
  704     if Writer is not None:
  705         Writer.close()
  706     
  707 def GeneratePlotOutfile(Outfile, Mol, MolNum, TorsionID, TorsionMatchNum, TorsionMols, TorsionEnergies, TorsionAngles):
  708     """Generate a plot corresponding to torsion angles and energies."""
  709 
  710     OutPlotParams = OptionsInfo["OutPlotParams"]
  711 
  712     # Initialize seaborn and matplotlib paramaters...
  713     if not OptionsInfo["OutPlotInitialized"]:
  714         OptionsInfo["OutPlotInitialized"] = True
  715         RCParams = {"figure.figsize":(OutPlotParams["Width"], OutPlotParams["Height"]),
  716                     "axes.titleweight": OutPlotParams["TitleWeight"],
  717                     "axes.labelweight": OutPlotParams["LabelWeight"]}
  718         sns.set(context = OutPlotParams["Context"], style = OutPlotParams["Style"], palette = OutPlotParams["Palette"], font = OutPlotParams["Font"], font_scale = OutPlotParams["FontScale"], rc = RCParams)
  719 
  720     # Create a new figure...
  721     plt.figure()
  722 
  723     # Draw plot...
  724     PlotType = OutPlotParams["Type"]
  725     if re.match("linepoint", PlotType, re.I):
  726         Axis = sns.lineplot(x = TorsionAngles, y = TorsionEnergies, marker = "o",  legend = False)
  727     elif re.match("scatter", PlotType, re.I):
  728         Axis = sns.scatterplot(x = TorsionAngles, y = TorsionEnergies, legend = False)
  729     elif re.match("line", PlotType, re.I):
  730         Axis = sns.lineplot(x = TorsionAngles, y = TorsionEnergies, legend = False)
  731     else:
  732         MiscUtil.PrintError("The value, %s, specified for \"type\" using option \"--outPlotParams\" is not supported. Valid plot types: linepoint, scatter or line" % (PlotType))
  733 
  734     # Setup title and labels...
  735     Title = OutPlotParams["Title"]
  736     if OptionsInfo["OutPlotTitleTorsionSpec"]:
  737         TorsionPattern = OptionsInfo["TorsionsPatternsInfo"]["Pattern"][TorsionID]
  738         Title = "%s: %s" % (OutPlotParams["Title"], TorsionPattern)
  739 
  740     # Set labels and title...
  741     Axis.set(xlabel = OutPlotParams["XLabel"], ylabel = OutPlotParams["YLabel"], title = Title)
  742     
  743     # Save figure...
  744     plt.savefig(Outfile)
  745 
  746     # Close the plot...
  747     plt.close()
  748 
  749 def SetupOutputFileNames(Mol, MolNum, TorsionID, TorsionMatchNum):
  750     """Setup names of output files. """
  751     
  752     FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
  753     MolName = GetOutputFileMolName(Mol, MolNum)
  754     
  755     OutfileRoot  = "%s_%s_Torsion%s_Match%s" % (FileName, MolName, TorsionID, TorsionMatchNum)
  756     
  757     StructureOutfile = "%s.%s" % (OutfileRoot, FileExt)
  758     EnergyTextOutfile = "%s_Energies.csv" % (OutfileRoot)
  759     PlotExt = OptionsInfo["OutPlotParams"]["OutExt"]
  760     PlotOutfile = "%s_Plot.%s" % (OutfileRoot, PlotExt)
  761 
  762     return (StructureOutfile, EnergyTextOutfile, PlotOutfile)
  763 
  764 def GetOutputFileMolName(Mol, MolNum):
  765     """Get output file prefix. """
  766     
  767     MolName = "Mol%s" % MolNum
  768     if OptionsInfo["OutfileMolName"]:
  769         MolName = re.sub("[^a-zA-Z0-9]", "_", RDKitUtil.GetMolName(Mol, MolNum), re.I)
  770 
  771     return MolName
  772 
  773 def AddHydrogens(Mol, AddCoords = True):
  774     """Check and add hydrogens. """
  775     
  776     if  not OptionsInfo["AddHydrogens"]:
  777         return Mol
  778 
  779     if OptionsInfo["HydrogensAdded"]:
  780         return Mol
  781 
  782     OptionsInfo["HydrogensAdded"] = True
  783     return Chem.AddHs(Mol, addCoords = AddCoords)
  784 
  785 def ProcessOptions():
  786     """Process and validate command line arguments and options."""
  787     
  788     MiscUtil.PrintInfo("Processing options...")
  789 
  790     # Validate options...
  791     ValidateOptions()
  792     
  793     OptionsInfo["ModeMols"] = Options["--modeMols"]
  794     OptionsInfo["FirstMolMode"] = True if re.match("^First$", Options["--modeMols"], re.I) else False
  795     
  796     OptionsInfo["ModeTorsions"] = Options["--modeTorsions"]
  797     OptionsInfo["FirstTorsionMode"] = True if re.match("^First$", Options["--modeTorsions"], re.I) else False
  798     
  799     OptionsInfo["Infile"] = Options["--infile"]
  800     OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"], Options["--infile"])
  801     OptionsInfo["Infile3D"] = True if re.match("^yes$", Options["--infile3D"], re.I) else False
  802     
  803     OptionsInfo["Outfile"] = Options["--outfile"]
  804     OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters("--outfileParams", Options["--outfileParams"])
  805 
  806     OptionsInfo["OutfileMolName"] = True if re.match("^yes$", Options["--outfileMolName"], re.I) else False
  807     
  808     OptionsInfo["OutPlotTitleTorsionSpec"] = True if re.match("^yes$", Options["--outPlotTitleTorsionSpec"], re.I) else False
  809     
  810     # The default width and height, 10.0 and 7.5, map to aspect raito of 16/9 (1.778)...
  811     DefaultValues = {'Type': 'linepoint', 'Width': 10.0, 'Height': 5.6, 'Title': 'Torsion Scan', 'XLabel': 'Torsion Angle (degrees)', 'YLabel': 'Energy (kcal/mol)'}
  812     OptionsInfo["OutPlotParams"] = MiscUtil.ProcessOptionSeabornPlotParameters("--outPlotParams", Options["--outPlotParams"], DefaultValues)
  813     if not re.match("^(linepoint|scatter|Line)$", OptionsInfo["OutPlotParams"]["Type"], re.I):
  814         MiscUtil.PrintError("The value, %s, specified for \"type\" using option \"--outPlotParams\" is not supported. Valid plot types: linepoint, scatter or line" % (OptionsInfo["OutPlotParams"]["Type"]))
  815     
  816     OptionsInfo["OutPlotInitialized"] = False
  817     
  818     # Procsss and validate specified SMILES/SMARTS torsion patterns...
  819     TorsionPatterns = Options["--torsions"]
  820     TorsionPatternsList = []
  821     for TorsionPattern in TorsionPatterns.split(","):
  822         TorsionPattern = TorsionPattern.strip()
  823         if not len(TorsionPattern):
  824             MiscUtil.PrintError("Empty value specified for SMILES/SMARTS pattern in  \"-t, --torsions\" option: %s" % TorsionPatterns)
  825         
  826         TorsionMol = Chem.MolFromSmarts(TorsionPattern)
  827         if TorsionMol is None:
  828             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))
  829         TorsionPatternsList.append(TorsionPattern)
  830     
  831     OptionsInfo["TorsionPatterns"] = TorsionPatterns
  832     OptionsInfo["TorsionPatternsList"] = TorsionPatternsList
  833     
  834     OptionsInfo["Overwrite"] = Options["--overwrite"]
  835     
  836     OptionsInfo["AddHydrogens"] = True if re.match("^yes$", Options["--addHydrogens"], re.I) else False
  837     OptionsInfo["HydrogensAdded"] = False
  838     
  839     if re.match("^ETDG$", Options["--conformerGenerator"], re.I):
  840         ConformerGenerator = "ETDG"
  841         UseExpTorsionAnglePrefs = True
  842         UseBasicKnowledge = False
  843     elif re.match("^KDG$", Options["--conformerGenerator"], re.I):
  844         ConformerGenerator = "KDG"
  845         UseExpTorsionAnglePrefs = False
  846         UseBasicKnowledge = True
  847     elif re.match("^ETKDG$", Options["--conformerGenerator"], re.I):
  848         ConformerGenerator = "ETKDG"
  849         UseExpTorsionAnglePrefs = True
  850         UseBasicKnowledge = True
  851     elif re.match("^SDG$", Options["--conformerGenerator"], re.I):
  852         ConformerGenerator = "SDG"
  853         UseExpTorsionAnglePrefs = False
  854         UseBasicKnowledge = False
  855     else:
  856         MiscUtil.PrintError("The value, %s, specified for option \"-c, --conformerGenerator\" is not supported." % (Options["--conformerGenerator"]))
  857     
  858     OptionsInfo["ConformerGenerator"] = ConformerGenerator
  859     OptionsInfo["UseExpTorsionAnglePrefs"] = UseExpTorsionAnglePrefs
  860     OptionsInfo["UseBasicKnowledge"] = UseBasicKnowledge
  861 
  862     if re.match("^UFF$", Options["--forceField"], re.I):
  863         ForceField = "UFF"
  864         UseUFF = True
  865         UseMMFF = False
  866     elif re.match("^MMFF$", Options["--forceField"], re.I):
  867         ForceField = "MMFF"
  868         UseUFF = False
  869         UseMMFF = True
  870     else:
  871         MiscUtil.PrintError("The value, %s, specified for \"--forceField\" is not supported." % (Options["--forceField"],))
  872     
  873     MMFFVariant = "MMFF94" if re.match("^MMFF94$", Options["--forceFieldMMFFVariant"], re.I) else "MMFF94s"
  874     
  875     OptionsInfo["ForceField"] = ForceField
  876     OptionsInfo["MMFFVariant"] = MMFFVariant
  877     OptionsInfo["UseMMFF"] = UseMMFF
  878     OptionsInfo["UseUFF"] = UseUFF
  879     
  880     if UseMMFF:
  881         OptionsInfo["EnergyLabel"] = "%s_Energy" % MMFFVariant
  882     else:
  883         OptionsInfo["EnergyLabel"] = "%s_Energy" % ForceField
  884     
  885     OptionsInfo["EnforceChirality"] = True if re.match("^yes$", Options["--enforceChirality"], re.I) else False
  886     
  887     OptionsInfo["MaxConfs"] = int(Options["--maxConfs"])
  888     OptionsInfo["MaxConfsTorsion"] = int(Options["--maxConfsTorsion"])
  889     OptionsInfo["MaxIters"] = int(Options["--maxIters"])
  890     
  891     OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
  892     OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
  893     
  894     OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
  895     
  896     RandomSeed = -1
  897     if not re.match("^auto$", Options["--randomSeed"], re.I):
  898         RandomSeed = int(Options["--randomSeed"])
  899     OptionsInfo["RandomSeed"] = RandomSeed
  900     
  901     OptionsInfo["RemoveHydrogens"] = True if re.match("^yes$", Options["--removeHydrogens"], re.I) else False
  902     
  903     OptionsInfo["TorsionMaxMatches"] = int(Options["--torsionMaxMatches"])
  904     OptionsInfo["TorsionMinimize"] = True if re.match("^yes$", Options["--torsionMinimize"], re.I) else False
  905 
  906     TorsionRange = Options["--torsionRange"]
  907     TorsionRangeWords = TorsionRange.split(",")
  908     
  909     TorsionStart = int(TorsionRangeWords[0])
  910     TorsionStop = int(TorsionRangeWords[1])
  911     TorsionStep = int(TorsionRangeWords[2])
  912     
  913     if TorsionStart >= TorsionStop:
  914         MiscUtil.PrintError("The start value, %d, specified for option \"--torsionRange\" in string \"%s\" must be less than stop value, %s." % (TorsionStart, Options["--torsionRange"], TorsionStop))
  915     if TorsionStep == 0:
  916         MiscUtil.PrintError("The step value, %d, specified for option \"--torsonRange\" in string \"%s\" must be > 0." % (TorsionStep, Options["--torsionRange"]))
  917     if TorsionStep >= (TorsionStop - TorsionStart):
  918         MiscUtil.PrintError("The step value, %d, specified for option \"--torsonRange\" in string \"%s\" must be less than, %s." % (TorsionStep, Options["--torsionRange"], (TorsionStop - TorsionStart)))
  919     
  920     if TorsionStart < 0:
  921         if TorsionStart < -180:
  922             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"]))
  923         if TorsionStop > 180:
  924             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"]))
  925     else:
  926         if TorsionStop > 360:
  927             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"]))
  928     
  929     OptionsInfo["TorsionRange"] = TorsionRange
  930     OptionsInfo["TorsionStart"] = TorsionStart
  931     OptionsInfo["TorsionStop"] = TorsionStop
  932     OptionsInfo["TorsionStep"] = TorsionStep
  933     
  934     OptionsInfo["UseTethers"] = True if re.match("^yes$", Options["--useTethers"], re.I) else False
  935     OptionsInfo["UseChirality"] = True if re.match("^yes$", Options["--useTethers"], re.I) else False
  936 
  937 def RetrieveOptions():
  938     """Retrieve command line arguments and options"""
  939     
  940     # Get options...
  941     global Options
  942     Options = docopt(_docoptUsage_)
  943     
  944     # Set current working directory to the specified directory...
  945     WorkingDir = Options["--workingdir"]
  946     if WorkingDir:
  947         os.chdir(WorkingDir)
  948     
  949     # Handle examples option...
  950     if "--examples" in Options and Options["--examples"]:
  951         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
  952         sys.exit(0)
  953 
  954 def ValidateOptions():
  955     """Validate option values"""
  956 
  957     MiscUtil.ValidateOptionTextValue("-a, --addHydrogens", Options["--addHydrogens"], "yes no")
  958     MiscUtil.ValidateOptionTextValue("-c, --conformerGenerator", Options["--conformerGenerator"], "SDG ETDG KDG ETKDG")
  959     
  960     MiscUtil.ValidateOptionTextValue("-f, --forceField", Options["--forceField"], "UFF MMFF")
  961     MiscUtil.ValidateOptionTextValue(" --forceFieldMMFFVariant", Options["--forceFieldMMFFVariant"], "MMFF94 MMFF94s")
  962     
  963     MiscUtil.ValidateOptionTextValue("--enforceChirality ", Options["--enforceChirality"], "yes no")
  964     
  965     MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
  966     MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi txt csv tsv")
  967     MiscUtil.ValidateOptionTextValue("--infile3D", Options["--infile3D"], "yes no")
  968 
  969     MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
  970     MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
  971     MiscUtil.ValidateOptionsDistinctFileNames("-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"])
  972     
  973     if not Options["--overwrite"]:
  974         FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
  975         FileNames = glob.glob("%s_*" % FileName)
  976         if len(FileNames):
  977             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"]))
  978     
  979     MiscUtil.ValidateOptionTextValue("--outPlotTitleTorsionSpec", Options["--outPlotTitleTorsionSpec"], "yes no")
  980     
  981     MiscUtil.ValidateOptionTextValue("--outfileMolName ", Options["--outfileMolName"], "yes no")
  982     
  983     MiscUtil.ValidateOptionTextValue("--modeMols", Options["--modeMols"], "First All")
  984     MiscUtil.ValidateOptionTextValue("--modeTorsions", Options["--modeTorsions"], "First All")
  985 
  986     MiscUtil.ValidateOptionIntegerValue("--maxConfs", Options["--maxConfs"], {">": 0})
  987     MiscUtil.ValidateOptionIntegerValue("--maxConfsTorsion", Options["--maxConfsTorsion"], {">": 0})
  988     MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0})
  989     
  990     MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
  991     MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
  992     
  993     if not re.match("^auto$", Options["--randomSeed"], re.I):
  994         MiscUtil.ValidateOptionIntegerValue("--randomSeed", Options["--randomSeed"], {})
  995     
  996     MiscUtil.ValidateOptionTextValue("-r, --removeHydrogens", Options["--removeHydrogens"], "yes no")
  997     
  998     MiscUtil.ValidateOptionIntegerValue("--torsionMaxMatches", Options["--torsionMaxMatches"], {">": 0})
  999     MiscUtil.ValidateOptionTextValue("--torsionMinimize", Options["--torsionMinimize"], "yes no")
 1000     MiscUtil.ValidateOptionNumberValues("--torsionRange", Options["--torsionRange"], 3, ",", "integer", {})
 1001     
 1002     MiscUtil.ValidateOptionTextValue("--useChirality", Options["--useChirality"], "yes no")
 1003     MiscUtil.ValidateOptionTextValue("--useTethers", Options["--useTethers"], "yes no")
 1004 
 1005 # Setup a usage string for docopt...
 1006 _docoptUsage_ = """
 1007 RDKitPerformTorsionScan.py - Perform torsion scan
 1008 
 1009 Usage:
 1010     RDKitPerformTorsionScan.py [--addHydrogens <yes or no>] [--conformerGenerator <SDG, ETDG, KDG, ETKDG>]
 1011                                [--forceField <UFF, or MMFF>] [--forceFieldMMFFVariant <MMFF94 or MMFF94s>]
 1012                                [--enforceChirality <yes or no>] [--infile3D <yes or no>] [--infileParams <Name,Value,...>]
 1013                                [--modeMols  <First or All>] [--modeTorsions  <First or All> ] [--maxConfs <number>]
 1014                                [--maxConfsTorsion <number>] [--maxIters <number>] [--mp <yes or no>] [--mpParams <Name.Value,...>]
 1015                                [--outfileMolName  <yes or no>] [--outfileParams <Name,Value,...>] [--outPlotParams <Name,Value,...>]
 1016                                [--outPlotTitleTorsionSpec <yes or no>] [--overwrite]  [--quiet <yes or no>] [--removeHydrogens <yes or no>]
 1017                                [--randomSeed <number>] [--torsionMaxMatches <number>] [--torsionMinimize <yes or no>]
 1018                                [--torsionRange <Start,Stop,Step>] [--useChirality <yes or no>] [--useTethers  <yes or no>]
 1019                                [-w <dir>] -t <torsions> -i <infile>  -o <outfile> 
 1020     RDKitPerformTorsionScan.py -h | --help | -e | --examples
 1021 
 1022 Description:
 1023     Perform torsion scan for molecules around torsion angles specified using
 1024     SMILES/SMARTS patterns. A molecule is optionally minimized before performing
 1025     a torsion scan. A set of initial 3D structures are generated for a molecule
 1026     by scanning the torsion angle across the specified range and updating the 3D
 1027     coordinates of the molecule. A conformation ensemble is optionally generated
 1028     for each 3D structure representing a specific torsion angle. The conformation
 1029     with the lowest energy is selected to represent the torsion angle. An option
 1030     is available to skip the generation of the conformation ensemble and simply
 1031     calculate the energy for the initial 3D structure for a specific torsion angle
 1032 
 1033     The torsions are specified using SMILES or SMARTS patterns. A substructure match
 1034     is performed to select torsion atoms in a molecule. The SMILES pattern match must
 1035     correspond to four torsion atoms. The SMARTS patterns containing atom indices may
 1036     match  more than four atoms. The atoms indices, however, must match exactly four
 1037     torsion atoms. For example: [s:1][c:2]([aX2,cH1])!@[CX3:3](O)=[O:4] for thiophene
 1038     esters and carboxylates as specified in Torsion Library (TorLib) [Ref 146].
 1039 
 1040     A set of four output files is generated for each torsion match in each
 1041     molecule. The names of the output files are generated using the root of
 1042     the specified output file. They may either contain sequential molecule
 1043     numbers or molecule names as shown below:
 1044         
 1045         <OutfileRoot>_Mol<Num>.sdf
 1046         <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>.sdf
 1047         <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Energies.csv
 1048         <OutfileRoot>_Mol<Num>_Torsion<Num>_Match<Num>_Plot.<ImgExt>
 1049         
 1050         or
 1051         
 1052         <OutfileRoot>_<MolName>.sdf
 1053         <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>.sdf
 1054         <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Energies.csv
 1055         <OutfileRoot>_<MolName>_Torsion<Num>_Match<Num>_Plot.<ImgExt>
 1056         
 1057     The supported input file formats are: Mol (.mol), SD (.sdf, .sd), .csv, .tsv .txt)
 1058 
 1059     The supported output file formats are: SD (.sdf, .sd)
 1060 
 1061 Options:
 1062     -a, --addHydrogens <yes or no>  [default: yes]
 1063         Add hydrogens before minimization.
 1064     -c, --conformerGenerator <SDG, ETDG, KDG, ETKDG>  [default: ETKDG]
 1065         Conformation generation methodology for generating initial 3D structure
 1066         of a molecule and conformation ensemble representing a specific torsion
 1067         angle. No conformation ensemble is generated for 'No' value of
 1068         '--torsionMinimize' option.
 1069         
 1070         Possible values: Standard Distance Geometry, (SDG), Experimental Torsion-angle
 1071         preference with Distance Geometry (ETDG), basic Knowledge-terms with Distance
 1072         Geometry (KDG),  and Experimental Torsion-angle preference along with basic
 1073         Knowledge-terms with Distance Geometry (ETKDG) [Ref 129] .
 1074     -f, --forceField <UFF, MMFF>  [default: MMFF]
 1075         Forcefield method to use for  energy minimization of initial 3D structure
 1076         of a molecule and conformation ensemble representing a specific torsion.
 1077         No conformation ensemble is generated during for 'No' value of '--torsionMinimze'
 1078         option and constrained energy minimization is not performed. Possible values:
 1079         Universal Force Field (UFF) [ Ref 81 ] or Merck Molecular Mechanics Force
 1080         Field [ Ref 83-87 ] .
 1081     --forceFieldMMFFVariant <MMFF94 or MMFF94s>  [default: MMFF94]
 1082         Variant of MMFF forcefield to use for energy minimization.
 1083     --enforceChirality <yes or no>  [default: Yes]
 1084         Enforce chirality for defined chiral centers during generation of conformers.
 1085     -e, --examples
 1086         Print examples.
 1087     -h, --help
 1088         Print this help message.
 1089     -i, --infile <infile>
 1090         Input file name.
 1091     --infile3D <yes or no>  [default: no]
 1092         Skip generation and minimization of initial 3D structures for molecules in
 1093         input file containing 3D coordinates.
 1094     --infileParams <Name,Value,...>  [default: auto]
 1095         A comma delimited list of parameter name and value pairs for reading
 1096         molecules from files. The supported parameter names for different file
 1097         formats, along with their default values, are shown below:
 1098             
 1099             SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes
 1100             
 1101             SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
 1102                 smilesTitleLine,auto,sanitize,yes
 1103             
 1104         Possible values for smilesDelimiter: space, comma or tab.
 1105     --modeMols <First or All>  [default: First]
 1106         Perform torsion scan for the first molecule or all molecules in input
 1107         file.
 1108     --modeTorsions <First or All>  [default: First]
 1109         Perform torsion scan for the first or all specified torsion pattern in
 1110         molecules up to a maximum number of matches for each torsion
 1111         specification as indicated by '--torsionMaxMatches' option. 
 1112     --maxConfs <number>  [default: 250]
 1113         Maximum number of conformations to generate for initial 3D structure of a
 1114         molecule. The lowest energy conformation is written to the output file.
 1115     --maxConfsTorsion <number>  [default: 50]
 1116         Maximum number of conformations to generate for conformation ensemble
 1117         representing a specific torsion. A constrained minimization is performed
 1118         using the coordinates of the specified torsion and the lowest energy
 1119         conformation is written to the output file.
 1120     --maxIters <number>  [default: 500]
 1121         Maximum number of iterations to perform for a molecule during minimization
 1122         to generation initial 3D structures. This option is ignored during 'yes' value
 1123         of  '--infile3D' option.
 1124     --mp <yes or no>  [default: no]
 1125         Use multiprocessing.
 1126          
 1127         By default, input data is retrieved in a lazy manner via mp.Pool.imap()
 1128         function employing lazy RDKit data iterable. This allows processing of
 1129         arbitrary large data sets without any additional requirements memory.
 1130         
 1131         All input data may be optionally loaded into memory by mp.Pool.map()
 1132         before starting worker processes in a process pool by setting the value
 1133         of 'inputDataMode' to 'InMemory' in '--mpParams' option.
 1134         
 1135         A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
 1136         data mode may adversely impact the performance. The '--mpParams' section
 1137         provides additional information to tune the value of 'chunkSize'.
 1138     --mpParams <Name,Value,...>  [default: auto]
 1139         A comma delimited list of parameter name and value pairs for to
 1140         configure multiprocessing.
 1141         
 1142         The supported parameter names along with their default and possible
 1143         values are shown below:
 1144         
 1145             chunkSize, auto
 1146             inputDataMode, Lazy   [ Possible values: InMemory or Lazy ]
 1147             numProcesses, auto   [ Default: mp.cpu_count() ]
 1148         
 1149         These parameters are used by the following functions to configure and
 1150         control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
 1151         mp.Pool.imap().
 1152         
 1153         The chunkSize determines chunks of input data passed to each worker
 1154         process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
 1155         The default value of chunkSize is dependent on the value of 'inputDataMode'.
 1156         
 1157         The mp.Pool.map() function, invoked during 'InMemory' input data mode,
 1158         automatically converts RDKit data iterable into a list, loads all data into
 1159         memory, and calculates the default chunkSize using the following method
 1160         as shown in its code:
 1161         
 1162             chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
 1163             if extra: chunkSize += 1
 1164         
 1165         For example, the default chunkSize will be 7 for a pool of 4 worker processes
 1166         and 100 data items.
 1167         
 1168         The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
 1169         'lazy' RDKit data iterable to retrieve data as needed, without loading all the
 1170         data into memory. Consequently, the size of input data is not known a priori.
 1171         It's not possible to estimate an optimal value for the chunkSize. The default 
 1172         chunkSize is set to 1.
 1173         
 1174         The default value for the chunkSize during 'Lazy' data mode may adversely
 1175         impact the performance due to the overhead associated with exchanging
 1176         small chunks of data. It is generally a good idea to explicitly set chunkSize to
 1177         a larger value during 'Lazy' input data mode, based on the size of your input
 1178         data and number of processes in the process pool.
 1179         
 1180         The mp.Pool.map() function waits for all worker processes to process all
 1181         the data and return the results. The mp.Pool.imap() function, however,
 1182         returns the the results obtained from worker processes as soon as the
 1183         results become available for specified chunks of data.
 1184         
 1185         The order of data in the results returned by both mp.Pool.map() and 
 1186         mp.Pool.imap() functions always corresponds to the input data.
 1187     -o, --outfile <outfile>
 1188         Output file name. The output file root is used for generating the names
 1189         of the output files corresponding to structures, energies, and plots during
 1190         the torsion scan.
 1191     --outfileMolName <yes or no>  [default: no]
 1192         Append molecule name to output file root during the generation of the names
 1193         for output files. The default is to use <MolNum>. The non alphabetical
 1194         characters in molecule names are replaced by underscores.
 1195     --outfileParams <Name,Value,...>  [default: auto]
 1196         A comma delimited list of parameter name and value pairs for writing
 1197         molecules to files. The supported parameter names for different file
 1198         formats, along with their default values, are shown below:
 1199             
 1200             SD: kekulize,no
 1201             
 1202     --outPlotParams <Name,Value,...>  [default: auto]
 1203         A comma delimited list of parameter name and value pairs for generating
 1204         plots using Seaborn module. The supported parameter names along with their
 1205         default values are shown below:
 1206             
 1207             type,linepoint,outExt,svg,width,10,height,5.6,
 1208             title,auto,xlabel,auto,ylabel,auto,titleWeight,bold,labelWeight,bold
 1209             style,darkgrid,palette,deep,font,sans-serif,fontScale,1,
 1210             context,notebook
 1211             
 1212         Possible values:
 1213             
 1214             type: linepoint, scatter, or line. Both points and lines are drawn
 1215                 for linepoint plot type.
 1216             outExt: Any valid format supported by Python module Matplotlib.
 1217                 For example: PDF (.pdf), PNG (.png), PS (.ps), SVG (.svg)
 1218             titleWeight, labelWeight: Font weight for title and axes labels.
 1219                 Any valid value.
 1220             style: darkgrid, whitegrid, dark, white, ticks
 1221             palette: deep, muted, pastel, dark, bright, colorblind
 1222             font: Any valid font name
 1223         
 1224      --outPlotTitleTorsionSpec <yes or no>  [default: yes]
 1225         Append torsion specification to the title of the torsion plot.
 1226     --overwrite
 1227         Overwrite existing files.
 1228     -q, --quiet <yes or no>  [default: no]
 1229         Use quiet mode. The warning and information messages will not be printed.
 1230     --randomSeed <number>  [default: auto]
 1231         Seed for the random number generator for generating initial 3D coordinates.
 1232         Default is to use a random seed.
 1233     --removeHydrogens <yes or no>  [default: Yes]
 1234         Remove hydrogens after minimization.
 1235     -t, --torsions <SMILES/SMARTS,...,...>
 1236         SMILES/SMARTS patterns corresponding to torsion specifications. It's a 
 1237         comma delimited list of valid SMILES/SMART patterns.
 1238         
 1239         A substructure match is performed to select torsion atoms in a molecule.
 1240         The SMILES pattern match must correspond to four torsion atoms. The
 1241         SMARTS patterns contain atom indices may match  more than four atoms.
 1242         The atoms indices, however, must match exactly four torsion atoms. For example:
 1243         [s:1][c:2]([aX2,cH1])!@[CX3:3](O)=[O:4] for thiophene esters and carboxylates
 1244         as specified in Torsion Library (TorLib) [Ref 146].
 1245     --torsionMaxMatches <number>  [default: 5]
 1246         Maximum number of torsions to match for each torsion specification in a
 1247         molecule.
 1248     --torsionMinimize <yes or no>  [default: no]
 1249         Perform constrained energy minimization on a conformation ensemble
 1250         for  a specific torsion angle and select the lowest energy conformation
 1251         representing the torsion angle.
 1252     --torsionRange <Start,Stop,Step>  [default: 0,360,5]
 1253         Start, stop, and step size angles in degrees for a torsion scan. In addition,
 1254         you may specify values using start and stop angles from -180 to 180.
 1255     --useChirality <yes or no>  [default: no]
 1256         Use chirrality during substructure matches for identification of torsions.
 1257      --useTethers <yes or no>  [default: yes]
 1258         Use tethers to optimize the final conformation by applying a series of extra forces
 1259         to align matching atoms to the positions of the core atoms. Otherwise, use simple
 1260         distance constraints during the optimization.
 1261     -w, --workingdir <dir>
 1262         Location of working directory which defaults to the current directory.
 1263 
 1264 Examples:
 1265     To perform a torsion scan on first molecule in a SMILES file using a minimum
 1266     energy structure of the molecule selected from an ensemble of conformations,
 1267     skipping generation of conformation ensembles for specific torsion angles and
 1268     constrained energy minimization of the ensemble, generate output files
 1269     corresponding to structure, energy and torsion plot, type:
 1270     
 1271         % RDKitPerformTorsionScan.py  -t "O=CNC" -i SampleSeriesD3R.smi 
 1272           -o SampleOut.sdf
 1273 
 1274     To run the previous example on all molecules in a SD file, type:
 1275     
 1276         % RDKitPerformTorsionScan.py  -t "O=CNC" --modeMols All
 1277           -i SampleSeriesD3R.sdf -o SampleOut.sdf
 1278 
 1279     To perform a torsion scan on first molecule in a SMILES file using a minimum
 1280     energy structure of the molecule selected from an ensemble of conformations,
 1281     generation of conformation ensembles for specific torsion angles and constrained
 1282     energy minimization of the ensemble, generate output files corresponding to
 1283     structure, energy and torsion plot, type:
 1284     
 1285         % RDKitPerformTorsionScan.py  -t "O=CNC" --torsionMinimize Yes
 1286            -i SampleSeriesD3R.smi -o SampleOut.sdf
 1287 
 1288     To run the previous example on all molecules in a SD file, type:
 1289     
 1290         % RDKitPerformTorsionScan.py  -t "O=CNC" --modeMols All
 1291            --torsionMinimize Yes -i SampleSeriesD3R.sdf -o SampleOut.sdf
 1292 
 1293     To run the previous example in multiprocessing mode on all available CPUs
 1294     without loading all data into memory and write out a SD file, type:
 1295 
 1296         % RDKitPerformTorsionScan.py  -t "O=CNC" -i SampleSeriesD3R.smi 
 1297           -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
 1298     
 1299     To run the previous example in multiprocessing mode on all available CPUs
 1300     by loading all data into memory and write out a SD file, type:
 1301 
 1302         % RDKitPerformTorsionScan.py  -t "O=CNC" -i SampleSeriesD3R.smi 
 1303           -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
 1304           --mpParams "inputDataMode,InMemory"
 1305     
 1306     To run the previous example in multiprocessing mode on specific number of
 1307     CPUs and chunk size without loading all data into memory and write out a SD file,
 1308     type:
 1309 
 1310         % RDKitPerformTorsionScan.py  -t "O=CNC" -i SampleSeriesD3R.smi 
 1311           -o SampleOut.sdf --modeMols All --torsionMinimize Yes --mp yes
 1312           --mpParams "inputDataMode,Lazy,numProcesses,4,chunkSize,8"
 1313 
 1314     To perform a torsion scan on first molecule in a SD file containing 3D coordinates,
 1315     skipping generation of conformation ensembles for specific torsion angles and
 1316     constrained energy minimization of the ensemble, generate output files
 1317     corresponding to structure, energy and torsion plot, type:
 1318     
 1319         % RDKitPerformTorsionScan.py  -t "O=CNC"  --infile3D yes
 1320           -i SampleSeriesD3R3D.sdf -o SampleOut.sdf
 1321 
 1322     To perform a torsion scan using multiple torsion specifications on all molecules in
 1323     a SD file containing 3D coordinates, generation of conformation ensembles for specific
 1324     torsion angles and constrained energy minimization of the ensemble, generate output files
 1325     corresponding to structure, energy and torsion plot, type:
 1326     
 1327         % RDKitPerformTorsionScan.py  -t "O=CNC,[O:1]=[C:2](c)[N:3][C:4]"
 1328           --infile3D yes --modeMols All  --modeTorsions All
 1329           --torsionMinimize Yes -i SampleSeriesD3R3D.sdf -o SampleOut.sdf
 1330 
 1331     To run the previous example using a specific torsion scan range, type:
 1332 
 1333         % RDKitPerformTorsionScan.py  -t "O=CNC,[O:1]=[C:2](c)[N:3][C:4]"
 1334           --infile3D yes --modeMols All --modeTorsions All --torsionMinimize
 1335           Yes --torsionRange 0,360,10 -i SampleSeriesD3R.smi -o SampleOut.sdf
 1336 
 1337 Author:
 1338     Manish Sud(msud@san.rr.com)
 1339 
 1340 See also:
 1341     RDKitCalculateRMSD.py, RDKitCalculateMolecularDescriptors.py, RDKitCompareMoleculeShapes.py,
 1342     RDKitConvertFileFormat.py, RDKitPerformConstrainedMinimization.py
 1343 
 1344 Copyright:
 1345     Copyright (C) 2020 Manish Sud. All rights reserved.
 1346 
 1347     The functionality available in this script is implemented using RDKit, an
 1348     open source toolkit for cheminformatics developed by Greg Landrum.
 1349 
 1350     This file is part of MayaChemTools.
 1351 
 1352     MayaChemTools is free software; you can redistribute it and/or modify it under
 1353     the terms of the GNU Lesser General Public License as published by the Free
 1354     Software Foundation; either version 3 of the License, or (at your option) any
 1355     later version.
 1356 
 1357 """
 1358 
 1359 if __name__ == "__main__":
 1360     main()