import os, sys, math
import re
from numpy import *
import subprocess
import GaussianRunPack.AtomInfo
import shutil
from rdkit import Chem
from rdkit.Chem import AllChem

class GaussianDFTRun:

    def __init__(self, functional, basis, nproc, value, solvent, in_file, error):

        self.in_file = in_file
        self.functional = functional.lower()
        self.basis = basis.lower()
        self.nproc = nproc
        self.value = value.lower()
        self.solvent = solvent.lower()
        self.error = error

        self.mem = ''

    def Extract_ExcitedState(self, lines):
        Egrd = 0.0
        Eext = 0.0
        Found = False
        WaveLength = []
        V_OS = []
        for line in lines:
            if line.find("SCF Done:  ") >=0:
                line_SCFEnergy = re.split("\s+", line)
                Egrd = float(line_SCFEnergy[5])
            if line.find("Total Energy, E(TD-HF/TD-DFT)") >=0:
                line_totalenergy = line.split('=')
                Eext = float(line_totalenergy[1])
            if line.find("Excitation energies and oscillator strengths:") >=0:
                 WaveLength = []
                 V_OS = []
            if line.find("Excited State  ") >=0:
                 line_StateInfo = line.split()
                 WaveLength.append(float(line_StateInfo[6]))
                 OS_info = line_StateInfo[8].split('=')
                 V_OS.append(float(OS_info[1]))
            if line.find("-- Stationary point found.") >=0:
                Found = True

        return Found, Egrd, Eext, WaveLength, V_OS


    def Extract_SCFEnergy(self, lines):

        Energy=[]

        for line in lines:
            if line.find("SCF Done:  ") >=0:
                line_StateInfo = line.split()
        #print (line_StateInfo[4])
                Energy.append(float(line_StateInfo[4]))

        return Energy[-1]

    def Extract_MO(self, lines):

            AlphaEigenVal = []
            BetaEigenVal = []
#            for line in lines:
            for line in lines:
                if line.find(" basis functions, ") >=0:
                    line_StateInfo = line.split()
            #        print (line_StateInfo[0])
                    
                    NumBasisFunc = int(line_StateInfo[0])
  
                if line.find(" alpha electrons ") >=0:
                    line_StateInfo = line.split()
            #        print (line_StateInfo[0])
            #        print (line_StateInfo[3])
            
                    NumAlphaElec = int(line_StateInfo[0])
                    NumBetaElec = int(line_StateInfo[3])

                if line.find("Alpha  occ. eigenvalues --") >=0:
  
                    if len(AlphaEigenVal) == NumBasisFunc:
                        AlphaEigenVal=[]
                    
                    line_removed = line.replace("Alpha  occ. eigenvalues --", " ")
                    line_StateInfo = line_removed.split()
                    for i in range(len(line_StateInfo)):
                        AlphaEigenVal.append(float(line_StateInfo[i]))
  
                if line.find("Alpha virt. eigenvalues --") >=0:
                    line_removed = line.replace("Alpha virt. eigenvalues --", " ")
                    line_StateInfo = line_removed.split()
                    for i in range(len(line_StateInfo)):
                        AlphaEigenVal.append(float(line_StateInfo[i]))
  
                if line.find("Beta  occ. eigenvalues --") >=0:
  
                    if len(BetaEigenVal) == NumBasisFunc:
                        BetaEigenVal=[]
                    
                    line_removed = line.replace("Beta  occ. eigenvalues --", " ")
                    line_StateInfo = line_removed.split()
                    for i in range(len(line_StateInfo)):
                        BetaEigenVal.append(float(line_StateInfo[i]))
  
                if line.find("Beta virt. eigenvalues --") >=0:
                    line_removed = line.replace("Beta virt. eigenvalues --", " ")
                    line_StateInfo = line_removed.split()
                    for i in range(len(line_StateInfo)):
                        BetaEigenVal.append(float(line_StateInfo[i]))

            return NumAlphaElec, NumBetaElec, AlphaEigenVal, BetaEigenVal 

    def Extract_Coordinate(self, lines):

        print ("Start finding coordinates...")

        Atom_index = []
        NumElement = []
        AtomicType = []
        X = []
        Y = []
        Z = []

        count = 0
        count_Station = 0
        Total_NumStation =0
        Index_Station = 0

######Counting Stationary points##############
        for line in lines:
            if line.find("-- Stationary point found.") >= 0:
                count_Station += 1

##############################################
        Total_NumStation = count_Station

        if Total_NumStation > 0:
            count_Station = 0
            print ("Total number of stationary points: ", Total_NumStation)
            for line in lines:
                if line.find("-- Stationary point found.") >= 0:
                    count_Station += 1
                    Index_Station += 1
                    print ("A stationary point is found! #", Index_Station)
                    continue
                if line.find("Standard orientation:") >=0 and count_Station > 0:
#                    print ("Standard orientaion was found")
#                    print ("Start reading coordinate")
                    count += 1
                    continue
                if count == 1:
                    Border_1 = line
#                    print (Border_1)
                    count += 1
                    continue
                if count == 2:
                    Index_1 = line
#                    print (Index_1)
                    count += 1
                    continue
                if count == 3:
                    Index_2 = line
#                    print (Index_2)
                    count += 1
                    continue
                if count == 4:
                    Border_2 = line
#                    print (Border_2)
                    count += 1
                    continue
                if count >= 5:
                    i_atom = line.split()
#                    print (i_atom)
                    if len(i_atom) == 6:
                        Atom_index.append(int(i_atom[0]))
                        NumElement.append(int(i_atom[1]))
                        AtomicType.append(int(i_atom[2]))
                        X.append(float(i_atom[3]))
                        Y.append(float(i_atom[3]))
                        Z.append(float(i_atom[3]))
                        count += 1
                        continue
                    else :
                        print ("Reading atom coordinates is finished...")
                        N = count-5
                        print ("Number of atoms: ", N)
                        count = 0
                        count_Station = 0
                    continue

        else:
            for line in lines:
                if line.find("Standard orientation:") >=0:
                    print ("Standard orientaion was found")
                    print ("Start reading coordinate")
                    count += 1
                    continue
                if count == 1:
                    Border_1 = line
#                    print (Border_1)
                    count += 1
                    continue
                if count == 2:
                    Index_1 = line
#                    print (Index_1)
                    count += 1
                    continue
                if count == 3:
                    Index_2 = line
#                    print (Index_2)
                    count += 1
                    continue
                if count == 4:
                    Border_2 = line
#                    print (Border_2)
                    count += 1
                    continue
                if count >= 5:
                    i_atom = line.split()
#                    print (i_atom)
                    if len(i_atom) == 6:
                        Atom_index.append(int(i_atom[0]))
                        NumElement.append(int(i_atom[1]))
                        AtomicType.append(int(i_atom[2]))
                        X.append(float(i_atom[3]))
                        Y.append(float(i_atom[3]))
                        Z.append(float(i_atom[3]))
                        count += 1
                        continue
                    else :
                        print ("Reading atom coordinates is finished...")
                        N = count-5
                        print ("Number of atoms: ", N)
                        count = 0
                    continue


#Translating  atomic number to element symbol
        Mol_atom = []
#        Mol_CartX = zeros(N)
#        Mol_CartY = zeros(N)
#        Mol_CartZ = zeros(N)

        for i in range(N):
            Mol_atom.append(GaussianRunPack.AtomInfo.AtomicNumElec(NumElement[i]))
#            print (Mol_atom[i])

#############
        del Atom_index[:]
        del NumElement[:]
        del AtomicType[:]
        del X[:]
        del Y[:]
        del Z[:]

        return Mol_atom


    def Extract_values(self, infilename,opt, nmr,uv,energy,ipe, eae, pne, nne, gap,dipole,deen, stable2o2, fluor,tadf):

        ifile = open(infilename,'r') #open file for reading
        lines = ifile.readlines()
        ifile.close()

        output = {}

        for line in lines:
            if line.find("Error termination") >=0:
                print("Gausssian is stopped due to some errors")
                return output

        print("Spliting links....")
        Links = self.SplitLinks(infilename)
        n = len(Links)
        print("The number of linkes = ", n)
  
        GS_lines = Links[1].splitlines()
        Links[1] = ""

        lines = ""

        if gap == 1:
  
            NumAlphaElec, NumBetaElec, AlphaEigenVal, BetaEigenVal = self.Extract_MO(GS_lines)

            if BetaEigenVal == []:
                Alpha_gap = 27.211*(AlphaEigenVal[NumAlphaElec]-AlphaEigenVal[NumAlphaElec-1])
                output["gap"] = Alpha_gap
            #    return Alpha_gap
            else:
                Alpha_gap = 27.211*(AlphaEigenVal[NumAlphaElec]-AlphaEigenVal[NumAlphaElec-1])
                Beta_gap = 27.211*(BetaEigenVal[NumBetaElec]-BetaEigenVal[NumBetaElec-1])
  
                output["gap"] = [Alpha_gap, Beta_gap]
            #    return Alpha_gap, Beta_gap
  
  
        if dipole == 1:
            Dipole_X = []
            Dipole_Y = []
            Dipole_Z = []
            Dipole_Total = []

            for line in GS_lines:
                if line.find(" X= ") >=0:
                    line_StateInfo = line.split()
                    #print (line_StateInfo[1])
                    Dipole_X.append(float(line_StateInfo[1]))
                    #print (line_StateInfo[3])
                    Dipole_Y.append(float(line_StateInfo[3]))
                    #print (line_StateInfo[5])
                    Dipole_Z.append(float(line_StateInfo[5]))
                    #print (line_StateInfo[7])
                    Dipole_Total.append(float(line_StateInfo[7]))
                
            output["dipole"] = [Dipole_X[-1], Dipole_Y[-1], Dipole_Z[-1], Dipole_Total[-1]]
        #    return Dipole_X[-1], Dipole_Y[-1], Dipole_Z[-1], Dipole_Total[-1]
  
        if energy == 1:

            output["Energy"] = self.Extract_SCFEnergy(GS_lines)


        if deen == 1:
            try:
                GS_Energy = output["Energy"]
            except KeyError:
                GS_Energy = self.Extract_SCFEnergy(GS_lines)
        
            Mol_atom = self.Extract_Coordinate(GS_lines)

#######Calculating Decomposed atoms total energy#######################
            decomposed_Energy = 0
#
            for i in range(len(Mol_atom)):    
                #print (Mol_atom[i], GaussianRunPack.AtomInfo.One_Atom_Energy(Mol_atom[i], self.functional, self.basis))
                decomposed_Energy += GaussianRunPack.AtomInfo.One_Atom_Energy(Mol_atom[i], self.functional, self.basis)
            print("Decomposed energy: ", decomposed_Energy)
#
#############################################################
  
        #    return deen 
            output["deen"] = GS_Energy - (decomposed_Energy)

        if stable2o2 == 1:
            NumAlphaElec, NumBetaElec, AlphaEigenVal, BetaEigenVal = self.Extract_MO(GS_lines)
            O2_SOMO, O2_LUMO = GaussianRunPack.AtomInfo.O2_MO_refer(self.functional, self.basis)

            if BetaEigenVal == []:
                
               # OxidizedbyO2  >  0 oxidation by O2 is hard to occure.
               # OxidizedbyO2 <=  0 oxidation by O2 is easy to occure.

               # ReducedbyO2  >  0 reduction by O2 is hard to occure.
               # ReducedbyO2 <=  0 reduction by O2 is easy to occure.

                OxidizedbyO2 = O2_LUMO -AlphaEigenVal[NumAlphaElec-1]

                ReducedbyO2 =  AlphaEigenVal[NumAlphaElec] - O2_SOMO

                output["stable2o2"] = [OxidizedbyO2,ReducedbyO2]

            else:

                OxidizedbyO2 = O2_LUMO - BetaEigenVal[NumBetaElec-1] 

                ReducedbyO2 = AlphaEigenVal[NumAlphaElec] - O2_SOMO

                output["stable2o2"] = [OxidizedbyO2,ReducedbyO2]
  
        ##########################################################################
        # Index of Links
        ##########################################################################
        # Index = 0           : (always blank)
        # Index = 1           : Ground state             [if opt==1]
        # Index = 1+nmr       : NMR chemical shift of S0       [if nmr==1]
        # Index = 1+nmr+ipe   : Ionization potential      [if ipe==1]
        # Index = 1+nmr+ipe+epe   : Electronic affinity      [if eae==1]
        # Index = 1+nmr+ipe+epe+pne   : neutrization energy from cation  [if pne==1]
        # Index = 1+nmr+ipe+epe+pne+nne  : neutrization energy from anion  [if nne==1]
        # Index = 1+nmr+ipe+epe+pne+nne+1 : Virtical excitation (S0 -> S1) [uv or fluor or tadf]
        # Index = 2+nmr+ipe+epe+pne+nne+1 : Optimization of S1             [fluor or tadf] 
        # Index = 3+nmr+ipe+epe+pne+nne+1 : Optimization of T1             [tadf]
        ##########################################################################
  

        if nmr == 1:
            Element = []
            ppm = []
  
            Index = 1+nmr
            nmr_lines = Links[Index].splitlines()
            Links[Index]=""

            for line in nmr_lines:
                if line.find("Isotropic =  ") >=0:
                    line_Info = line.split()
                    #print (line_Info[1])
                    Element.append(line_Info[1])
                    #print(line_Info[4])
                    ppm.append(float(line_Info[4]))

        # calculating chemical shift for H, C, or Si

            for i in range(len(Element)):
                if Element[i] =="H" or  Element[i] =="C" or  Element[i] =="Si":
                    ppm[i] = GaussianRunPack.AtomInfo.One_TMS_refer(Element[i],  self.functional, self.basis)-ppm[i]
  
        #    return Element, ppm 
            output["nmr"] = [Element, ppm]
  
        #print (output)

        if ipe == 1 or eae == 1:

            try:
                GS_Energy = output["Energy"]
            except KeyError:
                GS_Energy = self.Extract_SCFEnergy(GS_lines)

            print (GS_Energy)

            if ipe == 1:

                Index = 1+nmr+ipe
                IP_lines = Links[Index].splitlines()
                Links[Index]=""

                IP_Energy =  self.Extract_SCFEnergy(IP_lines)

                print (IP_Energy)

#                output["ipe"] = 27.211*(GS_Energy - IP_Energy)

######Normal ionization potential calculation####################
                output["ipe"] = 27.211*(IP_Energy - GS_Energy)
#################################################################

            if eae == 1:

                Index = 1+nmr+ipe+eae
                EA_lines = Links[Index].splitlines()
                Links[Index]=""

                EA_Energy =  self.Extract_SCFEnergy(EA_lines)

                print (EA_Energy)

#                output["eae"] = 27.211*(EA_Energy - GS_Energy)

######Normal electronic affinity calculation####################
                output["eae"] = 27.211*(GS_Energy - EA_Energy)
#################################################################


        if pne >= 1 or nne >= 1:

            if pne == 1:

                Index = 1+nmr+ipe+eae+pne
                PC_lines = Links[Index].splitlines()
                Links[Index]=""

                VNP_lines = Links[Index+1].splitlines()
                Links[Index+1]=""
                pne += 1

                PC_Energy =  self.Extract_SCFEnergy(PC_lines)
                VNP_Energy =  self.Extract_SCFEnergy(VNP_lines)

                print (PC_Energy)
                print (VNP_Energy)

                output["pne"] = 27.211*(PC_Energy - VNP_Energy)

            if nne == 1:

                Index = 1+nmr+ipe+eae+pne+nne 
                NC_lines = Links[Index].splitlines()
                Links[Index]=""

                VNN_lines = Links[Index+1].splitlines()
                Links[Index+1]=""
                nne +=1

                NC_Energy =  self.Extract_SCFEnergy(NC_lines)
                VNN_Energy =  self.Extract_SCFEnergy(VNN_lines)

                print (NC_Energy)
                print (VNN_Energy)

#                output["nne"] = 27.211*(NC_Energy - VNN_Energy)

######Normal electronic affinity calculation####################
                output["nne"] = 27.211*(VNN_Energy - NC_Energy)
#################################################################

  
        if uv == 1 or fluor == 1 or tadf == 1:
            Index = 1+nmr+ipe+eae+pne+nne+1 
            lines = "" if Index >= n else Links[Index].splitlines()
            Found, Egrd, Eext, WaveLength, V_OS = self.Extract_ExcitedState(lines)
#            print (Found)
#            print (Egrd)
#            print (Eext)
            output["uv"] = [WaveLength, V_OS]
            Links[Index]=""
  
        if fluor == 1 or tadf == 1:
            Index = 2+nmr+ipe+eae+pne+nne+1
            lines = "" if Index >= n else Links[Index].splitlines()
            S1_Found, S1_Egrd, S1_Eext, S1_WaveLength, S1_V_OS = self.Extract_ExcitedState(lines)
            output["S1 Total Energy"] = S1_Eext
            output["S1 Wavelength and Oscillator strengths"] = [S1_WaveLength, S1_V_OS]
            Links[Index]=""
  
        if tadf == 1:
            Index = 3+nmr+ipe+eae+pne+nne+1
            lines = "" if Index >= n else Links[Index].splitlines()
            T1_Found, T1_Egrd, T1_Eext, T1_WaveLength, T1_V_OS = self.Extract_ExcitedState(lines)
            output["T1 Total Energy"] = T1_Eext
            output["T1 Wavelength and Oscillator strengths"] = [T1_WaveLength, T1_V_OS]
            TADF_Eng = 0.0
            if S1_Found and T1_Found:
               TADF_Eng = S1_Eext - T1_Eext
            output["Energy difference (S1-T1)"] = TADF_Eng
            Links[Index]=""
  
       # Debug print for SplitLinks
       #      for i in range(0,n):
       #          f = open(infilename+"."+str(i), 'w')
       #          for s in Sections[i]:
       #              f.write(s)
       #          f.close() 
  
        return output

    def MakeSolventLine(self):

        s = ''
        s_solvent = ''

        try:
            float(self.solvent)
        except ValueError:
            print ('Solvent effect is included by PCM')
            s = s +'SCRF=(PCM, solvent=' + self.solvent + ')'
            s = s +'\n'

        else:
            if self.solvent == '0':
                return s, s_solvent

            else:
                print ('Solvent effect is included by PCM')
                s = s +'SCRF=(PCM, solvent=Generic, Read)'
                s = s +'\n'
                s_solvent = s_solvent+'EPS='+self.solvent
                s_solvent = s_solvent+'\n'
                s_solvent = s_solvent+'Radii=UA0'
                s_solvent = s_solvent+'\n'
                s_solvent = s_solvent+'\n'

        return s, s_solvent

    def MakeLinkTD(self, line_chk, line_method, State, SCRF, SCRF_read, Opt):

        self.MakeSolventLine()

        line_method_TD = 'TD(Nstate=10, '+ State + ')'
        line_readMOGeom = 'Geom=AllCheck Guess=Read'
        s = ''
        s = s + '--Link1--'
        s = s + '\n'
        if self.mem != '':
            line_mem = '%mem='+str(self.mem)
            s = s + line_mem      
            s = s + '\n'
        if self.nproc > 1 :
            line_proc = '%nproc='+str(self.nproc)
            s = s + line_proc
            s = s + '\n'
        s = s + line_chk
        s = s + '\n'
        s = s + line_method
        s = s + '\n'
        s = s + line_method_TD
        s = s + '\n'
        s = s + line_readMOGeom
        s = s + '\n'
        s = s + SCRF
        if Opt == True:
            s = s + " Opt"
            s = s + '\n'
        s = s + '\n' 

        s = s + SCRF_read

        return s


    def run_gaussian(self):
        
        infilename = self.in_file

        option_line = self.value    

        options = option_line.split()

        opt = 0
        nmr = 0
        uv = 0
        energy = 0
        gap = 0
        dipole = 0
        deen = 0
        stable2o2 = 0
        fluor = 0
        tadf = 0
        ipe = 0
        eae = 0
        pne = 0
        nne = 0

        for i in range(len(options)):
#            print(options[i])
            option = options[i]
            if option == 'opt':
                opt = 1
                print ('opt = ', opt)
            elif option == 'nmr':
                nmr = 1
                #print ('nmr')
            elif option == 'uv':
                uv = 1
                #print ('uv')
            elif option == 'energy':
                energy = 1
                #print ('energy')
            elif option == 'homolumo':
                gap = 1
                #print ('HOMO/LUMO')
            elif option == 'dipole':
                dipole = 1
                #print ('dipole')
            elif option == 'deen':
                deen = 1
                #print ('Decomposition energy')
            elif option == 'stable2o2':
                stable2o2 = 1
                #print ('Stability to O2')
            elif option == 'fluor':
                fluor = 1
                #print ('Fluorescence')
            elif option == 'tadf':
                tadf = 1
                #print ('Thermally Activated Delayed Fluorescence')
            elif option == 'ipe':
                ipe = 1
                #print ('Ionization potential')
            elif option == 'eae':
                eae = 1
                #print ('Electronic affinity')
            elif option == 'pne':
                pne = 1
                #print ('Neutraization energy from cation')
            elif option == 'nne':
                nne = 1
                #print ('Neutraization energy from anion')
            else:
                print('invalid option: ', option)

        PreGauInput = infilename.split('.')
    
        GauInputName = PreGauInput[0]+'.com'    

####File type of input?######################
        ReadFromchk = 0 
        ReadFromsdf = 0 

        if PreGauInput[1] == "sdf":
            ReadFromsdf = 1 
            Mol_atom, X, Y, Z, TotalCharge, SpinMulti = self.read_sdf()
        elif PreGauInput[1] == "chk":
            ReadFromchk = 1 
        else:
            print ("Invalid input file")
############################################

        line_chk = '%chk='+PreGauInput[0]
        line_oldchk = '%Oldchk='+PreGauInput[0]

        line_method = '#'+self.functional+'/'+self.basis

        line_comment = infilename

        ofile = open(GauInputName ,'w')

#####For reading geomerty and MO from checkpoint file###
        line_readMOGeom = 'Geom=AllCheck Guess=Read'
        line_readOnlyMOGeom = 'Geom=CheckPoint Guess=Read'
        line_readGeom = 'Geom=Checkpoint'
########################################################

        if self.mem != '':
            line_mem = '%mem='+str(self.mem)
            ofile.write(line_mem+'\n')

        if self.nproc > 1 :
            line_proc = '%nproc='+str(self.nproc)
            ofile.write(line_proc+'\n')
    
        ofile.write(line_chk+'\n')
        ofile.write(line_method+'\n')
        if opt == 1:
            ofile.write('Opt\n')


#######Solvent effect############################
        SCRF, SCRF_read = self.MakeSolventLine()
        ofile.write(SCRF)

#####Reading Geometry and MO from Checkpoint file
        if ReadFromchk == 1:
            ofile.write(line_readMOGeom+'\n')
#################################################

        ofile.write('\n')

#####Reading Geometry from sdf file#####

        if ReadFromsdf == 1:

            ofile.write(line_comment+'\n')
            ofile.write('\n')

            ofile.write('%5d %5d \n' % (TotalCharge, SpinMulti))

            for j in range(len(Mol_atom)):
                ofile.write('%-4s % 10.5f  % 10.5f  % 10.5f \n'
                 % (Mol_atom[j],X[j], Y[j], Z[j]))

            ofile.write('\n')
            ofile.write(SCRF_read)
#######################################

        if nmr == 1:

            ofile.write('--Link1--\n')

            if self.mem != '':
                line_mem = '%mem='+str(self.mem)
                ofile.write(line_mem+'\n')


            if self.nproc > 1 :
                line_proc = '%nproc='+str(self.nproc)
                ofile.write(line_proc+'\n')

            ofile.write(line_chk+'\n')
            ofile.write(line_method+'\n')

            line_method_nmr = 'NMR'
            ofile.write(line_method_nmr+'\n')
            ofile.write(SCRF)

            ofile.write(line_readMOGeom+'\n')
            ofile.write('\n')
            ofile.write(SCRF_read)

        if ipe == 1:

            ofile.write('--Link1--\n')

            if self.mem != '':
                line_mem = '%mem='+str(self.mem)
                ofile.write(line_mem+'\n')


            if self.nproc > 1 :
                line_proc = '%nproc='+str(self.nproc)
                ofile.write(line_proc+'\n')

            ofile.write(line_chk+'_IP\n')
            ofile.write(line_oldchk+'\n')
            ofile.write(line_method+'\n')
            ofile.write(SCRF)

            ofile.write(line_readOnlyMOGeom+'\n')
            ofile.write('\n')

            ofile.write('ionization potential calculation\n')
            ofile.write('\n')

            ofile.write('%5d %5d \n' % (TotalCharge+1, SpinMulti+1))
            ofile.write('\n')
            ofile.write(SCRF_read)

        if eae == 1:

            ofile.write('--Link1--\n')

            if self.mem != '':
                line_mem = '%mem='+str(self.mem)
                ofile.write(line_mem+'\n')


            if self.nproc > 1 :
                line_proc = '%nproc='+str(self.nproc)
                ofile.write(line_proc+'\n')

            ofile.write(line_oldchk+'\n')
            ofile.write(line_chk+'_EA\n')
            ofile.write(line_method+'\n')
            ofile.write(SCRF)

            ofile.write(line_readOnlyMOGeom+'\n')
            ofile.write('\n')

            ofile.write('electronic affinity calculation\n')
            ofile.write('\n')

            ofile.write('%5d %5d \n' % (TotalCharge-1, SpinMulti+1))
            ofile.write('\n')
            ofile.write(SCRF_read)

        if pne == 1:

            ofile.write('--Link1--\n')

            if self.mem != '':
                line_mem = '%mem='+str(self.mem)
                ofile.write(line_mem+'\n')


            if self.nproc > 1 :
                line_proc = '%nproc='+str(self.nproc)
                ofile.write(line_proc+'\n')

            ofile.write(line_chk+'_PC\n')
            ofile.write(line_oldchk+'\n')
            ofile.write(line_method+' opt'+'\n')
            ofile.write(SCRF)

            ofile.write(line_readOnlyMOGeom+'\n')
            ofile.write('\n')

            ofile.write('neutrization energy calculation from a cation\n')
            ofile.write('\n')

            ofile.write('%5d %5d \n' % (TotalCharge+1, SpinMulti+1))
            ofile.write('\n')
            ofile.write(SCRF_read)

            ofile.write('--Link1--\n')

            if self.mem != '':
                line_mem = '%mem='+str(self.mem)
                ofile.write(line_mem+'\n')

            if self.nproc > 1 :
                line_proc = '%nproc='+str(self.nproc)
                ofile.write(line_proc+'\n')

            ofile.write(line_chk+'_PC\n')
            ofile.write(line_oldchk+'\n')
            ofile.write(line_method+'\n')
            ofile.write(SCRF)

            ofile.write(line_readGeom+'\n')
            ofile.write('\n')

            ofile.write('neutrization energy calculation from a cation\n')
            ofile.write('\n')

            ofile.write('%5d %5d \n' % (TotalCharge, SpinMulti))
            ofile.write('\n')
            ofile.write(SCRF_read)


        if nne == 1:

            ofile.write('--Link1--\n')

            if self.mem != '':
                line_mem = '%mem='+str(self.mem)
                ofile.write(line_mem+'\n')


            if self.nproc > 1 :
                line_proc = '%nproc='+str(self.nproc)
                ofile.write(line_proc+'\n')

            ofile.write(line_oldchk+'\n')
            ofile.write(line_chk+'_NC\n')
            ofile.write(line_method+' opt'+'\n')
            ofile.write(SCRF)

            ofile.write(line_readOnlyMOGeom+'\n')
            ofile.write('\n')

            ofile.write('neutrization energy calculation from an anion\n')
            ofile.write('\n')

            ofile.write('%5d %5d \n' % (TotalCharge-1, SpinMulti+1))
            ofile.write('\n')
            ofile.write(SCRF_read)

            ofile.write('--Link1--\n')

            if self.mem != '':
                line_mem = '%mem='+str(self.mem)
                ofile.write(line_mem+'\n')

            if self.nproc > 1 :
                line_proc = '%nproc='+str(self.nproc)
                ofile.write(line_proc+'\n')

            ofile.write(line_oldchk+'\n')
            ofile.write(line_chk+'_NC\n')
            ofile.write(line_method+'\n')
            ofile.write(SCRF)

            ofile.write(line_readGeom+'\n')
            ofile.write('\n')

            ofile.write('neutrization energy calculation from an anion\n')
            ofile.write('\n')

            ofile.write('%5d %5d \n' % (TotalCharge, SpinMulti))
            ofile.write('\n')
            ofile.write(SCRF_read)


        if uv == 1 or fluor==1 or tadf == 1:
            sTD = self.MakeLinkTD(line_chk, line_method, "Singlet", SCRF, SCRF_read, False)
            ofile.write(sTD) 

        if fluor == 1 or tadf == 1:
            sTD = self.MakeLinkTD(line_chk, line_method, "Singlet",  SCRF, SCRF_read, True)
            ofile.write(sTD) 

        if tadf == 1:
            sTD = self.MakeLinkTD(line_chk, line_method, "Triplet",  SCRF, SCRF_read, True)
            ofile.write(sTD) 

        ofile.write('\n') 
        ofile.close()


############Run Gaussian##############################

        if os.path.isdir(PreGauInput[0]):
            shutil.rmtree(PreGauInput[0])
        os.mkdir(PreGauInput[0])
        shutil.move(GauInputName, PreGauInput[0])

        if ReadFromchk == 1 :
            shutil.move(infilename, PreGauInput[0]) 

        os.chdir(PreGauInput[0])

        tmp_test=open(GauInputName,'r')
        #tmp_test.close()
        #route = tmp_test.readline()

        #print (route)
        
        subprocess.call(["g16", PreGauInput[0]])

        logfile = PreGauInput[0]+'.log'

        output_dic = self.Extract_values(logfile,opt,nmr,uv,energy, ipe, eae, pne, nne, gap,dipole,deen,stable2o2, fluor,tadf)

        os.chdir("..")

        tmp_test.close()
        return(output_dic)

    def SplitLinks(self, logfile):
        f = open(logfile, 'r')
        lines = f.readlines()
        f.close()

        Links = []
        Link = ""

        for line in lines:
            if line.find("Initial command:")>0:
                Links.append(Link)
                Link = ""

            Link = Link + line

        Links.append(Link)
        return Links

#####################To get cartesian coordinates from sdf file###############
    def read_sdf(self):

        ifile = open(self.in_file, 'r')

        sdfinfo = Chem.SDMolSupplier(self.in_file, removeHs=False)

        for mol in sdfinfo:
 #           mol = Chem.AddHs(mol)
            N = mol.GetNumAtoms()
            N_Bond = mol.GetNumBonds()

     #   print ("Number of atoms = ", N)
     #   print ("Number of bonds = ", N_Bond)

        count = 0

        X = []
        Y = []
        Z = []
        element_symbol = []

        Bond_pair1 = []
        Bond_pair2 = []
        Bond_type = []

        TotalCharge = 0 
        CHG_atom = []
        CHG = []

        for line in ifile:
            if count == 0:
                Header1 = line

                count += 1
                continue

            if count == 1:
                Header2 = line

                count += 1
                continue

            if count == 2:
                Header3 = line

                count += 1
                continue

            if count == 3:
                a = line.split()
        #        N = int(a[0])
        #        N_Bond = int(a[1])

                count += 1
                continue

            if 3 < count <= N+4:
                i_atom = line.split()
                if len(i_atom) != 0:
                    X.append(float(i_atom[0]))
                    Y.append(float(i_atom[1]))
                    Z.append(float(i_atom[2]))
                    element_symbol.append(i_atom[3])
    
                count += 1
                continue

            if N+4 < count <= N+N_Bond+3 :
                bond_info = line.split()
                #print (bond_info)
                bond_info = line.split()
                Bond_pair1.append(int(bond_info[0]))
                Bond_pair2.append(int(bond_info[1]))
                Bond_type.append(int(bond_info[2]))

                count +=1
                continue

            if count > N+N_Bond+3:
                mol_info = line.split()
                #print (mol_info)
                if (mol_info[0] == "M"):
                    if (mol_info[1] == "END"):
                        break
                    if (mol_info[1] == "CHG"):
                        Num_CHGInfo = int(mol_info[2])
                        for k in range(Num_CHGInfo):
                            CHG_atom.append(int(mol_info[3+2*k]))
                            CHG.append(int(mol_info[4+2*k]))
                            TotalCharge += int(mol_info[4+2*k])
                else:
                    print("The sdf file is invalid!")
                    sys.exit()

                count +=1

#Copy to array of numpy###########################

        Mol_atom = []
        Mol_CartX = zeros(N)
        Mol_CartY = zeros(N)
        Mol_CartZ = zeros(N)

        CHG_atom = array(CHG_atom)
        CHG = array(CHG)

        for j in range(N):
            Mol_CartX[j] = X[j] 
            Mol_CartY[j] = Y[j] 
            Mol_CartZ[j] = Z[j] 
            Mol_atom.append(element_symbol[j])

    #del element_symbol[N:TotalStep]

        del element_symbol[:]
        del X[:]
        del Y[:]
        del Z[:]

###For debug#######################################
#
#        print (Mol_atom)
#        print (N)
#        print (len(Mol_CartX))
#
        print('Reading the sdf file has finished')

###Calculating the total number of electrons#################
        TotalNum_electron = 0 

        for j in range(N):
            if (len(CHG_atom) != 0):
                Judge = CHG_atom-j
                if (any(Judge) == 0):
                    TotalNum_electron += GaussianRunPack.AtomInfo.AtomicNumElec(Mol_atom[j])-CHG[where(Judge == 0)]
                else:
                    TotalNum_electron += GaussianRunPack.AtomInfo.AtomicNumElec(Mol_atom[j])
            else:
                TotalNum_electron += GaussianRunPack.AtomInfo.AtomicNumElec(Mol_atom[j])

        print('Total number of electron: %7d ' %  (TotalNum_electron))

        if (TotalNum_electron%2==0):
            print ("This system is a closed shell!")
            SpinMulti = 1
        else:
            print ("This system is a open shell!")
            SpinMulti = 2

        return Mol_atom, Mol_CartX, Mol_CartY, Mol_CartZ, TotalCharge, SpinMulti

        
