Superlattices in QSGW

This tutorial shows how to construct a superlattice using Questaal tools.

A short-period superlattice of InAs and GaSb is constructed on the [001] axis and a self-consistent QSGW potential is made for it. The superlattice is built with the superlattice editor.

The potential for longer period superlattices are assembled by using the sigma editor, embedding the QSGW self-energy into bulk regions away from the interface, calculated from QSGW calculations for bulk InAs and GaSb. This makes it possible to generate QSGW quality band structure for much large superlattices that would be possible by direct QSGW calculations.

A separate tutorial also shows how to build superlattices. Its purpose is to explain how to calculate transport in a metallic, magnetic trilayer using a Landauer-Buttiker framework.

Remark: This tutorial is written for the classic GW version of the code. It will be updated in due course.

Table of Contents

• Preliminaries
• Basic strategy
• Getting Started
  • Site files for zincblende InAs and GaSb
  • Site files for conventional cells of InAs and GaSb
  • Site, ctrl, and basp files using the (2,2) superlattice as template
  • Additional site files
• LDA + QSGW self-consistency, bulk systems (conventional cells)
  • Energy bands in conventional cells with spin-orbit coupling
• The InAs/GaSb [001] interface
• Site-resolved Partial and Integrated Densities of States
  • Partial DOS resolved by site
  • Integrated DOS to extract position of core levels
• Valence band offset
  • Reference points needed for superlattice band offsets
    • Dipole shifts by aligning d states in the superlattice to the bulk
    • Extracting electrostatic potentials from standard output
    • Aligning superlattice levels to reference
    • Interpretation of (6,6) superlattice
• Atomic Relaxation of Superlattices
• Longer period superlattices
• Footnotes and references

cp init.dat init.inas
blm -vinas=t --sort=no --shorten=no --wsitex,quad1 --findes,rmin=1.3,float,1spec --ehmx=-.3 --nit=20 --gw~gcutb=2.7~gcutx=2.3~pbtol=3e-4,3e-4,1e-3,1e-3,1e-2 --nk=6 --nkgw=4 --gmax=9.2 --eloc=-3.5 --ctrl=ctrl inas
mkdir inas-zb
cp ctrl.inas site.inas inas-zb

cp init.dat init.gasb
blm -vinas=f --wsitex,quad1 --findes,rmin=1.3,float,1spec --ctrl=ctrl gasb
mkdir gasb-zb
cp ctrl.gasb site.gasb gasb-zb

echo .5 .5 0 .5 -.5 0 0 0 1 | lmscell inas --wsite~short~fn=site
lmscell inas '--stack~sort@h3-z/1000~wsite@short@fn=site'
mkdir inas2gasb0
cp site.inas inas2gasb0/site.inasgasb

echo .5 .5 0 .5 -.5 0 0 0 1 | lmscell gasb --wsite~short~fn=site
lmscell ctrl.gasb -vz=11.5/11.43 '--stack~scale=1/z~stretch=z^3~sort@h3-z/1000~wsitex@short@fn=site'
mkdir inas0gasb2
cp site.gasb inas0gasb2/site.inasgasb

cp ./inas-zb/* .
cp inas2gasb0/site.inasgasb site.inas
cp inas0gasb2/site.inasgasb site.gasb
lmscell inas '--stack~show~file@site.gasb@dpos=0,0,-.032~show~wsite@short'
cp sites.inas sitein.inasgasb
blm --molstat --ewald=1d-10 --rdsite --omax=-.04 --addes,end --nosort --shorten=no --findes,rmin=1.3,float,1spec --ehmx=-.3 --nit=20 --gw~gcutb=2.7~gcutx=2.3~pbtol=3e-4,3e-4,1e-3,1e-3,1e-2 --nk=6,6,6 --nkgw=4,4,4 --gmax=9.2 --eloc=-3.5 --ctrl=ctrl inasgasb

lmfa -vles=1 --basfile=basp ctrl.inasgasb
lmfa ctrl.inasgasb
echo -vnk3=4 -vnkgw3=2 > switches-for-lm
mkdir -p inas2gasb2/relax0
cp ctrl.inasgasb basp.inasgasb site.inasgasb switches-for-lm inas2gasb2/relax0

echo -vnk3=4 -vnkgw3=4 > switches-for-lm
cp ctrl.inasgasb basp.inasgasb switches-for-lm inas2gasb0
cp ctrl.inasgasb basp.inasgasb switches-for-lm inas0gasb2

lmscell inas '--stack~show~file site.inas dup=1~file@site.gasb@dpos=0,0,-.032@dup=2~show~wsite@short'
echo -vnk3=2 -vnkgw3=1 > switches-for-lm
cp sites.inas site.inasgasb
mkdir -p inas4gasb4/relax0
cp ctrl.inasgasb basp.inasgasb site.inasgasb switches-for-lm inas4gasb4/relax0

lmscell inas '--stack~show~file site.inas dup=2~file@site.gasb@dpos=0,0,-.032@dup=3~show~wsite@short'
echo -vnk3=2 -vnkgw3=1 > switches-for-lm
cp sites.inas site.inasgasb
mkdir -p inas6gasb6/relax0
cp ctrl.inasgasb basp.inasgasb site.inasgasb switches-for-lm inas6gasb6/relax0

cp inas2gasb0/* .
lmfa ctrl.inasgasb
mpix -n 16 lmf ctrl.inasgasb

lmfgwd --job=-1 ctrl.inasgasb `cat switches-for-lm`
lmgwsc --wt --insul=28 --mpi=16,16 --tol=2e-5 --sym --getsigp --save= inas2gasb0 inasgasb > out.job
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm`; mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm` --rs=1,0 -vnit=1 >  out.lmf
cp rst.inasgasb  out.lmf inas2gasb0

lmf ctrl.inasgasb --wsig:fbz; cp sigm2.inasgasb sigm.inasgasb
head -2 ~/lm/misc/syml.fcc > syml.inasgasb
mpix -n 16 lmf ctrl.inasgasb -vso=1 `cat switches-for-lm` --rs=1,0 -vnit=1
mpix -n 16 lmf ctrl.inasgasb -vso=1 `cat switches-for-lm` --rs=1,0 -vnit=1 --band~fn=syml
echo -7,6 / | plbnds -nocol -fplot~sh~scl=.7 -ef=0 -scl=13.6 -lbl=L,G,X bnds.inasgasb
open fplot.ps

mpix -n 16 lmf ctrl.inasgasb -vso=1 `cat switches-for-lm` --rs=1,0 -vnit=1 --mixsig=0.8
mpix -n 16 lmf ctrl.inasgasb -vso=1 `cat switches-for-lm` --rs=1,0 -vnit=1 --band~fn=syml --mixsig=0.8
echo -7,6 / | plbnds -nocol -fplot~sh~scl=.7 -ef=0 -scl=13.6 -lbl=L,G,X bnds.inasgasb
open fplot.ps

cp $dir/* .
rm -f rst.inasgasb mixm.inasgasb
lmfa ctrl.inasgasb
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm`

lmfgwd --job=-1 ctrl.inasgasb `cat switches-for-lm`
lmgwsc --wt --insul=$nel --mpi=16,16 --tol=2e-5 --sym --getsigp --save= $dir inasgasb > out.job
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm`;  mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm` --rs=1,0 -vnit=1 >  out.lmf
cp rst.inasgasb out.lmf $dir

cp $dir/* . ; rm -f mixm.inasgasb; ln -s sigm sigm.inasgasb
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm` --rs=1,0 -vnit=1
lmf ctrl.inasgasb --wsig:fbz; cp sigm2.inasgasb sigm
mpix -n 16 lmf -vso=1 ctrl.inasgasb -vnk1=$nk1 -vnk2=$nk2 --rs=1,0 -vnit=1 --rsig:fbz --pdos~mode=0~sites=$sitelist --mixsig=.80 --pr35,40 > out.pdos
cp out.pdos $dir

lmdos -vso=1 ctrl.inasgasb -vnk1=$nk1 -vnk2=$nk2 --rs=1,0 -vnit=1 --rsig:fbz --pdos~mode=0~sites=$sitelist --mixsig=.80 --dos:npts=2001:window=-1,1
echo 0.999 4 -3 2 | pldos -fplot~sh~long~open~tmy=.25~dmin=0.40~xl=E -esclxy=13.6 -ef=0 -lst="$doschan" dos.inasgasb
cp dosp.dat $dir/pdos.80
open fplot.ps

lmdos -vso=1 ctrl.inasgasb -vnk1=$nk1 -vnk2=$nk2 --rs=1,0 -vnit=1 --rsig:fbz --pdos~mode=0~sites=$sitelist --mixsig=.80 --dos:npts=40001:window=-3,1:rdm:idos
cp dos.inasgasb $dir/idos.80
fplot -frme 0,.7,0,.4 -x -3,0 -y 0,7 -tmx .1:5 -colsy 1+1:1+4 -qr dos.inasgasb
open fplot.ps

fplot -frme 0,.7,0,.4 -x -1.2,-1 -y 0,5 -tmx .05 -1p -map -i 'x1>-1.2&x1<-1' 1:nr  -colsy 4:24:4 -lt 1,col=0,0,0 -qr dos.inasgasb
open fplot.ps

Preliminaries

This tutorial assumes ~/lm is the top-level directory for the Questaal repository, and that Questaal executables are in your path.

This tutorial uses mpix to execute MPI jobs. Typically mpirun is used for MPI jobs. We assume your system has 16 processors on a core, but you can use any number.

For drawing postscript files, this tutorial assumes you are use the Apple open. Substitute your postscript viewer of choice for open.

Basic strategy

QSGW can make very high quality band structures. In principle it can be done for many atoms, but in practice QSGW is limited to relative small sizes (up to about 50 atoms, depending on the context).

By using the “near sightedness” principle and cutting and pasting self-energies in various regions, we can using small calculations to construct a high-quality potential for much larger systems.

Here we use QSGW to perform calculations on :

1. Bulk InAs in the zincblende structure
2. Bulk GaSb in the a sheared zincblende structure, tailored to match InAs
3. Short-period superlattice of (InAs)_n(GaSb)_n along [001]

Finally the three calculations are merged to construct an approximate self-energy for a longer period superlattice.

There are many steps required to accomplish this, but the result is a very high quality band structure for the superlattice.

This tutorial is somewhat involved, and you may wish to go through some basic tutorial such as the PbTe tutorial, the Bi₂Te₃ tutorial which describes the basis set, and the introductory QSGW tutorial before embarking on this one.

Getting Started

Three files are required to make any calculation. All are autogenerated once supplied structural information.

• A ctrl file, which is the main input file. We will autogenerate this file using the blm utility.
• A site file with structural data. (The ctrl file can also hold structural information, but we keep it separate so only the site changes as the structure changes.) blm also makes this file.
• A basp file, which defines the basis set. The free-atom code lmfa makes this file.
• Along [001], the zincblende lattice has a repeat length of two unit cells. It is simpler (but not essential) to use a doubled unit cell as the building block for the superlattice. We will call this the “conventional cell.” ⁽¹⁾

One more user-supplied file is needed, switches-for-lm. It contains algebraic variables, enabling control of parameters such as the number of k points without altering the ctrl file.

For lmfa to make parameters for the basp file, it needs information about each kind of atom and the augmentation sphere radii. We want to use a fixed basis set for all calculations, which means we need to set up a minimal superlattice which contains all the atoms before doing any electronic structure calculations. This will be (InAs)₂(GaSb)₂ [001] superlattice.

The prescription is then:

• Use blm to build site files for bulk InAs and GaSb in the zincblende structure.
• Use these files and the lmscell utility to assemble a site file for the (InAs)₂(GaSb)₂ [001] superlattice.
• Construct the ctrl and basp files from this structure. They can remain (more or less) fixed for different structures; Only files site and switches-for-lm change between structures.

Site files for zincblende InAs and GaSb

This init file contains the requisite structural information for bulk InAs and GaSb. Copy the box below into file init.dat.

HEADER InAs and GaSb
LATTICE
% ifdef inas
% const a=11.43   # For InAs
% endif
% const a=11.50  # For GaSb
     ALAT={a} PLAT=    0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0

% ifdef inas
SITE
     ATOM=In    X= 0   0   0
     ATOM=As    X= 1/4 1/4 1/4
% else
SITE
     ATOM=Ga    X= 0   0   0
     ATOM=Sb    X= 1/4 1/4 1/4
% endif

Notes:

• init.dat is designed to work with blm to make site files for both InAs and GaSb.
• Both form in the zincblende structure, but InAs and GaSb are slightly lattice-mismatched. We assume here that InAs is thick and that the global lattice constant is that of InAs. To make GaSb conform to that lattice, it must be sheared by elongating it on the [001] direction. We make the simple (and reasonably accurate) assumption that the amount of shear is fixed by keeping the volume constant.
• Reliable treatment of the In and As d states will require local orbitals. blm will automatically include them; however it will be necessary to set the tolerances for identifying them a bit tighter than the default; see --eloc below.

Make the site file for zincblende GaSb:

cp init.dat init.gasb
blm -vinas=f --wsitex,quad1 --findes,rmin=1.3,float,1spec --ctrl=ctrl gasb
mkdir gasb-zb
cp ctrl.gasb site.gasb gasb-zb

Site files for conventional cells of InAs and GaSb

We define the conventional cell⁽¹⁾ to be a doubled unit cell along [001]. This will be the basic building block for (InAs)_n(GaSb)_m superlattices. The InAs and GaSb conventional cells are then the (InAs)₂(GaSb)₀ and The (InAs)₀(GaSb)₂ superlattices.

echo .5 .5 0 .5 -.5 0 0 0 1 | lmscell gasb --wsite~short~fn=site
lmscell ctrl.gasb -vz=11.5/11.43 '--stack~scale=1/z~stretch=z^3~sort@h3-z/1000~wsitex@short@fn=site'
mkdir inas0gasb2
cp site.gasb inas0gasb2/site.inasgasb

Make the corresponding site files for InAs. No shear is required in this case.

cp init.dat init.inas
blm -vinas=t --wsitex,quad1 --findes,rmin=1.3,float,1spec --ctrl=ctrl inas
mkdir inas-zb
cp ctrl.inas site.inas inas-zb

(The ctrl.ext file is retained only because it is needed to construct superlattices).

The conventional cell:

echo .5 .5 0 .5 -.5 0 0 0 1 | lmscell inas --wsite~short~fn=site
lmscell inas '--stack~sort@h3-z/1000~wsite@short@fn=site'
mkdir inas2gasb0
cp site.inas inas2gasb0/site.inasgasb

At the close of lmscell, the following table is printed:

 Found 2 q points in original cell that map to q=0:
             q (2*pi/a)                    q (multiples of qlat)
   0.000000   0.000000   1.000000      0.500000   0.500000   0.000000
   0.000000   0.000000   0.000000      0.000000   0.000000   0.000000

It tells you the other k point (an X point) that gets folded into the Γ point because of the artificial construction of the superlattice.

Site, ctrl, and basp files using the (2,2) superlattice as template

The following instructions build the ctrl file we will use for all calculations, and a site for the (InAs)₂(GaSb)₂ superlattice.

cp ./inas-zb/* .
cp inas2gasb0/site.inasgasb site.inas
cp inas0gasb2/site.inasgasb site.gasb
lmscell inas '--stack~show~file@site.gasb@dpos=0,0,-.032~show~wsite@short'
cp sites.inas sitein.inasgasb
blm --molstat --ewald=1d-10 --rdsite --omax=-.04 --addes,end --nosort --shorten=no --findes,rmin=1.3,float,1spec --ehmx=-.3 --nit=20 --gw~gcutb=2.7~gcutx=2.3~pbtol=3e-4,3e-4,1e-3,1e-3,1e-2 --nk=6,6,6 --nkgw=4,4,4 --gmax=9.2 --eloc=-3.5 --ctrl=ctrl inasgasb

Construct the basis set file, basp.inasgasb

lmfa -vles=1 --basfile=basp ctrl.inasgasb
lmfa ctrl.inasgasb
echo -vnk2=3 -vnkgw3=2 > switches-for-lm
mkdir -p inas2gasb2/relax0
cp ctrl.inasgasb basp.inasgasb site.inasgasb switches-for-lm inas2gasb2/relax0

lmfa creates the basis file basp.inas and the atomic density file atm.ext. It must be run twice, initially to determine what local orbitals are to be included. A second pass is necessary to properly partition the valence-core density. Arguments do the following:

^a This temporarily turns the floating orbitals into “empty spheres” with augmentation radii. Empty spheres are not used in the electronic structure calculations, but the radius is used to make the basis set.

The essential files needed for the superlattice are saved in inas2gasb2/relax0.

Inspect basp.inas. Note that In 4d is the normally augmented valence state; the 5d is treated as a local orbital. The corresponding As 3d state is deeper than the In 4d, so for As, the 3d is treated as a local orbital, and the 4d is the normally augmented valence state.

Copy necessary files to complete the set of required input files for basic building blocks

echo -vnk2=4 -vnkgw3=4 > switches-for-lm
cp ctrl.inasgasb basp.inasgasb switches-for-lm inas2gasb0
cp ctrl.inasgasb basp.inasgasb switches-for-lm inas0gasb2

Additional site files

The rigid shift (0,0,-.032). optimized for the (2,2) superlattice is nearly optimal for any (InAs)_n(GaSb)_m superlattice.

Make various (InAs)_n(GaSb)_n superlattices (with and without shifts)

lmscell inas '--stack~show~file site.inas dup=1~file@site.gasb@dpos=0,0,-.032@dup=2~show~wsite@short'
echo -vnk3=2 -vnkgw3=1 > switches-for-lm
cp sites.inas site.inasgasb
mkdir -p inas4gasb4/relax0
cp ctrl.inasgasb basp.inasgasb site.inasgasb switches-for-lm inas4gasb4/relax0

lmscell inas '--stack~show~file site.inas dup=2~file@site.gasb@dpos=0,0,-.032@dup=3~show~wsite@short'
echo -vnk3=2 -vnkgw3=1 > switches-for-lm
cp sites.inas site.inasgasb
mkdir -p inas6gasb6/relax0
cp ctrl.inasgasb basp.inasgasb site.inasgasb switches-for-lm inas6gasb6/relax0

LDA + QSGW self-consistency, bulk systems (conventional cells)

Make a self-consistent density and potential in the LDA:

cp inas2gasb0/* .
lmfa ctrl.inasgasb
mpix -n 16 lmf ctrl.inasgasb

Note that at the LDA level, there is no bandgap. You can see this in these lines from the standard output:

 BZWTS : --- Tetrahedron Integration ---
 ... only filled or empty bands encountered:  ev=0.100216  ec=0.100216
 VBmax = 0.100216  CBmin = 0.100216  gap = 0.000000 Ry = 0.00000 eV
 BZINTS: Fermi energy:      0.100216;  56.000000 electrons;  D(Ef):    0.000
         Sum occ. bands:  -72.6386525  incl. Bloechl correction:    0.000000

The following makes InAs self-consistent at the QSGW level and saves essential files in inas2gasb0

lmfgwd --job=-1 ctrl.inasgasb `cat switches-for-lm`
lmgwsc --wt --insul=28 --mpi=16,16 --tol=2e-5 --sym --getsigp --save= inas2gasb0 inasgasb > out.job
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm`
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm` --rs=1,0 -vnit=1 >  out.lmf
cp rst.inasgasb  out.lmf inas2gasb0

You should find lmgwsc converges to a high degree of self-consistency in 4 iterations. The first invocation of lmf isn’t necessary. It touches up the restart file making a self-consistent density for the given QSGW Σ⁰.

Before drawing energy bands, check the bandgap:

grep gap out.lmf

You should see that the bandgap is 0.76 eV, significantly larger than the experimental value of 0.42 eV. The error can be traced to several sources:

1. incomplete k convergence
2. neglect of spin-orbit coupling
3. neglect of ladder diagrams in the RPA dielectric function
4. neglect of phonon contributions to the self-energy

All of them are relatively small; but each reduces the gap by 0.05 eV to 0.1 eV.

Energy bands in conventional cells with spin-orbit coupling

Spin-orbit coupling plays a non-negligible role. It is readily taken into account because  ctrl.inas is already set up for a calculation with SO coupling.

However if you run lmf straight off it will fail with the message

 Exit -1 rdsigm unexpected value 21 for file sigm nqp ... expected 64

This is because the QSGW self-energy was constructed on the irreducible wedge of the zincblende lattice; however spin-orbit coupling reduces the symmetry. The simplest solution (though not most efficient) workaround is to remake the QSGW self energy file Σ⁰ at all k points

lmf ctrl.inasgasb --wsig:fbz
cp sigm2.inasgasb sigm.inasgasb

The following steps will:

• find the Fermi level with SO coupling included (only needed to tell the band plotter knows where to set the energy zero)
• generate the energy bands
• make a postscript file

We use a symmetry lines file for the fcc lattice from the Questaal distribution. since we have an artificial supercell, use only the first two symmetry lines :

head -2 ~/lm/misc/syml.fcc > syml.inasgasb
mpix -n 16 lmf ctrl.inasgasb -vso=1 `cat switches-for-lm` --rs=1,0 -vnit=1
mpix -n 16 lmf ctrl.inasgasb -vso=1 `cat switches-for-lm` --rs=1,0 -vnit=1 --band~fn=syml
echo -7,6 / | plbnds -nocol -fplot~sh~scl=.7 -ef=0 -scl=13.6 -lbl=L,G,X bnds.inasgasb
open fplot.ps

lmf is invoked initially to establish the Fermi level. This is not essential but it is convenient when drawing the energy bands.

Click here to view QSGW energy bands of InAs

Inspect the first row of bnd2.dat (the Γ point is the first point in the second panel of the energy bands). You should find the following:

The split-off band, at −0.36 eV, is split slightly less than the experimental splitting (−0.43 eV). The light-hole and heavy-hold bands should be degenerate; there is a slight splitting because the cubic symmetry was artificially broken⁽¹⁾ (the c axis is not equivalent, e.g. the k mesh is different). The bandgap, 0.66 eV, still significantly larger than the experimental value of 0.42 eV. Band gaps are overestimated, and indeed they should be because the RPA screened coulomb interaction W is not sufficiently screened. And also zero-point motion has been left out. Improvement of W by the addition of ladder diagrams indeed does improve the gaps, as Kutepov has shown. However, these higher order diagrams are expensive, and we adopt a simpler approach here. The dielectric constant is about 80% of experiment (it seems always to be the case for semiconductors and insulators), because ladder diagrams are missing in the RPA. We have found that scaling ε by 0.8 (assuming the full $\epsilon(\mathbf{r},\mathbf{r'},\omega)$ gets scaled by a fixed number), or alternatively using a hybrid of 80% QSGW and 20% LDA we can mimic the effect of the ladders. This eliminates most of the error.

mpix -n 16 lmf ctrl.inasgasb -vso=1 `cat switches-for-lm` --rs=1,0 -vnit=1 --mixsig=0.8
mpix -n 16 lmf ctrl.inasgasb -vso=1 `cat switches-for-lm` --rs=1,0 -vnit=1 --band~fn=syml --mixsig=0.8
echo -7,6 / | plbnds -nocol -fplot~sh~scl=.7 -ef=0 -scl=13.6 -lbl=L,G,X bnds.inasgasb

Click here to view the band structure of InAs and GaSb comparing QSGW with hybrid QSGW, $(\Sigma^0{\times}0.8{+}V_\mathrm{xc}^\mathrm{LDA}{\times}0.2)$. With the scaling, agreement with ellipsometry measurements is excellent (shown as circles). The fundamental gap drops to 0.42 eV in InAs.

GaSb proceeds in much the same way as InAs. Repeat the preceding LDA+QSGW instructions, substituting inas2gasb0 → inas0gasb2. The direct gap is 1.21 eV, larger than the experimental value 0.82 eV. With SO coupling and scaling it drops to 0.78 eV. Note that to get more perfect agreement with experiment, the scaling factor 0.8 can be modified slightly.

The InAs/GaSb [001] interface

The procedure for a superlattice is nearly identical to the bulk. We have already generated the necessary input files.

Set shell variable dir to inas2gasb2/relax0 or inas4gasb4/relax0, or another directory where a site file resides. The inas2gasb2 case is useful to construct the ctrl file and to determine the forces at the interface but the system is to small to be of relevance for electronic structure.

Start with LDA self-consistency:

cp $dir/* .
lmfa ctrl.inasgasb
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm`

You should find that the forces are all small. Forces on the frontier atoms will be modest but still less that 10 mRy/au. Lattice relaxations are discussed later.

For the QSGW calculation, do the following. The number of electrons (set here through shell variable nel) isn’t essential.

lmfgwd --job=-1 ctrl.inasgasb `cat switches-for-lm`
lmgwsc --wt --insul=$nel --mpi=16,16 --tol=2e-5 --sym --getsigp --save= $dir inasgasb > out.job
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm`; mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm` --rs=1,0 -vnit=1 >  out.lmf
cp rst.inasgasb  out.lmf $dir

Self-consistency should be reached in 4 iterations. (InAs)₂(GaSb)₂ takes about 30 minutes/iteration; (InAs)₄(GaSb)₄ about 12 hours/iteration; (InAs)₆(GaSb)₆ about two days.

Site-resolved Partial and Integrated Densities of States

If you are not familiar with Questaal’s implementation of Densities-of-States (DOS) resolved by site, it is recommended that you review the tutorial on partial DOS in Cr₃Si₆ before going through this section.

By projecting DOS onto sites, we can resolve the superlattice DOS layer by layer and compare to bulk InAs and GaSb. This is useful in determining the valence band offset. Here we set up the DOS for the reference bulk systems.

The instructions below refer to bulk InAs, bulk GaSb, and the (InAs)₄(GaSb)₄ superlattice. The quantities in the table below distinguish which system the instructions apply to

Assign shell variables to each of these quantities, e.g. in bash, dir=inas2gasb0 nk1=12 nk2=9 sitelist=1:8:2 doschan="1:4" for the InAs case.

To make the DOS and iDOS, a moments file must be created.

cp $dir/* . ; rm -f mixm.inasgasb
ln -s sigm sigm.inasgasb
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm` --rs=1,0 -vnit=1

The last step is not needed, but run it anyway to confirm that the RMS DQ is small and that the bandgap is as it should be (0.76 eV in InAs).

It is important to include spin-orbit coupling, so we need to render the QSGW  Σ⁰  suitable for reduced symmetry operations. Here we use the 80% hybrid to well reproduce the experimental band structure. Do:

lmf ctrl.inasgasb --wsig:fbz; cp sigm2.inasgasb sigm
mpix -n 16 lmf -vso=1 ctrl.inasgasb -vnk1=$nk1 -vnk2=$nk2 --rs=1,0 -vnit=1 --rsig:fbz --pdos~mode=0~sites=$sitelist --mixsig=.80 --pr35,40 > out.pdos
cp out.pdos $dir

Note the large number of k-divisions used. This is not necessary but a fine k mesh generates a smooth DOS. The --pdos switch is explained in some detail in the partial densities of states tutorial. See also documentation for the command-line switch for –pdos.

lmf generates file moms.inasgasb, which contains the weights necessary to generate site-resolved DOS. In what follows, DOS are generated in two forms: first to extract the partial DOS near the band edges, for drawing local DOS, and also the integrated DOS (iDOS) to extract the position of the deep d states, which will be used align reference levels in the bulk with the superlattice, to extract the valence band offset.

Partial DOS resolved by site

Partial DOS are made using the lmdos utility, once lmf has generated the appropriate file. This utility generates the projection of DOS into augmentation spheres², given a moms file generated by a band program, as documented in the partial densities of states tutorial, and also documentation for the --dos command-line argument. You must invoke lmdos with the same switches that were supplied to lmf, and also an additional argument, --dos.

lmdos -vso=1 ctrl.inasgasb -vnk1=$nk1 -vnk2=$nk2 --rs=1,0 -vnit=1 --rsig:fbz --pdos~mode=0~sites=$sitelist --mixsig=.80 --dos:npts=2001:window=-1,1

On exit, lmdos should print the following table (inas2gasb0 case)

Channels in dos file generated by LMDOS:
 site class label   spin-1                       spin-2
    1    1   In     1                            2
    3    2   As     3                            4
    5    4   In2    5                            6
    7    5   As2    7                            8

As the table shows, DOS are resolved for eight channels, interleaving (spin-up,spin-down). This system is nonmagnetic so the two are equal and there is no point in distinguishing them. Also the two In and two As should be equivalent. lmdos writes the DOS to file dos.inasgasb in a (not-particularly friendly) format. pldos reads files in this format to make figures, as documented in the pldos manual.

There are thus four channels per layer. Since we don’t need to resolve DOS by site yet, make a figure of the In+As DOS for a single layer by combining DOS for a single layer, i.e. channels 1:4. For the bulk InAs and GaSb cases, shell variable doschan should be set to “1:4”. Make a postscript file for the DOS with (instructions also apply to the superlattice, as described below)

echo 0.999 4 -3 2 | pldos -fplot~sh~long~open~tmy=.25~dmin=0.40~xl=E -esclxy=13.6 -ef=0 -lst="$doschan" dos.inasgasb
cp dosp.dat $dir/pdos.80
open fplot.ps

fplot~... tells pldos to create a DOS data file dosp.dat, formatted in Standard Questaal format for xy data. The second command saves this file for future use, to draw reference and superlattice DOS on the same figure. pldos also makes a script plot.dos for the fplot utility, which will in turn generate a postscript file, fplot.ps. (You can edit plot.dos to customize the figure as we will do when aligning superlattice and reference DOS). Modifiers ~long~open~tmy=.125~dmin=0.40~xl=E make the figure more aesthetic, while ~sh tells pldos to run fplot in a subshell, making fplot.ps directly. -esclxy=13.6 converts ordinate and abscissa from atomic units to eV and eV⁻¹, respectively, and the abscissa is shifted to place the Fermi level at zero (ef=0).

The arguments from standard input (0.999 5 -3 2) cap the DOS to just under 1, with a frame of about 5 cm, in an energy window (−3,2)  eV.

Superlattice DOS

For the (InAs)₄(GaSb)₄ or (InAs)₆(GaSb)₆ superlattice, follow essentially the same steps from the beginning of the pdos section, and it is instructive to resolve the DOS by layer. Now there are a multiplicity of layers (8 or 12) and 32 or 48 channels, four channels per layer; pldos will draw 8 or 12 panels of DOS, one panel for each layer. doschan is constructed to group the DOS for each layer into one panel (see sitelist and doschan shell variables explained).

Integrated DOS to extract position of core levels

Invoke lmdos as before, but change the --dos switch to --dos:npts=40001:window=-3,1:rdm:idos. This generates the integrated DOS on a very fine energy mesh over a wide energy window, writing data in Standard Questaal format for xy data.

lmdos -vso=1 ctrl.inasgasb -vnk1=$nk1 -vnk2=$nk2 --rs=1,0 -vnit=1 --rsig:fbz --pdos~mode=0~sites=$sitelist --mixsig=.80 --dos:npts=40001:window=-3,1:rdm:idos
cp dos.inasgasb $dir/idos.80

Inspect dos.inasgasb. Energy is stored in the first column and the  iDOS  starting in column 2. The second instruction saves the  iDOS  for later, when we align bands to extract the valence band maximum.

For the bulk cases, e.g. InAs, draw a picture of the integrated DOS over a broad energy window (−3,1) Ry:

fplot -frme 0,.7,0,.4 -x -3,0 -y 0,7 -tmx .1:5 -colsy 1+1:1+4 -qr dos.inasgasb
open fplot.ps

Four  iDOS  are drawn but spin-1 and spin-2 overlay each other. In one curve, at −2.8 Ry there are two closely spaced step-function changes in the  iDOS  to 2, and then 5 electrons. These are the five 5 As d states split into 2 + 3 levels by spin-orbit coupling. The other curve shows a similar shape near −1.1 . Above −0.8 Ry where the valence bands begin, the  iDOS  begin to rise again.

For the superlattices there are many levels. To see the evolution of In d levels in a (−1.2,−1) Ry energy window, try

fplot -frme 0,.7,0,.4 -x -1.2,-1 -y 0,5 -tmx .05 -1p -map -i 'x1>-1.2&x1<-1' 1:nr  -colsy 4:24:4 -lt 1,col=0,0,0 -qr dos.inasgasb
open fplot.ps

For the strained bulk GaSb or the superlattice, repeat the steps starting at Site-resolved Partial Densities of States, substituting the appropriate set of shell variables.

Valence band offset

Reference points needed for superlattice band offsets

Interfacial effects — apart from possible long range electrostatic shifts — die out quickly. Then the local valence band maximum (l-VBM), sufficiently far from an interface, will not depend on the interface. A good way to determine the alignment of the l-VBM on two sides of an interface, is to choose a reference energy or potential whose position relative to the l-VBM is constant far enough away from the interface; call it E_ref. Two good choices for E_ref are:

• the position of semicore states such as 3d or Sb 4d
• the potential on the augmentation sphere boundary

This construction enables us to determine a local valence band maximum (l-VBM) for any layer, but because of nonlocal quantum effects, its meaning is not precise close to the interface.

It is common practice to align bands from separate bulk calculations. This is how “natural band offsets” are determined, for example. Aligning bands for two independent materials relies on an assumption, namely that if an interface were formed between two materials, no change in the dipole contribution to the electrostatic potential would arise with respect superposition of the potentials from the independent materials. The “natural band offsets” scheme aligns bands from independent calculations of two surface structures. Other assumptions can be made, for example from the following. Imagine a lattice of spheres that fill space. If each sphere is neutral, a sensible definition of the electrostatic potential is to choose 0 at the sphere boundary. As a next step, transfer charge between spheres. The electrostatic potential (the Madelung potential) is now completely defined, since the reference was defined. In any case the assumption can be tested by comparing how the potential differs from the potential of an actual interface. Since in fact the offset depends on details of the relaxation of the interface, any assumption is likely to be suspect.

Absent additional dipoles forming at a junction relative to the bulk potentials as defined by the implementation, valence band maxima are merely those of the bulk systems.

This information is contained in the file out.pdos, of the bulk directory, e.g. for InAs (directory inas2gasb0) there should be this table:

 BZWTS : --- Tetrahedron Integration ---
 ... only filled or empty bands encountered:  ev=0.120480  ec=0.151794
 VBmax = 0.120480  CBmin = 0.151794  gap = 0.031314 Ry = 0.42587 eV

Extract the VBM for the two bulk systems:

grep -A 2 'VBmax' inas2gasb0/out.pdos | head -1 | awk '{print $3}' > vb1
grep -A 2 'VBmax' inas0gasb2/out.pdos | head -1 | awk '{print $3}' > vb2

If no additional dipole forms at the interface the band offset would be (0.120480-0.164545)×13.6 = −0.599 eV. That is, the VBM of InAs falls below that of GaSb by about 0.6 eV.

We need to add to this the dipole shift, which can be obtained by taking the difference in reference energies.

Dipole shifts by aligning d states in the superlattice to the bulk

We align the d states (16 of them in the (4,4) superlattice), against the corresponding reference states. Thus we need

1. Energies of the bulk InAs + bulk GaSb d states
2. Energies of the d states in the superlattice

For the former, figures were already drawn identifying the levels pictorially. There is some arbitrariness in how a numerical value is extracted, but the arbitrariness is not large except for a possible uniform constant which is irrelevant. For convenience let’s pick the energy where the electron count is exactly 4 for each lattice. This readily accomplished with the mcx calculator. For the bulk systems

mcx -qr inas2gasb0/idos.80 -at 4 x2  -w:nohead . > d1
mcx -qr inas2gasb0/idos.80 -at 4 x4  -w:nohead . >> d1
mcx -qr inas0gasb2/idos.80 -at 4 x2  -w:nohead . > d2
mcx -qr inas0gasb2/idos.80 -at 4 x4  -w:nohead . >> d2

-at val expr  finds the abscissa where the ordinate matches  expr. Note that the first and third channels (In spin-1 and As spin-1) lie in columns 2 and 4. Inspect files d1 and d2. You should get about −1.038 Ry for the In level and −2.673 Ry for the As level.

The d levels of a superlattice are a bit more complicated as there are many of them.
They are readily extracted in mcx using the looping construct [ ... ]. Set shell variable dir according to the table and invoke mcx for the (4,4) superlattice as follows:

mcx -qr $dir/idos.80 -a idos [ ix=1:31:2 idos -at 4 x\{ix+1} -tog -pop '?ix>1' -rcat ] > d3

For the (6,6) superlattice substitute 1:31:2 → 1:47:2.

Inspect file d3. You should see the 16 d levels of the superlattice. The difference in the bulk and superlattice data is the local shift in VBM. Compute this difference and convert to eV. For the (4,4) superlattice:

mcx [ ix=1:4 d1 '?ix>1' -rcat ] [ ix=1:4 d2 '?ix>1' -rcat ]  -rcat d3 -ccat -e2 x3 x1 -e1 x1-x2 -s13.6

For the (6,6) superlattice, substitute 1:4 → 1:6.

Inspect the output of this command. The (4,4) superlattice still has strong interface effects, so we discuss the (6,6) case instead. Around row 6 the shift is about 0.05 eV; this corresponds to the middle (InAs) layer. Around row 18 it is 0.01 eV. Note that all the elements can be shifted by a uniform constant, which arises from arbitrariness in the choice of the G=0 part of the potential. Thus net dipole shift is about 0.04 eV.

Extracting electrostatic potentials from standard output

An alternate way to obtain the band offset is to align bulk and superlattice electrostatic potentials on the augmentation sphere boundary. Apart from interfacial effects, it should yield the same result as the d level shifts.

lmf calculates the potential at each augmentation sphere boundary and sets the average to zero. This information is also contained in out.pdos; for example for InAs there should be a table

 Average es pot at MT boundaries after vconst shift
 Site    ves        Site    ves        Site    ves        Site    ves
   1  -0.033546  |    3   0.036403  |    5  -0.033548  |    7   0.036403
   2   0.714034  |    4   0.662523  |    6   0.714049  |    8   0.662523

Sites 2,4,6,8 are floating orbital sites which we will ignore for the purpose of aligning bulk and superlattice states.

For the (InAs)₄(GaSb)₄ and (InAs)₆(GaSb)₆ superlattices there are 32 or 48 potentials, of which the odd numbered sites correspond to the real atoms.

The following instructions extract these potentials from out.pdos for the two bulk systems and the (InAs)₄(GaSb)₄ superlattice, saving results in files v1, v2, v3. In each of the instructions grep, awk and sed are used to extract the table and render it into a form of xy data that mcx can read. Then mcx extracts the relevant potentials.

grep -A 2 'Average es pot at MT boundaries' inas2gasb0/out.pdos | awk '$1 ~ /^[0-9]+$/{print }' | sed 's/|//g' | mcx -nc=2 . -sort x1 -inc 'i%4==1|i%4==2' > v1
grep -A 2 'Average es pot at MT boundaries' inas0gasb2/out.pdos | awk '$1 ~ /^[0-9]+$/{print }' | sed 's/|//g' | mcx -nc=2 . -sort x1 -inc 'i%4==1|i%4==2' > v2
... for the (4,4) superlattice:
grep -A 9 'Average es pot at MT boundaries' inas4gasb4/relax0/out.pdos | awk '$1 ~ /^[0-9]+$/{print }' | sed 's/|//g' | mcx -nc=2 . -sort x1 -inc 'i%2==1' > v3
... for the (6,6) superlattice:
grep -A 13 'Average es pot at MT boundaries' inas6gasb6/relax0/out.pdos | awk '$1 ~ /^[0-9]+$/{print }' | sed 's/|//g' | mcx -nc=2 . -sort x1 -inc 'i%2==1' > v3

Aligning superlattice levels to reference

We now have files d1…d3… for the d levels and files v1…v3… for the electrostatic potentials.

It is convenient to tabulate the levels as a function of height from the basal plane. This is easily accomplished with the lmplan utility. For one of the superlattices, e.g. the (6,6) superlattice, set shell variables according to the table and do

cp $dir/* .
echo q | lmplan ctrl.inasgasb

lmplan prints out a table of positions of each site its projection h on the axis perpendicular to the basal plane:

 Gtplan: 24 different planes normal to (    0.000000    0.000000    1.000000)
 Projection of normal onto P1, P2, P3:      0.000000    0.000000    6.055456
   ib      Class   Plane  PL       x-x.h     y-y.h     z-z.h      h       del h

    1       1:In      1    0      0.00000   0.00000   0.00000   0.00000   0.28662
    2       5:E       1    0      0.50000   0.00000   0.00000   0.00000   0.00000
    3       2:As      2    0      0.25000   0.25000   0.00000   0.25000   0.25000
    4      10:E2      2    0      0.25000  -0.25000   0.00000   0.25000   0.00000
    5       6:In2     3    0      0.50000   0.00000   0.00000   0.50000   0.25000
    6      15:E3      3    0      0.00000   0.00000   0.00000   0.50000   0.00000
    7       7:As2     4    0      0.25000  -0.25000   0.00000   0.75000   0.25000
    8      20:E4      4    0      0.25000   0.25000   0.00000   0.75000   0.00000

Extract h from this table and store in file h, excluding floating sites:

echo q | lmplan ctrl.inasgasb | grep -A 100 Plane | grep -v :E | awk '$1 ~ /^[0-9]+$/{print $8,"#",$2}' > h

Make a single file with the superlattice and bulk electrostatic potentials, and another with the corresponding d levels. Note the conversion to electron volts. Also plot these changes as a function of h, solid lines for the electrostatic potential and dashed for the d levels. The following applies for the (4,4) case; for the (6,6) superlattice, substitute 1:4 → 1:6.

mcx v3 [ ix=1:4 v1 '?ix>1' -rcat ] [ ix=1:4 v2 '?ix>1' -rcat ] -rcat -ccat -e2 x2 x4 -s13.6 -w:nohead . | pr -m -t - h > estat
mcx d3 [ ix=1:4 d1 '?ix>1' -rcat ] [ ix=1:4 d2 '?ix>1' -rcat ] -rcat -ccat -e2 x1 x3 -s13.6 -w:nohead . | pr -m -t - h > dlevel
fplot -frme 0,.7,0,.7 -ab x3 -ord x1-x2 estat -ab x3 -ord x1-x2 -lt 2 dlevel

In the figure, the vertical dot-dashed lines show the center point of the (InAs) and (GaSb) layers, respectively. Note the following points:

1. The two references are slightly inconsistent, especially at the boundary layers. These are consequences of finite-size effects.
   Even where the potential is most bulk-like there is still a nearly constant 0.02 eV discrepancy.
2. A global constant can be added to either of these shifts, since they are determined by the arbitrary definition of the G=0 potential.
   Taking this into account, the relevant dipole shift is the change at the midpoints (compare blue and red lines). The shift is approximately 0.05 eV.
3. The potential shift is not quite flat even in the middle regions. This is not an artifact of the small cell: if a still longer cell were calculated it would become clear that there is a linear term, which if many layers are present, would cause a dielectric breakdown. This is to be expected because the [001] interface is polar. In the real world the crystal will reconstruct somehow, possibly by introducing charged defects at the interface.
4. The dipole adds to the “natural” (or reference) valence band offset which is −0.599 eV. Thus the band offset is 0.55 eV.

Interpretation of (6,6) superlattice

With this information in hand let’s revisit the partial DOS of the (6,6) superlattice and observe the following

1. DOS in the frontier layers (InAs)₁, (InAs)₆, (GaSb)₁ (GaSb)₆ differs significantly from the bulk counterpart. GaSb valence bands just near the valence band maximum tunnel into the InAs (see (InAs)₁ near the VBM); similarly InAs conduction band states near the conduction band minimum tunnel into the GaSb (see (GaSb)₁ near the CBM). These states decay exponentially away from the interface.
2. Even while the band offset exceeds the bulk InAs bandgap (0.43 eV), the gap does not close, even in the (6,6) superlattice. Look for the gap in inas6gasb6/relax0/out.pdos: it is 0.28 eV. This appears to originate from confinement effects: the GaSb acts like an effective quantum well confining the (highly extended) states in the InAs conduction band.

To more clearly see this, make and inspect the energy bands along the c axis around a (−1,1) eV window around the Fermi energy

echo   11   0  0  0     0 0 .5  > symlz.inasgasb
mpix -n 11 lmf -vso=1 ctrl.inasgasb -vnk1=12 -vnk3=2 --rs=1,0 -vnit=1 --rsig:fbz --mixsig=.80 --band~mq~col=1:30,dup=936,seq=79:391:78~col2=40:69,dup=936,seq=118:430:78~fn=symlz
echo -1,1,5,10 | plbnds -fplot~scl=.7~sh -ef=0 -scl=13.6 -dat=green -lt=1,bold=3,col=0,0,1,colw=1,0,0,colw2=0,1,0 -lbl=G,X  bnds.inasgasb

Bands colored blue are of GaSb character. There is a nearly dispersionless state at the VBM (E=0). To the extent a band is red, it is of InAs character. The lowest conduction band has significant, but by no means all InAs character: InAs has a very small effective mass, and this state it is higher than it would be in a longer period superlattice because of confinement effects. The extended range of this state is also evident by its significant dispersion. States colored green are of As character. The highest valence-band state of nearly pure As character lies at around −0.8 eV. Those between −0.8 eV and −0.1 eV are of mixed GaSb/As character, with increasing GaSb approaching the E=0.

Atomic Relaxation of Superlattices

This section explains how a limited relaxation of the (2,2) superlattice is accomplished. As noted when it was constructed (see lmscell and blm arguments explained The (2,2) superlattice was used as a template to construct ctrl.inasgasb, and a the (GaSb) layers were stacked on with a particular offset (0,0,−0.032) that minimized the total energy.

Here we explain how that shift was determined by parameterizing the shift as (0,0,delta) and examining the total energy and forces as a function of delta

Set shell variable delta to (0,0,−0.032), e.g. in the bash shell, delta=-.032, and revisit the generation of the (InAs)₂(GaSb)₂ site file

cp ./inas-zb/* .
cp inas2gasb0/site.inasgasb site.inas
cp inas0gasb2/site.inasgasb site.gasb
lmscell inas "--stack~show~file@site.gasb@dpos=0,0,$delta~show~wsite@short"
cp inas2gasb2/relax0/* .
cp sites.inas site.inasgasb

You should find, for delta=−0.032, that site.inasgasb and inas2gasb2/relax0/site.inasgasb are identical, but change with delta.

Repeat the instructions above for delta=0 and make the calculation self-consistent

rm -f mixm.inasgasb rst.inasgasb
lmfa ctrl.inasgasb
mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm` -vdelta=$delta

Forces are printed in a table from stdout:

 Forces, with eigenvalue correction
  ib           estatic                  eigval                    total
   1    0.00    0.00  -13.34     0.00    0.00   60.85     0.00    0.00   47.51
...
   7    0.00    0.00   18.44     0.00    0.00   17.52     0.00    0.00   35.96
   9    0.00    0.00   26.55     0.00    0.00  -56.90     0.00    0.00  -30.34
...
  15    0.00    0.00   -7.94     0.00    0.00  -45.71     0.00    0.00  -53.65

Forces are large and opposite at the (1,15) interface, and similarly so at the (7,9) interface By rigidly translating the GaSb block all of the forces should drop.

Repeat the setup and relaxation for delta=−0.01 and the maximum force should drop from 53.7 to 38.9. Also you should find that the total energy drops by about 9 mRy. Repeat the steps for the sequence delta=0, −0.01, −0.02, −0.03, −0.04 and you should see the following emerge

The mean absolute value of the force is minimum around delta=−0.034, as is the minimum in total energy. For the SL calculations we selected delta=−0.032 — not quite optimum because a slightly cruder minimization process was used. At the relaxed point the maximum force is 9.6 (delta=−0.032) and 7.2 (delta=−0.034).

The lattice should properly be fully relaxed. This is readily accomplished through the DYN category. Adding -vminx=t includes this category in the input:

mpix -n 16 lmf ctrl.inasgasb `cat switches-for-lm` -vminx=t > out &

but the forces are already small at the first stage and the additional relaxation hardly affects the band structure.

Longer period superlattices

… To be finished

Footnotes and references

¹ The zincblende lattice doubled along the z is a conventional cell in the sense that all primitive lattice vectors are orthogonal. It differs from the proper conventional for the zincblende lattice (a cube, with 8 atoms), and is not a true conventional cell because it reduces the symmetry from cubic to tetragonal symmetry.

² The In+As DOS do not sum to the total DOS because there is a missing interstitial contribution. This is unimportant for the present purposes.