# Converting a periodic lattice into a trilayer geometry (ASA)

This tutorial shows how render a self-consistent calculation for crystal lattice into a trilayer geometry suitable for Landauer-Buttiker transport through an active region sanwiched between semi-infinite leads. This tutorial works with the cubic compound CsPbI3, but lattice and basis vectors are deformed slightly to simulate the effect of finite temperature.

As a first step a self-consistent calculation is performed with the band code lm, and also with the crystal Green’s function code lmgf, in a periodic structure containing one formula unit. It is shown that the two codes generate the same densities.

Next, the single formula unit is doubled, using the superlattice editor. The system is still treated as periodic, but with a doubled unit cell. It is shown that the density made self-consistent for the single formula unit remains so in the doubled cell, for both lm and lmgf.

Finally two-unit superlattice is converted into a trilayer geometry suitable for the layer Green’s function code lmpg. The first unit becomes both the first layer of the active region and the the semi-infinite left lead, and the second becomes a second layer in the active region and the semi-infinite right lead. It is shown that the entire trilayer (consisting of the same formula unit in each layer) continues to be self-consistent.

Single Formula Unit

cp sitein.cspbi
blm --asa --gf --nk=8,8,8 --nl=3 --omax=0.16,0.21,0.25 --wsrmax=3.3 --rdsite --noshorten --nfile --wsite cspbi
lmscell ctrl.cspbi --stack~sort@h3~show@planes~wsite@short@fn=site1
lm -vfile=1 ctrl.cspbi -vnit=0
lmstr -vfile=1 ctrl.cspbi
mpix -n 16 lm -vnit=100 -vfile=1 -vmet=2 -vccor=f -vtwoc=t ctrl.cspbi --rs=0,1 > out
cp rsta.cspbi rsta1.cspbi [for backup only]
lmgf -vgamma=t -vgf=1 -vnit=-1 -vfile=1 -vmet=2 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho


Doubled Formula Unit

lmscell -vfile=1 ctrl.cspbi cspbi --stack~rsasa@fn=rsta.cspbi@rdim~file@site1.cspbi@rsasa@rdim@fn=rsta.cspbi~wsite:short:fn=site2:rsasa=rsta2
cp rsta2 rsta.cspbi
lmstr -vfile=2 ctrl.cspbi
lm -vfile=2 ctrl.cspbi -vnit=0 --rs=1,0
mpix -n 16 lm -vnit=1 -vfile=2 -vmet=2 -vccor=f -vtwoc=t ctrl.cspbi --quit=rho
mpix -n 16 lmgf -vgamma=t -vgf=1 -vnit=-1 -vfile=2 -vmet=2 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho


Conversion from Crystal to Trilayer Geometry

lmscell -vfile=2 -vpgf=2 cspbi --stack~rsasa@fn=rsta2~setpl:sort:range=9:shorps=0,0,1~wsite:fn=site3:rsasa=rsta3 >> out.lmscell.cspbi3
cp rsta.cspbi rsta2.cspbi [for backup only]
cp rsta3 rsta.cspbi
lmpg -vgamma=t -vpgf=1 -vnit=0 -vfile=3 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho
lmstr -vgamma=t -vpgf=2 -vnit=1 -vfile=3 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho
lmpg -vgamma=t -vpgf=2 -vnit=1 -vfile=3 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho
lmpg -vgamma=t -vpgf=1 -vnit=-1 -vfile=3 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho


### 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, or run serially by removing the prefix mpix -n 16 where it appears.

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

### Single Formula Unit

In this section the ctrl and site files are prepared.

Copy the contents of the box below to file site.cspbi. It is a site file for a pseudocubic structure of CsPbI3. (The cubic structure has higher symmetry, but this structure provides a reasonable representation of CsPbI3 at room temperature.) Empty spheres have already been added to satisfy the space-filling requirement of the ASA.

% site-data vn=3.0 fast io=15 nbas=16 alat=11.82540057 plat= 1.00496473 0.0 0.0 0.00024653 0.9953719 0.0  0.0223297 -0.0000584 1.0183501
#                        pos
Cs        0.5238550   0.4975648   0.5289579
Pb       -0.0336829  -0.0002003   0.9843686
I        -0.1007007   0.4975195   0.0129042
I        -0.0804758  -0.0001690   0.4726700
I         0.4549649  -0.0001821   0.9282127
E         0.4217061  -0.0000008   0.4664323
E        -0.0160774   0.4975539   0.4836709
E         0.4451329   0.4976335   0.9997768
E1       -0.3508811   0.2180172   0.2044577
E1       -0.3509509  -0.2183507   0.2045600
E1        0.2004869   0.2553894   0.2279732
E1        0.2002804  -0.2556073   0.2282405
E1        0.1967869  -0.2731194   0.7477251
E1        0.1969437   0.2729484   0.7477949
E1       -0.2615046  -0.2467878   0.7555042
E1       -0.2614808   0.2464068   0.7555879


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

# Autogenerated from init.cspbi using:
# blm --asa --gf --nk=8,8,8 --nl=3 --omax=0.16,0.21,0.25 --wsrmax=3.3 --rdsite --noshorten --nfile --wsite cspbi

# Variables entering into expressions parsed by input
% const nit=10 asa=t
% const met=(asa?1:5)
% const so=0 nsp=so?2:1
% const lxcf=2 lxcf1=0 lxcf2=0     # for PBE use: lxcf=0 lxcf1=101 lxcf2=130
% const ccor=T twoc=F sx=0 gamma=sx scr=4 # ASA-specific
% const gf=1 nz=16 ef=-.238 c3=t pgf=0 sparse=0 # lmgf-specific variables
% const nk1=8 nk2=nk1 nk3=nk2 beta=.3

VERS  LM:7 ASA:7 # FP:7
IO    SHOW=f HELP=f IACTIV=f VERBOS=35,35  OUTPUT=*
EXPRESS
# Lattice vectors and site positions
%ifdef file>=1
file=   site{file}
%endif
file=   site

# Self-consistency
nit=    {abs(nit)}               # Maximum number of iterations
mix=    B2,b={beta},k=7          # Charge density mixing parameters
conv=   1e-5                     # Convergence tolerance (energy)
convc=  3e-5                     # tolerance in RMS (output-input) density

# Brillouin zone
nkabc=  {nk1},{nk2},{nk3}        # 1 to 3 values
metal=  {met}                    # Management of k-point integration weights in metals

# Potential
nspin=  {nsp}                    # 2 for spin polarized calculations
so=     {so}                     # 1 turns on spin-orbit coupling
xcfun=  {lxcf},{lxcf1},{lxcf2}   # set lxcf=0 for libxc functionals

SYMGRP  find

OPTIONS
ASA[ CCOR={ccor} TWOC={twoc} ADNF=F GAMMA={gamma}] SCR={scr} SX={sx}

HAM     PMIN=-1                      # Constrain Pnu > P_free
#       RDVEXT=2

GF    MODE={gf} GFOPTS={?~c3~p3;~p2;}
BZ    EMESH={nz},10,{ef}-1,{ef},.5,.3  # nz-pts;contour mode;emin;emax;ecc;bunching
PGF   MODE={pgf} SPARSE={sparse} PLATL=0.02232969 -0.00005842 1.01835015 PLATR=0.02232969 -0.00005842 1.01835015

SPEC
# SCLWSR=11 WSRMAX=3.3 OMAX1=0.16 0.21 0.25
ATOM=Cs       Z= 55  R= 3.300000  LMX=2
ATOM=Pb       Z= 82  R= 3.300000  LMX=2 IDMOD=0,0,1
ATOM=I        Z= 53  R= 3.300000  LMX=2
ATOM=E        Z=  0  R= 3.237117  LMX=2
ATOM=E1       Z=  0  R= 2.471028  LMX=1

START CNTROL={nit<=0} BEGMOM={nit<=0}


The main part the ctrl was autogenerated using the blm utility. There are some modifications (added by hand) to the autogenerated template, as the box below explains.

blm was used to generate an initial template ctrl file from site.cspbi

$blm --asa --gf --nk=8,8,8 --nl=3 --omax=0.16,0.21,0.25 --wsrmax=3.3 --rdsite --noshorten --nfile --wsite cspbi  blm’s command-line arguments are summarized here. • –asa and –gftarget the ctrl for the ASA codes lm, lmgf, lmpg • –nl=3, –omax=0.16,0.21,0.25, –wsrmax=3.3 affect parameters in augmentation spheres. These values were selected to allow spheres to fill space, while not letting any sphere get too large. • –rdsite, –noshorten, and –wsitecause blm to read data from sitein.cspbi, and write the structural data without modifying it. blm creates a template actrl.cspbi. Some modifications were necessary; the following diff command reveals the modifications of the template: $ diff actrl.cspbi ctrl.cspbi


The output of diff is shown below, with annotations:

< % const ccor=T sx=0 gamma=sx scr=4 # ASA-specific
< % const gfmode=1 nz=16 ef=0 c3=t   # lmgf-specific variables
---
> % const ccor=T twoc=F sx=0 gamma=sx scr=4 # ASA-specific
> % const gf=1 nz=16 ef=-.238 c3=t pgf=0 sparse=0 # lmgf-specific variables

<       ASA[ CCOR={ccor} TWOC=F ADNF=F GAMMA={gamma}] SCR={scr} SX={sx}
---
>       ASA[ CCOR={ccor} TWOC={twoc} ADNF=F GAMMA={gamma}] SCR={scr} SX={sx}


Variable twoc is added to enable control of token TWOC Reasons for this are explained in the main text.

Variable gfmod was shortened to gf. It is use to control program flow in lmgf (MODE={gf})

Variable ef is used by lmgf to set the energy integration window. An estimate of it was obtained by running lm to self-consistency (see below).

Variables pgf and sparse are used by lmpg to control program flow there;

< GF    MODE={gfmode} GFOPTS={?~c3~p3;~p2;}
< BZ    EMESH={nz},10,-1,{ef},.5,.3  # nz-pts;contour mode;emin;emax;ecc;bunching
---
> HAM     PMIN=-1                      # Constrain Pnu > P_free
> #       RDVEXT=2
>
> GF    MODE={gf} GFOPTS={?~c3~p3;~p2;}
BZ>     EMESH={nz},10,{ef}-1,{ef},.5,.3  # nz-pts;contour mode;emin;emax;ecc;bunching
> PGF   MODE={pgf} SPARSE={sparse} PLATL=0.02232969 -0.00005842 1.01835015 PLATR=0.02232969 -0.00005842 1.01835015

• A new HAM_PMIN= was added for the same reason IDMOD was added (see below)
• The energy window in BZ_EMESH was modified
• PLATL and PLATR are tokens required by lmpg to define the left- and right- semi-infinite regions; PGF_GFOPTS carries switches specific to lmpg. Note that they are the same vectors as the third vector in PLAT (inspect the site file).
<   ATOM=Pb       Z= 82  R= 3.300000  LMX=2
---
>   ATOM=Pb       Z= 82  R= 3.300000  LMX=2 IDMOD=0,0,1


The ASA can be a bit tricky when floating the band center-of-gravity for localized orbitals. The Pb d orbital is such a case, and IDMOD=0,0,1 freezes its logarithmic derivative at its starting point (which is a safe value). HAM_PMIN=-1 (above) also restricts logarithmic derivatives in all channels to remain above the free-electron value

51c56
<   ATOM=E1       Z=  0  R= 2.471028  LMX=2
---
>   ATOM=E1       Z=  0  R= 2.471028  LMX=1


The radius of the second empty sphere (E1) is small enough to reduce the basis to sp orbitals on that site. This significantly reduces the rank of the hamiltonaian at minimal loss of accuracy.

When later we convert to a trilayer geometry, the trilayer axis will be parallel to PLATL and PLATR (approximately along z). Definition of the trilayer requires partitioning of the sites into principal layers.

lmscell can make this construction automatically using its superlattice editor, but it requires that the sites be ordered by increasing distance along the trilayer axis. It is advisable at this stage to rearrange the sites to simplify synchronization of the crystal and trilayer geometries.

The following instruction

$lmscell ctrl.cspbi --stack~sort@h3~show@planes~wsite@short@fn=site1  • sorts the basis by h3 (h3 is the trilayer axis), • displays the projection of (now reordered) site positions onto h3 (it is called h in the printout) • writes the reordered structural information to file site1.cspbi. site.cspbi and site1.cspbi contain the same structural information but order the basis differently. The rest of the tutorial will use site1.cspbi (or later site2.cspbi); which file is used will be controlled from the command line. Editor instructions are documented here. #### Self-consistent density through the band code In the following we use the band code lm to a self-consistent ASA density $ lmstr -vfile=1 ctrl.cspbi
$lm -vfile=1 ctrl.cspbi -vnit=0$ mpix -n 16 lm -vnit=100 -vfile=1 -vmet=2 -vccor=f -vtwoc=t ctrl.cspbi --rs=0,1 > out
$cp rsta.cspbi rsta1.cspbi  The first two steps create the structure matrix and a (crude) trial density, and the next step converges the density to self-consistency. The (last) copy command is not necessary but it is useful if you want to redo the tutorial starting from this point. (Just copy rsta1.cspbi to rsta.cspbi). In the iterations to reach self-consistency, there are several points to note: • At the end of the file you should see  it 17 of100 ehf= -99965.003750 ehk= -99965.003749 From last iter ehf= -99965.003748 ehk= -99965.003748 diffe(q)= -0.000001 (0.000027) tol= 0.000010 (0.000030) more=F c nit=100 file=1 met=2 ccor=0 twoc=1 ehf=-99965.0037495 ehk=-99965.0037494 cpudel ... Time this iter: time(s): 1.58 total: 28.1s Jolly good show ! You converged to rms DQ= 0.000061  Self-consistency was reached in 17 iterations. Note that the HK and HF functionals are nearly identical. • See what Fermi level lm finds: $ grep Fermi out


The Fermi energy should be about −0.238 Ry. This is basis for selecting variable ef=−.238 in the ctrl file. lm doesn’t need this tag, but an energy window (part of BZ_EMESH) are parameters lmgf and lmpg require to perform a contour integration.

• Note the Hamitonian rank
$lmgf -vgamma=t -vgf=1 -vfile=1 -vmet=2 -vccor=f -vtwoc=t ctrl.cspbi -vnit=-1 --rs=1,0 --quit=rho > out  • -vgamma=t tells lmgf to work in the “gamma representation.” On the real axis, the poles are equivalent, but there can be false poles as well. They are pushed far from the real axis in this representation. Using this switch can be especially important when third-order (i.e. three-center) potential functions are used. • -vgf=1 selects the mode lmgf operates in (self-consistency cycle) • vnit=-1 makes the START category appear as START CNTROL=1 BEGMOM=1  which causes lmgf to begin with energy moments $Q_0$, $Q_1$, and $Q_2$ of the density, which are sufficient to generate a density and corresponding potential. The potential is used in turn to construct the hamiltonian for the Green’s function. • --rs=1,0 --quit=rho causes lmgf to read $Q_0$, $Q_1$, and $Q_2$ from rsta.cspbi, and also to stop after the charge density is generated. rsta.cspbi stores this information by site. You should find something the following at the end of the output:  PQMIX: read 0 iter from file mixm. RMS DQ=1.45e-5 AMIX: nmix=0 mmix=8 nelts=160 beta=0.3 tm=10 rmsdel=1.45e-5 h gamma=1 gf=1 file=1 met=2 ccor=0 twoc=1 nit=-1 ehf=-99965.0036905 ehk=-99965.0036905  Deviation from exact self-consistency (1.45e-5) is small, establishing that the same density is self-consistent in lm and lmgf. Midway through the output you should see  --- GFASA (gamma) : make qnu, idos GFASA: integrated quantities to efermi = -0.238 PL D(Ef) N(Ef) E_band 2nd mom Q-Z 1.923116 31.999960 -20.311026 15.739134 -0.000040 deviation from charge neutrality: -0.00004  It shows that the Fermi level (−0.238 Ry) found by lm also corresponds to the charge-neutrality point in lmgf. Also inspect the potential shift file vshft.cspbi. The first line should read something similar to  ef=-0.2380000 vconst=-0.0033597  lmgf shifted the crystal potential to keep the charge neutrality point fixed at −0.238 Ry. The shift is small (about 3 mRy) and causes the Fermi level to fall just above the valence band maximum. The system is an insulator without a density-of-states at the Fermi level. Thus, the Fermi level (or the shift) is very “soft” and is sensitive to details, e.g. the number of points on the contour energy integration. ### Doubled Formula Unit The trilayer geometry requires both left semi-infinite and right semi-infinite regions. The original unit cell can serve as both, if we duplicate the unit and double the cell. Then the first unit will correspond to the left region and the second to the right. To prepare for the trilayer geometry, construct a doubled cell with periodic boundary conditions using the superlattice editor $ lmscell -vfile=1 ctrl.cspbi cspbi --stack~rsasa@fn=rsta.cspbi@rdim~file@site1.cspbi@rsasa@rdim@fn=rsta.cspbi~wsite:short:fn=site2:rsasa=rsta2


See the user’s manual for the editor instructions.

This invocation of lmscell doubles the unit cell along the third axis, creating both a site file for structural data (file site2.cspbi) and an rsta file for sphere moments rsta2. Thus it automatically sets up files for a self-consistent potential in the doubled cell.

To confirm that this is the case, do the following.

• copy rsta2 to the appropriately named file
• make the structure constants for the doubled cell
• make the potential from the moments in rsta.cspbi
• run lm one pass to confirm it is self-consistent in the doubled cell
• run lmgf one pass to confirm it is self-consistent in the doubled cell
$cp rsta2 rsta.cspbi$ lmstr -vfile=2 ctrl.cspbi
$lm -vfile=2 ctrl.cspbi -vnit=0 --rs=1,0$ mpix -n 16 lm -vnit=1 -vfile=2 -vmet=2 -vccor=f -vtwoc=t ctrl.cspbi --quit=rho -vnk3=4
$mpix -n 16 lmgf -vgamma=t -vgf=1 -vnit=-1 -vfile=2 -vmet=2 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho -vnk3=4  A new argument -vnk3=4 halves the number of k divisions on the third axis, while keeping fixed the first two axes (see the ctrl file). Since the unit cell was doubled along the third axis, the number of k divisions can be halved on that axis to recover the same mesh. Search for RMS DQ in the output and note that it remains less than 10−5 for both lm and lmgf. Note also that vshft.cspbi is essentially the same as before. ### Conversion from Crystal to Trilayer Geometry Having established a self-consistent density with two formula units for the periodic case, we can proceed to converting the geometry to a trilayer form. In this tutorial we make the active region two layers: layer 0 will be taken from the first unit of the preceding section (as will the left semi-infinite layers); layer 1 will be taken from the second unit of the preceding section (as will the right semi-infinite layers). lmpg requires that the active region be divided into principal layers. lmscell has the ability to make the conversion to a trilayer geometry, including partitioning the active region into principal layers. It requires as input PLATL and PLATR. They are already included in the ctrl file; note that are identical to the third axis in the original site file (site1.cspbi). The active region will consist of the same lattice vectors as in the periodic system, only now the third lattice vector does not repeat periodically. Instead, attached on the left are an infinite stack of repeating layers each with the third vector PLATL; similarly stacked are repeating layers on the right with the third vector PLATR. The assignment of principal layers proceeds as follows. 1. Recall that atoms in site1.cspbi are ordered by projection along the third axis. All sites with projection between 0 and the height corresponding to PLATL will form the one layer in the left semi-infinite region; similarly all sites between the last site and sites below it with projection PLATR form a layer in the right semi-infinite region. 2. Once the end layers are established (the are taken from the active region and shifted left and right by PLATL and PLATR respectively), lmscell will start from the left layer, form the first layer (called layer 0) by finding all sites with a connecting vector that touch any site in the left semi-infinite layer, the layer 1 by finding all sites that touch any site in layer 0, and so on. This is accomplished with the following instruction: $ lmscell -vfile=2 -vpgf=2 cspbi --stack~rsasa@fn=rsta2~setpl:sort:range=9:shorps=0,0,1~wsite:fn=site3:rsasa=rsta3 >> out.lmscell.cspbi3


Switch -vfile=2 tells lmscell (and later lmpg) to read structural information from site2.cspbi. In this tutorial, the active region has only these two layers, so the site2.cspbi can be used both for the periodic case and the trilayer geometry; only the third lattice vector has different meanings in two two contexts, as noted above.

Switch -pgf=2 gives token PGF_MODE the value 2. When it is nonzero, lmscell will try to read PLATL and PLATR.

Breaking down the syntax of the --stack command (see also the instruction manual):

• ~rsasa@fn=rsta2~ tells lmscell to read energy moments of the current (2 formula-unit) structure into memory. lmscell has not yet been told to switch to a trilayer geometry.

• ~setpl causes lmscell to convert the periodic geometry to a trilayer form, and assign a principal layer index to each site following the prescription noted above.
The modifiers serve the following functions:

• :sort re-orders sites by increasing height along the PL axis
• :range=9 is required. This supplies to lmscell the length of the longest connecting vector It is defined by the range of the structure constants, and you should pick a number consistent with this range.
• :shorps=0,0,1 places site in 1st quadrant on the third axis
• ~wsite:fn=site3:rsasa=rsta3 causes lmscell to write a new site file, and a moments file rsta3 which contains moments for all sites in the active region, and also sites in the left- and right- principal layers.

Compare site3.cspbi to site2.cspbi. They share in common the same site and lattice vectors, but site3.cspbi also includes information about principal layers.

Compare rsta3 to rsta2. rsta3 has moments for 64 sites, while rsta3 has moments for 32 sites. The left and right layers contain 16 sites each (this information was copied from the first and last layer in the active region)

#### Self-consistent density in the trilayer geometry

This section demonstrates that the density in the trilayer geometry is self-consistent. This should be the case since all the layers are identical, but it is a nontrivial check because lmpg does not make use of periodic boundary conditions on the third axis.

There are in fact three regions that can be checked: the left- and right- semi-infinite regions, and the active region.

The following step is not needed, but you might want to save the existing rsta file so you can repeat earlier steps without starting from the beginning:

$cp rsta.cspbi rsta2.cspbi  • copy rsta3 to the appropriately named file • make the potential from the moments in rsta.cspbi for all atoms in the trilayer geometry • make the structure constants for the trilayer geometry • run lmpg in mode 2, to converge the left- and right- regions to self-consistency. By construction, the potential should already be self-consistent: running lmpg confirms that it is. • run lmpg in mode 1, to converge the active region to self-consistency. Similarly in this case the potential should already be self-consistent. $ cp rsta3 rsta.cspbi
$lmpg -vgamma=t -vpgf=1 -vnit=0 -vfile=3 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho$ lmstr -vgamma=t -vpgf=2 -vnit=1 -vfile=3 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho
$lmpg -vgamma=t -vpgf=2 -vnit=1 -vfile=3 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho$ lmpg -vgamma=t -vpgf=1 -vnit=-1 -vfile=3 -vccor=f -vtwoc=t ctrl.cspbi --rs=1,0 --quit=rho


You should find that the end layers are converged to a precision of order 10−5, and the trilayer converged to a precision of order 10−4. Both of these are acceptable.

Inspect vshft.cspbi now. You should see something similar to the following:

 ef=-0.2380000  vconst=0.0152048 vconst(L)=-0.0017305 vconst(R)=-0.0017305


Potential shifts for the end layers are generated in the self-consistency cycle for the end layers (PGF_MODE=2), and shift for the active layer (vconst=) in the self-consistency cycle for the trilayer (PGF_MODE=1).

### 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).

Implemention 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).