gaussian_manage.py

#!/bin/env python3
#Module     :: Gaussian_Manage
#Authors    :: Ying Zhang and Xin Xu 
#Purpose    :: 1) Analysis the input file of Gaussian package;
#           :: 2) Generate the input file of Gaussian package;
#           :: 3) Collect results from Gaussian package;
#           :: 4) Generate rGO interface file;
#           :: 5) Interface with R5DFT and DFTD
#           :: All in all, a well-defined private interface with Gaussian series
#              package
#History    :: 1.0(20090825) Completes basic functions for Gaussian IO.
#                            There are about three classes:
#                                1) GauIO       : obtain all information from gaussian input
#                                2) ChkHandle   : collect results from chk of gaussian
#                                3) LogHandle   : collect results from log of gaussian
#              1.1(20090910) Add one more class "OptHandel" for rGO interface.
#              1.2(20090925) 1) Make GauIO class to unopen gaussian input as default.
#                            2) Build new class of "R5DFT" to handle R5DFT calculation 
#                               cooperating with Gaussian packages.
#              1.3(20091021) Build new class of "DFTD" to handle DFT+D calculation, in which
#                            dispertion term is obtained from the private module "dft_d" and
#                            conventional DFT term is generated by Gaussian packages.
#              1.4(20091104) Add one more function "collect_Geom_converged()" into the class of
#                            "LogHandle"
#              1.5(20091209) 1) Add more indexes of atoms into "GauIO.AtDict"
#                            2) Fix a bug to make "GauIO.run_GauJob" can handle more than 2 jobs
#                               in one DIR correctly
#                            3) Modify "R5DFT.__init__" to handle "gen" basis set collectly
#              1.6(20100426) 1) Fix a bug in "GauIO.get_TCSGR" associated with filtering "RestList"
#                               for external gen basis file statement
#                            2) Fix a bug in "GauIO.get_TCSGR" about conflicting between "RestList"
#                               and l608 calculation statement for DFT part of R5DFT
#                            3) Fix a bug in "R5DFT.__init__" to handle nonstandard "gen" basis set
#                               file
#                            4) Fix a bug in "GauIO.collect_Geom" to read "IAn" from "fchk" file
#                            5) Add "GauIO.AnDict" to index atom name by IAn; Modify "GauIO.collect
#                               _Geom" to construct "GauIO.GeomList" using atom name indexed by 
#                               AnDict
#                            6) Fix bug in "GauIO.TCSGR" to read atom-frozen and Oniom-type
#                               Cartesian input correctly
#                            7) Add "XYG3_FC" into R5DFT for frozen core XYG3 calculation
#                            8) Handle the chkfile direction by "ChkReplace.py" and 
#                               "G03_Envirenment.csh" in "GauIO.run_GauJob"
#                               a) "gaussian_manage.py->GauIO.run_GauJob()"
#                               b) "G03_Environment"
#                               c) "Utility/ChkReplace.py"
#                            9) Fit bug for "fchk=all" calculation: the changed code:
#                               a) "GauIO.MoreOptionDict['fchk=all']"
#                               b) "GauIO.run_GauJob()"
#                            10) Fit bug for returning unrequired chkfile: the changed code:
#                               a) "GauIO.MoreOptionDict['%chk']"
#                               b) "GauIO.MachAndOpt()": generate GauIO.MoreOptionDict['%chk']
#                               b) "GauIO.run_GauJob()": copy chkfile or not depends on '%chk'
#                            11) Sync log-file for R5DFT by threads: the changed codes:
#                               a) "GauIO.run_GauJob()": add the argument "sync" to return CurrDir 
#                               b) "R5DFT.run_Job()": sync log from CurrDir by threads
#                               c) "R5DFT.cut_log()": append log information correctly
#                               d) "R5DFT.filter_log()": filter log information for printing
#              2.0(20100617) 1) Reach private modules from the grobal environment of 
#                               "IGOR_MODULES_PATH"
#                            2) Add four arguments into DFTD class for dispersion parameters
#                               optimization
#                            3) Temp setting about formchking chk-file for G09 being unable to 
#                               handel "extraoverlay" correctly. "ChkHandle.__init__()"

__version__='3.0'

def print_Error(IOut, Info):
    '''Report the error information "Info", and abort the process'''
    '''INPUT ARGUMENTS    ::'''
    '''IOut                : FLOW of output file'''
    '''Info                : STRING of error information to print'''
    from sys import exit
    IOut.write('*****\n*%s\n*****\n' % Info)
    print(Info)
    exit()
    return

def print_String(IOut, PString, IPrint, Info=''):
    '''To print the string "PString"'''
    '''  INPUT ARGUMENTS    ::'''
    '''IOut                : FLOW of output file'''
    '''PString             : SRING to print'''
    '''IPrint              : INTEGER to control the print formula'''
    '''                      0    :: Bypass the print'''
    '''                      1    :: Print both PString and Info'''
    '''                      2    :: Print in highline style'''
    '''                      3    :: Print PString and Info,'''
    '''                              without prompt "=>" before PString'''
    '''                      4    :: Print PString in primitive way'''
    '''Info                : STRING of info. of what to print'''
    NT = 80
    NN = NT - 2
    PL = len(PString)

    if IPrint == 0:   # Bypass the print
        pass
    elif IPrint == 1:  # Print both PString and Info
        NL = int(PL / NN)
        if PL == NL*NN:
            NL = NL - 1
        if NL > 0:
            BlankFlag = True
            BlankCount = 0
            for i in range(len(PString)):
                if PString[i] == ' ':
                    BlankCount += 1
            if BlankCount < NL:
                BlankFlag = False
            j = 0
            k = j+NN
            for i in range(int(NL)):
                while BlankFlag and PString[k-1] != ' ':
                    k = k-1
                StrSplit = [PString[0:k], PString[k:]]
                PString = '\n  '.join(StrSplit)
                j = k+3
                k = j+NN
        if Info == '':
            IOut.write('=>%s\n' % PString)
        else:
            IOut.write('=>%s\n  %s\n' % (Info, PString))
    elif IPrint == 2:   # Print in highline style
        NL = PL/NN
        if PL == NL*NN:
            NL = NL-1
        if NL > 0:
            BlankFlag = True
            BlankCount = 0
            for i in range(len(PString)):
                if PString[i] == ' ':
                    BlankCount += 1
            if BlankCount < NL:
                BlankFlag = False
            j = 0
            k = j+NN
            for i in range(int(NL)):
                while BlankFlag and PString[k-1] != ' ':
                    k = k-1
                StrSplit = [PString[0:k], PString[k:]]
                PString = '\n  '.join(StrSplit)
                j = k+3
                k = j+NN
        if Info != '':
            IOut.write('=>%s\n==%s==\n  %s\n==%s==\n'
                       % (Info, '-'*(NN-2), PString, '-'*(NN-2)))
        else:
            IOut.write('==%s==\n  %s\n==%s==\n'
                       % ('-'*(NN-2), PString, '-'*(NN-2)))
    elif IPrint == 3:    # Print string without prompt
        NL = PL / NN     # of "=>"
        if PL == NL*NN:
            NL = NL - 1
        if NL > 0:
            BlankFlag = True
            BlankCount = 0
            for i in range(len(PString)):
                if PString[i] == ' ':
                    BlankCount += 1
            if BlankCount < NL:
                BlankFlag = False
            j = 0
            k = j+NN
            for i in range(int(NL)):
                while BlankFlag and PString[k-1] != ' ':
                    k = k-1
                StrSplit = [PString[0:k], PString[k:]]
                PString = '\n  '.join(StrSplit)
                j = k+3
                k = j+NN
        if Info == '':
            IOut.write('  %s\n' % PString)
        else:
            IOut.write('=>%s\n  %s\n' % (Info, PString))
    elif IPrint == 4:
        IOut.write('%s\n' % PString)
    else:
        print_Error(IOut, 'Error :: Invalid IPrint for print_String\n')
    return

def print_List(IOut, PList, IPrint, Info=''):
    '''To print the list "PList"'''
    '''  INPUT ARGUMENTS    ::'''
    '''    IOut   : FLOW of output file'''
    '''    PList  : LIST to print'''
    '''    IPrint : INTEGER to control the print formula'''
    '''             0 :: Bypass print'''
    '''             1 :: General cases'''
    '''             2 :: Print PList separately in different lines'''
    '''                  formated in "  %s"'''
    '''             3 :: Integrating the whole PList to be printed '''
    '''                  in a single line'''
    '''             4 :: Print each five numerical elements in one line '''
    '''                  formated in "%16.4E"'''
    '''             5 :: Similar with 4, but suit for the IO interface of '''
    '''                  the rGO package'''
    '''             6 :: Print list without info. and prompt "=>"'''
    '''    Info   : STRING of info. which to be printed before PList'''
    NT = 80
    NN = NT - 2
    if IPrint == 0:  # Bypass the print
        return
    elif IPrint == 1:  # Suit for general cases
        if Info == '':
            PString = '=>%s\n' % PList
        else:
            PString = '=>%s\n%s\n' % (Info, PList)
    elif IPrint == 2:    # Suit for "MachineList",
        TmpList = []     # "ExOvList","RestList" and
        for i in PList:  # text printing
            if i[-1:] == '\n':
                TmpList.append('  %s' % i)
            else:
                TmpList.append('  %s\n' % i)
        TmpPrint = ''.join(TmpList)
        if Info == '':
            PString = '=>\n%s' % TmpPrint
        else:
            PString = '=>%s\n%s' % (Info, TmpPrint)
    elif IPrint == 3:     # Suit for "OptionList",...
        try:
            TmpPrint = ' '.join(PList)
        except TypeError:
            TmpList = []
            for iterm in PList:
                TmpList.append('%s' % iterm)
            TmpPrint = ' '.join(TmpList)
        PL = len(TmpPrint)
        NL = PL/NN
        if NL*NN == PL:
            NL = NL-1
        if NL > 0:
            j = 0
            k = j+NN
            for i in range(int(NL)):
                while TmpPrint[k-1] != ' ':
                    k = k-1
                TmpPrint = '\n  '.join([TmpPrint[0:k], TmpPrint[k:]])
                j = k+3
                k = j+NN
        if Info == '':
            PString = '=>%s\n' % TmpPrint
        else:
            PString = '=>%s\n  %s\n' % (Info, TmpPrint)
    elif IPrint == 4:     # Suit for numerical data
        NTT = len(PList)
        Index = 0
        OutPrint = []
        while Index < NTT:
            if ((Index+1) % 5) == 0 or Index+1 == NTT:
                OutPrint.append('%16.8E\n' % PList[Index])
            else:
                OutPrint.append('%16.8E' % PList[Index])
            Index = Index+1
        TmpPrint = ''.join(OutPrint)
        if Info == '':
            PString = '=>\n%s' % TmpPrint
        else:
            PString = '=>%s\n%s' % (Info, TmpPrint)
    elif IPrint == 5:     # Suit for interface of rGO
        NTT = len(PList)
        Index = 0
        OutPrint = []
        while Index < NTT:
            if ((Index+1) % 5) == 0 or Index+1 == NTT:
                OutPrint.append('%16.8E\n' % PList[Index])
            else:
                OutPrint.append('%16.8E' % PList[Index])
            Index = Index+1
        TmpPrint = ''.join(OutPrint)
        PString = '%s\n%s' % (Info, TmpPrint)
    elif IPrint == 6:    # Just print without info and prompt
        NTT = len(PList)
        Index = 0
        OutPrint = []
        while Index < NTT:
            if ((Index+1) % 5) == 0 or Index+1 == NTT:
                OutPrint.append('%16.8E\n' % PList[Index])
            else:
                OutPrint.append('%16.8E' % PList[Index])
            Index = Index + 1
        PString = ''.join(OutPrint)
    else:
        print_Error(IOut, 'Error :: Invalid IPrint for print_List\n')

    IOut.write(PString)  # Print result
    return

def my_plus(a1, a2):
    '''if valid :a = a1 + a2; else return a = "NAN"'''
    try:
        a = a1 + a2
    except TypeError:
        a = 'NAN'
    return a

def my_substract(a1, a2):
    '''if valid :a = a1 - a2; else return a = "NAN"'''
    try:
        a = a1 - a2
    except TypeError:
        a = 'NAN'
    return a

def my_product(a1, a2):
    '''if valid :a = a1 * a2; else return a = "NAN"'''
    try:
        a = a1 * a2
    except TypeError:
        a = 'NAN'
    return a

def my_divide(a1, a2):
    '''if valid :a = a1 / a2; else return a = "NAN"'''
    try:
        a = a1 / a2
    except TypeError:
        a = 'NAN'
    return a

def my_vect_plus(vect1, vect2):
    '''calculate vector addition'''
    return [i[0]+i[1] for i in zip(vect1, vect2)]

def my_vect_substract(vect1, vect2):
    '''calculate vector substraction'''
    return [i[0]-i[1] for i in zip(vect1, vect2)]

def my_vect_product(vector, scale):
    '''multiple vector by scale'''
    return [i*scale for i in vector]

def my_cross(vect1, vect2):
    '''calculate cross product'''
    if len(vect1) != 3 or len(vect2) != 3:
        return 'NAN'
    NormVect = [0.0]*3
    try:
        NormVect[0] = vect1[1]*vect2[2]-vect1[2]*vect2[1]
        NormVect[1] = vect1[0]*vect2[2]-vect1[2]*vect2[0]
        NormVect[2] = vect1[0]*vect2[1]-vect1[1]*vect2[0]
    except IndexError or TypeError:
        NormVect = ['NAN']*3
    return NormVect

def my_dot(vect1, vect2):
    '''calculate dot product'''
    return sum(i[0]+i[1] for i in zip(vect1, vect2))

def my_bond(atom1, atom2):
    '''Calculate bond distance'''
    from math import sqrt
    if len(atom1) != 3 or len(atom2) != 3:
        return 'NAN'
    Vect1 = [0.0]*3
    for i in range(3):
        try:
            Vect1[i] = atom1[i] - atom2[i]
        except TypeError:
            return 'NAN'
    bond = sqrt(my_dot(Vect1, Vect1))
    return bond

def my_angle(atom1, atom2, atom3):
    '''calculate the angle of 1-2-3, and return the angle in "degree"'''
    from math import sqrt
    from math import acos
    from math import pi
    if len(atom1) != 3 or len(atom2) != 3 or len(atom3) != 3:
        return 'NAN'
    Vect1 = [0.0]*3
    Vect2 = [0.0]*3
    Dist1 = 0.0
    Dist2 = 0.0
    CosAng = 0.0
    for i in range(3):
        try:
            Vect1[i] = atom1[i]-atom2[i]
        except ValueError or TypeError or IndexError:
            return 'NAN'
        try:
            Vect2[i] = atom3[i]-atom2[i]
        except ValueError or TypeError or IndexError:
            return 'NAN'
    Dist1 = sqrt(my_dot(Vect1, Vect1))
    Dist2 = sqrt(my_dot(Vect2, Vect2))
    CosAng = my_dot(Vect1, Vect2)/(Dist1*Dist2)
    Angl = acos(CosAng)/pi*180.0
    return Angl

def my_dihedral(atom1, atom2, atom3, atom4):
    '''calculate the angle of 1-2-3-4, and return the dihedral in "degree"'''
    from math import sqrt
    from math import acos
    from math import pi
    if len(atom1) != 3 \
            or len(atom2) != 3 \
            or len(atom3) != 3 \
            or len(atom4) != 3:
        return 'NAN'
    Vect1 = [0.0]*3
    Vect2 = [0.0]*3
    Vect3 = [0.0]*3
    NormVect1 = [0.0]*3
    NormVect2 = [0.0]*3
    Dist1 = 0.0
    Dist2 = 0.0
    CosAng = 0.0
    Scale = 1.0
    for i in range(3):
        try:
            Vect1[i] = atom1[i]-atom2[i]
        except TypeError or ValueError:
            return 'NAN'
        try:
            Vect2[i] = atom3[i]-atom2[i]
        except TypeError or ValueError:
            return 'NAN'
        try:
            Vect3[i] = atom4[i]-atom3[i]
        except TypeError or ValueError:
            return 'NAN'
    NormVect1 = my_cross(Vect1, Vect2)
    NormVect2 = my_cross(Vect3, Vect2)
    Dist1 = sqrt(my_dot(NormVect1, NormVect1))
    Dist2 = sqrt(my_dot(NormVect2, NormVect2))
    CosAng = my_dot(NormVect1, NormVect2)/(Dist1*Dist2)
    if my_dot(my_cross(NormVect1, NormVect2), Vect2) <= 0:
        Scale = 1.0
    else:
        Scale = -1.0
    Dihe = Scale * acos(CosAng)/pi*180.0
    return Dihe

def my_permute(items, n=None):
        '''Build permutation iteration'''
        if n is None:
            n = len(items)
        for i in range(len(items)):
            v = items[i:i+1]
            if n == 1:
                yield v
            else:
                rest = items[:i] + items[i+1:]
                for p in my_permute(rest, n-1):
                    yield v + p


class GauIO:
    '''\
 Manage the input file of the Gaussian package.\n\
   INPUT VARIABLES    ::\n\
 iout                : FLOW of output file\n\
 fn                  : STRING of the input file name\n\
                       fn = None -> default, bypass the action to open the file\n\
 bugctrl             : INTEGER to control the bebug information print\n\
                       0: default\n\
                       1: more results \n\
                       2: more results add detail debugging info.\n\
    MANAGE VARIABLES   :: \n\
       (***Note***: if need, all the strings following have to be specified in lower style\n\
 self.IOut           : FLOW of output file\n\
 self.FileName       : STRING of the input file name\n\
 self.ModuDir        : STRING of the path name to reach the gaussian environment
 self.JobName        : STRING of this job name\n\
 self.ChkName        : STRING, the name of the check file\n\
 self.CartesianFlag  : LOGICAL to state geomerty input is Cartesian coordinate or not\n\
 self.IPrint         : INTEGER of print level\n\
 self.Charge         : INTEGER of input Chage   \n\
 self.Spin           : INTEGER of input Spin    \n\
 self.NAtom          : INTEGER of atoms number \n\
 self.MachineList    : LIST of machine commands\n\
 self.OptionList     : LIST of options\n\
 self.TitleList      : LIST of this job title \n\
 self.GeomList       : LIST of input geometry   \n\
 self.IAn            : LIST of Atom index\n\
 self.AtLabel        : LIST of Atom label\n\
 self.CList          : LIST of geometry coordinate\n\
 self.ZList          : LIST of geometry Z-matrix\n\
 self.ZListR         : LIST of Z-matrix parameters\n\
 self.RestList       : LIST of rest content after geometry\n\
 --------------------------Results Collected--------------------------------------\n\
 self.EngyReal       : FLOAT, the total energy of the quesion\n\
 self.ForcList       : LIST, Force of the quesion\n\
 self.HessList       : LIST, Hessian of the quesion\n\
 self.DipoList       : List, Dipole\n\
 self.DpDvList       : List, Dipole derivatives\n\
 self.PolaList       : List, Polarizability\n\
 --------------------------Several varibles for options control--------------------------\n\
 self.KickOptionList : DICTIONARY of initial disable options\n\
              dict.keys() = []\n\
 self.MoreOptionDict : DICTIONARY of options which need more detailed handle\n\
              dict.keys() = ['checkpoint','allcheck','extraoverlay','fchk=all','%chk']\n\
 self.ExOvList       : LIST of IOPs for the option "extraoverlay"\n\
 ----------------------------------------------------------------------------------------\
 '''

    AtDict    = {\
   'x':0  ,
   'h':1  , 'he':2  ,\
  'li':3  , 'be':4  ,  'b':5  ,  'c':6  ,  'n':7  ,  'o':8  ,  'f':9  , 'ne':10 ,\
  'na':11 , 'mg':12 , 'al':13 , 'si':14 ,  'p':15 ,  's':16 , 'cl':17 , 'ar':18 ,\
   'k':19 , 'ca':20 , 'ga':31 , 'ge':32 , 'as':33 , 'se':34 , 'br':35 , 'kr':36 ,\
  'sc':21 , 'ti':22 ,  'v':23 , 'cr':24 , 'mn':25 ,\
  'fe':26 , 'co':27 , 'ni':28 , 'cu':29 , 'zn':30 ,\
  'rb':37 , 'sr':38 , 'in':49 , 'sn':50 , 'sb':51 , 'te':52 ,  'i':53 , 'xe':54 ,\
   'y':39 , 'zr':40 , 'nb':41 , 'mo':42 , 'tc':43 ,\
  'ru':44 , 'rh':45 , 'pd':46 , 'ag':47 , 'cd':48 ,\
  'cs':55 , 'ba':56 , 'tl':81 , 'pb':82 , 'bi':83 , 'po':84 , 'at':85 , 'rn':86 ,\
  'la':57 , 'hf':72 , 'ta':73 ,  'w':74 , 're':75 ,\
  'os':76 , 'ir':77 , 'pt':78 , 'au':79 , 'hg':80 \
                  }
    AnDict    = {\
     0:'X' ,
     1:'H' ,   2:'He',\
     3:'Li',   4:'Be',   5:'B' ,   6:'C' ,  7 :'N' ,   8:'O' ,   9:'F' ,   10:'Ne',\
    11:'Na',  12:'Mg',  13:'Al',  14:'Si',  15:'P' ,  16:'S' ,  17:'Cl',   18:'Ar',\
    19:'K' ,  20:'Ca',  31:'Ga',  32:'Ge',  33:'As',  34:'Se',  35:'Br',   36:'Kr',\
    21:'Sc',  22:'Ti',  23:'V' ,  24:'Cr',  25:'Mn',\
    26:'Fe',  27:'Co',  28:'Ni',  29:'Cu',  30:'Zn',\
    37:'Rb',  38:'Sr',  49:'In',  50:'Sn',  51:'Sb',  52:'Te',  53:'I',    54:'Xe',\
    39:'Y' ,  40:'Zr',  41:'Nb',  42:'Mo',  43:'Tc',\
    44:'Ru',  45:'Rh',  46:'Pd',  47:'Ag',  48:'Cd',\
	55:'Cs',  56:'Ba',  81:'Tl',  82:'Pb',  83:'Bi',  84:'Po',  85:'At',   86:'Rn',\
	57:'La',  72:'Hf',  73:'Ta',  74:'W',   75:'Re',\
	76:'Os',  77:'Ir',  78:'Pt',  79:'Au',  80:'Hg' \
                  }
    def __init__(self,iout,fn=None,bugctrl=0):
        '''\
        Initialize variables belonged to GauIO\
        '''
        import sys
        import os
        import os.path
        from os      import getcwd
        from os      import getenv
        from os.path import isfile

        self.IOut        = iout                                       # Flow of the output file
        self.IPrint        = bugctrl                                    # to control the printing out
        self.FileName    = fn                                         # Name of the input file

        self.WorkDir    = getcwd().strip()                           # STRING, current DIR 
        self.HomeDir    = getenv('HOME')                             # STRING, Home DIR
        if isfile('%s/.xdh_modules_path' %self.HomeDir):           # Load Private Modules DIR
            tmpf    = open('%s/.xdh_modules_path'\
                %self.HomeDir,'r')
            self.ModuDir=tmpf.readline().strip()                              # STRING, PATH of my modules
            sys.path.append(self.ModuDir)                                     # Append it into "sys.path"
            tmpf.close()
        else:
            print(('Error for loading \"$HOME/.xdh_modules_path\" \n'+\
                'which contains absolute path of relevant python modules'))
            sys.exit(1)

        global __version__
        if isfile('%s/version.txt' %self.ModuDir):                   # Load Private Modules DIR
            tmpf    = open('%s/version.txt'\
                %self.ModuDir,'r')
            __version__=tmpf.readline().strip()                      
        else:
            __version__='no version info.'
         
        if self.FileName==None:
            self.f  = 'None'
            self.JobName  = 'TmpName'
            if self.IPrint>=1:
                print_String(self.IOut,
                    'Do not open input file for GauIO class',1)
        else:
            try:
                self.f     = open(self.FileName,'r')                    # Open it
            except IOError:
                self.f  = open('Error_%s' %self.FileName,'w')
            if self.IPrint>=1:
                print_String(self.IOut,
                    'Open Gau-Input file "%s" for GauIO class'
                    % self.FileName,1)
            #
            #Generating "self.JobName" which is head of "self.FileName"
            #For example, if    self.FileName  = 'g03_1.gjf'
            #             then  self.JobName   = 'g03_1'
            path, filename  = \
                os.path.split(os.path.abspath(self.FileName))
            name, extension = os.path.splitext(filename)
            self.JobName= name                                       # Name of this job

            #if bugctrl>=1:
            #    print_String(iout,'Enter the job : \"%s\"'
            #        % self.JobName,2)


        self.MachineList=[]                                          # List, machine commands
        self.OptionList= []                                          # List, options
        self.KickOptionList=['nonstd']                               # Default disable options
        self.MoreOptionDict={'checkpoint':0,'check':0,'allcheck':0,
                'scrf':0,'fchk=all':0,'extraoverlay':0,'%chk': 0}    # Dict., options complicated
        self.ExOvList    = []
        self.TitleList    = []                                         # List of this job title
        self.Charge    = ''                                             # Input Chage
        self.Spin    = ''                                             # Input Spin
        self.GeomList    = []                                         # Input Geometry
        self.AtLabel    = []                                         # List, Atom labels
        self.IAn    = []                                             # List, Atom indexs
        self.CList    = []                                             # List, Atom coordinates
        self.CartesianFlag=False                                     # "True"  : Cartesian
        self.ZList  = []                                             # List, Atom Z-matrix
        self.ZListR = []                                             # List, Z-matrix parameters
                                                                     # "False" : Z-Matrix;
        self.NAtom    = 0                                              # Number of atoms
        self.RestList    = []                                         # Rest content after Geometry

        self.EngyReal    = 0.0                                        # REAL, the total energy
        self.ForcList    = []                                         # LIST, Force 
        self.HessList    = []                                         # LIST, Hessian
        self.DipoList    = []                                         # List, Dipole
        self.DpDvList    = []                                         # List, Dipole derivatives 
        self.PolaList    = []                                         # List, Polarizability
        return
    def __del__(self):
        '''Close the document flow of input file'''
        if self.f=='None':
            if self.IPrint>=2:
                print_String(self.IOut,
                    'Do not close input file for GauIO class',1)
        else:
            if self.IPrint>=2:
                print_String(self.IOut,
                    'Close Gau-Input file "%s" for GauIO class'
                    % self.FileName,1)
            self.f.close()                                           # Close input file
        return
    def get_MachAndOpt(self):
        '''Loading machine commands and initial options from input file\n\
  Note: the file flow locates in the next blank line below options\n\
 '''
        line=self.f.readline().strip()
        if len(line)==0:
            print_Error(self.IOut,
                'Error occurs in reading option keywords from' +\
                ' "%s". Please make it valid ' % self.FileName)
        while line[0]=='%':                                          # Machine commands first
            self.MachineList.append(line)
            line=self.f.readline().strip()
            if len(line)==0:
                print_Error(self.IOut,
                'Error occurs in reading option keywords from' +\
                ' "%s". Please make it valid ' % self.FileName)
        if line[0]=='#':                                             # Then for options
            TmpOptions=[]
            TmpOptions=line[1:].strip().split()
            for option in TmpOptions:
                self.OptionList.append(option)
            line=self.f.readline().strip()
            while len(line)!=0:
                TmpOptions=[]
                TmpOptions=line.strip().split()
                for option in TmpOptions:
                    self.OptionList.append(option)
                line=self.f.readline().strip()
        else:
            print_Error(self.IOut,
                'Error occurs in reading option keywords from' +\
                ' "%s". Please make it valid ' % self.FileName)

        for option in self.MachineList:                              # Get the name of chkfile
            if option.lower().find('%chk')!=-1:
                TmpList        =\
                    option.strip().split('=')[1].split('.')
                if len(TmpList)==1 or len(TmpList)==2:
                    self.ChkName= TmpList[0]
                else:
                    self.ChkName= '.'.join(TmpList[0:-1])
                break
        else:
            self.MoreOptionDict['%chk'] = 1                          # Avoid unrequired chkfile
                                                                     # copy back
            self.ChkName    = self.JobName
            self.MachineList.insert(0,
                '='.join(['%chk','%s.chk' % self.ChkName]))
        if self.IPrint>=2:
            print_String(self.IOut,
                'ChkName is %s.chk' % self.ChkName,1)
            if len(self.MachineList)>0:
                print_List(self.IOut,self.MachineList,2,
                    'Machine commands(%s) :' % len(self.MachineList))
            print_List(self.IOut,self.OptionList,3, 
                'Job options(%s) :' % len(self.OptionList))
        return
    def ctrl_Option(self):
        '''\
        Analysis the Options.\n\
            1) "self.KickOptionList" shall be specified before this function\n\
            2) "self.MoreOptionList" will be specified in this function\n\
            3) "self.ExOvLayList"    will be loaded, if the option of "extraoverlay" exists\
        '''
        #
        #First to kick out some specifical option which could not be
        #handled in this version
        #
        if self.IPrint>=2:
            if len(self.KickOptionList)>0:
                print_List(self.IOut,self.KickOptionList,3,
                    'KickOptionList(%s) :' % len(self.KickOptionList))
        if len(self.KickOptionList)>0:
            for option in self.OptionList:
                for ivdmd in self.KickOptionList:
                    tmp1=option.lower()
                    if tmp1.find(ivdmd)!=-1:
                        print_Error(self.IOut,
                            'Invalid option(%s) is found in gausian'\
                            %ivdmd + ' input file.')
        #
        #Then to determine the complicated option initialization
        #
        tmpKeys = list(self.MoreOptionDict.keys())
        for key in tmpKeys:
            self.MoreOptionDict[key]=0
        for option in self.OptionList:
            tmpKey=option.strip().lower().split('=')[-1]
            if tmpKey in tmpKeys:
                self.MoreOptionDict[tmpKey]=1
                continue
            tmpKey=option.strip().lower().split('(')[0]
            if tmpKey in tmpKeys:
                self.MoreOptionDict[tmpKey]=1
        if self.MoreOptionDict['checkpoint']==1:                     # To make sure that: 
            for option in self.OptionList:                           # Check "geom=checkpoint"
                tmpList = option.strip().lower().split('=')
                if len(tmpList)!=2: continue
                if tmpList[-1]=='checkpoint' and tmpList[0]=='geom':
                    self.MoreOptionDict['checkpoint']=1
                    break
            else:
                self.MoreOptionDict['checkpoint']=0
        if self.MoreOptionDict['check']==1:                          # To make sure that: 
            for option in self.OptionList:                           # Check "geom=check"
                tmpList = option.strip().lower().split('=')
                if len(tmpList)!=2: continue
                if tmpList[-1]=='check' and tmpList[0]=='geom':
                    self.MoreOptionDict['checkpoint']=1
                    break
        if self.MoreOptionDict['scrf']==1:
            for option in self.OptionList:
                if option.lower().find('scrf')!=-1:
                    self.MoreOptionDict['scrf']=option
                    break
            else:
                self.MoreOptionDict['scrf']=='unknown solvation'
        #
        #Now get the addtion iop command by option "extraoverlay"
        #
        if self.MoreOptionDict['extraoverlay']==1:
            tmpExOvLay=self.f.readline().strip()
            while len(tmpExOvLay)!=0:
                self.ExOvList.append(tmpExOvLay)
                tmpExOvLay=self.f.readline().strip()
        if self.IPrint>=2:
            if len(self.ExOvList)>0:
                print_List(self.IOut,self.ExOvList,3,
                    "ExtraOverlay IOPs(%s) :" % len(self.ExOvList))
        return
    def get_TCSGR(self):
        '''Get the (T)itle, (C)harge, (S)pin, (G)eom and (R)est content from input file'''
        from re import compile
        #
        # If option of "allcheck" is stated,
        #   bypass (T)(C)(S)(G) and load (R)
        #
        if self.MoreOptionDict['allcheck']==1:
            self.RestList=self.f.readlines()
            if self.IPrint>=1:
                print_String(self.IOut,
                    '\"AllCheck\" is stated, then title,charge, '+\
                    'spin and geometry are specified in CheckFile',1)
        else:
            #
            # First to get "self.TitleList" (T)
            #
            self.TitleList.append(self.f.readline().strip())
            if len(self.TitleList[0])==0:
                print_Error(self.IOut, 
                    'Missing job title in %s ' % self.FileName +\
                    '"GauIO.get_TCSGR"')
            addTitle=self.f.readline().strip()
            while len(addTitle)!=0:                                  # Read the title more row
                self.TitleList.append(addTitle)
                addTitle=self.f.readline().strip()
            if self.IPrint>=2:
                if len(self.TitleList)==1:
                    print_String(self.IOut,
                        'Title : "%s"' %self.TitleList[0],1)
                else:
                    print_List(self.IOut,self.TitleList,2,
                        'Title of job:')
            #
            # Then to get "self.Charge" (C) and "self.Spin" (S)
            #
            p1      = compile(' +|, *')
            ChargeSpin=self.f.readline().strip()
            if len(ChargeSpin)==0:
                print_Error(self.IOut, 
                    'Missing Charge and Spin in %s' % self.FileName)
            TmpList = p1.split(ChargeSpin)
            TmpList = [x.strip() for x in TmpList]
            if TmpList.count('')!=0:
                for i in range(TmpList.count('')):
                    TmpList.remove('')
            try:
                self.Charge=int(TmpList[0])
            except ValueError:
                print_Error(self.IOut, 
                    'Invalid Charge in %s' % self.FileName)
            try:
                self.Spin=int(TmpList[1])
            except ValueError:
                print_Error(self.IOut,
                    'Invalid Spin in %s' % self.FileName)
            if self.IPrint>=2:                                       # Debugging
                print_String(self.IOut, 
                    'Charge and Spin : (%d, %d)'
                    % (self.Charge, self.Spin),1)
            #
            # If option of "checkpoint" is stated,
            #  bypass "self.GeomList" (G) and load "self.RestList" (R)
            #
            if self.MoreOptionDict['checkpoint']==1:
                #self.f.readline()
                self.RestList=self.f.readlines()
                if self.IPrint>=1:
                    print_String(self.IOut,
                        'Molecular geometry will be loaded'+\
                        ' from the CheckFile',1)
                if self.IPrint>=2:
                    print_List(self.IOut,self.RestList,1)
            else:
                #
                #Then to get "self.GeomList" (G) and "self.NAtom"
                #
                self.NAtom=0
                p1      = compile(' {1,}|, {0,}')
                tmpGeom = self.f.readline().strip().replace('\t',' ')
                tmpList = [x.strip() for x in p1.split(tmpGeom)]
                for i in range(tmpList.count('')):
                    tmpList.remove('')
                if self.IPrint>=2:
                    print_List(self.IOut,tmpList,1)
                if len(tmpList)==4:                                  # To get Cartesian coordinate
                    self.CartesianFlag=True
                    if self.IPrint>=1:
                        print_String(self.IOut,
                            'Loading Cartesian coordinates',1)
                    while len(tmpList)==4:
                        self.GeomList.append(tmpGeom)
                        self.NAtom += 1
                        self.AtLabel.append(tmpList[0])
                        try:
                            self.CList.append([float(tmpList[1]),
                                float(tmpList[2]),float(tmpList[3])])
                        except ValueError:
                            print_String(self.IOut,
                            'Warning: could not get CList "GauIO.get_TCSGR"',1)
                        tmpGeom=self.f.readline().strip().replace('\t',' ')
                        tmpList     =\
                            [x.strip() for x in p1.split(tmpGeom)]
                        for i in range(tmpList.count('')):
                            tmpList.remove('')
                    else:
                        print_List(self.IOut,self.GeomList,2,
                            Info='In job=%s' %self.JobName)
                elif len(tmpList)>=5:                                # "5" : for atom fixing
                                                                     # "6" : for Onion-type input
                    self.CartesianFlag=True
                    if self.IPrint>=2:
                        print_String(self.IOut,
                            'Loading Cartesian coordinates',1)
                    while len(tmpList)>=5:
                        self.GeomList.append(tmpGeom)
                        self.NAtom += 1
                        self.AtLabel.append(tmpList[0])
                        try:
                            self.CList.append([float(tmpList[2]),
                                float(tmpList[3]),float(tmpList[4])])
                        except ValueError:
                            print_String(self.IOut,
                            'Warning: could not get CList "GauIO.get_TCSGR"',1)
                        tmpGeom=self.f.readline().strip().replace('\t',' ')
                        tmpList     =\
                            [x.strip() for x in p1.split(tmpGeom)]
                        for i in range(tmpList.count('')):
                            tmpList.remove('')
                elif len(tmpList)==1:                                # To get Z-Matrix coordinate
                    if self.IPrint>=1:
                        print_String(self.IOut,
                                'Loading Z-Matrix coordinates',1)
                    self.CartesianFlag=False
                    tmpLength=1
                    while len(tmpList)==tmpLength or \
                        len(tmpList) == 8:
                        self.GeomList.append(tmpGeom)
                        self.NAtom=self.NAtom+1
                        self.AtLabel.append(tmpList[0])
                        self.ZList.append(tmpList)
                        tmpGeom=self.f.readline().strip().replace('\t',' ')
                        tmpList     =\
                            [x.strip() for x in p1.split(tmpGeom)]
                        for i in range(tmpList.count('')):
                            tmpList.remove('')
                        if self.IPrint>=2:
                            print_List(self.IOut, tmpList,
                                3, 'Igor Debugging')
                            print_String(self.IOut, 
                                'Len is %s, meeting %s'\
                                %(tmpLength,tmpList), 2)
                        if tmpLength<7:
                            tmpLength = tmpLength+2

                    LocPos  = self.f.tell()                          # get Z-Matrix parameter
                    p1      = compile('= {0,}| {1,}|, {0,}')
                    tmpGeom=self.f.readline().strip().replace('\t',' ')
                    tmpList =\
                        [x.strip() for x in p1.split(tmpGeom)]
                    for i in range(tmpList.count('')):
                        tmpList.remove('')
                    if len(tmpList)==2:
                        tmpItem = tmpList[0].strip()
                        tmpFlag = False
                        for i in self.GeomList:
                            if i.find(tmpItem)!=-1:
                                tmpFlag = True
                                break
                        for j in self.AtLabel:
                            if j.lower().find(tmpItem.lower())!=-1:
                                tmpFlag = False
                                break
                        if tmpFlag:
                            if self.IPrint>=1:
                                print_String(self.IOut,
                                'Loading Z-Matrix parameters',
                                1)
                            while len(tmpList)==2:
                                self.ZListR.append(tmpGeom)
                                tmpGeom = self.f.readline().strip().replace('\t',' ')
                                tmpList =\
                                    [x.strip()\
                                    for x in p1.split(tmpGeom)]
                                for i in range(tmpList.count('')):
                                    tmpList.remove('')
                        else:
                            self.f.seek(LocPos)
                    else:
                        self.f.seek(LocPos)
                else:
                    print_Error(self.IOut,'Error in Geom. Info.')
  
                p1    = compile('\d+')                                 # Gen. IAn from AtLabel
                for TmpLabel in self.AtLabel:
                    if p1.match(TmpLabel[0:2]):
                        try:
                            self.IAn.append(int(TmpLabel))
                        except ValueError:
                            print_Error(self.IOut, 
                                "Error in IAn identification")
                    else:
                        TmpRe    = TmpLabel.lower()
                        for key in sorted(GauIO.AtDict.keys()):
                            if TmpRe==key:
                                self.IAn.append(GauIO.AtDict[key])
                                break
                        else:
                            print_Error(self.IOut,\
                               'Error in IAn identification')
                if self.IPrint>=2:
                    print_List(self.IOut,self.IAn,3,
                         'IAn List(%d) :' % len(self.IAn))

                self.RestList       = self.f.readlines()             # Load "self.RestList"
                if self.IPrint>=2:                                   # Debugging
                    if self.CartesianFlag:
                        print_List(self.IOut,self.GeomList,2,
                           'Cartesian Coordinate(%s)' 
                           % len(self.GeomList))
                    else:
                        print_List(self.IOut,self.GeomList,2,
                            'Z-Matrix Coordinate(%s)' 
                            % len(self.GeomList))
                        if len(self.ZListR)>0:
                            print_List(self.IOut,self.ZListR,2,
                                'Z-Matrix Parameters(%s)'
                                % len(self.ZListR))

        for i in range(len(self.RestList)):                          # filter the RestList
            self.RestList[i]=self.RestList[i].strip()
        if len(self.RestList)==0:
            pass
        elif len(self.RestList)==1:
            if len(self.RestList[0])==0:
                self.RestList   = []
        else:
            TmpIndex    = 0
            if self.RestList[0]=='' and self.RestList[1]=='':
                self.RestList = []
            else:
                while TmpIndex <= len(self.RestList)-1:
                    if self.RestList[TmpIndex]=='' and TmpIndex==0:
                        self.RestList   = self.RestList[1:]
                        TmpIndex = TmpIndex - 1
                    elif self.RestList[TmpIndex]=='' and TmpIndex!=0:
                        if self.RestList[TmpIndex-1]=='':
                            self.RestList   =\
                                self.RestList[:TmpIndex-1]
                            break
                    TmpIndex    = TmpIndex + 1
                # Fix bug for conflicting between RestList and 
                # l608 statement for R5DFT
                if self.RestList[-1].strip()=='':
                        self.RestList =\
                            self.RestList[:-1]
        if self.IPrint>=2:
            if len(self.RestList)!=0:
                print_List(self.IOut,self.RestList,2,
                    'RestList(%s):' % len(self.RestList))
                print_List(self.IOut,self.RestList,1)
        return 
    def form_Inp(self):
        '''Form gaussian input file for given statements'''
        #
        # Now creat a new gaussian input file
        #
        if (self.FileName!=None):
            wf=open('Job_%s.com' % self.JobName,'w')
        else:
            wf=open('%s.com' % self.ChkName,'w')
        # Form machine commands
        for line in self.MachineList:
            wf.write('%s\n' % line)
        # Input options
        self.OptionList = [x.lower() for x in self.OptionList]
        self.OptionList.insert(0,'p')
        for item in self.OptionList:
            for i in range(self.OptionList.count(item)-1):
                self.OptionList.remove(item)
                if self.IPrint>=2:
                    print_String(self.IOut,
                        'Remove option "%s"' %item,1)
        if len(self.OptionList)>10:
             self.OptionList[9]="%s\n" % self.OptionList[9]
        if len(self.OptionList)>20:
             self.OptionList[19]="%s\n" % self.OptionList[19]
        Option=' '.join(self.OptionList)
        wf.write('#%s\n' % Option)
        wf.write('\n')
        # Fill extraoverlay IOps
        if self.MoreOptionDict['extraoverlay']==1:
            for line in self.ExOvList:
                wf.write('%s\n' %line)
            wf.write('\n')
        # Write Title, Charge, Spin, Geom and Rest, if required
        if self.MoreOptionDict['allcheck']==1:                       # "Geom=AllCheck"
            for line in self.RestList:
                wf.write('%s\n' %line) 
            wf.write(' \n')
            wf.write(' \n')
        else:                                                        # Normal case and CheckPoint
            # Write the job title
            for title in self.TitleList:
                wf.write('%s\n' % title)
            wf.write('\n')
            # Write the Charge and Spin
            wf.write('%d%5d\n' %(self.Charge,self.Spin))
            if self.MoreOptionDict['checkpoint']==0:                 # For normal case
                # Write geometry info.
                for geom in self.GeomList:
                    wf.write('%s\n' %geom)
                wf.write(' \n')
                if not self.CartesianFlag and len(self.ZListR)!=0:
                    if self.IPrint>=2:
                        print_String(self.IOut,
                            'Now Loading ZListR',2)
                    for item in self.ZListR:
                        wf.write('%s\n' %item)
                    wf.write('  \n')
            else:
                wf.write(' \n')
            for line in self.RestList:
                wf.write('%s\n' % line)
            wf.write(' \n')
            wf.write(' \n')
        # Now complete the input file
        wf.close()
        return
    def update_CList(self):
        '''Update CList from GeomList'''
        from re     import compile
        self.CList    = []
        p1          = compile(' {1,}|, {0,}')
        for i in range(self.NAtom):                                  # Update CList
            TmpList =\
                [x.strip() for x in p1.split(self.GeomList[i])]
            for i in range(TmpList.count('')):
                TmpList.remove('')
            if len(TmpList)==4:
                self.CList.append([float(TmpList[1]),
                    float(TmpList[2]),float(TmpList[3])])
            elif len(TmpList)==5:
                self.CList.append([float(TmpList[2]),
                    float(TmpList[3]),float(TmpList[4])])
            else:
                print_Error(self.IOut,"Error in CList updating")
        return
    def run_GauJob(self,iop=0,sync=None):
        '''run Gaussian Job through a specified Cshell named \'Gau_Environment\'
        '''
        from os import system
        from os import mkdir
        from os import removedirs
        from os import remove
        from os import chdir
        from os import listdir
        from os import rename
        from os.path    import exists
        from os.path    import isfile
        from os import getpid                                        # To make more then 2 jobs
        import shutil
        import subprocess
                                                                     # running in one DIR
        CurrGID = getpid()
        CurrScr = '%s/.script%s' %(self.WorkDir,CurrGID)
        if not exists(CurrScr):
            mkdir(CurrScr)                                    # Creat the current work dir
        src = 'Job_%s.com' % self.JobName
        dst = '%s/Job_%s.com' % (CurrScr, self.JobName) 
        rename(src, dst)
        src = '%s/ChkReplace.py' % self.ModuDir
        dst = '%s/ChkReplace.py' % self.WorkDir
        shutil.copy2(src, dst) # ChkReplace.py for chkpath

        #print('Script is %s' %CurrScr)
        #print('self.ModuDir %s' %self.ModuDir)
        if sync!=None:                                               # Message CurrScr between 
            sync.value  = '%s' % CurrScr                             #   threads

        dst = '%s/Job_%s.com' % (CurrScr, self.JobName) 
        if iop==0:                                                   # normal gaussian calculation
            #chdir(CurrScr)
            exe = '%s/Gau_Environment' % self.ModuDir
            #chdir(self.WorkDir)
        elif iop==3:                                                 # just g03, for R5DFT class
            #chdir(CurrScr)
            exe = '%s/G03_Environment' % self.ModuDir
            #chdir(self.WorkDir)
        elif iop==9:                                                 # just g09, for R5DFT class
            #chdir(CurrScr)
            exe = '%s/G09_Environment' % self.ModuDir
        elif iop==16:                                                # just g16, for R5DFT class
            #chdir(CurrScr)
            exe = '%s/G16_Environment' % self.ModuDir
        else:
            raise RuntimeError('invalid Gaussian version %s' % iop)
        subprocess.call([exe, dst])

        shutil.copy2('%s/Job_%s.log' % (CurrScr, self.JobName),
                'Job_%s.log' % self.JobName)
        #if self.MoreOptionDict['%chk']==1:                            # Do not save chk
        #    remove('%s.chk' % self.ChkName)
        remove('ChkReplace.py')
        #if self.MoreOptionDict['fchk=all']==1:                       # save Test.FChk for fchk=all
        #    system('mv %s/Test.FChk Test.FChk' % CurrScr)
        #if self.IPrint<2:
        #    for tmpFile in listdir('%s' % CurrScr):
        #        remove('%s/%s' % (CurrScr,tmpFile))
        #    removedirs('%s' %CurrScr)
        return
    def copy_DataIO(self,DataIO):
        '''Copy all the IO data from DataIO'''
        from copy   import deepcopy
        # Initial parameters
        self.IOut           = DataIO.IOut                            # Flow
        #self.FileName       = DataIO.FileName                       # Strings
        self.ModuDir        = DataIO.ModuDir
        self.JobName        = DataIO.JobName
        self.ChkName        = DataIO.ChkName
        self.CartesianFlag  = DataIO.CartesianFlag                   # Logic
        self.IPrint         = DataIO.IPrint                          # Integers
        self.Charge         = DataIO.Charge
        self.Spin           = DataIO.Spin
        self.NAtom          = DataIO.NAtom
        self.MachineList    = DataIO.MachineList[:]                  # Lists
        self.OptionList     = DataIO.OptionList[:]
        self.TitleList      = DataIO.TitleList[:]
        self.GeomList       = DataIO.GeomList[:]
        self.IAn            = DataIO.IAn[:]
        self.AtLabel        = DataIO.AtLabel[:]
        self.CList          = DataIO.CList[:]
        self.ZList          = DataIO.ZList[:]
        self.ZListR         = DataIO.ZListR[:]
        self.RestList       = DataIO.RestList[:]
        self.MoreOptionDict = deepcopy(DataIO.MoreOptionDict)        # Dictionaries  
        return

class OptHandle:
    '''\
    Handle Geometry Optimization\
    '''
    def __init__(self,iout,IOClass,bugctrl=0):
        '''\
        Initialize variables belonged to OptHandle\
        '''
        import sys
        from os        import getcwd
        from os        import getenv
        from os.path   import isfile
        from os.path   import exists
        from os        import mkdir
        from re        import compile

        self.IOut       = iout
        self.IPrint     = bugctrl
        self.IOClass    = IOClass

        self.WorkDir    = getcwd().strip()                           # STRING, Working DIR 
        self.HomeDir    = getenv('HOME')                             # STRING, HOME DIR
        if isfile('%s/.xdh_modules_path' %self.HomeDir):           # Load Private Modules DIR
            tmpf    = open('%s/.xdh_modules_path'\
                %self.HomeDir,'r')
            self.ModuDir=tmpf.readline().strip()                     # STRING, PATH of my modules
            ModuDir=self.ModuDir
            sys.path.append(ModuDir)                                 # Append it into "sys.path"
            tmpf.close()
        else:
            print(('Error for loading \"$HOME/.xdh_modules_path\" \n'+\
                'which contains absolute path of relevant python modules'))
            sys.exit(1)
         
        self.Opt        = False                                      # LOGIC for geom. opt. task
        self.Conv       = False                                      # LOGIC for opt convergency
        self.Hessian    = False                                      # LOGIC for freq calc. task
        self.MaxStep    = 30                                         # INTEGER, maximum step size
        self.NP            = 1                                          # INTEGER, opt points
        self.OptType    = 0                                          # INTEGER :: 
                                                                     #  0 for Estimated Hessian;
                                                                     #  1 for "calcfc";
                                                                     #  2 for "calcall"

        self.OptionList = []                                         # LIST, options prepered
        self.Energy     = self.IOClass.EngyReal                      # REAL, Total energy
        self.ForcList   = self.IOClass.ForcList                      # LIST, Force in question
        self.HessList   = self.IOClass.HessList                      # LIST, Hessian in question
        for option in self.IOClass.OptionList:
            tmpOpt=option.lower().strip()
            if tmpOpt.find('opt')!=-1:
                self.Opt= True
                print_String(self.IOut,
                    'Geometry optimization turns on',1)
                p1    = compile('[(),]')
                for tmp1 in p1.split(tmpOpt):
                    if tmp1.find('maxstep')!=-1:
                        self.MaxStep    = int(tmp1.split('=')[1])
                        tmpStep        = float(self.MaxStep)*0.01
                        print_String(self.IOut,
                            'Maximum step size (%8.4f) in Opt.'
                            % tmpStep, 1)
                    if tmp1.find('calcfc')!=-1:
                        self.OptType=1
                        print_String(self.IOut,
                            'OptType: \"CalcFC\" turns on', 1)
                    elif tmp1.find('calcall')!=-1:
                        self.OptType=2
                        print_String(self.IOut,
                            'OptType: \"CalcAll\" turns on', 1)
            elif tmpOpt.find('freq')!=-1:
                self.Hessian=True
                print_String(self.IOut,
                    'Frequence calculation turns on', 1)
            elif tmpOpt=='force':
                pass
            else:
                self.OptionList.append(option)
       #self.IOClass.OptionList    = self.OptionList[:]                 # recopy Options into IOClass
        if self.IPrint >= 2:
            print_List(self.IOut,self.OptionList,3,
                'Filtering OptionList by "OptHandle" '+\
                'for geomerty optimization')
        if self.Opt:
            if not exists('%s/GeomOpt' % self.WorkDir):
                mkdir('%s/GeomOpt' % self.WorkDir)
            tmpf=open('%s/GeomOpt/geom' % self.WorkDir,'w')
            tmpf.write('%5i\n' % self.IOClass.NAtom)
            for i in self.IOClass.GeomList:
                tmpf.write('%s\n' %i)
            tmpf.close()
        return
    def __del__(self):
        return