Questaal Home
Navigation

Transmission for a model step potential (ASA)

This tutorial uses empty spheres in the ASA to construct a simple model of a free electrons encountering a step potential, and uses lmpg’s implementation Landauer-Buttiker theory to calculate the transmission through it. It checks how well the ASA does for this analytically known result. The crystal structure is modeled by a bcc lattice of empty sites.

Also demonstrated here is lmpg’s use of normal modes to generate k(E)k(E) in the left end layer. This layer is a flat potential, and the bands should be simple parabolas. The machinery of the ASA indeed recovers this simple result, showing that the tight-binding hamiltonian accurately describes free electron bands.

lmpg is documented on this web page.


Table of Contents

Command Summary

To see structure

lmplan step -vnit=0 -vnk1=4 -vpgf=1 --pledit~q

Setup potential for and and execute transmission calculation

lmpg step -vnit=0 -vnk1=4 -vpgf=1
lmstr step -vnit=1 -vnk1=4 -vpgf=5 --no-iactiv > /dev/null
lmpg step -vnit=1 -vnk1=4 -vpgf=5

Make a picture of the transmission for the first k point

$ mcx jzk.step -inc x2==1 -e2 x1 x3 > dat
$ fplot dat
$ 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.

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


Getting started

This tutorial starts from input file ctrl.step. It consists of a flat potential with a constant barrier of 0.1 Ry in the middle. There are 3+4+2 layers:

layers potential
1..3 0
4..7 0.1
8..9 0

Copy the contents of the box below to file ctrl.step. It is the input (ctrl) file for this system.

Click here to view the ctrl file
# Free Electron Barrier: 0 0 .1 .1 .1 .1 0 0 0
# Adjust width of potential step (variable barrier) to see effect on Fano resonances
% const nk1=4 ef0=0 nz=10 pgf=5 nit=1 barrier=.1
PGF MODE={pgf} PLATL=0 0 2*ha PLATR=0 0 2*ha GFOPTS= p3
VERS LM:7 ASA:7
IO SHOW=f HELP=F VERBOS=41 20 WKP=f IACTIV=f
ITER MIX=B6,w=2,1,b=.2,n=4;A6,w=1,2,n=4 XIPMX=t BETV=.03 NIT={nit} CONVC=1D-6 CONV=0
CONST a=5.44 # lattice constant
pi4b3=4*pi/3 # to construct sphere radius Ra
ha=.5 npl=9 # plane spacing, number of principal layers
vola=a^3*ha Ra=(vola/pi4b3)^(1/3) # sphere radius
bzj=11 # k mesh will be offset around gamma
OPTIONS ASA[ ]
EWALD TOL=1d-12 NKDMX=1500 NKRMX=1500
SYMGRP MX MY R4Z
BZ NKABC={nk1} {nk1} 1 BZJOB=bzj
# for DOS, PGF
% if pgf==1
EMESH=400/1 1 0 0.5 0.0000001
% endif
# for CURRENT, PGF
% if pgf==5
EMESH=200/1 1 0.0 0.5 0.00001
% endif
# for bands, PGF
% if pgf==3
EMESH=400 0 0 1 0
% endif
# for integrated properties, PGF
EMESH={nz} 10 -.8 {ef0} .5
STR RMAX=2.7 MODE=0 SHOW=T EQUIV=t
% const nl=3
STRUC NBAS=npl*2 NSPEC=1 NL={nl}
ALAT=a PLAT= 1 0 0 0 1 0 0 0 10.0
SPEC
# We are dealing with free electrons so
# XA is just an empty sphere with Z=0
ATOM=XA Z=0 IDMOD=0 0 0 R=Ra NR=601 A=.02 GRP2=1
SITE
ATOM=XA POS= 0/2 0/2 0/2 PL=0
ATOM=XA POS= 1/2 1/2 1/2 PL=0
ATOM=XA POS= 0/2 0/2 2/2 PL=1
ATOM=XA POS= 1/2 1/2 3/2 PL=1
ATOM=XA POS= 0/2 0/2 4/2 PL=2
ATOM=XA POS= 1/2 1/2 5/2 PL=2
ATOM=XA POS= 0/2 0/2 6/2 PL=3
ATOM=XA POS= 1/2 1/2 7/2 PL=3
ATOM=XA POS= 0/2 0/2 8/2 PL=4
ATOM=XA POS= 1/2 1/2 9/2 PL=4
ATOM=XA POS= 0/2 0/2 10/2 PL=5
ATOM=XA POS= 1/2 1/2 11/2 PL=5
ATOM=XA POS= 0/2 0/2 12/2 PL=6
ATOM=XA POS= 1/2 1/2 13/2 PL=6
ATOM=XA POS= 0/2 0/2 14/2 PL=7
ATOM=XA POS= 1/2 1/2 15/2 PL=7
ATOM=XA POS= 0/2 0/2 16/2 PL=8
ATOM=XA POS= 1/2 1/2 17/2 PL=8
START NIT={nit} FREE=F BEGMOM={nit==0} CNTROL=T CNVG=1D-6 RDVES=T
ATOM=XA V=0
P= 1.7132985 2.3670680 3.1798926
Q= 0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
0.0000000 0.0000000 0.0000000
ATOM=XA2 V=0
ATOM=XA3 V=0
ATOM=XA4 V=0
ATOM=XA5 V=0
ATOM=XA6 V=0
ATOM=XA7 V={barrier}
ATOM=XA8 V={barrier}
ATOM=XA9 V={barrier}
ATOM=XA10 V={barrier}
ATOM=XA11 V={barrier}
ATOM=XA12 V={barrier}
ATOM=XA13 V={barrier}
ATOM=XA14 V={barrier}
ATOM=XA15 V=0
ATOM=XA16 V=0
ATOM=XA17 V=0
ATOM=XA18 V=0
ATOM=XA19 V=0
ATOM=XA20 V=0
ATOM=XA21 V=0
ATOM=XA22 V=0

Click here for an interpretation of the ctrl file
  • Lines beginning with  #  are comments.
  • Lines beginning with  %  are not part of the input file, but are directives to the file preprocessor. They do not become part of the preprocessed input file, but can modify its contents.
  • Otherwise lines beginning with a nonblank character denote the start of a category. Categories organize tags into groups, as described in the input file documentation. In particular  %const  declares preprocessor variables that may be parsed later in expressions in curly brackets, e.g.  {barrier}. Lines like the following
    % if pgf==1

    % endif
    will add lines between  % if  and  % endif  to the input stream if the conditional expression ( pgf==1  in this case) evaluates to nonzero.
  • %const nk1=4 … barrier=.1  creates variable nk1 and assigns it to 4; it is later parsed as data associated with  NKABC  in the  BZ   category. Its value can be superseded by a command-line argument,  -vnk1=number. Similarly  barrier  is parsed by  V  in the  START  category. It controls the barrier height in this file.

Category  PGF  supplies information needed specific to the layer code lmpg, in particular the two vector defining the cladding layer.  MODE  specifies broadly what lmpg will calculate.

CONST  supplies other variables, similar to the  %const  preprocessor variables. The difference is that they get assigned after the preprocessor completes, and you can use them in expressions without curly brackets  {…} .

VERS,  IO,  ITER,  OPTIONS,  SYMGRP, are all described in the input file documentation, but are not of great importance in this context.

BZ  controls the mesh of k-points in the plane of the interface; see below

STRUC  defines the lattice structure; it is explained below.

SPEC  specifies the chemical species. In this case there is only one, an empty site with no charge in it.

First, see what structure the ctrl corresponds to:

$ lmplan step -vnit=0 -vnk1=4 -vpgf=1 -vsparse=0 --pledit~q

At the top you should see the following:

LMPLAN: nbas = 18 nspec = 1 verb 41,20
PGFSET: 9 principal layers, 18+2+2 sites. Normal: 0 0 1
PL | 0 1 2 3 4 5 6 7 8
size | 2 2 2 2 2 2 2 2 2
PLV | 0 0 0 0 0 0 0 0 0

lmplan says that there are 18 atoms in the active region (supplied by  STRUC_NBAS), cladded by left and right end regions with two atoms in each. There are 9 principal layers in the active region, each with two atoms.

The lattice vectors (taken from  STRUC  in this file) are printed out a little farther down.

LATTIC: Padding Plat(3) with end principal layers:
0.000000 0.000000 10.000000 Plat(3) as input
0.000000 0.000000 1.000000 PlatL: angle (deg) with Plat(3) = 0
0.000000 0.000000 1.000000 PlatR: angle (deg) with Plat(3) = 0
0.000000 0.000000 14.000000 Plat(3) including double padding
Plat Qlat
1.000000 0.000000 0.000000 1.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000 1.000000 0.000000
0.000000 0.000000 14.000000 0.000000 0.000000 0.071429

The first two lattice vectors specify the plane defined by periodic boundary conditions. They are the Cartesian axes x^\mathbf{\hat{x}} and y^\mathbf{\hat{y}} in this case.

The third lattice vector, Plat(3), is of length 10, and is cladded by a left layer of length 1 and a right layer of length 1. A second cladding layer is added to converge the electrostatics. (lmpg uses periodic boundary conditions for now to compute the electrostatic Madelung potential. This will change in future.) The active region with doubly cladded layers combine to make a the third lattice vector of the entire cell, with length 14.

Next follows a list of all the atoms and their projection on the axis normal to the periodic plane.

Gtplan: 22 different planes normal to ( 0.000000 0.000000 1.000000)
Projection of normal onto P1, P2, P3: 0.000000 0.000000 14.000000
ib Class Plane PL x-x.h y-y.h z-z.h h del h
19(L) 19:XA19 1 0 0.00000 0.00000 0.00000 -1.00000 3.50000
20(L) 20:XA20 2 0 0.50000 0.50000 0.00000 -0.50000 0.50000
1 1:XA 3 1 0.00000 0.00000 0.00000 0.00000 0.50000
2 2:XA2 4 1 0.50000 0.50000 0.00000 0.50000 0.50000
3 3:XA3 5 2 0.00000 0.00000 0.00000 1.00000 0.50000
4 4:XA4 6 2 0.50000 0.50000 0.00000 1.50000 0.50000
5 5:XA5 7 3 0.00000 0.00000 0.00000 2.00000 0.50000
6 6:XA6 8 3 0.50000 0.50000 0.00000 2.50000 0.50000
7 7:XA7 9 4 0.00000 0.00000 0.00000 3.00000 0.50000
8 8:XA8 10 4 0.50000 0.50000 0.00000 3.50000 0.50000
9 9:XA9 11 5 0.00000 0.00000 0.00000 4.00000 0.50000
10 10:XA10 12 5 0.50000 0.50000 0.00000 4.50000 0.50000
11 11:XA11 13 6 0.00000 0.00000 0.00000 5.00000 0.50000
12 12:XA12 14 6 0.50000 0.50000 0.00000 5.50000 0.50000
13 13:XA13 15 7 0.00000 0.00000 0.00000 6.00000 0.50000
14 14:XA14 16 7 0.50000 0.50000 0.00000 6.50000 0.50000
15 15:XA15 17 8 0.00000 0.00000 0.00000 7.00000 0.50000
16 16:XA16 18 8 0.50000 0.50000 0.00000 7.50000 0.50000
17 17:XA17 19 9 0.00000 0.00000 0.00000 8.00000 0.50000
18 18:XA18 20 9 0.50000 0.50000 0.00000 8.50000 0.50000
21(R) 21:XA21 21 10 0.00000 0.00000 0.00000 9.00000 0.50000
22(R) 22:XA22 22 10 0.50000 0.50000 0.00000 9.50000 0.50000

There is exactly one atom/plane, and the hamiltonian is short ranged enough that only two planes are needed to make a principal layer. Atoms 19-20 are the left cladding layer; atoms 21-22 are the right cladding layer.


Transmission

First set up the potential and tight-binding structure constants. lmpg uses mode 1 for self-consistent calculations. Here we don’t need self-consistency but we do need to set up the potential parameters, even for this flat potential.

The reduced Green’s function is PSP-S where SS are the structure constants. PP can be build from the potential parameters. By running lmpg with 0 iterations it will make the potential parameters without doing any Green’s function calculations. lmstr makes the structure constants.

$ lmpg step -vnit=0 -vnk1=4 -vpgf=1 -vsparse=0
$ lmstr step -vnit=1 -vnk1=4 -vpgf=1 -vsparse=0 --no-iactiv > /dev/null

Finally, calculate the transmission. lmpg calculates transmission in mode 5.

$ lmpg step -vnit=1 -vnk1=4 -vpgf=5 -vsparse=0

The transmission probabilities are written to files jzk.step (resolved by k) and jz.step (integrated over k)

$ fplot -inc x2==1 -ord x3 jz.step
$ open fplot.ps

Three Fano resonances are seen, corresponding to kinetic energy where each of the three k points crosses the step height, 0.1 Ry.

It is perhaps more instructive to look at the transmission from one k point. Use the mcx calculator to extract two columns of data (energy, transmission) from just the first k point and plot it:

$ mcx jzk.step -inc x2==1 -e2 x1 x3 > dat
$ fplot dat
$ open fplot.ps

The transmission switches from 0 to 1 in a smooth way, around 0.15 Ry. The oscillations are the well known Fano resonances.

View The Figure

Square well barrier transmission

The k point mesh

The k mesh was controlled by tags  BZ_NKABC, which controlled the number of divisions along each of the three axes of the lattice vectors (note for lmpg the third number should always be 1); and by  BZJOB, which tells lmpg whether to include Γ as a mesh point, or defined the grid so that it straddles Γ.

The k-mesh is not printed out unless you use a high verbosity. To see the mesh we used in this example, try the following:

$ lmpg step -vnk1=4 --quit=ham --pr60,20

You should see in the output

BZMESH: 3 irreducible QP from 16 ( 4 4 1 ) shift= T T F
Qx Qy Qz Multiplicity Weight
1 0.125000 0.125000 0.000000 4 0.500000
2 0.375000 0.125000 0.000000 8 1.000000
3 0.375000 0.375000 0.000000 4 0.500000

The mesh has 16 points; all but three are symmetry-equivalent to other points. Note there is no point at k=0.

Try another mesh, e.g.

$ lmpg step -vnk1=5 --quit=ham --pr60,20
$ lmpg step -vnk1=4 -vbzj=0 --quit=ham --pr60,20

The first line also generates 16 points but one of them lies at Γ. In this case there are 6 inequivalent points.

The second line generates 25 points (6 inequivalent ones).


Energy bands of left cladding layer

lmpg has a special mode (PGF_MODE=3) to find k(E)k_\parallel(E) for the left cladding layer, as if it were fully periodic (as opposed to being semi-infinite) on the third axis. It does this by solving a quadratic eigenvalue problem in the principal layers, which give the normal modes for a given energy. This is the inverse of the usual eigenvalue problem, where EE is calculated for a given k\mathbf{k}. On the third axis periodic boundary conditions are not imposed, so kk_\parallel can be complex. In mode 3, lmpg calculates for a particular energy all of the kk_\parallel for a given EE, and writes out those for which Imk\mathrm{Im}k_\parallel is small.

In the present test case the energy bands are free electrons, with dispersion 2k2/(2m)=k2\hbar^2 k^2/(2m) = k^2 since =2m=1\hbar=2m=1 in atomic Rydberg units. We can see how well the ASA recovers the exact result.

In the transmission calculation of the preceding section, we singled out the first k-point and will do the same here. The first point is (0.125,0.125,0) in units 2π/a2\pi/a, so this mode should trace out kk_\parallel for which E(k)=k2+k2E(\mathbf{k}) = k^2_\parallel + k^2_\perp.

Note that ctrl.step has the lines

% if pgf==3
EMESH=200 0 0 .5 0
% endif

When lmpg is run in mode 3, it uses the energy mesh specified by this tag which is a uniformly spaced set of 200 points between 0 and 0.5 Ry on the real axis.

Run the calculation as follows:

lmpg step -vnit=1 -vnk1=4 -vpgf=3

lmpg It will print to stdout and also file bnds.step, for each energy on the mesh values of k(E)k_\parallel(E) for which kk_\parallel is real. Note that values it prints out are in units of 2π/a2\pi/a.

Inspect ctrl.step to see that aa is 5.44 a.u. The analytic values of kk it should find should correspond to E(k)=(k2+k2)2π/a=k2×2π/a+0.0416881E(k_\parallel) = (k^2_\parallel + k^2_\perp)2\pi/a = k^2_\parallel \times 2\pi/a + 0.0416881.

Make a figure comparing the contents of bnds.step to this analytic formula:

$ fplot -s circ:.6 -lt 0 bnds.step -ord '(x1*2*pi/5.44)^2'+0.0416881 -lt 1,bold=3,col=1,0,0 -tp -.5:.5:.01
$ open fplot.ps

It plots circles for k(E)k_\parallel(E) in the file, against a continuous red line for the analytic result. They should agree almost exactly. Note that there are higher lying parabolas as well. These correspond to free electron energy bands in the folded zone.

View The Figure

Free electron energy bands


Things to consider

  1. The transmission for the first k point “turns on” at about 0.14 Ry. If the Γ point had been the first mesh point (see k point mesh). The transmission would have turned on at the barrier height, 0.1 Ry. Explain why the transmission turns on at 0.14 Ry in this test.

  2. For energy EE above a barrier of length LL and height V0V_0, the exact result for the transmission probability is given by

    T=[1+2mV02sin2(kL)4E2k2]1T = {\left[ {1 + \frac{2mV_0^2\sin^2(k'L)}{4E\hbar^2{k'}^2}} \right]^{- 1}}

    where (kL)2=2mL2(EV0)/2(k'L)^2 = 2mL^2(E - V_0)/\hbar^2. Work out numerically what it should be for this case and compare to what was calculated numerically.


Footnotes and references

Implementation of Green’s functions in the ASA context, and its connection to the LMTO-ASA hamiltonian is described in detail in O. K. Andersen, Z. Pawlowska, and O. Jepsen, Phys. Rev. B 34, 5253 (1986).

Implementation of the layer Green’s function technique, including the non-equilibrium capability can be found in S. V. Faleev, et al, Phys. Rev. B71, 195422 (2005).

The original derivation of normal modes used to find the energy band structure of the left cladding layer can be found in Phys. Rev. B39, 923 (1989).