Calculate effective mass - k relation
from band structure

Download script from .\EffectiveMass.py


import sys
import numpy as np
from numpy import sqrt, exp, sin, cos, tan, pi
import numpy.linalg as LA
import csv
from matplotlib import pyplot as plt


"""
Calculate effective mass - k relation
from band structure
"""



#===================================
# physical constants
#===================================
pi = 3.14159265358979323846
pi2 = 2.0 * pi
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
e0 = 8.854418782e-12; # C2N-1m-2";
kB = 1.380658e-23 # JK-1";
me = 9.1093897e-31 # kg";
R = 8.314462618 # J/K/mol
a0 = 5.29177e-11 # m";


#===================================
# parameters
#===================================
a = 4.0e-10 #m

infile = 'band.csv'
cutline = 1

#===================================
# figure configuration
#===================================
fontsize = 12
legend_fontsize = 8


#=============================
# 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(fname):
    x = []
    y = []
    with open(fname) as f:
        fin = csv.reader(f)
        labels = next(fin)
        xlabel = labels[0]
        ylabel = labels[1]
        for row in fin:
            try:
                x.append(float(row[0]))
                y.append(float(row[1]))
            except:
                print("Warning: Invalid float data [{}] or [{}]".format(row[0], row[1]))

    return xlabel, ylabel, x, y


def main():
    klabel, Elabel, k, E = read_csv(infile)
    print("k=", k)
    print("E=", E)

    nk = len(k)
    dk = k[1] - k[0]

    km = hbar * hbar * (pi2 / a)**2.0
    nskip = 1
    xk = []
    ymc = []
    signprev = None
    print("")
    print("Differentiate for nskip = ", nskip)
    print("i E[i-1] E[i] E[i+1] 1.0/m*")
    for i in range(nskip, nk - nskip, nskip):
        d2Edk2c = (E[i+nskip] + E[i-nskip] - 2 * E[i]) * e / pow(nskip * dk, 2.0)
        minv = d2Edk2c / km
        print(i, E[i-1], E[i], E[i+1], minv)
        if abs(minv) <= 1.0e20: # << 1.0/me ~ 1e30
            if cutline:
                xk.append(k[i])
                ymc.append(None)
            signprev = -signprev
            continue
        else:
            m = km / d2Edk2c

        if signprev is None:
            signprev = m
        elif signprev * m < 0.0:
            if cutline:
                xk.append(k[i])
                ymc.append(None)
            signprev = m

        xk.append(k[i])
        ymc.append(m / me)

    nskip = 4
    xk2 = []
    ymc2 = []
    signprev = None
    print("")
    print("Differentiate for nskip = ", nskip)
    print("i E[i-1] E[i] E[i+1] 1.0/m*")
    for i in range(nskip, nk - nskip, nskip):
        d2Edk2c = (E[i+nskip] + E[i-nskip] - 2 * E[i]) * e / pow(nskip * dk, 2.0)
        minv = d2Edk2c / km
        print(i, E[i-1], E[i], E[i+1], minv)
        if abs(minv) <= 1.0e20: # << 1.0/me ~ 1e30
            if cutline:
                xk2.append(k[i])
                ymc2.append(None)
            signprev = -signprev
            continue
        else:
            m = km / d2Edk2c

        if signprev is None:
            signprev = m
        elif signprev * m < 0.0:
            if cutline:
                xk2.append(k[i])
                ymc2.append(None)
            signprev = m

        xk2.append(k[i])
        ymc2.append(m / me)

    plt.plot(xk, ymc, linewidth = 0.5, marker = 'o', markersize = 1.0, label = 'nskip = 1')
    plt.plot(xk2, ymc2, linewidth = 0.5, marker = 'o', markersize = 1.0, label = 'nskip = 2')
    plt.xlabel(klabel)
    plt.ylabel("m$_e$ / m$_e^0$")
    plt.xlim([-0.5, 0.5])
#    plt.ylim([-0.5, 0.5])
    plt.tight_layout()

    plt.pause(0.1)
    print("Press ENTER to exit>>", end = '')
    input()


if __name__ == "__main__":
    main()