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 CsPbI₃, 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.

Table of Contents

• Preliminaries
• Single Formula Unit
  • Self-consistent density through the band code
  • Self-consistent density through the Green’s function code
• Doubled Formula Unit
• Conversion from Crystal to Trilayer Geometry
  • Self-consistent density in the trilayer geometry
• Footnotes and references

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

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

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 CsPbI₃. (The cubic structure has higher symmetry, but this structure provides a reasonable representation of CsPbI₃ 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
         GFOPTS=p3;specfr;declead;padtol=1.001

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 --asa --gf --nk=8,8,8 --nl=3 --omax=0.16,0.21,0.25 --wsrmax=3.3 --rdsite --noshorten --nfile --wsite cspbi

$ diff actrl.cspbi ctrl.cspbi

< % 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}

< 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
>          GFOPTS=p3;specfr;declead;padtol=1.001

<   ATOM=Pb       Z= 82  R= 3.300000  LMX=2
---
>   ATOM=Pb       Z= 82  R= 3.300000  LMX=2 IDMOD=0,0,1

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

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

  $ grep Makidx out
  

  The Hamiltonian consists of a basis of 104 orbitals. There are 40 orbitals omitted because of the basis for empty sphere E1 was reduced (LMX=1).

• Note the command line argument -vfile=1. When the command-line switch blm ... --nfile ... is used, the following lines are inserted into the ctrl file:

  %ifdef file>=1
    file=   site{file}
  %endif
  

  If variable file is defined and its value exceeds unity, the preprocessed ctrl file picks up an extra line, which if file=1, reads file= site1 , thus causing structural data to read from this file.

• Variables ... -vccor=f -vtwoc=t affect tags OPTIONS_ASA_CCOR and OPTIONS_ASA_TWOC. They turn off lm’s “combined correction term”, and also reduce lm’s standardd 3-center hamiltonian with a simpler two-center one. Neither of these approximations is necessary (both reduce the accuracy of the bands) but in making them the hamiltonian generated by lm is formally identical to the Green’s function lmgf makes. lmgf has no analogous combined correction term, and while it does have a third order hamiltonian, the third order addition is slightly different from that of lm. Thus we employ these approximations to demonstrate that lm, lmgf, and lmpg are essentially equivalent.

• Switch --rs=0,1 causes lm to write a restart file rsta.cspbi (ASA format).

Self-consistent density through the Green’s function code

Starting from a self-consistent density generated by lm, we should find that the density generated by lmgf is the same. This is not trivial, since lmgf requires a contour integration that must encompass all the poles of the occupied valence bands. The integration must be performed on the real axis below the deepest pole to the Fermi level.

To confirm that lmgf indeed generates a self-consistent density, try

$ rm -f mixm.cspbi
$ 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⁻⁵ 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⁻⁵, and the trilayer converged to a precision of order 10⁻⁴. 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).