import os
import sys
from math import exp, sqrt
import numpy as np
from math import log, exp
from numpy import arange
from scipy import integrate         # 数値積分関数 integrateを読み込む
from scipy import optimize          # newton関数はscipy.optimizeモジュールに入っている
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt


from tklib.tkfile import tkFile
import tklib.tkre as tkre
from tklib.tkutils import IsDir, IsFile, SplitFilePath
from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg
from tklib.tksci.tksci import Reduce01, Round
from tklib.tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3
from tklib.tkcrystal.tkcif import tkCIF, tkCIFData
from tklib.tkcrystal.tkcrystal import tkCrystal
from tklib.tkcrystal.tkvasp import tkVASP
from tklib.tksci import tkequation
import tklib.tkcsv


"""
Calculate carrier densities and Fermi level based on DOSCAR of VASP.
"""


# constants
pi   = 3.14159265358979323846
h    = 6.6260755e-34    # Js";
hbar = 1.05459e-34      # "Js";
c    = 2.99792458e8     # m/s";
e    = 1.60218e-19      # C";
kB   = 1.380658e-23     # JK<sup>-1</sup>";
me   = 9.1093897e-31    # kg";


#==================================
# global variables
#==================================
Debug = 0

#mode: 'EF', 'T'
mode = 'T'

# Sample information (DOS.data is converted to states/eV/Vcell)
# read from DOSCAR
CAR_path = '.'
DOSCAR_path = "TotalDOS-SnSe.dat"
#DOSCAR_path = "DOS.dat"

UseEF0 = 0
T0   = 300       # K
EF0  = 0.2    # EF at charge neutral
Emid = EF0    # E in bandgap, used to search EC and EV
Egth = 0.01   # Critera DOS value to search band edges

# Other levels
dEA =  0.3    # Acceptor level measured from EV, eV
NA  =  0.0    # 1.0e14
dED = -0.28   # Donor level measured from EC, eV
ND  =  0.0    # 0.5e17
EA  = None    # Acceptor states are not considered if EA = None
ED  = None    # Donor states are not considered if ED = None


# EF range for mode = 'EF'
dEFmin = -0.2  # measured from EV, eV
dEFmax =  0.2  # measured from EC, eV
nEF    = 50

# Temperature range for mode = 'T'
Tmin  = 300  # K
Tmax  = 600
nT = 11
Tstep = None

# integration parameters
# Integration range w.r.t. kBT
Einteg0 = -3.0
Einteg1 =  3.0
nrange  =  6.0

# Relative accuracy of the quad functin
eps = 1.0e-5
nmaxdiv = 150

# Bisection method parameters
# Initial search range
dEbisec  = 0.4
eps_bisec = 1.0e-6
nmaxiter_bisection = 200

# Newton method parmeters
dump_newton = 0.5
eps_newton  = 1.0e-6
h_newton    = 1.0e-6
nmaxiter_newton   = 30

iprintinterval = 1

# Global variables
Ne0K = None

Vcell   = None  # A^3
E_raw   = None
dos_raw = None
nDOS    = None

Emin    = None
Emax    = None
Estep   = None

EV      = None
EC      = None

EFmin = None
EFmax = None
EFstep = None


#=============================
# Treat argments
#=============================
def pfloat(str):
    try:
        return float(str)
    except:
        return None

def pint(str):
    try:
        return int(str)
    except:
        return None

def getarg(position, defval = None):
    try:
        return sys.argv[position]
    except:
        return defval

def getfloatarg(position, defval = None):
    return pfloat(getarg(position, defval))

def getintarg(position, defval = None):
    return pint(getarg(position, defval))

def usage():
    global mode
    global CAR_path, DOSCAR_path, Vcell
    global E_raw, dos_raw, nDOS, Emin, Emax, Estep
    global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
    global T0, EF0

    if CAR_path == '':
        CAR_path = '.'
    
    print("")
    print("Usage: Variables in () are optional")
    print(" python {} mode (file args)".format(sys.argv[0], ))
    print(" (i) python {} T CAR_path DOSCAR_path Tmin Tmax nT)".format(sys.argv[0]))
    print("     ex: python {} T {} {} {} {} {}".format(sys.argv[0], CAR_path, DOSCAR_path, Tmin, Tmax, nT))
    print(" (ii) python {} EF CAR_path DOSCAR_path nEF EF0 T0)".format(sys.argv[0]))
    print("     ex: python {} EF {} {} {} {} {}".format(sys.argv[0], CAR_path, DOSCAR_path, nEF, EF0, T0))

def updatevars():
    global mode
    global CAR_path, DOSCAR_path, Vcell
    global E_raw, dos_raw, nDOS, Emin, Emax, Estep
    global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
    global T0, EF0
    global Tmin, Tmax, nT, Tstep

    argv = sys.argv
    if len(argv) == 1:
        usage()
        exit()

    mode = getarg( 1, mode)
    CAR_path    = getarg( 2, CAR_path)
    DOSCAR_path = getarg( 3, DOSCAR_path)
    if mode == 'EF':
        nEF = getintarg(4, nEF)
    if mode == 'T':
        Tmin = getintarg  (4, Tmin)
        Tmax = getintarg  (5, Tmax)
        nT   = getintarg  (6, nT)
        EF0  = getfloatarg(7, EF0)
        T0   = getfloatarg(8, T0)
        Tstep = (Tmax - Tmin) / (nT - 1)


#=============================
# other functions
#=============================
def savecsv(outfile, header, datalist):
    try: 
        print("Write to [{}]".format(outfile))
        f = open(outfile, 'w')
    except:
#    except IOError:
        print("Error: Can not write to [{}]".format(outfile))
    else:
        fout = csv.writer(f, lineterminator='\n')
        fout.writerow(header)
#        fout.writerows(data)
        for i in range(0, len(datalist[0])):
            a = []
            for j in range(len(datalist)):
                a.append(datalist[j][i])
            fout.writerow(a)
        f.close()

def read_csv(infile, xmin = None, xmax = None, delimiter = ','):
    print("xrange=", xmin, xmax)
    data = []
    try:
        infp = open(infile, "r")
        f = csv.reader(infp, delimiter = delimiter)
        header = next(f)
        print("header=", header)
        for j in range(len(header)):
            data.append([])

        for row in f:
            x = pfloat(row[0])
            if xmin is not None and xmin <= x <= xmax:
                y = pfloat(row[1])
                data[0].append(x)
                data[1].append(y)
    except:
        print("Error: Can not read [{}]".format(infile))
        exit()
    return header, data[0], data[1]


def FindBandEdges(E, DOS, Emidgap, DOSth):
    EV = 1.0e30
    EC = 1.0e30
    i0 = -1
    for i in range(len(E)):
        if E[i] >= Emidgap:
            i0 = i
            break

    for i in range(i0, 0, -1):
        if DOS[i] > DOSth:
            EV = E[i+1]
            break

    for i in range(i0, len(E), +1):
        if DOS[i] > DOSth:
            EC = E[i-1]
            break

    return EV, EC


def rieman(x0, dx, y, xmin, xmax):
    S = 0.0
    for i in range(len(y)):
        x = x0 + i * dx
        if x < xmin or xmax < x:
            continue
        S += y[i]
    return S * dx

# define the DOS function
def DOS(E):
    global E_raw, dos_raw
    global Emin, Emax, Estep, nDOS

    if E < Emin or Emax < E:
        print("Error in f_dos: Given E={} exceeds the DOS E range [{}, {}]".format(E, Emin, Emax))
        exit()

    fdos = interp1d(E_raw, dos_raw, kind = 'cubic')
    return fdos(E)

# Fermi-Dirac distribution function
def fe(E, T, EF):
    global e, kB

    if T == 0.0:
        if E < EF:
            return 1.0
        else:
            return 0.0
    return 1.0 / (exp((E - EF) * e / kB / T) + 1.0)

def fh(E, T, EF):
    return 1.0 - fe(E, T, EF)

# Define the function to be integrated
def DOSfe(E, T, EF):
    return DOS(E) * fe(E, T, EF)

def DOSfh(E, T, EF):
    return DOS(E) * fh(E, T, EF)

def integrate(func, E0, E1, h):
    n = int((E1 - E0) / h + 1.000001)
    h = (E1 - E0) / n
    y = [func(E0 + i * h) for i in range(n+1)]
    
    S = 0.5 * (y[0] + y[n]) + sum(y[1:n])
    return [h * S, -1.0]

# Calculate the number of electrons from E0 to E1 with the Fermi level EF at T
def Ne(T, EF, E0, E1):
    global Estep
#    N = integrate.quad(lambda E: DOSfe(E, T, EF), E0, E1, limit = nmaxdiv, epsrel = eps)
    N = integrate(lambda E: DOSfe(E, T, EF), E0, E1, Estep)
    if Debug:
        print("  Ne integ error: {}".format(N[1]))
    return N[0]

def Nh(T, EF, E0, E1):
    global Estep
#    N = integrate.quad(lambda E: DOSfh(E, T, EF), E0, E1, limit = nmaxdiv, epsrel = eps)
    N = integrate(lambda E: DOSfh(E, T, EF), E0, E1, Estep)
    if Debug:
        print("  Nh integ error: {}".format(N[1]))
    return N[0]

# ionized donor density
def NDp(EF, T):
    global ND, ED, kB
    return ND * (1.0 - fe(ED, T, EF))

# ionized acceptor density
def NAm(EF, T):
    global NA, EA, kB
    n = NA * fe(EA, T, EF)
#    print("n=", EF, T, NA, EA, n)
    return n
    
# Calculate EF base on 
# the charge neutrality (electron number) condition, Ne(T, EF, dE) - N = 0,
# using the Newton method with the initial value EF0
def CalEF(T, EF0, E0, E1, totNe):
    EF = optimize.newton(lambda EF: Ne(T, EF, E0, E1) - totNe, EF0)
    return EF

def cal_densities(T, EF, ECmin, ECmax, EVmin, EVmax):
    Ne1 = Ne(T, EF, ECmin, ECmax)
    Nh1 = Nh(T, EF, EVmin, EVmax)
    NDp1 = NDp(EF, T)
    NAm1 = NAm(EF, T)
    
    return Ne1, Nh1, NDp1, NAm1, Ne1 + NAm1 - Nh1 - NDp1

def dQ(T, EF, ECmin, ECmax, EVmin, EVmax, N0):
    Ne1, Nh1, NDp1, NAm1, dN = cal_densities(T, EF, ECmin, ECmax, EVmin, EVmax)
    return dN - N0

def diff(h, T, EF, ECmin, ECmax, EVmin, EVmax, N0):
    return (dQ(T, EF + h, ECmin, ECmax, EVmin, EVmax, N0) - dQ(T, EF - h, ECmin, ECmax, EVmin, EVmax, N0)) / 2.0 / h
    
def callbackfunc(obj):
    if obj.xa is None:
        print("Iter {:5d}: x: {:>16.12f}, dx = {:>10.4g}"
            .format(obj.iter, obj.x, obj.dx))
    else:
        print("Iter {:5d}: x: {:>16.12f} in [{:>16.12f}, {:>16.12f}], dx = {:>10.4g}"
            .format(obj.iter, obj.x, obj.xa, obj.xb, obj.dx))
    return 1


def exec_T():
    global mode, T0, EF0, Ne0K
    global file, E_raw, dos_raw
    global nDOS, Emin, Emax, Estep
    global EV, EC
    global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
    global dEA, NA, dED, ND, EA, ND

    xT = []
    xTinv = []
    yEF = []
    yEint = []
    yEentropy = []
    yEFa = []
    yNe = []
    yNh = []
    yNDp = []
    yNAm = []
    yRH  = []
    EF = EF0
    print ("%5s\t%8s\t%14s" % ("T(K)", "EF(eV)", "Ne(0 K) (check)"))
    for iT in range(nT):
        T = Tmin + iT * Tstep

        print("")
        print("iT {:03d}: {} K".format(iT, T))

        kBTe = kB * T / e
        dE = nrange * kBTe

        if EV - dE < Emin or Emax < EC + dE:
            print("*** Integration range error: [{}, {}] exceeds DOS E range [{}, {}]".format(E0, E1, Emin, Emax))
            exit()

# Calculate excess electron density N0 if UseEF0 != 0
        Ne0, Nh0, NDp0, NAm0, dQ0 = cal_densities(T0, EF0, EC - Estep, EC + dE, EV - dE, EV + Estep)
        if UseEF0:
            N0 = Ne0 - Nh0
            print("Ne0 for EF={:8.4f} eV at {:8.4f} K = {:12.6g} cm^-3 (={:12.4g} - {:12.4g})".format(EF0, T0, N0, Ne0, Nh0))
        else:
            N0 = 0

        if iT == 0:
#        if iT == 0 or iT > 0:
            EFmin = EF - dEbisec
            EFmax = EF + dEbisec
            solver = tkequation.Equation(
                debug_explicit = 1,
                method    = 'bisection',
#               method    = 'newton', 
                func      = lambda EF: dQ(T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep, N0),
                xa        = EFmin,
                xb        = EFmax,
                nmaxiter  = nmaxiter_bisection,
                eps       = eps_bisec,
                callback  = callbackfunc,
                isprint   = 0
                )
        else:
                   solver = tkequation.Equation(
                debug_explicit = 1,
                method    = 'newton', 
                func      = lambda EF: dQ(T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep, N0),
                diff1func = lambda EF: diff(h_newton, T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep, N0), 
                x0        = EF,
                dump      = dump_newton,
                nmaxiter  = nmaxiter_newton,
                eps       = eps_newton,
                callback  = callbackfunc,
                isprint   = 0
                )

        x = solver.solve()
        if solver.iter >= 0:
            print("  Success: Convergence reached at iter = {}, x = {:10.6g}, f = {:8.4g}"
                .format(solver.iter, solver.x, solver.f))
        else:
            print("  Failed: Convergence did not reach")
            return 0

        EF  = solver.x
        dQh = solver.f
        Neh, Nhh, NDph, NAmh, dQh = cal_densities(T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep)
        RH = 1.0 / e * (Nhh - Neh) / pow(Nhh + Neh, 2.0)

        xT.append(T)
        xTinv.append(1000.0 / T)
        yEF.append(EF)
        yNe.append(Neh)
        yNh.append(Nhh)
        yNDp.append(NDph)
        yNAm.append(NAmh)
        yRH.append(RH)
        print("  ***(T,EF,Ne,Nh,NDp,NAm,RH)=({:6.4f},{:6.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g})"
                    .format(T, EF, Neh, Nhh, NDph, NAmh, RH))
        print ("%5.0f\t%12.6g" % (T, EF))


    print("")
    yEaNe   = []
    yEaNeTh = []
    yEaNh   = []
    yEaNhTh = []
    print ("%8s %10s %8s %8s %8s %8s %10s %10s %12s %12s %12s" 
            % ("T(K)", "EF(eV)", "Ea(Ne,meV)", "Ea(Ne,T1/2)", "Ea(Nh,meV)", "Ea(Nh,T1/2)", "Ne", "Nh", "ND+", "NA-", "RH"))
    for i in range(nT):
        if i == 0:
            i1 = i
            i2 = i+1
        elif i == nT - 1:
            i1 = i-1
            i2 = i
        else:
            i1 = i-1
            i2 = i+1

        Th1 = sqrt(xT[i1])
        Th2 = sqrt(xT[i2])
        dT = 1.0 / (xTinv[i2] - xTinv[i1])

        EaNe   = (log(yNe[i2]+1.0e-10)     - log(yNe[i1]+1.0e-10))     * dT # meV
        EaNeTh = (log(yNe[i2]/Th2+1.0e-10) - log(yNe[i1]/Th1+1.0e-10)) * dT # meV
        EaNh   = (log(yNh[i2]+1.0e-10)     - log(yNh[i1]+1.0e-10))     * dT # meV
        EaNhTh = (log(yNh[i2]/Th2+1.0e-10) - log(yNh[i1]/Th1+1.0e-10)) * dT # meV
        yEaNe.append  (-1000.0 * 1000.0 * kB / e * EaNe)
        yEaNeTh.append(-1000.0 * 1000.0 * kB / e * EaNeTh)
        yEaNh.append  (-1000.0 * 1000.0 * kB / e * EaNh)
        yEaNhTh.append(-1000.0 * 1000.0 * kB / e * EaNhTh)
        print("%8.3f %10.6f %10.6g %10.6g %10.6g %10.6g " % (xT[i], yEF[i], yEaNe[i], yEaNeTh[i], yEaNh[i], yEaNhTh[i])
            + "%12g %12g %12g %12g %12g" % (yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i]))


#=============================
# Plot graphs
#=============================
    fig = plt.figure(figsize = (12, 8))

    ax1 = fig.add_subplot(2, 3, 1)
    ax2 = fig.add_subplot(2, 3, 2)
    ax3 = fig.add_subplot(2, 3, 4)
    ax4 = fig.add_subplot(2, 3, 3)
    ax5 = fig.add_subplot(2, 3, 6)
    ax1.plot(xT, yEF, label = '$E_F$', color = 'black')
    ax1.set_xlabel("T (K)")
    ax1.set_ylabel("E$_F$ (eV)")
    ax1.legend()
    ax2.plot(E_raw, dos_raw, label = 'DOS', color = 'black')
    ax2.plot([EF0, EF0], [min(dos_raw), max(dos_raw)], color = 'red', linestyle = 'dashed')
    ax2.set_xlabel("$E$ (eV)")
    ax2.set_ylabel("DOS (states/cm$^3$)")
    ax2.legend()
    ax3.plot(xT, yRH, label = '$R_H$', color = 'black', marker = 'o')
    ax3.set_xlabel("$T$ (K)")
    ax3.set_ylabel("$R_H$ (C$^{-1}$cm$^{-3}$)")
    ax3.legend()
    ax4.plot(xTinv, yNe,  label = '$N_e$',   color = 'black', marker = 'o')
    ax4.plot(xTinv, yNh,  label = '$N_h$',   color = 'red',   marker = 'o')
    ax4.plot(xTinv, yNDp, label = '$N_D^+$', color = 'blue',  marker = 'x')
    ax4.plot(xTinv, yNAm, label = '$N_A^-$', color = 'green', marker = 'x')
    ax4.set_yscale('log')
#    ax4.plot([0, 0], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
    ax4.set_xlabel("1000/$T$ (K$^{-1}$)")
    ax4.set_ylabel("$N$ (cm$^{-3}$)")
    ax4.legend()
    ax5.plot(xTinv, yEaNe,   label = '$E_a$($N_e$)',          color = 'black', marker = 'o')
    ax5.plot(xTinv, yEaNeTh, label = '$E_a$($N_eT^{-0.5}$)', color = 'black', marker = 'o')
    ax5.plot(xTinv, yEaNh,   label = '$E_a$($N_h)',           color = 'red',   marker = 'o')
    ax5.plot(xTinv, yEaNhTh, label = '$E_a$($N_hT^{-0.5})',  color = 'red',   marker = 'o')
    ax5.set_xlabel("1000/$T$ (K$^{-1}$)")
    ax5.set_ylabel("$E_a$ (meV)")
    ax5.legend()
# Rearange the graph axes so that they are not overlapped
    plt.tight_layout()

    plt.pause(0.1)

    print("")
    usage()
    print("")

    print("Press ENTER to exit>>", end = '')
    input()

def exec_EF():
    global mode, T0, EF0
    global file, E_raw, dos_raw
    global nDOS, Emin, Emax, Estep
    global EV, EC
    global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
    global dEA, NA, dED, ND, EA, ND

    kBTe = kB * T0 / e
    dE = nrange * kBTe

    xEF = []
    yNe = []
    yNh = []
    yNDp = []
    yNAm = []
    yRH  = []
    EFprev = EF0
    print("at {} K".format(T0))
    print ("%5s\t%8s\t%14s\t%14s\t%14s\t%14s" % ("EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH"))
    for iEF in range(nEF):
        EF = EFmin + iEF * EFstep
        dE = nrange * kBTe

        Ne1, Nh1, NDp1, NAm1, dQ = cal_densities(T0, EF, EC - Estep, EC + dE, EV - dE, EV + Estep)
        RH = 1.0 / e * (Nh1 - Ne1) / pow(Nh1 + Ne1, 2.0)

        xEF.append(EF)
        yNe.append(Ne1)
        yNh.append(Nh1)
        yNDp.append(NDp1)
        yNAm.append(NAm1)
        yRH.append(RH)
        print ("%10.4f\t%12.6g\t%12.6g\t%12.6g\t%12.6g\t%12.6g" % (EF, Ne1, Nh1, NDp1, NAm1, RH))

    yT0Ne = []
    yT0Nh = []
    print("")
    print ("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%8s\t%8s" % ("EF(eV)", "Ne", "Nh", "ND+", "NA-", "RH", "T0(Ne)", "T0(Nh)"))
    for i in range(nEF):
        if i == 0:
            i1 = i
            i2 = i+1
        elif i == nEF - 1:
            i1 = i-1
            i2 = i
        else:
            i1 = i-1
            i2 = i+1

        dEF = xEF[i2] - xEF[i1]
        dlnNedEF = (log(yNe[i2]) - log(yNe[i1])) / dEF
        dlnNhdEF = (log(yNh[i2]) - log(yNh[i1])) / dEF
        T0Ne =  1.0 / kB / dlnNedEF * e
        T0Nh = -1.0 / kB / dlnNhdEF * e
        yT0Ne.append(T0Ne)
        yT0Nh.append(T0Nh)
        print ("%10.4f\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%10.4g\t%8.4g\t%8.4g" % (xEF[i], yNe[i], yNh[i], yNDp[i], yNAm[i], yRH[i], T0Ne, T0Nh))

    print("")
    Emid = (EV + EC) / 2.0
    print("Effective density of states at the mid gap {} eV:".format(Emid))
    fNe = interp1d(xEF, yNe, kind = 'cubic')
    fNh = interp1d(xEF, yNh, kind = 'cubic')
    dlnNedEF = (log(fNe(Emid + Estep)) - log(fNe(Emid - Estep))) / 2.0 / Estep
    dlnNhdEF = (log(fNh(Emid + Estep)) - log(fNh(Emid - Estep))) / 2.0 / Estep
    T0Ne =  1.0 / kB / dlnNedEF * e
    T0Nh = -1.0 / kB / dlnNhdEF * e
    
    NC = fNe(Emid) * exp((EC - Emid) * e / kB / T0Ne)
    NV = fNh(Emid) * exp((Emid - EV) * e / kB / T0Nh)
    print("   NC = {:12.6g} cm-3  (T0 = {:12.6g} K)".format(NC, T0Nh))
    print("   NV = {:12.6g} cm-3  (T0 = {:12.6g} K)".format(NV, T0Nh))
    
    xEFapprox = [EFmin + i * EFstep for i in range(nEF)]
    yNeapprox = [NC * exp(-(EC - xEFapprox[i]) * e / kB / T0) for i in range(nEF)]
    yNhapprox = [NV * exp(-(xEFapprox[i] - EV) * e / kB / T0) for i in range(nEF)]


#=============================
# Plot graphs
#=============================
    fig = plt.figure(figsize = (8, 8))

    ax2 = fig.add_subplot(2, 2, 1)
    ax3 = fig.add_subplot(2, 2, 2)
    ax4 = fig.add_subplot(2, 2, 3)
    ax1 = fig.add_subplot(2, 2, 4)
    ax2.plot(E_raw, dos_raw, label = 'DOS', color = 'black')
    ax2.plot([EF0, EF0], [min(dos_raw), max(dos_raw)], color = 'red', linestyle = 'dashed')
    ax2.set_xlabel("$E$ (eV)")
    ax2.set_ylabel("DOS (states/cm$^3$)")
#    ax2.legend()
    ax3.plot(xEF, yNe,  label = '$N_e$',   color = 'black')   #, marker = 'o')
    ax3.plot(xEF, yNh,  label = '$N_h$',   color = 'red')     #, marker = 'x')
    ax3.plot(xEF, yNDp, label = '$N_D^+$', color = 'blue')    #, marker = 'x')
    ax3.plot(xEF, yNAm, label = '$N_A^-$', color = 'green')   #, marker = 'x')
    ax3.plot(xEFapprox, yNeapprox, label = '$N_e$(Boltz)', color = 'gray',   linewidth = 0.5) #, marker = 'x')
    ax3.plot(xEFapprox, yNhapprox, label = '$N_h$(Boltz)', color = 'purple', linewidth = 0.5) #, marker = 'x')
    ax3.plot([EV, EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
    ax3.plot([EC, EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
    ax3.set_yscale('log')
    ax3.set_xlabel("$E_F$ (eV)")
    ax3.set_ylabel("N (cm$^{-3}$)")
    ax3.legend()
    ax4.plot(xEF, yRH, label = 'R$_H$', color = 'black', marker = 'o', markersize = 1.5)
    ax4.set_xlabel("$E_F$ (eV)")
    ax4.set_ylabel("$R_H$ (C$^{-1}$cm$^{-3}$)")
#    ax4.legend()
    ax1.plot(xEF, yT0Ne, label = '$T_0$($N_e$)', color = 'black')   #, marker = 'o')
    ax1.plot(xEF, yT0Nh, label = '$T_0$($N_h$)', color = 'red')     #, marker = 'x')
    ax1.plot([EV, EV], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
    ax1.plot([EC, EC], ax3.get_ylim(), color = 'red', linestyle = 'dashed', linewidth = 0.5)
    ax1.set_ylim([0.0, 2.0*T0])
    ax1.set_xlabel("$E_F$ (eV)")
    ax1.set_ylabel("$T_0$ (K)")
    ax1.legend()
# Rearange the graph axes so that they are not overlapped
    plt.tight_layout()

    plt.pause(0.1)

    print("")
    usage()
    print("")

    print("Press ENTER to exit>>", end = '')
    input()


def main():
    global mode, T0, EF0, Ne0K
    global CAR_path, DOSCAR_path, Vcell
    global E_raw, dos_raw, nDOS, Emin, Emax, Estep
    global EV, EC
    global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep
    global dEA, NA, dED, ND, EA, ED
    global Einteg0

    updatevars()

    print("")
    print("=======================================================")
    print(" Calculate EF/T dependence of semiconductor properties")
    print("=======================================================")
    print("mode: ", mode)

    print("*** Read VASP files:")
    vasp = tkVASP()
    CAR_path = vasp.getdir(CAR_path)
    INCAR_path   = vasp.get_INCAR(CAR_path)
    POSCAR_path  = vasp.get_POSCAR(CAR_path)
    CONTCAR_path = vasp.get_CONTCAR(CAR_path)
    OUTCAR_path  = vasp.get_OUTCAR(CAR_path)
    print("  CAR dir: ", CAR_path)
    print("  INCAR  : ", INCAR_path)
    print("  POSCAR : ", POSCAR_path)
    print("  CONTCAR: ", CONTCAR_path)
    print("  OUTCAR : ", OUTCAR_path)

    print("")
    print("*** Read crystal structure from [{}]".format(POSCAR_path))
    cry = vasp.read_poscar(POSCAR_path)
#    cry.PrintInf("cell")
    a, b, c, alpha, beta, gamm = cry.LatticeParameters()
    Vcell = cry.Volume()
    print("  cell: {:12.8f} {:12.8f} {:12.8f} A   {:10.6f} {:10.6f} {:10.6f}".format(a, b, c, alpha, beta, gamm))
    print("  volume: {:12.6f} A^-3".format(Vcell))

    print("")
    print("*** Read Total DOS from [{}]".format(DOSCAR_path))
    E_raw, dos_raw = np.genfromtxt(DOSCAR_path, skip_header=1, usecols=(0,1), unpack=True)
    nDOS = len(E_raw)
    for i in range(len(E_raw)):
        dos_raw[i] *= 1.0 / (Vcell * 1.0e-24)
    Emin = E_raw[0]
    Emax = E_raw[nDOS-1]
    Estep = E_raw[1] - E_raw[0]
    print("  E range: {} - {}, {} eV step".format(Emin, Emax, Estep))

    print("*** Search EV and EC from the midgap energy {} eV with threshold DOS value of {}".format(Emid, Egth))  
    EV, EC = FindBandEdges(E_raw, dos_raw, Emid, Egth)
    Eg = EC - EV
    print("  Band edges: EV={:10.6g}  EC={:10.6g}  Eg={:10.6g} eV".format(EV, EC, Eg))

    EFmin = EV + dEFmin
    EFmax = EC + dEFmax
    EFstep = (EFmax - EFmin) / (nEF - 1)

    print("")
    print("*** Additional levels")
    ED = EC + dED
    EA = EV + dEA
    print("  Donor level   : {} eV, {} cm^-3".format(ED, ND))
    print("  Acceptor level: {} eV, {} cm^-3".format(EA, NA))

    print("")
    print("EF dependence")
    print("  EF range: {} - {}, {} eV step".format(EFmin, EFmax, EFstep))

    print("")
    print("T dependence")
    print("  Temperature range: {} - {}, {} K step".format(Tmin, Tmax, Tstep))

    print("")
    print("Integration configuration")
    print("  Integration E range: {} - {} eV, or EF+-{}*kBT eV".format(Einteg0, Einteg1, nrange))
    print("  quad() parameters")
    print("    releps=", eps)
    print("    nmaxdiv=", nmaxdiv)

    print("")
    print("EF at {} K = {} eV".format(T0, EF0))
    Einteg0 = -5.0
    Ne0K = Ne(0.0, EF0, Einteg0, EF0)
    print("  Ne at {} K from {} eV = {:12.6g} cm^-3.".format(T0, Einteg0, Ne0K))
    print("")

    if mode == 'T':
        exec_T()
    elif mode == 'EF':
        exec_EF()
    else:
        terminate("Error: Invalid mode [{}]".format(mode), usage = usage)


if __name__ == '__main__':
    main()