class1_read.py

# In[1]:

import numpy as np
import os
import sys
import re
import shutil
from pathlib import Path
import pymatgen as pmg
from pymatgen.core import Structure 
from pymatgen.io.vasp.outputs import BSVasprun, Vasprun, Eigenval
from pymatgen.electronic_structure.plotter import BSPlotter, BSDOSPlotter, DosPlotter

from collections import Counter
#sys.path.append(os.environ['SCRIPT'])
from class0_functions1 import read_incar, find_files, parent_folder

# read_file_values: Read data
# rotate_read: Read the same data in many folders

# In[2]:

class read_file_values:
    # No POSCAR here. USE pymatgen to read POSCAR, KPOINTS!
    def __init__(self, folder):
        #if folder[-1] == '/':
        #    self.folder = folder
        #else:
        #    self.folder = folder + '/'
        #assert folder[-1] == '/', 'folder has to end with /'
        self.folder = str(Path(folder)) + '/'
        #hierdirecs = parent_folder(self.folder, hier=10)  # to find 1st-class materials folder
        #assert len(hierdirecs) <= 4
        #self.hierdirecs = hierdirecs
        #self.hier1direc = self.hierdirecs[0]    # 1st-class folder Eg. '/home/yubi/work/berylliumO/'
        #self.hier2direc = self.hierdirecs[1]    # 2nd-class folder Eg. '/home/yubi/work/berylliumO/wurtzite_2_Vo_HSE/'
        #self.sources_direc= self.hier1direc+'sources/' # Eg. '/home/yubi/work/berylliumO/sources/'

    def incar(self, para, var_format=float):
        # read the values of certain variables (came as a list)
        # return a list
        assert isinstance(para, list), 'read a list of variables: should keep the variable a list'
        dictionary = read_incar(self.folder)
        values=np.zeros(len(para))
        for j in range(len(para)):
            values[j]=var_format(dictionary[para[j]])
        return values

    def oszicar(self, folder=''):
        '''
        read energy from oszicar
        folder=self.folder by default
        '''
        energy=False
        if folder == '':
            folder = self.folder
        with open(folder+'OSZICAR', "r") as f:
            lines = f.readlines()
        for i in range(len(lines)-1, 0, -1):
            ll=lines[i].split()
            if 'E0=' in ll:
                energy = float(ll[ll.index('E0=')+1])
                break
        assert energy!=False,'This folder might not finish calculation because OSZICAR file does not have energy.'
        return energy

    def bandgap(self, fullinfo=False):
        run = BSVasprun(self.folder+"vasprun.xml", parse_projected_eigen=True)
        bs = run.get_band_structure(self.folder+"KPOINTS")
        bandgapinfo=bs.get_band_gap() # {'direct': True, 'energy': 11.126, 'transition': '(0.000,0.000,0.000)-(0.000,0.000,0.000)'}
        # bs.get_direct_band_gap()
        # bs.get_direct_band_gap_dict()
        energy=np.round(bandgapinfo['energy'],4) # keyword: energy
        print('direct bandgap=%s '%(bandgapinfo['direct']))
        if fullinfo:
            return bandgapinfo # bandgapinfo is a dictionary
        else:
            return energy 

    def cbm_vbm(self, c_v=0 , fullinfo=False):
        # if c_v = 0 ,return cbm; if c_v=1, return vbm
        run = BSVasprun(self.folder+"vasprun.xml", parse_projected_eigen=True)
        bs = run.get_band_structure(self.folder+"KPOINTS")
        if c_v == 0:
            info=bs.get_cbm() # there is the kpoint related to CBM
            # {'band_index': defaultdict(<class 'list'>, {<Spin.up: 1>: [192], <Spin.down: -1>: [192]}), 'kpoint_index': [0], 'kpoint': <pymatgen.electronic_structure.bandstructure.Kpoint>, 'energy': 13.286, 'projections': {}}
        elif c_v == 1:
            info=bs.get_vbm() # Valence band maximum
        else:
            print('Error! choose CBM (c_v=0) or VBM (c_v=1)')
        energy=np.round(info['energy'],4) # keyword: energy
        print('%senergy=%s'%(['CBM','VBM'][c_v], energy))
        if fullinfo:
            return info # bandgapinfo is a dictionary
        else:
            return energy 

    def fermienergy(self):
        run = BSVasprun(self.folder+"vasprun.xml", parse_projected_eigen=True)
        bs = run.get_band_structure(self.folder+"KPOINTS")
        return bs.efermi
    
    def read_single_eigenE(self,eigenvaltype='ho',spinpolarized=False):
        run = Eigenval(self.folder+'EIGENVAL')
        #run.eigenvalues #  dict of {(spin): np.ndarray(shape=(nkpt, nbands, 2))}.
        assert spinpolarized == False, 'Only read non-polarized spin values for now'
        from pymatgen import Spin
        spinupEnergy= run.eigenvalues[Spin.up] # format is (nkpt, nbands, 2) # len(2) is energy eigenvalue + occupation
        if eigenvaltype=='ho':
            # only look at spin up
            ## spinupEnergy[0] eigenvalues of gamma point, the first K point
            ## spinupEnergy[0,:,1] = [Occband1, Occband2,...]; 1 refers to occupations at gamma point
            ## spinupEnergy[0,:,0] = [Eband1, Eband2,...]; 0 refers to energies at gamma point
            numocc=sum(spinupEnergy[0,:,1]>0.1) # this number of bands are occupied
            print('The high occupied state at Gamma is ', spinupEnergy[0,numocc-1])
            # spinupEnergy[0,numocc-1,0] will be the energy of the highest occupied state
            return spinupEnergy[0,numocc-1,0] #np.round(,2)
            

    def lattice_para(self, structure=''):
        # return a list of unique lattice parameters
        # for certain structures, only need some parameters
        stru=Structure.from_file(self.folder+'CONTCAR')
        latt=stru.lattice
        a,b,c=latt.abc
        a,b,c = np.round([a,b,c], 4)
        print('From CONTCAR: a=%s,b=%s,c=%s' % (a,b,c))
        if structure=='wurtzite':
            assert np.round(a,5)==np.round(b,5), 'Error! not %s structure' % (structure)
            return [a,c], ['a','c']
        elif structure=='rock-salt':
            assert np.round(a,5)==np.round(b,5) and np.round(a,5)==np.round(c,5), 'Error! not %s structure' % (structure)
            return [a], ['a']
        else: # no symmetry, return all lattice parameters
            return [a,b,c], ['a','b','c']

    def get_u(self):
        u0=0.375
        stru=Structure.from_file(self.folder+'CONTCAR')
        speci=stru.species
        num_atoms=len(speci)
        #choose two atoms: atom 0 and the one in its z direction
        atom0=stru[0]
        for i in range(1,num_atoms):
            atom1=stru[i]
            diff_corrds= (atom0.coords - atom1.coords)[:2] # [x,y,z] -> [x,y]
            if np.all(np.round(diff_corrds,5)==[0.,0.]): # [x,y]==[0,0] in the z direction of atom 0
                u= round(stru.get_distance(0,i)/stru.lattice.c,7)
                print('From CONTCAR: atomX=%s \n\t\tatomY=%s'%(atom0,atom1))
                return u,(u-u0)/u0
                break


    def outcar(self, var='NELECT'):
        # Can use pymatgen to read it!!!
        with open(self.folder+'OUTCAR', 'r') as f:
            lines= f.readlines()
            found=0
            reNELECT = re.compile(var)
            for i in range(len(lines)):
                if reNELECT.search(lines[i]):
                    # if NELECT doesn't exist, it will give an error
                    nelect = float(lines[i].split('=')[1].split()[0]) # number of electrons
                    found=1
            assert found ==1, 'Parameter %s doesn\'t exist in OUTCAR' % var
        return nelect

    def outcar_dielec_const(self, tensor=False):
        '''
        This is to read dielectric constants from OUTCAR in unit cell with different parameters
        Need to save all the data in a single file to be used in other situations
        '''
        print('read OUTCAR from a 4-th class subfolder: dielec_eps/')
        filetoread='dielec_elec-part_eps/OUTCAR'
        if tensor:
            filetoread='OUTCAR'
        with open(self.folder+filetoread, 'r') as f:
            lines= f.readlines()
            found=0
            string='MACROSCOPIC STATIC '
            restring = re.compile(string)
            for i in range(len(lines)):
                if restring.search(lines[i]):
                    # if dielectric tensor doesn't exist, it will give an error
                    tensorstr=lines[i+2:i+5]
                    print('The dielectric tensor is %s' % tensorstr )
                    found=1
            assert found ==1, 'DIELECTRIC TENSOR doesn\'t exist in OUTCAR'
        diagonal=np.zeros(3)
        for j in range(len(tensorstr)):
            diagonal[j]=np.round(float(tensorstr[j].split()[j]), 6)
        print('diagonal eps matrix:', diagonal)
        if tensor:
            return diagonal # a list of three diagonal values
        else:
            eps=np.average(diagonal)
            return eps


    def read_formation_enthalpy(self,compound_folder, return_unit_formula=True):
        '''
        formation enthalpy for the unit of compound
        eg. BeO rock-salt has 4Be and 4O. The calculation should also return the energy of 4BeO
        '''
        compound_folder = str(Path(compound_folder))+'/'
        #if Path(compound_folder) == Path(self.folder):
        #    native = True ## compound_folder = self.folder # use the native compound like BeO by default. 
        #else:
        #    native = False # If another formation enthalpy, say Li2O, enter the folder to that compound
        # compound energy
        read_fil=read_file_values(compound_folder)
        compound_energy=read_fil.oszicar() #/formula[-1]

        # read compound species and formula to determine the folders of elemental phases
        struc=Structure.from_file(compound_folder+'POSCAR')
        speci=struc.species # speci[i] is Element Be; speci[i].name is 'Be'
        unique_speci=list(set(speci))
        coun=Counter(speci) # coun[speci[i]]=2 is the number of element atom
        num_elements=len(coun) # the number of species
        energies = np.zeros(num_elements) # initialize element energy list
        print(' compound formula=%s energy=%s' %(struc.formula,compound_energy))
        folder_header='energyf_'
        #if native:
        #    folder_header = 'enthalpyf_' # want ele_folder = enthalpy_Be/enthalpy_O
        #else:
        #    # split folder name to get impurity specie name 'energyf_Li_O'.split('_') = ['energyf', 'Li', 'O'], read Li
        #    impurity_specie = compound_folder.split('_')[1] # 'Li'
        #    folder_header = 'energyf_' # want ele_folder = energyf_Li/energyf_O

        # element energies
        formula=[]
        dict_elementformula= {}
        unique_speci_name = []
        for i in range(num_elements): # go over unique_speci
            # get the number of element atom from compound formula
            unique_speci_name.append(unique_speci[i].name)
            formula.append(int(coun[unique_speci[i]])) # coun[speci[i]]=2 is the number of element atom; formula=[2,2] or [4,4]
            dict_elementformula[unique_speci[i].name] = formula[-1] # compound formula: ith-species is unique_speci[i].name
            # element calculation
            ele_folder=self.folder+folder_header+ unique_speci[i].name+'/' # enthalpyf_Be/, enthalpyf_O/ #elements[i]
            ele_struc=Structure.from_file(ele_folder+'POSCAR')
            ele_formula=ele_struc.formula # like 'Be2','O2'; formula is to find the number of atoms in element[i] of enthalpyf_element[i]
            ele_formula_num=float(ele_formula.split(unique_speci[i].name)[-1]) # the number of atoms is 'Be2'.split('Be')[-1]
            read_fil=read_file_values(ele_folder) # initialize a class
            energies[i]= np.round(read_fil.oszicar() / ele_formula_num * formula[i],3) # make the number of element atoms agree with the compound formula
            print('  In %s, element %s with %s atoms has energy=%s' % (ele_folder.split('/')[-2],unique_speci[i].name, formula[i], energies[i]))
        enthalpy=np.round(compound_energy - np.sum(energies),3) # round to 3 decimal number
        if return_unit_formula:
            # return formation enthalpy per formula unit
            from class0_functions2 import hcf, get_printformula
            # calculate highest common factor and divide enthalpy by hcf
            formula_hcf=hcf(*formula) # hcf=highest common factor
            print('enthalpyf=%s for formula=%s and hcf=%s' % (enthalpy, get_printformula(dict_elementformula), formula_hcf))
            enthalpyf_unitformula = np.round(enthalpy / formula_hcf,3) # round to 3 decimals
            formula_unit =np.array(formula) / formula_hcf
            # edit number of atoms in dictionary
            elements = dict_elementformula.keys()
            for ele in elements:
                dict_elementformula[ele]=dict_elementformula[ele]/formula_hcf # make sure the formula is an integer
            print('enthalpy of formation per unit formula is %s eV\n' % (enthalpyf_unitformula) )
            return enthalpyf_unitformula, formula_unit,unique_speci_name, dict_elementformula
        else:
            from class0_functions2 import get_printformula
            print('enthalpy of formation for formula=%s is %s eV\n' % (get_printformula(dict_elementformula), enthalpy))
            formula = np.array(formula)
            return enthalpy, formula,unique_speci_name, dict_elementformula

    def read_delta_miu4energyf(self, condition,impurity_atomnames=[]):
        '''
        get delta miu that is used for formation energy
            condition = 'O-rich' # the condition in formation energy
            impurity_atomnames = ['Li','F','H']  # a list of impurity
        '''
        assert condition[-5:] == '-rich', 'assume condition must have one atom is rich'
        
        # (1): read formation enthalpy for native atoms' compound
        compound_native_folder = self.folder # read native material formation enthalpy by default
        enthalpyf_native,formula_native,speci_name_native,dict_elementformula_native = self.read_formation_enthalpy(compound_native_folder, return_unit_formula=True)
        formula_native = np.array(formula_native) # [1,1] for beryllium oxide with formula one Be and one O
        speci_name_native = np.array(speci_name_native) # ['O','Be']

        # (2): calculate delta miu (chemical potential) for native atoms # see tutorial_python_solveEq.py for how to solve equations
        # a. analyze condition to get which atom is rich 
        richatom=condition.split('-rich')[0] # get 'O' atom is rich
        # b. form a binary linear equation to be solve by numpy like x=np.linalg.solve(A,b) is the solution to Ax=b where A is coefficient matrix, x b are column vectors
        # b.1 get the second row in coefficient matrix, given by the rich atom condition
        rich_condition = np.zeros(len(speci_name_native))
        rich_condition[speci_name_native == richatom] = 1 # get the rich condition like 1*delta_O + 0*delta_Be = 0 
        # which is (1 0)*(delta_O delta_Be) = 0 here (1 0) is the second row in the coefficient matrix
        # b.2 form coefficient matrix and vector b
        A = np.stack([formula_native,rich_condition]) # the first row in A is the (1 1)*(delta_O delta_Be) = enthalpyf, the second row is the rich atom condition
        b = np.array([enthalpyf_native,0]).T
        # c. solve delta_O and delta_Be
        delta_miu_native = np.linalg.solve(A,b)
        # d. form a dictionary for delta miu
        dict_delta_miu = {}
        for i in range(len(speci_name_native)):
            dict_delta_miu[speci_name_native[i]] = delta_miu_native[i]
        print('Under condition %s, the native atoms have delta miu per atom = %s\n\n' % (condition, dict_delta_miu))

        # (3): solve delta miu for each impurity atom
        if len(impurity_atomnames)>=1:
            for impurity_atomname in impurity_atomnames:
                delta_mius_impurity = [0] # delta miu should not be positive. If positive, take it to 0
                compound_impurity_folder_header = 'energyf_' + impurity_atomname + '_' # 'energyf_Li_'
                compound_impurity_folders_all = find_files(self.folder, header=compound_impurity_folder_header, var='', remove=False)
                # like [energy_Li_O/, energy_Li_Be/]
                for compound_impurity_folder in compound_impurity_folders_all:
                    ## a. get constrains for impurity atom
                    #compound_impurity_folder = self.folder+'energyf_' + impurity_atomname + '_' + native_atom_i # assume only a two-atom compound is formed like Li2O, but not Li2BeO2
                    native_atom_i = compound_impurity_folder[:-1].split(compound_impurity_folder_header)[1].split('_')[0] #[:-1] is to remove '/'
                    compound_impurity_folder = self.folder+compound_impurity_folder
                    if os.path.isdir(compound_impurity_folder):
                        #read_fil = read_file_values(pwd)
                        enthalpyf_impurity,formula_impurity,speci_name_impurity,dict_elementformula_impurity = self.read_formation_enthalpy(compound_impurity_folder, return_unit_formula=True)
                        # b. calculate delta miu by constrain
                        # the condition is (1 2)*(delta_O delta_Li)=enthalpyf_impurity, delta_Li = (enthalpyf_impurity - delta_O*formula_O)/formula_Li
                        delta_mius_impurity.append( np.round((enthalpyf_impurity - dict_delta_miu[native_atom_i] * dict_elementformula_impurity[native_atom_i])/dict_elementformula_impurity[impurity_atomname],3) )
                    else:
                        print('folder %s does not exist' % (compound_impurity_folder)) # this line is no longer needed because folder names come from find_files now
                # c. select the smallest delta miu (chemical potential) for impurity atoms
                print('All constrains for impurity atom %s gives delta miu=%s' %(impurity_atomname, delta_mius_impurity) )
                dict_delta_miu[impurity_atomname] = np.round(np.min(delta_mius_impurity),3)
                print('impurity atom %s has delta miu=%s\n' %(impurity_atomname,dict_delta_miu[impurity_atomname]) )
        print('\tUnder condition %s, All delta miu = %s' % (condition,dict_delta_miu))
        return dict_delta_miu


    def eigenval_ho_lu(self,ho_lu=0):
        '''
        Read EIGENVAL to get the energies of the highest occupied (ho) states
        return the highest occupied energy and nband
        energyHOLO is used in Generalized-Koopmans' condition
        nband refers to the band where that electron occupied -> used in PCHARG calculation
        ho_lu whether we want highest occupied (ho) or lowest unoccupied (lu) energies: 0 for ho, 1 for lu
        '''
        ho_lu_energy = False
        with open(self.folder+'EIGENVAL', 'r') as f:
            lines=f.readlines()
            lines=lines[8:]
            for i in range(len(lines)):
                splited= np.array(lines[i].split()).astype(float)
                if np.any(splited[-2:] != np.array([1.,1.])):
                    #print('\n'+''.join(lines[i-1:i+1])[:-1]) # [:-1] is to remove '\n' at the end
                    # needs to modify for more generalized case
                    # so far only works when neutral has double electrons
                    if ho_lu==0: # highest occupied
                        # redefine splited to be the previous band with both electrons occupied
                        splited= np.array(lines[i-1].split()).astype(float) # First get the highested energy in splited. Then modify it if there is even higher energy occupied.
                        ho_lu_energy=max(splited[1],splited[2]) # splited[0] is line number, splited[1] is the first energy, splited[2] is the second energy
                        nband=int(splited[0]) 
                        splitted2=np.array(lines[i].split()).astype(float)
                        # if any state of the next band is occupied
                        # modify the energy to the correct position in the next band
                        if splitted2[-2]==1.: # find one higher occupied state 'lines[i]=193 6.397799 10.874288 1.000000 0.000000\n'
                            ho_lu_energy=splitted2[1] # energy should be 6.397799 
                            nband=int(splitted2[0])
                        elif splitted2[-1]==1:
                            ho_lu_energy=splitted2[2]
                            nband=int(splitted2[0])
                    elif ho_lu==1: # lowest unoccupied
                        splited= np.array(lines[i].split()).astype(float)
                        if splited[-2]==0. and splited[-1]==0.: # find one higher occupied state 'lines[i]=193 6.397799 10.874288 1.000000 0.000000\n'
                            ho_lu_energy=min(splited[1],splited[2])
                        elif splited[-1]==0.: # energy should be 6.397799   ho_lu_energy=splitted2[1]
                            ho_lu_energy=splited[2]
                        elif splited[-2]==0.: # in case the state energies are not in order
                            ho_lu_energy=splited[1]
                        else:
                            print('Error! Not finding the highest-lowest intersection')
                        nband=int(splited[0])
                    break
        assert ho_lu_energy != False, 'This folder might not finish calculation because EIGENVAL file does not have eigenenergies.'
        return [ho_lu_energy,nband]


    def locpot2freysoldt_correction(self,i, fol1,fol2, bulk=True,x_str='AEXX'):
        '''
        # the folders with different defects: fol1 is the reference, should be neutral; fol2 is positive
        Put LOCPOTs in the same folder and do Freysoldt correction
        1-generate a folder (to store LOCPOTS and run Freysoldt)
        2-copy LOCPOT from perfect_supercell
        3-obtain eps value from unit cell calculations
        4-run shell script 'sxdefectalign ..' to calculate Freysoldt correction
        '''
        #assert len(self.hierdirecs) == 3, 'You need to be in a 3rd class folder'
        #the LOCPOT files
        f1=fol1+'LOCPOT'
        f2=fol2+'LOCPOT'
        # generate a folder to run freysoldt calculation
        if bulk:
            new_folder = self.folder+'freysoldt%s_correction_ref-bulk/'% (i)
        else:
            new_folder = self.folder+'freysoldt%s_correction_ref-defect/'% (i)
        if os.path.isdir(new_folder): # remove old calculations
            shutil.rmtree(new_folder)
        os.mkdir(new_folder)
        os.system('cp %s %s' % (f1, new_folder+'neutralrefLOCPOT'))
        os.system('cp %s %s' % (f2, new_folder+'positiveLOCPOT'))
        # step 3
        dic=read_incar(fol2)

        hierdirecs = parent_folder(self.folder, hier=10)  # to find 1st-class materials folder
        self.hierdirecs = hierdirecs
        self.hier1direc = self.hierdirecs[0]    # 1st-class folder Eg. '/home/yubi/work/berylliumO/'
        self.hier2direc = self.hierdirecs[1]    # 2nd-class folder Eg. '/home/yubi/work/berylliumO/wurtzite_2_Vo_HSE/'
        self.sources_direc= self.hier1direc+'sources/' # Eg. '/home/yubi/work/berylliumO/sources/'

        sys.path.append(self.hier1direc)
        #import time
        #time.sleep(1) # pause to see the comment above
        if x_str == 'AEXX':
            # from 'AEXX_EPS.py' import y(eps values)
            from AEXX_EPS import x,y
        elif x_str == 'HFSCREEN':
            # from 'HFSCREEN_EPS.py' import y (=eps)
            from HFSCREEN_EPS import x,y
        else:
            print('Error! x_str is not recognized to find eps values. Should be either AEXX or HFSCREEN')
        print('\nAssume AEXX_EPS or HFSCREEN_EPS has been updated!!! \nEspecially HFSCREEN_EPS should come from the right AEXX file!!!')
        
        x_list=list(x)
        eps_list=list(y)
        this_x=float(dic[x_str])
        this_x=np.round(this_x,6)
        if this_x in x_list:
            this_eps=eps_list[x_list.index(this_x)]
        else:
            paras=np.polyfit(x,y,2) # paras=[a1,a2,a3]
            func = lambda a1,a2,a3,x: a1*(x**2) + a2*(x) + a3
            this_eps = np.round(func(*paras, this_x),6)
            print('Use fitting to get eps of corresponding AEXX')
        print('current %s=%s eps=%s' % (x_str,this_x, this_eps))

        # step 4 get defect position and run freysoldt correction
        defectcenter=read_incar(fol2, incar='DEFECT')['CENTER']
        os.chdir(new_folder)
        command='sxdefectalign --ecut 30 --charge -1 --eps %s --center %s --relative --vdef positiveLOCPOT --vref neutralrefLOCPOT --vasp ' % (this_eps, defectcenter)

        os.system(command+' > sx1.fc')
        with open('sx.sh', 'w') as f:
            f.write(command+'-C  > sx2.fc\ncat sx2.fc')
        os.system('chmod +x sx.sh')
        #os.system('cat sx1.fc')

    def get_DEFECT_info(self, para, var_format=float):
        # read the values of certain variables (came as a list)
        # can read 'EPS', 'FREYCORR', 'EPS_ELEC', 'CENTER'
        # return a list
        assert isinstance(para, list), 'read a list of variables: should keep the variable a list'
        dictionary = read_incar(self.folder, incar='DEFECT')
        values=np.zeros(len(para))
        for j in range(len(para)):
            values[j]=var_format(dictionary[para[j]])
        return values

    def poscar(self, variable):
        print('read CONTCAR')
        from pymatgen import Structure
        mypos = Structure.from_file(self.folder+'CONTCAR')
        if variable == 'abc' or variable == 'a':
            # return lattice constant a
            # abs: assume lattice constant is the same for a,b,c. Return any of it
            return mypos.lattice.a #'a'
        elif variable == 'bond_length':
            return mypos.get_distance(0,1) # the bond length between atom 0 and atom 1
        elif variable == 'volume':
            return np.round(mypos.lattice.volume,3)
        else:
            print('Error! variable %s in POSCAR not recognized' % (variable) )

class rotate_read:
    def __init__(self, folder):
        '''
        self.folder: the working folder
        '''
        #if folder[-1] != '/':
        #    folder=folder+'/' # make sure the folder has correct form: ....../INCAR
        #self.folder=folder
        self.folder = str(Path(folder)) + '/'
        #hierdirecs = parent_folder(self.folder, hier=10)  # to find 1st-class materials folder
        #assert len(hierdirecs) <= 3, 'You are in the son folder. Go back to the parent folder! Otherwise another son folder will be created'
        #self.hierdirecs = hierdirecs
        #self.hier1direc = self.hierdirecs[0]    # 1st-class folder Eg. '/home/yubi/work/berylliumO/'
        #self.hier2direc = self.hierdirecs[1]    # 2nd-class folder Eg. '/home/yubi/work/berylliumO/wurtzite_2_Vo_HSE/'
        #self.sources_direc= self.hier1direc+'sources/' # Eg. '/home/yubi/work/berylliumO/sources/'

    def rotate(self, x_str, y_str, header,var,middle='',avoid=''):
        '''
        x_str = 'AEXX', 'ENCUT', 'HFSCREEN'
        y_str = 'energyHO','energyLU', 'eps', 'equil_energy', 'run_freysoldt', 'freysoldt_corr'
        '''
        existed_folders=find_files(self.folder, header=header, var=var,middle=middle, remove=False,avoid=avoid)
        existed_folders = sorted(existed_folders)
        print('The folders satisfied the pattern: \n',existed_folders)
        # rotate between folders
        #existed_folders=['defect9_AEXX/', 'defect10_AEXX/', 'defect11_AEXX/', 'defect12_AEXX/']# currently only have this data!!!
        xx=np.zeros(len(existed_folders))
        yy=[]
        for i in range(len(existed_folders)):
            # name of each sub-folder
            fol_name= self.folder + existed_folders[i]
            os.chdir(fol_name)
            read_fil=read_file_values(fol_name)
            if x_str == 'AEXX' or x_str == 'ENCUT' or x_str == 'HFSCREEN':
                xx[i]=read_fil.incar([x_str])
            if y_str == 'energyHO':
                # the energies of the highest occupied states
                #returns a list yy[i] = [energyHOLU, nband]
                charge=int(existed_folders[i].split('e_')[0][-1])
                yy.append(read_fil.eigenval_ho_lu(ho_lu=0))
            elif y_str == 'energyLU':
                # the energies of the highest occupied states
                #returns a list yy[i] = [energyHOLU, nband]
                charge=int(existed_folders[i].split('e_')[0][-1])
                yy.append(read_fil.eigenval_ho_lu(ho_lu=1))
            elif y_str == 'eps':
                yy.append(read_fil.outcar_dielec_const())
            #elif y_str == 'koopmansEdiff':
            #    yy.append(read_fil.koopmans_energylevel())
            elif y_str == 'equil_energy': # equilibrium energy
                yy.append(read_fil.oszicar())
            elif y_str == 'run_freysoldt':
                read_fil=read_file_values(self.folder) # create folders in self.folder, not in the subfolder: fol_name 
                ind=int(existed_folders[i].split('koopmans')[1].split('_charge')[0])
                bulk=False
                if bulk:
                    fol1=self.folder+'perfect_supercell/'
                else:
                    assert '1e_' in fol_name, 'The folder name is not right. It should contain 1e_'
                    fol1=fol_name.replace('1e_','0e_')
                #read_fil.locpot2freysoldt_correction(ind, fol1, fol_name,bulk=False)
                read_fil.locpot2freysoldt_correction(ind, fol1, fol_name,bulk=bulk,x_str=x_str)
                yy.append(0.0)
            elif y_str == 'freysoldt_corr':
                yy.append( read_fil.get_DEFECT_info(['FREYCORR'], var_format=float)[0] )
            print('%s=%f %s=%s' % (x_str,xx[i], y_str, yy[i]))
        yy=np.array(yy)
        xx=np.round(xx,6)
        yy=np.round(yy,6)
        os.chdir(self.folder) # go back to the original folder
        return xx,yy