In [1]:
import numpy as np
import ipywidgets as ipw

from bokeh.io import push_notebook, output_notebook, show
from bokeh.layouts import row, column
from bokeh.plotting import figure
from bokeh.models import ColumnDataSource
output_notebook()

from ipywidgets import interact

from collections import OrderedDict
old_settings = np.seterr(over = 'ignore') #Ignore warnings about overflow data points
Loading BokehJS ...

This notebook shows evolution of density of states g(E) in 2D by considering all allowed stated in an electron shell between consequtive kx values (n and n+1). The three graphs presented are electron energy dispersion (left), allowed k states for a particular k shell (center) and the density of states g(E) (right).¶

In [2]:
# Set up constants and parameters
m = 9.109e-31; #Define electron mass
hbar = 6.626e-34/2.0/np.pi; #Reduced planck constant
size = 1e-9; #Define the width of the potential well in real space
Emax = 2.0; # Pick the highest energy (eV) the electron could have by distribution
k0 = np.pi/size; #Lowest allowed k in reciprocal space

# Set up the variable space
kF = np.sqrt(Emax*1.602e-19*2*m)/hbar; # Calculate the Fermi wavevector
N = int(np.pi*kF**2/k0**2); # calculate the number of allowed states by the energy constraint
x = np.linspace(-N*k0,N*k0,2*N,endpoint=False); # Double the variable space by symmetry
y = np.linspace(-N*k0,N*k0,2*N,endpoint=False); # Double the variable space by symmetry
kx,ky = np.meshgrid(x,y);
En = (hbar*x*1e9)**2/(2*m);  #Energy for a free electron in eV
In [3]:
# Plot the dispersion relation of the free electrons
freeE = figure(height=300, width=320, title="Dispersion relationship of a free electron",
              tools="pan,reset,save,wheel_zoom", x_range = [x.min(), x.max()],
              y_range = [0,1.1*En.max()]);
freeE.xaxis[0].axis_label='k (nm-1)';
freeE.xaxis[0].formatter.precision = 1;
freeE.xaxis[0].ticker.desired_num_ticks = 4;
freeE.yaxis[0].axis_label='Energy (eV)';
freeE.xaxis[0].formatter.precision = 1;
source0 = ColumnDataSource(data = {'xVal':kx, 'yVal':En});
freeE.scatter('xVal', 'yVal', source = source0, line_width=1, line_alpha=0.2);
sourceC = ColumnDataSource(data = {'x':[x.min(), x.max()],'y':[En[N-5], En[N+5]], 'label':[f'Ef = {En[N-5]:.2f} eV']*2});
freeE.line(x = 'x', y = 'y', legend_field = 'label', source = sourceC, line_color = 'red', line_alpha = 1.0);
#show(freeE)
In [4]:
# Plot the states per energy interval
states = figure(height=300, width=350, title="Allowed states in k-space",
              tools="reset,pan,wheel_zoom", x_range = [1.1*x.min(), 1.1*x.max()],
              y_range = [1.1*y.min(),1.1*y.max()]);
states.xaxis[0].axis_label='kx';
states.xaxis[0].formatter.precision = 1;
states.xaxis[0].ticker.desired_num_ticks = 3;
states.yaxis[0].axis_label='ky';
states.yaxis[0].formatter.precision = 1;
source1 = ColumnDataSource(data = {'x':kx,'y':ky});
x0 = np.linspace(x[N-5],x[N+5],11,endpoint = True);
y0 = np.sqrt(abs(x0.max()**2-x0**2));
x0 = np.concatenate((x0,x0));
y0 = np.concatenate((y0,-y0));
sourceC1 = ColumnDataSource(data = {'r':[x[N+5]-k0/2],'rr':[x[N+5]+k0/2]})
sourceC2 = ColumnDataSource(data = {'x':x0, 'y':y0});
states.rect('x', 'y', width = k0, height = k0, source = source1, line_color = 'black', line_alpha = 1.0, fill_alpha = 0.0);
states.annulus(x = 0, y = 0, inner_radius = 'r', outer_radius = 'rr', source = sourceC1, line_color = 'red', line_alpha = 1.0, fill_alpha = 0.0);
states.rect(x = 'x', y = 'y', width = k0, height = k0, source = sourceC2, line_alpha = 0.0, fill_alpha = 1.0);
#show(states)
In [5]:
#calculate Density of states
#Env = En/6.242e+18;
gE = [m/hbar**2/np.pi]*x.size;
gEplot = figure(height=300, width=320, title="Density of States (2D)",
              tools="pan,reset,save,wheel_zoom", x_range=[En.min(), En.max()], y_range=[0, 1.5*np.amax(gE)])
gEplot.xaxis[0].axis_label='Energy (eV)';
gEplot.xaxis[0].formatter.precision = 1;
gEplot.xaxis[0].ticker.desired_num_ticks = 5;
gEplot.yaxis[0].axis_label='Density of States';
gEplot.yaxis[0].formatter.precision = 1;
source2 = ColumnDataSource(data = {'xVal': En, 'yVal': gE});
sourceC3 = ColumnDataSource(data = {'x':[En[N+5]],'y':[gE[1]], 'label':[f'Ef = {En[N-5]:.2f} eV']});
gEplot.line('xVal', 'yVal', source = source2, line_width=2, line_alpha=1.0);
gEplot.scatter('x', 'y', legend_field = 'label', source = sourceC3, fill_color = 'red', fill_alpha = 1.0, line_width=1.5, line_alpha=1.0, line_color = 'red');
#show(row(freeE, states, gEplot))
In [6]:
# Set up callbacks to live update the plots
def update_data(n):
    # Generate the new curve
    currenE = [f'Ef = {En[N+n]:.2f} eV'];
    sourceC.data = {'x':[x.min(), x.max()],'y':[En[N-n], En[N+n]], 'label':currenE*2};
    
    x0 = np.linspace(x[N-n],x[N+n],2*n+1,endpoint = True);
    y0 = np.sqrt(abs(x0.max()**2-x0**2));
    x0 = np.concatenate((x0,x0));
    y0 = np.concatenate((y0,-y0)); 
    sourceC1.data = {'r':[x[N+n]+k0/2], 'rr': [x[N+n]-k0/2]}
    sourceC2.data = {'x':x0, 'y':y0};
    
    sourceC3.data = {'x':[En[N+n]],'y':[gE[1]], 'label':currenE};
    push_notebook(handle = handle);
In [7]:
handle = show(row(freeE,states,gEplot), notebook_handle = True);
In [8]:
interact(update_data, n = ipw.IntSlider(min = 5, max = 14, step = 1, value = 5, description = 'k/k0'));

interactive(children=(IntSlider(value=5, description='k/k0', max=14, min=5), Output()), _dom_classes=('widget-…
In [ ]: