Aestimo Tutorial #3

In this tutorial, we will look at modelling a structure with a parabolic potential rather than the stepped potentials that aestimo usually considers. In order to do this, we will need to understand how aestimo's code is structured in more depth.

This example is a replication of that given in Harrison's book[1] in Ch.3 sec.5.

First, we import the libraries that we need as usual.

In [3]:
import aestimo.aestimo as solver
import aestimo.config as ac
ac.messagesoff = True # turn off logging in order to keep notebook from being flooded with messages.
import aestimo.database as adatabase
import numpy as np
import matplotlib.pyplot as plt
import copy
from pprint import pprint

Then, we redefine some of the database entries since Harrison's initial examples don't take account of different effective masses in the different materials.

In [4]:
meV2J = solver.meV2J # conversion factor
q = solver.q # electron charge

adatabase.materialproperty = {
        'Eg':0.0, #1.426, # set the energy scale origin to be at the GaAs condution band
        'm_e':0.067, # normally Harrison would be using 0.15
        'Eg':2.673-1.426, #2.673, # set the energy scale origin to be at the GaAs condution band

adatabase.alloyproperty = {

Here we define the system's properties. However, we can't create our 'Structure' object in the usual way using the StructureFrom class, instead we will work directly with its parent class Structure. We willh ave to create some of the structure object's attributes ourselves, specifically the arrays that describe the system w.r.t. position.

The Structure class is used by Aestimo to hold all of the details about the system that is being modelled. In essence, it is a collection of attributes hanging off a class instance plus a method for calculating the effective mass (array) with respect to energy for those times that we ask non-parabolicity to be included in the simulation. These attributes are either scalars or they are 1d arrays of parameters w.r.t. position across the structure.

We normally use the subclass StructureFrom instead of the Structure class since StructureFrom knows how to turn the input files into a Structure object.The StructureFrom class uses the method create_structure_arrays() to calculate the attribute arrays from the material list structure that we normally use in our input files. So we also need to run this method if we ever update any of the model's parameters during a simulation so that the Structure object's arrays are up to date.

First, here are the parameters that StructureFrom needs:

  • inputfile object

    • material - list of layers describing structure. Each layer is a list of [Thickness (nm) | Material | Alloy fraction | Doping(cm^-3) | n or p type ]
    • T - (float) Temperature (K)
    • gridfactor - grid step size (m)
    • computation_scheme - (int) computing scheme - see example files for details.
    • meff_method - chooses effective mass energy dependent function for times that nonparabolicity is included in simulation - see example files for details. (not needed for the aestimo_eh solver
    • fermi_np_scheme - (bool) chooses whether to calculate fermi-level including non-parabolicity terms to in-plane dispersion (not needed for the aestimo_eh solver)
    • subnumber_e - number of subbands to look for.
    • subnumber_h - number of valence subbands to look for. (only needed for aestimo_numpy_h solver)
    • Fapplied - applied field (Vm**-1)
    • maxgridpoints - maximum number of grid points allowed (used only to prevent simulation of structures that are accidentally too large)
  • database object/module - we give the reference explicitly in case we want to supply a customised object. This works because modules, classes and namedtuples can often substitute for each other in python.

However, if we want to create our own version of the Structure object for aestimo to use, we will need to define

Custom Structure Class

  • value attributes

    • T - Temperature (K)
    • dx - grid step size (m) (gridfactor in input files)
    • comp_scheme - (int) computing scheme
    • meff_method - (int) effective mass function (see above) (only aestimo not aestimo_eh)
    • fermi_np_scheme - (bool) (see above) (only aestimo not aestimo_eh)
    • subnumber_e - number of subbands to look for.
    • subnumber_h - number of valence bands (only aestimo_numpy_h)
    • Fapp - applied field (Vm**-1) (Fapplied in input files)
    • n_max - number of grid points
    • x_max - total thickness of structure (m) (currently optional)
  • array attributes

    • fi - Bandstructure potential (J) (array, len n_max)
    • eps - dielectric constant (including eps0) (array, len n_max)
    • dop - doping distribution (m**-3) (array, len n_max)
    • cb_meff - conduction band effective mass (kg) (array, len n_max)
  • method attributes

    • cb_meff_E(E,fi) - given the energy and the effective potential array, returns an array for the effective mass.
  • array attributes 2 (non-parabolicity) - These attributes are only necessary if you are using a non-parabolic model for the effective mass of conduction band elections.

    • cb_meff_alpha - non-parabolicity constant. Used for Nelson's empirical 2-band model[2])
    • Eg - band gap energy (Used for k.p model found in Vurgaftman[3])
    • Ep - (Used for k.p model found in Vurgaftman[3])
    • F - (Used for k.p model found in Vurgaftman[3])
    • delta_S0 - spin split-off energy (Used for k.p model found in Vurgaftman[3])
  • array attributes 3 (valence band modelling) - In aestimo_numpy_h there are many more arrays that should be defined but I will not discuss this any further since things are getting complicated enough.

In [9]:
s0 = {} # this will be our datastructure

s0['T'] = 300.0 #Kelvin

# 0: Schrodinger
# 1: Schrodinger + nonparabolicity
# 2: Schrodinger-Poisson
# 3: Schrodinger-Poisson with nonparabolicity
# 4: Schrodinger-Exchange interaction
# 5: Schrodinger-Poisson + Exchange interaction
# 6: Schrodinger-Poisson + Exchange interaction with nonparabolicity
s0['comp_scheme'] = 0

# Non-parabolic effective mass function
# 0: no energy dependence
# 1: Nelson's effective 2-band model
# 2: k.p model from Vurgaftman's 2001 paper
s0['meff_method'] = 0 

# Non-parabolic Dispersion Calculations for Fermi-Dirac
s0['fermi_np_scheme'] = True

# Total subband number to be calculated for electrons
s0['subnumber_e'] = 10

s0['Fapp'] = 0.00 # (V/m)

# For 1D, z-axis is choosen
gridfactor= 0.01 #nm
s0['dx'] = gridfactor*1e-9
maxgridpoints = 200000 #for controlling the size

# a finite parabolic well between two barriers
a = 10.0 #nm #maximum width of parabolic well (at the barrier's energy level)
b = 10.0 #nm #width of the first barrier layer
b2 = 5.0 #nm #width of the first barrier layer
xmin = 0.0 #minimum Al alloy in structure
xmax = 10.0 #maximum Al alloy in structure

def alloy_profile(z):
    """function of alloy profile for finite parabolic well"""
    well = xmin + (z - (b+a/2.0))**2 / (a/2.0)**2 * (xmax - xmin)    
    return np.where(well<xmax,well,xmax) #cuts off parabolic well at barrier height

def bandstructure_profile(x):
    mat1 = adatabase.materialproperty['AlAs']
    mat2 = adatabase.materialproperty['GaAs']
    alloy = adatabase.alloyproperty['AlGaAs']
    return alloy['Band_offset']*(x*mat1['Eg'] + (1-x)* mat2['Eg']-alloy['Bowing_param']*x*(1-x))*solver.q # for electron. Joule
    #return alloy['Band_offset']*(x*mat1['Eg'] + (1-x)* mat2['Eg'])*aestimo.q # for electron. Joule

# Nb. we could have avoided using the database module at all since we are 
# defining our own Structure object and calculating its attribute arrays 
# ourselves but I wanted to show how easy it is to customise the database's values.

## Create structure arrays -----------------------------------------------------

# Calculate the required number of grid points
s0['x_max'] = sum([b,a,b2])*1e-9 #total thickness (m)
s0['n_max'] = n_max = int(s0['x_max']/s0['dx'])
# Check on n_max
if s0['n_max'] > maxgridpoints:
    solver.logger(" Grid number is exceeding the max number of %d", maxgridpoints)
meff = adatabase.materialproperty['GaAs']['m_e']*solver.m_e
epsGaAs = adatabase.materialproperty['GaAs']['epsilonStatic']

z = np.arange(0.0,s0['x_max'],s0['dx'])*1e9 #nm

s0['cb_meff'] = np.ones(n_max)*meff	#conduction band effective mass
s0['cb_meff_alpha'] = np.zeros(n_max)   #non-parabolicity constant.
s0['eps'] = np.ones(n_max)*epsGaAs	#dielectric constant
s0['dop'] = np.ones(n_max)*0.0          #doping
s0['fi'] = bandstructure_profile(alloy_profile(z))   #Bandstructure potential

# Initialise structure class
model = solver.Structure(**s0)
In [10]:
#config.d_E = 1e-5*meV2J
#config.Estate_convergence_test = 3e-10*meV2J
result= solver.Poisson_Schrodinger(model)

#Plot QW representation
%matplotlib inline
ac.wavefunction_scalefactor = 5000
calculation time  2.03413 s
INFO:aestimo:calculation time  2.03413 s

We can see that (as expected) the parabolic energy levels are equally spaced, at least for the lower levels which are not so sensitive to the finite height of the barrier potentials.

In [11]:
Energies = np.array(result.E_state)
Energies[1:] - Energies[:-1]
array([ 871.88179039,  871.87321772,  871.86216698,  871.82993394,
        871.66001195,  870.77565   ,  866.77248077,  850.23643253,
  1. Paul Harrison, Quantum Wells, Wires and Dots, 3rd. Ed. (Wiley, 2010).
  2. D.F. Nelson et al., Phys. Rev. B 35, 7770 (1987)
  3. I. Vurgaftman, J.R. Meyer, JR, L.R. Ram-Mohan, Band parameters for III–V compound semiconductors and their alloys, Journal of applied physics, 89(11), 5815–5875 (2001).