# Band offset tutorial

## Requirements

• nanotools

• RESCU+

• ASE

### Pseudopotentials

I will need the following pseudopotentials.

• pseudos/band_alignment/Al_DZDP.mat

• pseudos/band_alignment/As_DZP.mat

• pseudos/band_alignment/Ga_DZP.mat

Let’s copy them from the pseudo database to the current working directory and export RESCUPLUS_PSEUDO=$PWD. ### References ## Briefing In this tutorial, I explain how to calculate the band offsets at a GaAs/AlAs interface. ## Explanations The calculations are summarized in the following figure from Dandrea et al. The macroscopic potential ($V_{macro}$) is the necessary heterostructure property. It is essentially an averaged version of the effective potential. It should recover it’s bulk value near the center of each material composing the heterostructure. One can then use it to position the valence band maximum ($\tilde E_V$) of each material, a quantity obtainable from bulk calcualtions. The conduction band minimum are obtained by using the experimental values of the band gaps. One can then readily derive the valence and conduction band offsets ($\Delta E_V$and$\Delta E_C$) I shall go through the following steps: 1. Accuracy check for AlAs and GaAs. 2. Calculation of the lattice parameter (a) of zinc-blende AlAs and GaAs. 3. Calculation of the lattice parameter (c) of tetragonal AlAs and GaAs under epitaxial constaints. 4. Calculation of the valence band maximum ($\tilde E_V$) with respect to the macroscopic potential in tetragonal GaAs and AlAs. 5. Calculation of the macroscopic potential ($V_{macro}$) in the relaxed heterostructure. 6. Put everything together and calculate the band offsets. ### Accuracy control In this section, I converge the total energy of AlAs and GaAs with respect to resolution and k-sampling. I shall only present the scripts briefly, for more details read this tutorial. I begin with executing the following Python script for AlAs. from nanotools import Atoms, Cell, System, TotalEnergy from nanotools.checkconv import CheckPrecision a = 5.7 b = a / 2. cell = Cell(avec=[[0.,b,b],[b,0.,b],[b,b, 0.]], resolution=0.20) fxyz = [[0. , 0. , 0. ],[0.25, 0.25, 0.25]] atoms = Atoms(fractional_positions=fxyz, formula="AlAs") sys = System(cell=cell, atoms=atoms) ecalc = TotalEnergy(sys) ecalc.solver.set_mpi_command("mpiexec --bind-to core -n 16") ecalc.solver.set_stdout("resculog.out") calc = CheckPrecision(ecalc, parameter="resolution", etol=1.e-3) calc.solve()  The script for GaAs is the same except that a = 5.6 atoms = Atoms(fractional_positions=fxyz, formula="GaAs")  A resolution of 0.12 Ang yields a total energy error of less than 1 meV per atom for both materials. I then use that resolution in both materials and setup a k-sampling convergence object changing CheckPrecision as follows calc = CheckPrecision(ecalc, parameter="k-sampling", etol=1.e-3)  to find a k-sampling of$8\times 8\times 8$. ### Lattice constants In this section, I show how to obtain the lattice constant of AlAs and GaAs. I use the EquationOfStates calculator to do so. For more details on the EquationOfStates calculator, read this tutorial. The following script sets up 5 calculations for primitive cells of varying volumes. from nanotools import Atoms, Cell, System, TotalEnergy from nanotools.eos import EquationOfState as EOS import numpy as np a = 5.7 b = a / 2. cell = Cell(avec=[[0.,b,b],[b,0.,b],[b,b, 0.]], resolution=0.12) fxyz = [[0. , 0. , 0. ],[0.25, 0.25, 0.25]] atoms = Atoms(fractional_positions=fxyz, formula="AlAs") sys = System(cell=cell, atoms=atoms) sys.kpoint.set_grid([8,8,8]) ecalc = TotalEnergy(sys) ecalc.solver.set_mpi_command("mpiexec --bind-to core -n 16") ecalc.solver.set_stdout("resculog.out") calc = EOS.from_totalenergy(ecalc, relative_volumes=np.linspace(0.9,1.1,5)) calc.solve()  Once the calculations are completed, I use the plot_eos function to show an equation of state fit to the energy and volume data (per atom). from nanotools.eos import EquationOfState as EOS calc = EOS.read("nano_eos_out.json") calc.plot_eos(filename='AlAs-eos.png', show=False)  It will show the following This gives me the equilibrium volume of the FCC primitive cell. To get the equilibrium lattice constant, I multiply by 2 (atoms) and by 4, the FCC to CUB factor, and take the cubic root. I obtain 5.642 Ang. Repeating the procedure for GaAs I find a lattice parameter of 5.588 Ang. ### Epitaxial lattice constants In this section, I show how to obtain the lattice constant of epitaxial AlAs and GaAs under epitaxial constraint. In the previous sections, I found the lattice constant of AlAs to be 5.642 Ang and that of GaAs to be 5.588 Ang. Down the line, I will calculate the band offsets for the (001) interface between AlAs and GaAs. I now make the assumption that each material will take roughly 50% of the strain, which means that the heterojunction will have a lattice constant of 5.615 = (5.642 + 5.588) / 2 Ang in the x- and y-directions. The stress will dissipate in the z-direction resulting in a tetragonal geometry. I will thus recalculate the lattice constant of AlAs in the z-direction using a cubic fit. The following script sets up 5 calculations for tetragonal cells of varying lengths. from nanotools import Atoms, Cell, System, TotalEnergy import numpy as np a = 5.615 b = a / np.sqrt(2) fxyz = [[0. , 0. , 0. ],[0.5 , 0. , 0.25], [0.5 , 0.5 , 0.5 ],[0. , 0.5 , 0.75]] atoms = Atoms(fractional_positions=fxyz, formula="AlAsAlAs") cell = Cell(avec=np.diag([b,b,a]), resolution=0.12) sys = System(cell=cell, atoms=atoms) sys.kpoint.set_grid([8,8,6]) ecalc = TotalEnergy(sys) ecalc.solver.set_mpi_command("mpiexec --bind-to core -n 16") ecalc.solver.set_stdout("resculog.out") vols = []; etots = []; n = 0 for f in np.linspace(0.9,1.1,5): cell = Cell(avec=np.diag([b,b,a*f]), resolution=0.12) ecalc.system.set_cell(cell) ecalc.solver.restart.clear_paths() ecalc.solve(output=f"nano_scf_out_{n}") vols.append(a*f) etots.append(ecalc.get_total_energy_per_atom().m) np.savez("data",vols=vols,etots=etots) n += 1  The energy and volume data is saved to data.npz as it is computed. I then use numpy to perform a cubic fit to get the equilibrium lattice constant: 5.674 Ang. For GaAs, the same procedure yields: 5.574 Ang ### Valence band edge of GaAs I now seek the position of the valence band edge ($\tilde{E_V}$) of AlAs and GaAs with respect to the macroscopic potential ($V_{macro}$). This is a bulk property which I obtain from a band structure calculation of the tetragonal cell. For detailed instructions regarding band structure calculations, refer to tutorial_bs.md. I use the following script to get the ground state potential and band structure of AlAs. import numpy as np from nanotools import Atoms, Cell, System, TotalEnergy from nanotools.bandstructure import BandStructure as BS a = 5.615 b = a / np.sqrt(2) c = 5.674 fxyz = [[0. , 0. , 0. ],[0.5 , 0. , 0.25], [0.5 , 0.5 , 0.5 ],[0. , 0.5 , 0.75]] atoms = Atoms(fractional_positions=fxyz, formula="AlAsAlAs") cell = Cell(avec=np.diag([b,b,c]), resolution=0.12) sys = System(cell=cell, atoms=atoms) sys.kpoint.set_grid([8,8,6]) ecalc = TotalEnergy(sys) ecalc.solver.set_mpi_command("mpiexec --bind-to core -n 16") ecalc.solver.set_stdout("resculog.out") ecalc.solve() bcalc = BS.from_totalenergy(ecalc) bcalc.solve()  I now locate the valence band edge with respect to the macroscopic potential. from nanotools import TotalEnergy, BandStructure as BS import numpy as np ecalc = TotalEnergy.read("nano_scf_out.json") bcalc = BS.read("nano_bs_out.json") veff = ecalc.get_field("potential/effective") vmac = np.mean(veff) vbm = bcalc.energy.get_vbm() print("%20s = %f" % ("vmacro (eV)", vmac.m)) print("%20s = %f" % ("vbm (eV)", vbm.m)) print("%20s = %f" % ("vbm - vmacro (eV)", vbm.m - vmac.m))  I read the ground state results using TotalEnergy.read and load the effective potential ($v_{pp} + v_H + v_{xc}) using the get_field function. The units are automatically converted from Hartree (in the HDF5 file) to eV. Next, I use the Energy class method get_vbm which returns the maximal eigenvalue under the Fermi level. Finally I print the information to screen and get  vmacro (eV) = -20.098377 vbm (eV) = -6.625077 vbm - vmacro (eV) = 13.473300  The same procedure for GaAs yields  vmacro (eV) = -20.566220 vbm (eV) = -6.022645 vbm - vmacro (eV) = 14.543575  ### Relaxation of the heterostructure In this section, I build and relax the heterostructure. When the two materials come into contact, they will exchange electrons near the interface, leading the a so-called interface dipole. The charge transfer induces interface-localized stress, which leads to non-negligible forces on the atoms near the interface. Relaxation of the atoms usually result in an additional dipole which must be taken into account in the band alignment. I thus build and relax the heterostructure as follows. from ase.build import bulk, stack, make_supercell from ase.optimize import BFGS from nanotools.ase.calculators.rescuplus import Rescuplus import numpy as np import os a = 5.615 c = 5.674 alas = bulk('AlAs', crystalstructure='zincblende', a=a, orthorhombic=True) alas.set_cell(alas.cell * np.array([1.,1.,c/a]), scale_atoms=True) alas = make_supercell(alas, [[1,0,0],[0,1,0],[0,0,4]]) alas.center() a = 5.615 c = 5.574 gaas = bulk('GaAs', crystalstructure='zincblende', a=a, orthorhombic=True) gaas.set_cell(gaas.cell * np.array([1.,1.,c/a]), scale_atoms=True) gaas = make_supercell(gaas, [[1,0,0],[0,1,0],[0,0,4]]) gaas.center() slab = stack(alas, gaas) inp = {} inp["system"] = { "cell" : {"resolution": 0.12}, "kpoint" : {"grid":[8,8,1]}} inp["solver"] = {"mix": {"metric": "ker", "precond": "ker"}} inp["solver"]["eig"] = {"lwork" : 2**24, "algo" : "mrrr"} inp["solver"]["restart"] = {"DMkPath": "nano_scf_out.h5"} inp["energy"] = {"forces_return": True} mpicmd = "mpiexec -n 36 " rsccmd = mpicmd+" rescuplus -i PREFIX.rsi && mv nano_scf_out.json PREFIX.rso" calc = Rescuplus(command=rsccmd, input_data=inp) slab.calc = calc os.system("rm -f nano_scf_out.h5") opt = BFGS(slab, trajectory="hetero.traj", logfile="relax.out") opt.run(fmax=0.01)  I create two1 \times 1 \times 4$supercells from the tetragonal cells of GaAs and AlAs (orthorhombic option), then stack them using the stack function of ase.build. I then get a 32-monolayer heterostructure which I can visualize using the view function of the ase.visualize module. I also set solver.restart.DMRPath to nano_scf_out.h5 to employ real space density matrix restart during the relaxation procedure. This will typically save at least a few self-consistent iterations in the computation of the ground state of each atomic configuration. Atomic relaxation is usually most significant within the first one to three monolayers. However, one should be careful to converge his results with respect to heterostructure thickness. During the relaxation, the atomic configurations are saved to hetero.traj, as commanded by the parameter trajectory of the BFGS solver. You can visualize the relaxation process typing ase gui hetero.traj  You may think that atoms move very slightly, but the effect on the interface dipole is actually significant. ### Band diagram of the heterostructure In this section, I create the band diagram from the results obtained in the previous sections. I first concentrate on massaging the effective potential of the heterostructure into the macroscopic potential ($V_{macro}$). I first obtain the average potential ($V_{avg}$) by taking the average in the directions transverse to the heterostructure (here the x- and y-directions). from nanotools import TotalEnergy ecalc = TotalEnergy.read("nano_scf_out.json") fig = ecalc.plot_field("potential/effective", axis=2)  I begin by reading the results with TotalEnergy.read. Then I plot the effective potential$V_{avg}along the z-direction (axis=2) using the method plot_field. You should get the following plot. It is difficult to align the valence bands on the averaged potential since it still shows large oscillations due to the atomic potentials. I thus use smooth_field which will smooth the average field by convolution with an erf pulse. import os os.chdir("/home/vincentm/repos/1-rescuf08/examples/band_offset/relax") from nanotools import TotalEnergy ecalc = TotalEnergy.read("nano_scf_out.json") veff = ecalc.smooth_field("potential/effective", axis=2) fig = ecalc.plot_field(veff, axis=2)  I then get a figure like It is now apparent that the macroscopic potential reaches some bulk value in each material away from the interfaces. It is those values that I use to align the valence band minimum (using\tilde E_V\$ calculated above).

The final step is to summarize everything in a figure and calculate the band offsets. This is done by the following Python script (assuming vmac is defined as in the above script).

import numpy as np
import matplotlib.pyplot as plt
from nanotools import TotalEnergy
vmac = ecalc.smooth_field("potential/effective", axis=2)
z = ecalc.system.cell.get_grid(axis=2)
lz = ecalc.system.cell.get_length(axis=2)

# vmacro-alas
gapa = 2.23
ca = 5.674; lima = np.array([ca, 3*ca])
vmac_alas = np.mean(vmac[np.logical_and(z > lima[0], z < lima[1])])
vbm_alas = vmac_alas + 13.473300
cbm_alas = vbm_alas + gapa
# vmacro-gaas
gapb = 1.52
cb = 5.574; limb = np.array([cb, 3*cb]) + 4*ca
vmac_gaas = np.mean(vmac[np.logical_and(z > limb[0], z < limb[1])])
vbm_gaas = vmac_gaas + 14.543575
cbm_gaas = vbm_gaas + gapb
# offsets
vbm_off = vbm_gaas - vbm_alas
cbm_off = -(cbm_gaas - cbm_alas)
print("vbm offset (+gaas higher): %f (eV)" % vbm_off)
print("cbm offset (+alas higher): %f (eV)" % cbm_off)


There are still some wiggles in the macroscopic potential, but I need one number for each material. I thus average the macroscopic potential over the two central unit cells. I then add the valence band maxima to the macroscopic potentials and subtract them to obtain the valence band offsets. The conduction band offsets are obtained by adding the experimental gap values to the valence band maxima and subtracting again. I finally find

vbm offset (+gaas higher): 0.499948 (eV)
cbm offset (+alas higher): 0.210052 (eV)


The valence band offset is about 50 meV smaller than the experimental value of 0.55 eV. The conduction band offset on the other hand is identical at 0.20 eV. This is partly luck as I have neglected many effects (e.g. spin-orbit coupling). But you can now go through the band alignment steps by yourself for various material combinations.

The band diagram may be plotted as follows

fig = plt.figure()
plt.plot(z, vmac, label='vmacro')
plt.plot(lima, np.ones(2)*(vbm_alas), 'b', label='vbm_alas')
plt.plot(lima, np.ones(2)*(cbm_alas), 'r', label='cbm_alas')
plt.plot(limb, np.ones(2)*(vbm_gaas), 'b', label='vbm_gaas')
plt.plot(limb, np.ones(2)*(cbm_gaas), 'r', label='cbm_gaas')
plt.xticks([0, ca, 2*ca, 3*ca, 4*ca, cb+4*ca, 2*cb+4*ca, 3*cb+4*ca, 4*cb+4*ca, lz])
plt.grid(axis='x')
plt.xlabel("z (ang)")
plt.ylabel("energy (eV)")
plt.show()
fig.savefig("hetero_gaas_alas_diag.png", dpi=600)


which looks as follows