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

Download script from .\EF-T-DOS.py

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


"""
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-1";
me = 9.1093897e-31 # kg";


Debug = 0


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

# EF range for mode = 'EF'
T0 = 300 # K
dEFmin = -0.2 # measured from EV, eV
dEFmax = 0.2 # measured from EC, eV
nEF = 50
EFmin = None
EFmax = None
EFstep = None

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

# Sample information (DOS.data is converted to states/eV/Vcell)
Vcell = 212.309 # A^3

# read from DOSCAR
file = "DOS.dat"
E_raw = None
dos_raw = None
nDOS = None
Emin = None
Emax = None
Estep = None
EV = None
EC = None

dEA = 0.32 # Acceptor level measured from EV, eV
NA = 0.8e17
dED = -0.28 # Donor level measured from EC, eV
ND = 0.5e17
EA = None # Acceptor states are not considered if EA = None
ED = None # Donor states are not considered if ED = None

# other parameters
EF0 = 0.4 # EF at charge neutral
Egth = 0.01 # Critera DOS value to search band edges

# 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
epsbisec = 1.0e-7
# max iteration number
nmaxiter = 200
# Output cycle
iprintiterval = 1

# graph configuration
fontsize = 16
legend_fontsize = 12


#=============================
# 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 Vcell, T0, EF0
    global file, E_raw, dos_raw, nDOS, Emin, Emax, Estep
    global dEFmin, dEFmax, EFmin, EFmax, nEF, EFstep

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

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

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

    mode = getarg( 1, mode)
    file = getarg( 2, file)
    if mode == 'EF':
        nEF = getintarg(3, nEF)
    if mode == 'T':
        Tmin = getintarg(3, Tmin)
        Tmax = getintarg(4, Tmax)
        nT = getintarg(5, nT)
        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, EF0, Egth):
    EV = 1.0e30
    EC = 1.0e30
    i0 = -1
    for i in range(len(E)):
        if E[i] >= EF0:
            i0 = i
            break

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

    for i in range(i0, len(E), +1):
        if DOS[i] > Egth:
            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 - 1)
    y = [func(E0 + i * h) for i in range(n)]

    S = 0.5 * (y[0] + y[n-1]) + sum(y[1:n-1])
    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 exec_T():
    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

    xT = []
    xTinv = []
    yEF = []
    yEint = []
    yEentropy = []
    yEFa = []
    yNe = []
    yNh = []
    yNDp = []
    yNAm = []
    yRH = []
    EFprev = EF0
#初期範囲
    EFmin = EF0 - dEbisec
    EFmax = EF0 + dEbisec
    dEFbisec = 0.1
    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

        E0 = EV - dE
        E1 = EC + dE

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

        Ne0K = Ne(0.0, EF0, E0, EF0)

# まず、EFmin,EFmaxにおけるΔQを計算し、それぞれが正・負あるいは負・生となっていることを確認する
        print("E range:", EF0, E1)
#        Nemin = Ne(T, EFmin, EF0, E1)
        Nemin = Ne(T, EFmin, EC - Estep, E1)
#        Nhmin = Nh(T, EFmin, E0, EF0)
        Nhmin = Nh(T, EFmin, E0, EV + Estep)
        NAmmin = NAm(EFmin, T)
        NDpmin = NDp(EFmin, T)
        dQmin = Nemin + NAmmin - Nhmin - NDpmin
#        Nemax = Ne(T, EFmax, EF0, E1)
        Nemax = Ne(T, EFmax, EC - Estep, E1)
#        Nhmax = Nh(T, EFmax, E0, EF0)
        Nhmax = Nh(T, EFmax, E0, EV + Estep)
        NAmmax = NAm(EFmax, T)
        NDpmax = NDp(EFmax, T)
        dQmax = Nemax + NAmmax - Nhmax - NDpmax
        print(" min (EF,Ne,Nh,NDp,NAm,dQ)=({:6.4f},{:10.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g})"
                    .format(iT, EFmin, Nemin, Nhmin, NDpmin, NAmmin, dQmin))
        print(" max (EF,Ne,Nh,NDp,NAm,dQ)=({:6.4f},{:10.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g})"
                    .format(EFmax, Nemax, Nhmax, NDpmax, NAmmax, dQmax))
#        dQmin = Ne(EFmin, T) + NAm(EFmin, T) - Nh(EFmin, T) - NDp(EFmin, T)
#        dQmax = Ne(EFmax, T) + NAm(EFmax, T) - Nh(EFmax, T) - NDp(EFmax, T)
#        print(" EFmin = {:12.8f} dQmin = {:12.4g}".format(EFmin, dQmin))
#        print(" EFmax = {:12.8f} dQmax = {:12.4g}".format(EFmax, dQmax))
        print(" dQ=", dQmin, dQmax)
        print(" dQmin*dQmax=", dQmin * dQmax)
        if dQmin * dQmax > 0.0:
            print("Error: Initial Emin and Emax should be chosen as dQmin * dQmax < 0")
            return 0

# 2分法開始
        for i in range(nmaxiter):
            EFhalf = (EFmin + EFmax) / 2.0
#            Neh = Ne(T, EFhalf, EF0, E1)
            Neh = Ne(T, EFhalf, EC - Estep, E1)
#            Nhh = Nh(T, EFhalf, E0, EF0)
            Nhh = Nh(T, EFhalf, E0, EV + Estep)
            NAmh = NAm(EFhalf, T)
            NDph = NDp(EFhalf, T)
            dQhalf = Neh + NAmh - Nhh - NDph
            print(" iter {:03d}: half (EF,Ne,Nh,NDp,NAm,dQ)=({:6.4f},{:10.4g},{:10.4g},{:10.4g},{:10.4g},{:10.4g})"
                    .format(i, EFhalf, Neh, Nhh, NDph, NAmh, dQhalf))
# EFの精度がepsより小さくなったら計算終了
            if abs(EFmin - EFhalf) < epsbisec and abs(EFmax - EFhalf) < epsbisec:
#                    print(" Success: Convergence reached at EF = {}".format(EFhalf))
                    break
            if dQmin * dQhalf < 0.0:
                EFmax = EFhalf
                dQmax = dQhalf
            else:
                EFmin = EFhalf
                dQmin = dQhalf
        else:
            print(" Failed: Convergence did not reach")
            return 0

        EF = EFhalf

        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))

        EFprev = EF
        EFmin = EF - dEFbisec
        EFmax = EF + dEFbisec
        print ("%5.0f\t%12.6g" % (T, EF))

    print("")
    yEaNe = []
    yEaNeTh = []
    yEaNh = []
    yEaNhTh = []
    print ("%8s %8s %8s %8s %8s %8s %12s %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]) - log(yNe[i1])) * dT # meV
        EaNeTh = (log(yNe[i2]/Th2) - log(yNe[i1]/Th1)) * dT # meV
        EaNh = (log(yNh[i2]) - log(yNh[i1])) * dT # meV
        EaNhTh = (log(yNh[i2]/Th2) - log(yNh[i1]/Th1)) * 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("%5.3f %8.3f %8.3g %8.3g %8.3g %8.3g " % (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 = 'EF', color = 'black')
    ax1.set_xlabel("T (K)", fontsize = fontsize)
    ax1.set_ylabel("E$_F$ (eV)", fontsize = fontsize)
    ax1.legend(fontsize = legend_fontsize)
    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)", fontsize = fontsize)
    ax2.set_ylabel("DOS (states/cm$^3$)", fontsize = fontsize)
    ax2.legend(fontsize = legend_fontsize)
    ax3.plot(xT, yRH, label = 'RH', color = 'black', marker = 'o')
    ax3.set_xlabel("T (K)", fontsize = fontsize)
    ax3.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)", fontsize = fontsize)
    ax3.legend(fontsize = legend_fontsize)
    ax4.plot(xTinv, yNe, label = 'Ne', color = 'black', marker = 'o')
    ax4.plot(xTinv, yNh, label = 'Nh', color = 'red', marker = 'o')
    ax4.plot(xTinv, yNDp, label = 'ND$^+$', color = 'blue', marker = 'x')
    ax4.plot(xTinv, yNAm, label = 'NA$^-$', 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}$)", fontsize = fontsize)
    ax4.set_ylabel("N (cm$^{-3}$)", fontsize = fontsize)
    ax4.legend(fontsize = legend_fontsize)
    ax5.plot(xTinv, yEaNe, label = 'Ea(Ne)', color = 'black', marker = 'o')
    ax5.plot(xTinv, yEaNeTh, label = 'Ea(Ne*T$^{-0.5}$)', color = 'black', marker = 'o')
    ax5.plot(xTinv, yEaNh, label = 'Ea(Nh)', color = 'red', marker = 'o')
    ax5.plot(xTinv, yEaNhTh, label = 'Ea(Nh*T$^{-0.5})', color = 'red', marker = 'o')
    ax5.set_xlabel("1000/T (K$^{-1}$)", fontsize = fontsize)
    ax5.set_ylabel("Ea (meV)", fontsize = fontsize)
    ax5.legend(fontsize = legend_fontsize)
    ax1.tick_params(labelsize = fontsize)
    ax2.tick_params(labelsize = fontsize)
    ax3.tick_params(labelsize = fontsize)
    ax4.tick_params(labelsize = fontsize)
    ax5.tick_params(labelsize = fontsize)
# Rearange the graph axes so that they are not overlapped
    plt.tight_layout()

#    plt.pause(0.1)
    print("Close the graph window to terminate this program")
    plt.show()

    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
    E0 = EV - dE
    E1 = EC + dE

    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 = Ne(T0, EF, EF0, E1)
        Ne1 = Ne(T0, EF, EC - Estep, E1)
#        Nh1 = Nh(T0, EF, E0, EF0)
        Nh1 = Nh(T0, EF, E0, EV + Estep)
        NDp1 = NDp(EF, T0)
        NAm1 = NAm(EF, T0)
        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 ("%5s\t%8s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s" % ("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%10.4g\t%10.4g" % (EF, Ne1, Nh1, NDp1, NAm1, RH, 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 = {} cm-3 (T0 = {} K)".format(NC, T0Nh))
    print(" NV = {} cm-3 (T0 = {} 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)]

# calculate DC(E)*fe(E), DV(E)*(1-fe(E))
    print("")
    print("Calculate fe(E), fh(E)=1-fe(E), DC(E)*fe(E), DV(E)*fh(E) at EF={} eV, T={} K".format(EF0, T0))
    yfe = []
    yfh = []
    yDOSfe = []
    yDOSfh = []
    for i in range(len(E_raw)):
        E = E_raw[i]
        yfe.append(fe(E, T0, EF0))
        yfh.append(fh(E, T0, EF0))
        if E >= EC:
            yDOSfe.append(DOSfe(E, T0, EF0))
        else:
            yDOSfe.append(0.0)
        if E <= EV:
            yDOSfh.append(DOSfh(E, T0, EF0))
        else:
            yDOSfh.append(0.0)

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

    ax2 = fig.add_subplot(2, 2, 1)
    ax2r = ax2.twinx()
    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(E_raw, yDOSfe, label = 'D(E)fe', color = 'blue', linewidth = 0.5)
    ax2.plot(E_raw, yDOSfh, label = 'D(E)fh', color = 'red', linewidth = 0.5)
    ax2.plot([EF0, EF0], [min(dos_raw), max(dos_raw)], color = 'red', linestyle = 'dashed')
    ax2r.plot(E_raw, yfe, label = 'fe', color = 'blue', linewidth = 0.5, linestyle = 'dashed')
    ax2r.plot(E_raw, yfh, label = 'fh', color = 'red', linewidth = 0.5, linestyle = 'dashed')
    ax2.set_ylim([0.0, max(dos_raw) * 1.1])
    ax2r.set_ylim([0.0, 2.0])
    ax2.set_xlabel("E (eV)", fontsize = fontsize)
    ax2.set_ylabel("DOS (states/cm$^3$)", fontsize = fontsize)
    ax2r.set_ylabel("fe, fh", fontsize = fontsize)
#    ax2.legend()
    h2l, l2l = ax2.get_legend_handles_labels()
    h2r, l2r = ax2r.get_legend_handles_labels()
    ax2.legend(h2l + h2r, l2l + l2r, fontsize = legend_fontsize)
    ax3.plot(xEF, yNe, label = 'Ne', color = 'black') #, marker = 'o')
    ax3.plot(xEF, yNh, label = 'Nh', color = 'red') #, marker = 'x')
    ax3.plot(xEF, yNDp, label = 'ND$^+$', color = 'blue') #, marker = 'x')
    ax3.plot(xEF, yNAm, label = 'NA$^-$', color = 'green') #, marker = 'x')
    ax3.plot(xEFapprox, yNeapprox, label = 'Ne(Boltz)', color = 'gray', linewidth = 0.5) #, marker = 'x')
    ax3.plot(xEFapprox, yNhapprox, label = 'Nh(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("EF (eV)", fontsize = fontsize)
    ax3.set_ylabel("N (cm$^{-3}$)", fontsize = fontsize)
    ax3.legend(fontsize = legend_fontsize)
    ax4.plot(xEF, yRH, label = 'R$_H$', color = 'black', marker = 'o')
    ax4.set_xlabel("EF (eV)", fontsize = fontsize)
    ax4.set_ylabel("R$_H$ (C$^{-1}$cm$^{-3}$)", fontsize = fontsize)
    ax4.legend(fontsize = legend_fontsize)
    ax1.plot(xEF, yT0Ne, label = 'T$_0$(Ne)', color = 'black') #, marker = 'o')
    ax1.plot(xEF, yT0Nh, label = 'T$_0$(Nh)', 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("EF (eV)", fontsize = fontsize)
    ax1.set_ylabel("T$_0$ (K)", fontsize = fontsize)
    ax1.legend(fontsize = legend_fontsize)
    ax1.tick_params(labelsize = fontsize)
    ax2.tick_params(labelsize = fontsize)
    ax2r.tick_params(labelsize = fontsize)
    ax3.tick_params(labelsize = fontsize)
    ax4.tick_params(labelsize = fontsize)
# Rearange the graph axes so that they are not overlapped
    plt.tight_layout()

#    plt.pause(0.1)
    print("Close the graph window to terminate this program")
    plt.show()

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

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


def main():
    global mode, T0, EF0
    global Vcell, 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, ED

    updatevars()

    print("mode: ", mode)

    print("DOS configuration")
    print(" Read from [{}]".format(file))
    E_raw, dos_raw = np.genfromtxt(file, 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("")
    EV, EC = FindBandEdges(E_raw, dos_raw, EF0, Egth)
    Eg = EC - EV
    print(" Band edge: EV={} EC={} Eg={} eV".format(EV, EC, Eg))
    EFmin = EV + dEFmin
    EFmax = EC + dEFmax
    EFstep = (EFmax - EFmin) / (nEF - 1)
    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 0 K = ", EF0, " eV")
    Ne0K = Ne(0.0, EF0, -5.0, EF0)
#    Ne0K = Ne(0.0, EF0, Einteg0, EF0)
    print(" Ne at 0 K = ", Ne0K, " /u.c.")
    print("")

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


if __name__ == '__main__':
    main()