Questaal Home
Navigation

Nb/Ni/Nb(110) Metallic Trilayer, Landauer-Buttiker Transport and Andreev levels

Note: this tutorial is under construction.

This tutorial’s final objective is to compute the Andreev levels induced by a thin ferromagnetic Ni layer sandwiched between superconducting Nb leads. It uses electronic structure codes that make the Atomic Spheres Approximation. It uses a Landauer-Buttiker technique to compute the scattering matrix SN in the normal state, for a trilayer oriented along the [110] direction. Fron SN it computes the Andreev levels which appear when Nb makes a transition to the supercondcting state.

There is a related tutorial concerning Landauer-Buttiker transmission through a Nb/Fe/Nb(001) trilayer; the setup of the trilayer is a bit simpler to manage in this case since we handle the lattice mismatch between Nb and Ni in a less elaborate manner than was done for the Nb/Fe case.

This tutorial uses a relaxed Nb/Ni/Nb superlattice, where the atoms were relaxed in a constrained manner, allowing only the Ni and Nb frontier atoms to movel. In this tutorial the relaxations are taken as an input; how the relaxation is to be accomplished will be explained elsewhere.


Table of Contents

Sections build on prior sections, though any section can be treated independently by copying all the files indicated in the Setup instructions.

Section 1. Nb and Ni reoriented along the [110] direction and self-consistent ASA calculation

Creates files related to elemental Nb site4.nb, rsta4.nb, site2a.nb and site2b.nb, and files related to elemental Ni: site4.ni, rsta4.ni, site2a.ni and site2b.ni.

Requires files init.nb and init.ni. Setup:

touch init.nb; rm *.nb *.ni; cp ~/lmf/pgf/test/nbni/init.{nb,ni}  .

Create site4.nb, rsta4.nb

blm --nfile --nl=3 --nk=12,12,12 --asa nb > out.blm
cat actrl.nb | sed 's/LMX=2/LMX=2  GRP2=1/'  | awk '{ if ($1 == "OPTIONS") {print "HAM   PMIN=-1"}; print}' > ctrl.nb
lmchk ctrl.nb --wspec > out.lmchk.nb
lmstr nb >/dev/null
lm ctrl.nb -vnit=0 > /dev/null
lm -vmet=3 ctrl.nb > out.lm.nb
lmscell --rot~x:-pi/4 --plx~m~0,1,1,0,-1,1,2,1,1  nb --wsite~short~fn=site4 --sort:'x3 x2' > out.lmscell.nb
lmstr nb -vfile=4 >/dev/null
lm ctrl.nb -vmet=3 -vfile=4 -vnk1=12 -vnk2=8 -vnk3=8 --rs=0,1 > out.lm.nb
lmscell ctrl.nb -vfile=4 --stack~rsasa~wsite@nosite@rsasa=rsta4.nb >> out.lmscell.nb

Create site2a.nb and site2b.nb.

lmscell -vfile=4 ctrl.nb '--stack~rmsite@targ=3,4~pad@0,0,-1/sqrt(2)~wsite@fn=site2a' >> out.lmscell.nb
lmscell -vfile=4 ctrl.nb '--stack~rmsite@targ=1,2~addpos@dx=0,0,-1/sqrt(2)~pad@0,0,-1/sqrt(2)~wsite@fn=site2b' >> out.lmscell.nb

Create site4.ni, rsta4.ni

blm --nfile --nl=3 --nk=12,12,12 --asa --mag ni > out.blm.ni
cat actrl.ni | sed 's/\(MMOM.*\)/\1 GRP2=1/' | awk '{ if ($1 == "OPTIONS") {print "HAM   PMIN=-1"; print "BZ    INVIT=F"}; print}' > ctrl.ni
lmchk ctrl.ni --wspec > out.lmchk.ni
lmscell --rot~x:-pi/4 --plx~m~1,1,-1,2,-2,0,0,0,1  ni --wsite~short~fn=site4  --sort:'x3 x2' > out.lmscell.ni
lmstr ni -vfile=4 >/dev/null
lm ctrl.ni -vbeta=.1 -vmet=3 -vfile=4 -vnit=-30 -vnk1=12 -vnk2=8 -vnk3=17 --rs=0,1 > out.lm.ni
lmscell ctrl.ni -vfile=4 --stack~rsasa~wsite@nosite@rsasa=rsta4.ni >> out.lmscell.ni

Create site2a.ni and site2b.ni.

lmscell -vfile=4 ctrl.ni '--stack~rmsite@targ=3,4~pad@0,0,-1/2/sqrt(2)~wsite@fn=site2a' >> out.lmscell.ni
lmscell -vfile=4 ctrl.ni '--stack~rmsite@targ=1,2~addpos@dx=0,0,-1/2/sqrt(2)~pad@0,0,-1/2/sqrt(2)~wsite@fn=site2b' >> out.lmscell.ni

Section 2. ctrl and template site files for (Nb(6)/Ni(8)/Nb(6) superlattice, and ES insertion files

Make ctrl.nbni and site4.nbni, and ES site files esitee1.nbni and esitee2.nbni needed for the two frontier layers.

Requires files site8.fp, dposx.6+8+6, spec.nb and spec.ni. Setup:

cp ~/lmf/pgf/test/nbni/site8.fp site1.nbni; cp ~/lmf/pgf/test/nbni/dposx.6+8+6 dposx.nbni
cp ~/lmf/pgf/test/nbni/spec.{nb,ni} .

Make ctrl.nbni and site4.nbni

blm --nl=3 --nfile --ewald=1d-12 --rdsite=site1 --rdpos=dposx --rdspec~fn=spec.nb~fn=spec.ni --usefiler --nosort --shorten=no --nit=20 --nk=12,8,2 --mag --mixs=B4,b={beta},k=4 nbni --findes~rmin=.9~shorten=0,0,1~omax=0.12~sclwsr=20~spec~nspmx=4~nesmx=36 --omax=.18,.26,.26 --asa~pfree~rmaxs=7.2 >> out.blm
cp actrl.nbni ctrl.nbni
lmchk nbni -vfile=0 --mino~z --wsite~short~fn=site2 >> out.lmchk
blm --nfile --ewald=1d-12 --rdsite=site2 --rdspec~fn=ctrl --usefiler~targ=1:4 --nosort --shorten=no --nit=20 --nk=12,8,2 --mag --mixs=B4,b={beta},k=4 nbni --findes~rmin=.9~shorten=0,0,1~omax=0.12~sclwsr=20~spec~nspmx=4~nesmx=36 --omax=.18,.26,.26 --asa~pfree~rmaxs=7.2  --ctrl=ctrl >> out.blm
lmchk nbni -vfile=0 --mino~z --wsite~short~fn=site3 >> out.lmchk
blm --nfile --ewald=1d-12 --rdsite=site3 --rdspec~fn=ctrl --usefiler~targ=1,4 --nosort --shorten=no --nit=20 --nk=12,8,2 --mag --mixs=B4,b=\{beta},k=4 nbni --findes~rmin=1.0~shorten=0,0,1~omax=0.12~sclwsr=20~spec~1spec --omax=.19,.23,.26 --asa~pfree~rmaxs=7.2 '--pgf~platl=0,0,sqrt(2)/2~platr=0,0,sqrt(2)/2' --gf --ctrl=actrl --convp >> out.blm
cp actrl.nbni ctrl.nbni
sed -i.bak -e 's/LMX=1/LMX=2 IDMOD=0,0,1/' ctrl.nbni
lmchk nbni -vfile=0 --wsite~fn=site4 >> out.lmchk

Extract ES sites at interfaces 1 and 2 to make esitee1.nbni and esitee2.nbni

lmscell ctrl.nbni -vfile=4 --stack~sort@x3~rmsite@targ=1:5,7,12:30~wsite@fn=esitee1 >> out.lmscell
lmscell ctrl.nbni -vfile=4 --stack~sort@x3~rmsite@targ=1:19,24,26:30~wsite@fn=esitee2 >> out.lmscell

Section 3. Site files for even-Ni-plane superlattices and trilayers

Makes ASA site files site{n}.nbni for n=8,12,16,…32.

Requires files ctrl.6+8+6, site8.fp, dposx.6+8+6, spec.{nb,ni}, esitee[12].nbni, site4.ni Setup:

cp ~/lmf/pgf/test/nbni/site8.fp site1.nbni
cp ~/lmf/pgf/test/nbni/ctrl.6+8+6 ctrl.nbni
cp ~/lmf/pgf/test/nbni/dposx.6+8+6 dposx.nbni
cp ~/lmf/pgf/test/nbni/{spec.{nb,ni},esitee[12].nbni,site4.ni} .

Make first in series, site8.nbni

lmscell -vfile=1 '-vdy=1/2/sqrt(2)' ctrl.nbni --rdpos=dposx --stack~file@esitee2@insert=15@scale3=0~file@esitee1@insert=7@scale3=0~initpl@left=1:2@right=nbas-1:nbas~wsite@fn=site8 >> out.lmscell

Make remainder by induction (tcsh shell)

foreach j ( 8 12 16 20 24 28)
  @ i = $j + 4
  lmscell -vdy='1/2/sqrt(2)' -vfile=$j ctrl.nbni --stack~file@site4.ni@dpos=0,dy,'5.5*dy'@xinsert=16~setrlx@targ=16:20@ijk=0~initpl@left=1:2@right=nbas-1:nbas~show@planes~wsite@fn=site$i >> out.lmscell
end

Find and add principal layer indices (tcsh shell)

foreach j ( 8 12 16 20 24 28 32)
  lmscell -vpgf=1 -vfile=$j ctrl.nbni --stack~setpl:sort:unsort:withinit:rangea=1.1~wsite@fn=site$j >> out.lmscell
end

Section 4. Self-consistent ASA band calculation of the Nb(6)/Ni(20)/Nb(6) superlattice

Requires files ctrl.6+8+6 site20.nbni, nb.nb, ni.ni. Setup:

cp ~/lmf/pgf/test/nbni/ctrl.6+8+6 ctrl.nbni; cp ~/lmf/pgf/test/nbni/site20.nbni .
cp ~/lmf/pgf/test/nbni/nb.nb nb.nbni; cp ~/lmf/pgf/test/nbni/nb.nb nbx.nbni
cp ~/lmf/pgf/test/nbni/ni.ni ni.nbni; cp ~/lmf/pgf/test/nbni/ni.ni nix.nbni

Make psta.nbni

sed -i.bak -e 's/B4,b={beta},k=4/B8,b={beta},k=8/' ctrl.nbni
rm -f mixm.nbni psta.nbni sv.nbni
lmstr -vfile=20 nbni > /dev/null
lm -vfile=20 nbni -vnit=0 --rs=0,1 >> out.lm
lm -vfile=20 nbni -vnit=-1 --pr30,20 --rs=1,0 -vscr=11 -vccor=f >> out.lm

Make the superlattice self-consistent and remake psta with the self-consistent density

lm -vbeta=.1 -vconvc=1e-6 -vmet=3 -vfile=20 nbni -vccor=f -vnit=50 -vscr=2 --zerq~all~qout --pr40,21 --rs=0,1 >> out.lm
lm -vfile=20 nbni --rsedit~rs~save@fn=rsta20~q >> out.lm
lm -vfile=20 nbni -vnit=-1 --pr30,20 --rs=1,0 -vscr=11 -vccor=f >> out.lm

Section 5. Self-consistent Green’s function calculation of the Nb(6)/Ni(20)/Nb(6) trilayer

Requires files ctrl.6+8+6 site20.nbni, psta.20, rsta20.lm. Setup:

cp ~/lmf/pgf/test/nbni/ctrl.6+8+6 ctrl.nbni; cp ~/lmf/pgf/test/nbni/psta.20 psta.nbni
cp ~/lmf/pgf/test/nbni/{site20.nbni,rsta20.lm} .
sed -i.bak -e 's/B4,b={beta},k=4/B8,b={beta},k=8/' ctrl.nbni
touch nb.nbni; rm -f nb*.nbni ni*.nbni e*.nbni *bc*.nbni* mixm.nbni vshft.nbni sv.nbni

Convert rsta20.lm to trilayer form (file rsta.nbni) suitable for lmpg, and make initial potential parameters

rm -f out.lmpg
lmpg -vpgf=1 -vfile=20 ctrl.nbni -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rsedit~rs@3d@fn=rsta20.lm~save~q >> out.lmpg
lmpg -vpgf=1 -vfile=20 ctrl.nbni -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rs=1,0 >> out.lmpg

Make end layers and active region each self-consistent

lmstr -vpgf=1 -vfile=20 ctrl.nbni -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t > /dev/null
lmpg -vpgf=2 -vconvc=1e-6 -vfile=20 ctrl.nbni -vccor=f -vnit=30 -vscr=0 --pr40,21 -vgamma=t >> out.lmpg
lmpg -vpgf=1 -vconvc=1e-6 -vfile=20 ctrl.nbni -vccor=f -vnit=50 -vscr=2 --pr40,21 -vgamma=t --rs=0,1 >> out.lmpg

Section 6. Self-consistent layer Green’s function calculation of the Nb(6)/Ni(24)/Nb(6) trilayer

Requires files ctrl.6+8+6,site24.nbni,psta.20,rsta20.lmpg,rsta4.ni, vshft.20, batch20. Setup:

cp ~/lmf/pgf/test/nbni/ctrl.6+8+6 ctrl.nbni; cp ~/lmf/pgf/test/nbni/rsta20.lmpg rsta.nbni
cp ~/lmf/pgf/test/nbni/{site24.nbni,rsta4.ni,psta.20,batch20} .
cp ~/lmf/pgf/test/nbni/vshft.20 vshft.nbni
sed -i.bak -e 's/B4,b={beta},k=4/B8,b={beta},k=8/' ctrl.nbni

Estimate psta.nbni from psta.20

set nni=20
mcx -f12f12.6 -vnni=20 -vm=nni/2+10 '-vn=66+3*nni' psta.20 -a p0 p0 -sub n/2+1,n/2+12,1,n/2 -a i1 p0 -sub 1,n/2,n/2+1,n/2+12 -a i2 p0 -sub n/2+1,n/2+12,n/2+1,n -a i3 p0 -sub n/2+1,n,n/2+1,n/2+12 -a i4  p0 -sub n/2+1,n/2+12,n/2+1,n/2+12 -a i5 p0 -split p 1,n/2+1,n+1 1,n/2+1,n+1 p11 i2 p12 i1 i5 i3 p21 i4 p22 '-sarray[3,3]' > psta.24
cp psta.24 psta.nbni

Estimate rsta.nbni from rsta20.lmpg

lmscell -vpgf=1 ctrl.nbni -vfile=24 -vnb=file+18 --stack~batch=batch20 >> out.lmscell
cp rsta2 rsta.nbni

Further initialization

touch nb.nbni ; rm -f nb*.nbni ni*.nbni e*.nbni *bc*.nbni* mixm.nbni sv.nbni
lmstr -vpgf=1 -vfile=24 ctrl.nbni -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t > /dev/null

Make end layers and active region self-consistent

lmpg -vpgf=1 -vfile=24 ctrl.nbni -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rs=1,0 >> out.lmpg
lmpg -vpgf=2 -vconvc=1e-6 -vfile=24 ctrl.nbni -vccor=f -vnit=30 -vscr=0 --pr40,21 -vgamma=t >> out.lmpg
lmpg -vbeta=.1 -vpgf=1 -vconvc=1e-6 -vfile=24 ctrl.nbni -vccor=f -vnit=100 -vscr=2 --pr40,21 -vgamma=t --rs=0,1 >> out.lmpg

Section 7. Transmission with normal and scaled moments

Requires files ctrl.6+8+6,site24.nbni,psta.20,rsta20.lmpg,rsta4.ni, vshft.20, batch20. Setup:

cp ~/lmf/pgf/test/nbni/ctrl.6+8+6 ctrl.nbni; cp ~/lmf/pgf/test/nbni/rsta20.lmpg rsta.nbni
cp ~/lmf/pgf/test/nbni/{site24.nbni,rsta4.ni,psta.20,batch20} .
cp ~/lmf/pgf/test/nbni/vshft.20 vshft.nbni
sed -i.bak -e 's/B4,b={beta},k=4/B8,b={beta},k=8/' ctrl.nbni

Estimate psta.nbni from psta.20

set nni=20
mcx -f12f12.6 -vnni=20 -vm=nni/2+10 '-vn=66+3*nni' psta.20 -a p0 p0 -sub n/2+1,n/2+12,1,n/2 -a i1 p0 -sub 1,n/2,n/2+1,n/2+12 -a i2 p0 -sub n/2+1,n/2+12,n/2+1,n -a i3 p0 -sub n/2+1,n,n/2+1,n/2+12 -a i4  p0 -sub n/2+1,n/2+12,n/2+1,n/2+12 -a i5 p0 -split p 1,n/2+1,n+1 1,n/2+1,n+1 p11 i2 p12 i1 i5 i3 p21 i4 p22 '-sarray[3,3]' > psta.24
cp psta.24 psta.nbni

Confirm potential is self-consistent with with full magnetic moment

touch nb.nbni ; rm -f nb*.nbni ni*.nbni e*.nbni *bc*.nbni* mixm.nbni sv.nbni
lmstr -vpgf=1 -vfile=24 ctrl.nbni -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t > /dev/null
lmpg -vpgf=1 -vfile=24 ctrl.nbni -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rs=1,0 >> out.lmpg
lmpg -vpgf=2 -vconvc=1e-6 -vfile=24 ctrl.nbni -vccor=f -vnit=30 -vscr=0 --pr40,21 -vgamma=t >> out.lmpg
lmpg -vbeta=.1 -vpgf=1 -vconvc=1e-6 -vfile=24 ctrl.nbni -vccor=f -vnit=100 -vscr=2 --pr40,21 -vgamma=t --quit=rho >> out.lmpg

… not yet done: transmission for this case

Repeat test, scale magnetic moments by 1/2. Set up potential and find vshft.nbni.

lmpg -vpgf=1 -vfile=24 ctrl.nbni '--rsedit~rs~set all 1 zers~read~set all 2 zerq~add all 1 .5 1:nbas 1:nbas~save~q' >> out.lmpg2
lmpg -vpgf=1 -vfile=24 ctrl.nbni -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rs=1,0 >> out.lmpg2
lmpg -vbeta=.1 -vpgf=1 -vconvc=1e-6 -vfile=24 ctrl.nbni -vccor=f -vnit=100 -vscr=2 --pr40,21 -vgamma=t --quit=rho >> out.lmpg2
lmpg -vbeta=.1 -vpgf=1 -vconvc=1e-6 -vfile=24 ctrl.nbni -vccor=f -vnit=100 -vscr=2 --pr40,21 -vgamma=t --quit=rho >> out.lmpg2

… not yet done: transmission for this case

Preliminaries

The tutorial has several sections. Intermediate sections may be skipped: you can start the tutorial at any of Sections 1,2,3,4, provided you copy the necessary starting files from your Questaal source directory (where you cloned the Questaal package). For purposes of this tutorial we will assume it is ~/lmf. A condensed form of the entire tutorial has been bundled into a standard test case. You can invoke ~/lmf/pgf/test/test.nbni from any empty directory and let the script perform the steps for you. It’s an easy way to see what you should have been able to produce froom following this tutorial, should have produced, if you have problems.

This tutorial makes extensive use of many of the Questaal packages. including lmscell, lmplan, lm, lmgf, lmpg, and miscellaneous utilities such as vextract. They are assumed to be in your path.

Questaal invokes mpix to execute MPI jobs. mpix is a proxy for mpirun, the default command that launches MPI jobs. We assume your system has 16 processors on a core, but you can use any number, or run in serial.

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

Required Input files

This tutorial uses the following files that can be found in directory ~/lmf/pgf/test/nbni.

  1. init.ni : an init file for fcc Ni, with its equilibrium lattice constant a = ALAT = 6.643 a.u.

  2. init.nb : an init file for bcc Nb, with the lattice constant matched to Ni, ALAT=6.643 a.u. Nb is a larger atom than Ni, but because it forms in the bcc structure (2 atoms/a3) while Ni forms in the fcc (4 atoms/a3), Nb’s equilibrium lattice constant is a = ALAT = 6.226 a.u. – about 6.7% smaller than Ni’s. A more careful calculation would build a complex interface to account for the lattice mismatch, as was done in the Nb/Fe/Nb trilayer tutorial. Here we just stretch the Nb lattice constant to match that of Ni.

  3. site8.fp : a site file for an Nb(6)/ .. bunrelaxed structure, taken from the FP setup

  4. dposx.6+8+6 the displacements of frontier planes obtained by relaxing the structure using the lmf code. A separate tutorial will explain how to make this relaxation using lmf.

Introduction and Basic Strategy

This tutorial’s final objective is to compute the Andreev levels induced by a thin ferromagnetic Ni layer sandwiched between superconducting Nb leads. We will use electronic structure codes that make the Atomic Spheres Approximation: the the band programs lm and the principal layer Green’s function code lmpg, which can calculate Landauer-Buttiker transport; in this case the scattering matrix S.

The final result is reached through the following steps:

  1. The electronic structure of elemental Nb and Ni are computed in small supercells of 4 atoms oriented along the [110] direction. It is convenient to build stacks of Nb and Ni planes out of these supercells because the coordinate system becomes orthogonal. Self-consistency is accomplished with the ASA band code, lm.

  2. An Nb(6)/Ni(8)/Nb(6) superlattice is built from the cells of elemental Nb and Ni. Such stacks are unrealistic in pristine form because stacking the bcc and fcc crystal structures of Nb and Ni result in unrealistically small Nb-Ni bond lengths. To accomodate this, the superlattice has relaxed using the lmf code, where only the frontier Nb and Ni atoms are allowed to move. That relaxation is not performed here; instead the displacements obtained by the relaxing the frontier layers is taken as input.

    At the Nb/Ni interface, sphere packing become poor, even with relaxation of the atoms there. the ASA becomes problematic when the sphere packing is not good, as is also seen in the Nb/Fe/Nb trilayer. The solution is to add empty spheres at the interface between Nb and Ni where the bcc and fcc lattice mismatch gives rise to large interstitial spaces. This step of the tutorial explains how to add empty spheres to restore a reasonably close packed geometry at the Nb/Ni boundary.

  3. lm

  4. Finally a structure is set up for trilayer geometry (… Nb / Fe / Nb …), which is related to but distinct from the Nb/Fe superlattice, as it is not periodic along the stacking axis. Transmission calculations through the trilayer are performed using lmpg.

In the Nb/Fe superlattice tutorial a superlattice of Nb and Fe was constructed; structural information is contained in site4.nbfe. Starting from it this tutorial will do:

Note: that the two codes, lm and lmpg handle the boundary conditions differently: lm uses periodic boundary conditions, so the Nb(1)/Ni/Nb(2) unit cell is the same as an infinite sequence of Nb(1+2)/Ni layers. lmpg on the other hand does not have periodic boundary conditions except in the plane normal to the interface, but rather a semi-infinite sequence of Nb layers on the left, a semi-infinite sequence of Nb layers on the right, with an active region (mostly Ni in this case) sandwiched in the middle. See the description in the lmpg documentation. lm and lmpg can make use of the same site file. But for lm (when boundary conditions are periodic) the site file corresponds to a superlattice, while for lmpg it corresponds to a trilayer.

1. Self-consistent ASA-LDA energy bands for Nb and Ni with [110] rotated to the z axis

Since we will be working with superlattices and trilayers oriented along the [110] direction, it is convenient to rotate the crystal lattices of elemental Nb and ni, orienting the [110] direction along z. Doing this makes the stacking of Nb and Ni planes easier. A 90&‌deg; rotation around x orients the [110] becomes along z axis and the [110] along y. We will use that rotation in building up site files below. It is convenient to work in these orthogonal coordinates, but we pay a price in that the unit cell contains two atoms/plane, and moreover it has two planes, staggered in a AB pattern. Thus the basic unit cell of both Ni and Nb have 4 atoms, each with 2 planes. Smaller “half-cells” with one plane are built using the lmscell utility; these will form the basic building blocks of Nb/Ni/Nb stacks with any number of planes.

The instructions below copy the Nb init file from a test directory in the git repository (where all steps of this tutorial are packaged into a shell script), builds the input files bcc Nb, rotates the bcc lattice to a 4-atom cell as noted above.

Files generated for elemental Nb

Use the blm utility to make template ctrl and site files for elemental Nb with one atom/cell.

touch init.nb; rm *.nb *.ni
cp ~/lmf/pgf/test/nbni/init.{nb,ni}  .
blm --nfile --nl=3 --nk=12,12,12 --asa nb > out.blm
cat actrl.nb | sed 's/LMX=2/LMX=2  GRP2=1/'  | awk '{ if ($1 == "OPTIONS") {print "HAM   PMIN=-1"}; print}' > ctrl.nb

Command line arguments to blm are documented here. In particular --nfile tells blm to add lines to the ctrl file so a particular site file can be specified via a command-line argument, a feature we will use extensively here.

The last step modifies the template input file by adding tags HAM_PMIN and SPEC_GRP2. HAM_PMIN tells the ASA codes not to allow the linearization energy of any partial wave to float below the free-electronic value, by constraining the continuous principal quantum number P. SPEC_GRP2 tells it to average P and moments Q for different classes of the same species in the self-consistency cycle (see here, footnote 2). We use it because but the 4-atom cell loses some symmetry, so that different sites belong to different classes. and ths tag will force all the Nb to be equivalent. This isn’t necessary, and the effect is small, but we do it for consistency so the different nb classes are identical in this artificial supercell.

This next instruction creates a special file spec containing the SPEC category for Nb.

lmchk ctrl.nb --wspec > out.lmchk.nb

spec.nb is not needed in this step, but it will be used in Step 3, when blm is used to construct ctrl file for the trilayer.

Run lm to make the Nb density self-consistent

lmstr nb >/dev/null
lm ctrl.nb -vnit=0 > /dev/null
lm -vmet=3 ctrl.nb > out.lm.nb

The second line creates an initial atom file, nb.nb, with a (crude) starting guess for P and Q. The last step runs through the self-consistency cycle; after completion nb.nb contains the self-consistent values for P and Q.

Next, generate the the 4-atom supercell with z along [110] as noted at the beginning of this section.

lmscell --rot~x:-pi/4 --plx~m~0,1,1,0,-1,1,2,1,1  nb --wsite~short~fn=site4 --sort:'x3 x2' > out.lmscell.nb

The supercell utility in lmscell is documented on this page Make the 4-atom supercell self-consistent (noter that it reads nb.nb created earlier so that it is already essentially self-consistent. It is not exactly so, because of incomplete k convergence. The finite k mesh corresponds to a different grid in the supercell.

lmstr nb -vfile=4 >/dev/null
lm ctrl.nb -vmet=3 -vfile=4 -vnk1=12 -vnk2=8 -vnk3=8 --rs=0,1 > out.lm.nb

Switch -vfile=4 tells lmstr and lm to read structural data from site4.nb. This is because blm inserted extra lines into ctrl.nb (which you can see by inspecting the file).

The second step makes a self-consistent density and writes it to a restart file, which contains P and Q for each site. The next step rewrites the data in condensed form, in a file rsta4.nb.

lmscell ctrl.nb -vfile=4 --stack~rsasa~wsite@nosite@rsasa=rsta4.nb >> out.lmscell.nb

This will be used in Step 3 to assemble trial P and Q for superlattices from self-consistent densities of elements Nb and Ni.

Finally, for stacking planes Nb and Ni planes on top of each other, we need single planes “half cells” of the Nb. The first two atoms in site4.nb comprise the first plane and the second two the second. Create site2a.nb and site2b.nb, each containing a single plane.

lmscell -vfile=4 ctrl.nb '--stack~rmsite@targ=3,4~pad@0,0,-1/sqrt(2)~wsite@fn=site2a' >> out.lmscell.nb
lmscell -vfile=4 ctrl.nb '--stack~rmsite@targ=1,2~addpos@dx=0,0,-1/sqrt(2)~pad@0,0,-1/sqrt(2)~wsite@fn=site2b' >> out.lmscell.nb

See this page superlattice utility in lmscell. Sites are removed from teh 4-atom cell and the 3rd lattice vector is scaled by 1/2. Note that these are pseudo-site files, as essential information about the AB stacking is lost. Note also the positions in site2b.nb are shifted so that the plane passes through z=0.

Files generated for elemental Ni

The commands for Ni are very similar to those for Nb. A couple of extra tags are needed, notably a trial magnetic moment must be specified. Also the superlattice vectors (in units of multiples of the primitive vectors) are different for the fcc case.

blm --nfile --nl=3 --nk=12,12,12 --asa --mag ni > out.blm.ni
cat actrl.ni | sed 's/\(MMOM.*\)/\1 GRP2=1/' | awk '{ if ($1 == "OPTIONS") {print "HAM   PMIN=-1"; print "BZ    INVIT=F"}; print}' > ctrl.ni
lmchk ctrl.ni --wspec > out.lmchk.ni
lmscell --rot~x:-pi/4 --plx~m~1,1,-1,2,-2,0,0,0,1  ni --wsite~short~fn=site4  --sort:'x3 x2' > out.lmscell.ni
lmstr ni -vfile=4 >/dev/null
lm ctrl.ni -vbeta=.1 -vmet=3 -vfile=4 -vnit=-30 -vnk1=12 -vnk2=8 -vnk3=17 --rs=0,1 > out.lm.ni
lmscell ctrl.ni -vfile=4 --stack~rsasa~wsite@nosite@rsasa=rsta4.ni >> out.lmscell.ni

We also need the half-cells. Note that the spacing beween planes is half as large as for Nb.

lmscell -vfile=4 ctrl.ni '--stack~rmsite@targ=3,4~pad@0,0,-1/2/sqrt(2)~wsite@fn=site2a' >> out.lmscell.ni
lmscell -vfile=4 ctrl.ni '--stack~rmsite@targ=1,2~addpos@dx=0,0,-1/2/sqrt(2)~pad@0,0,-1/2/sqrt(2)~wsite@fn=site2b' >> out.lmscell.ni

2. Setup for for even-Ni-plane superlattices

Required inputs for step 2:

  • site8.fp: a site file for the unrelaxed structure, taken from the FP setup
  • dposx.6+8+6 : the displacements obtained by relaxing the structure using lmf
  • spec.{nb,ni} (these were generated in step 1)

In this section, a ctrl file for a 6(Nb)/8(Ni)/6(Nb) superlattice is made, i.e. one with 3 (Nb), 4(Ni), 3(Nb) layers, respectively. It also serves as a template to made 6(Nb)/4n(Ni)/6(Nb) superlattices and trilayers for any n. Note that we must treat superlattices with odd numbers of Ni planes differently from those with even numbers, because of the AB stacking of the Ni planes. The even numbered case is simpler and we will consider only that case in this tutorial. A presciption has also been worked out for the odd-numbered case; you can try it yourself starting from ~/lmf/pgf/test/nbni/site6.fp (the analog of ~/lmf/pgf/test/nbni/site8.fp used here.) Alternatively, you can run the script ~/lmf/pgf/test/test.nbni 2 3 4 and the steps will be performed for you.

We would like to build a ctrl file that is suitable for any Nb/Ni(4n)/Nb, modifying only the site file. Further we would like to build the site in a lego-like fashion, layer by layer, so that construction is automatic. We can do this (it is the reason why single-layer site files were made in step 1), but but there are several factors we must take into account.

  • We must incorporate relaxation the frontier atoms
  • We must add empty spheres there since frontier sphere radii must be reduced from their bulk value
  • Sphere radii must satisfy the constraints.

CHECK hyperlink

These steps make the procedure a bit complicated, but blm and lmchk have been designed to accomplish them in a semi-automated way. Copy the following files from the git repository:

cp ~/lmf/pgf/test/nbni/site8.fp ~/lmf/pgf/test/nbni/dposx.6+8+6 .

The first is a site file for an uneconstructed superlattice, where the site positions are ordered by increasing height. Inspect the file and note that it has 4 Nb atoms (2 layers) at both top and bottom ends, a pair of (Nbx, Nix) atoms on each frontier layer, and finally 4 Ni atoms in the middle. The frontier atoms are rendered as different species because they get squeezed by atoms across the frontier. We must allow their sphere radii to be smaller than the bulk, and can to that by defining them as different species. The sphere radii have essentially the same constraints as described in the Nb/Fe/Nb tutorial [check hyperlink]: overlaps cannot be too large even while the sum of sphere radii must reach the cell volume; and we want both end layers and Ni atoms in the middle region to be defined in the same way as elemental Nb and Ni, respectively, as these regions are bulk like.

The second contains the lattice relaxation of the frontier atoms for this particular cell. Note that only frontier atoms have any displacement. This is an approximation, but it is a reasonable one. In any case the real material will be highly non-ideal, which is very difficult to model, so there is little point in trying to allow more general relaxations. The main thing is that the frontier atoms are not pressed too close together or kept too far apart. The fact that the displacements are limited only to the interface greatly simplifies the adoption the lego-like, layer-by-layer construction, since each layer is essentially independent of the others.

The tricky part is identifying placements and radii of the fictitous empty spheres. blm does have the capacity to build ctrl files while find places to insert ES; but the algorithm is not that sophisticated. Some iteration is needed, using lmchk to improve on blm’s guess at ES placement. We follow this prescription:

  1. Use blm to make a first cut at a ctrl file, with sphere radii chosen for the relaxed superlattice and ES added. Now we use the extension nbni
  2. Use lmchk with --mino~z to improve on ES placements by minimize overlap between ES and other sites
  3. Use blm a second time to make a ctrl file, with sphere radii chosen for the improved ES positions.
  4. Use blm a fnal time, with slightly relaxed constraints to fulfill the constraint that sphere volumes sum to the cell volume.
    This provides a working ctrl file for the 6,8,6 superlattice, and by simple extension 6,4n,6 superlattices and trilayers.
  5. The lm
blm --nl=3 --nfile --ewald=1d-12 --rdsite=site1 --rdpos=dposx --rdspec~fn=spec.nb~fn=spec.ni --usefiler --nosort --shorten=no --nit=20 --nk=12,8,2 --mag --mixs=B4,b={beta},k=4 nbni --findes~rmin=.9~shorten=0,0,1~omax=0.12~sclwsr=20~spec~nspmx=4~nesmx=36 --omax=.18,.26,.26 --asa~pfree~rmaxs=7.2 >> out.blm
cp actrl.nbni ctrl.nbni
lmchk nbni -vfile=0 --mino~z --wsite~short~fn=site2 >> out.lmchk
blm --nfile --ewald=1d-12 --rdsite=site2 --rdspec~fn=ctrl --usefiler~targ=1:4 --nosort --shorten=no --nit=20 --nk=12,8,2 --mag --mixs=B4,b={beta},k=4 nbni --findes~rmin=.9~shorten=0,0,1~omax=0.12~sclwsr=20~spec~nspmx=4~nesmx=36 --omax=.18,.26,.26 --asa~pfree~rmaxs=7.2  --ctrl=ctrl >> out.blm
lmchk nbni -vfile=0 --mino~z --wsite~short~fn=site3 >> out.lmchk
blm --nfile --ewald=1d-12 --rdsite=site3 --rdspec~fn=ctrl --usefiler~targ=1,4 --nosort --shorten=no --nit=20 --nk=12,8,2 --mag --mixs=B4,b=\{beta},k=4 nbni --findes~rmin=1.0~shorten=0,0,1~omax=0.12~sclwsr=20~spec~1spec --omax=.19,.23,.26 --asa~pfree~rmaxs=7.2 '--pgf~platl=0,0,sqrt(2)/2~platr=0,0,sqrt(2)/2' --gf --ctrl=actrl --convp >> out.blm
cp actrl.nbni ctrl.nbni
lmchk nbni -vfile=0 --wsite~fn=site4 >> out.lmchk

The last instruction does merely rewrites the site file in a long format, with the name site4.nbni. It is not necesssary here, but will be later when we use the same file for lmpg, where the principal layers will need to be specified.

The following step is also not necessary for the present purposes, but in preparation for lmpg we modify the template by setting LMX on the E sites to 2, and freeze the log derivative parameter on the ES d orbital.

sed -i.bak -e 's/LMX=1/LMX=2 IDMOD=0,0,1/' ctrl.nbni

Finally, to prepare for making the site files for a general Ni(4n) structure, we extract the ES sites in the left and right frontier layers and park them in respective site files. In that way they can be seamlessly folded into a superlattice site file.

lmscell ctrl.nbni -vfile=4 --stack~sort@x3~rmsite@targ=1:5,7,12:30~wsite@fn=esitee1 >> out.lmscell
lmscell ctrl.nbni -vfile=4 --stack~sort@x3~rmsite@targ=1:19,24,26:30~wsite@fn=esitee2 >> out.lmscell

3. Site files for even-plane superlattice of variable length

Required inputs to be copied from ~/lmf/pgf/test/nbni/dposx.6+8+6, or generated by a prior step

  • ctrl.6+8+6 Either use ctrl.nbni from step 2 or copy this file into ctrl.nbni.
  • site8.fp site file for the unrelaxed structure
  • dposx.6+8+6 the displacements obtained by relaxing the structure using lmf
  • spec.{nb,ni} ASA species categories for Nb and Ni (step 1)
  • esitee[12] site files at the ES first and second frontier planes (see step 2)
  • site4.ni site files of a 4-atom Ni supercell (see step 1)

The ASA is more difficult to set up because of the constraint that the spheres fill space: the sum-of-sphere volumes must equal the volume of the unit cell. Also In the end regions we want to keep the spheres the same as the bulk Nb. At the interface there is considerable reconstruction which makes good sphere packing difficult.

We require the structure generated by the Nb/Fe superlattice tutorial. If you did that tutorial, copy the structure file of the relaxed interface, site4.nbfe, to file site9.asa to your working directory.

Alternatively you can begin immediately starting with Section 3 by doing the following:

touch site9.asa; rm -f *.asa ; cp /h/ms4/lmf/b/./testing/nbfe/{site9.asa,rsta.nbasa,rsta.feasa} .

Note: a script has been written to perform all these steps automatically (part of the standard Questaal tests). The steps are a bit involved, and the script is useful as a sanity check to check that you’ve done all the steps in the proper order. To perform all the steps of this section automatically, do the following:

~/lmf/testing/test.nbfeasa --quiet 1

A working ctrl file is needed to use Questaal’s tools. The ctrl file in the Nb/Fe superlattice tutorial was designed for lmf, and it is a bit different from the ASA.

We make an temporary ctrl file, taking minimum information from that tutorial and build ctrl.asa from scratch. Make a preliminary file with

blm --nfile --rdsite=site9 --nosort --shorten=no --ctrl=ctrl asa

To manage sphere packing in the ASA we must distinguish sites immediately at the interface from those in the bulk. Recall the superlattice geometry:

  Layer    Nb 1      Nb 2      Fe 1       Fe 2     Fe 3      Nb 3      Nb 4
  sites    1-3       4-6       7-10       11-14    15-18     19-21     22-24

The frontier sites are 4-6 and 19-21 (Nb) and 7-10 and 15-18 (Fe). Use the lmscell, utility to assign species Nbx and Fex to these files, or alternatively just use a text editor and do it by hand. The following creates site8.asa

lmscell ctrl.asa --stack~setlbl@spid=Nbx@targ=4:6,19:21~setlbl@spid=Fex@targ=7:10,15:18~wsite@short@fn=site8

site8.asa differs from site9.asa only in that some labels changed.

To obtain a site file with adequate sphere packing, empty spheres will need to be added. This is not trivial. The automatic empty sphere finder in the blm utility greatly facilitates this step, and even with blm it takes some tinkering. Invoking blm with this set of switches was found to work:

blm --pr35 --nl=3 --gf --nfile --rdsite=site8 --nosort --shorten=no --nk=9,4,2 --mag --asa~pfree~rmaxs=7.2 --wsite@fn=site7 asa --findes~rmin=.9~shorten=0,0,1~omax=0~sclwsr=20~spec~nspmx=4~nesmx=36 --omax=.18,.26,.26
cp actrl.asa ctrl.asa

(This page documents switches for blm.) In construction of the ctrl and site files, attention needs to be paid to the following:

  1. The total sphere volume equals the unit cell volume. Confirm that the autogenerated files meet this requirement:
    lmchk ctrl.asa -vfile=7 | grep 'Sum of sphere'
    
  2. No sphere overlap is too large. –omax=.18,.26,.26 imposes constraints on the atom-atom, atom-ES, and ES-ES overlaps, respectively. To see what these overlaps actually are, do:
    lmchk ctrl.asa -vfile=7 | grep 'max ovlp'
    
  3. The Nb sphere radius should be the same as the bulk value. This isn’t strictly necessary, but if this condition is not met, the density inside the Nb spheres will not be exactly the same as the bulk, so we should impose that radius as a constraint. If the Nb end spheres are shrunk, the extra volume will have to come from the other atoms. lmchk can do this; by telling it to resize the spheres (tag SPEC_SCLWSR=11), while also constraining the Nb sphere radius to that of bulk Nb, R=3.065511 (see ctrl.nbasa in Section 1). Temporarily the adjustments ctrl.asa, e.g.
    sed -i.bak -e 's/R= 3.[0-9]*/R= 3.065511 CSTRMX=1/' -e 's/# SCLWSR/  SCLWSR/' ctrl.asa
    rm -f ctrl.asa.bak
    lmchk ctrl.asa -vfile=7 --wspec       ← --wspec causes lmchk to write a template SPEC category to file spec
    

    You should see the following table:

    spec  name        old rmax    new rmax     ratio
      1   Nb          3.065511    3.065511    1.000000
      2   Nbx         2.741494    2.741494    1.000000
      3   Fex         2.500719    2.500719    1.000000
      4   Fe          2.720097    2.720097    1.000000
      5   E           1.629572    1.828081    1.121817
      6   E1          1.495031    1.677150    1.121817
      7   E2          1.448620    1.625086    1.121817
      8   E3          1.053943    1.182331    1.121817
    

    The E spheres were resized, but the others were not. You can use your text editor to substitute sphere radii in the table above for the original ones. Or simply rerun blm with –rdspec, which uses the SPEC template written by lmchk, and –usefiler, to tells blm to not attempt to resize any spheres:

    blm --pr35 --nl=3 --gf --nfile --rdsite=site7 --rdspec --usefiler --nosort --shorten=no --nk=9,4,2 --mag --asa~pfree~rmaxs=7.2 --wsite@fn=null asa --omax=.18,.26,.26
    rm null.asa
    

    Note: --gf isn’t needed here but it is included to anticipate its use in Section 5.

    Finally lmchk has an algorithm to shift locations of E spheres in order to reduce their overlap with neighboring sites. Do this and create a new site6.asa

    lmchk ctrl.asa -vfile=7 --mino~z --wsite@short@fn=site6   ← --mino~z triggers the E position relocator
    

    Check once more the total sphere volume and the overlaps

    lmchk ctrl.asa -vfile=6 | grep 'Sum of sphere'
    lmchk ctrl.asa -vfile=6 | grep 'max ovlp'
    

    You can confirm that the maximum overlap between real atoms remains at 18%, while the overlaps with E sites increased. This the best we can do without adding more very small E spheres and it is not necessary.

  4. This last step is not essential, but makes matters simpler, if the sites are ordered by increasing projection along the line normal to the interface. lmscell can do this automatically:
    lmscell -vfile=6 ctrl.asa --stack~sort@h3~initpl@left=1,2,3@right=58,59,60~wsite@fn=site5
    

    Note that h3 evaluates to the projection of a given site normal to the plane given by PLAT1×PLAT2. The rearranged sites in site5.asa clarifies where the E spheres actually go.

    Note: initpl is not needed here, but it will be needed for the trilayer geometry, Section 6. The superlattice (stack) editor is documented here.

    This completes development of the structure file. Rename site5.asa:

    cp site5.asa site.asa
    

    As a sanity check, verify that the final site file matches the one in the repository:

    mcx -cmpf~fn1=site5.asa~fn2=`echo ~/lmf/testing/nbfe/site5.asa`
    

The remainder of this page should be ignored …

4. Self-Consistent ASA calculation of the Nb/Fe superlattice

The ASA band program lm. is used to make a self-consistent calculation of the Nb/Fe superlattice, starting from the structure in site5.asa generated in the previous section, and restart files for elemental Nb and Fe.

Files site5.asa and rsta.nbasa and rsta.feasa are created in Sections 3, 1, and 2 respectively. If you want to proceed directly without working through preceding sections, do the following

~/lmf/testing/test.nbfeasa --quiet 1

Trial density

You can use lm’s restart editor to set up a trial density from P, Q taken from elemental Nb and elemental Fe for their respective nuclei, and no charge on E sites. (boundary conditions P and energy moments Q are sufficient to completely determine the density).

Do the following:
a. Use lm to create an initial density. The restart editor copies rsta files for bulk Nb and Fe into the appropriate slots of the superlattice
b. Use restart file editor in lm to install moments into Nb and Fe sites of the superlattice (output = rsta4)

cp site5.asa site.asa                 ← site5 was generated in the preceding section
touch nb.asa; rm {rsta,nb,fe,e}*.asa  ← Make a fresh start
lm asa -vnit=0 --pr30,20 --rs=0,1     ← Moments for a "free sphere" (crude starting guess for density)
cp rsta.nbasa rsta2.asa               ← lm reads this file with P,Q for elemental Nb and copies to all Nb in SL
cp rsta.feasa rsta3.asa               ← lm reads this file with P,Q for elemental Fe and copies to all Fe in SL
lm asa --rsedit~rs~rs@ib=1@src=1:3@fn=rsta2~rs@ib=14@src=1:3@fn=rsta2~~rs@ib=44@src=1:3@fn=rsta2~rs@ib=58@src=1:3@fn=rsta2~rs@ib=27@src=1:4@fn=rsta3~rs@ib=31@src=1:4@fn=rsta3~rs@ib=35@src=1:4@fn=rsta3~save@fn=rsta4~q
cp rsta4.asa rsta.asa                 ← This restart makes a good initial guess for a starting density

You could start use a density generated from the 3rd line alone. You should be able to converge to self-consistency from such a trial density, but the starting guess is a pretty crude one.

«–rsedit» reads Nb P, Q from rsta2.asa and Fe P, Q from rsta3.asa to create a new rsta.asa. To know what sites correspond to Nb and to Fe, do

lmchk asa | egrep Class\|Nb | head -20
lmchk asa | egrep Class\|Fe | head -20

Sites 1:3 are the leftmost Nb layer, 14:16 the next. On the right are sites 44:46 and 58:60. Sites 27:30, 31:34 and 35:38 are the three Fe layers.

The restart file editor is self-documenting. To see its instruction set, invoke lm --rsedit and type ? <enter>

Most the editors can be run interactively, or in batch mode. For documentation of this editor’s instructions, see the Appendix.

If character following «–rsedit» is non-blank the editor starts in batch mode using that character as delimiter («~» in this case). Thus the batch commands are:

  • rs
    Reads rsta.asa into the current density

  • rs@ib=1@src=1:3@fn=rsta2
    Reads P, Q from sites 1:3 in file rsta2.asa and copies them to sites 1,2,3 in the current density.

  • rs@ib=14@src=1:3@fn=rsta2~rs@ib=44@src=1:3@fn=rsta2~rs@ib=58@src=1:3@fn=rsta2
    Reads P, Q from sites 1:3 in file rsta2.asa and copies them to sites 14,15,16 in the current density. The process is repeated for sites 44,45,46 and sites 58,59,60.

  • rs@ib=27@src=1:4@fn=rsta3~rs@ib=31@src=1:4@fn=rsta3~rs@ib=35@src=1:4@fn=rsta3
    Reads P, Q from sites 1:4 in file rsta3.asa and copies them to sites 27:30, 31:34 and 35:38, respectively, in the current density.

  • save@fn=rsta4
    Saves the current density in file rsta4.asa

  • q
    Exits the editor

Self-consistent density with in the ASA band program lm

Because of the long unit cell, care must be taken to converge to self-consistency. If the input and output density are mixed in a simple, linear manner, i.e. , convergence to self-consistency is very slow. (Mixing parameter  b={beta}) in the ctrl file corresponds to .) This is for two reasons. First, the density-of states of this system is very high, allowing density changes by many small, low-energy fluctuations in state occupations near the Fermi level. Second, and very important, the cell is long and it is necessary to damp out against low-energy, long-wave charge fluctuations along the long axis. That can readily be seen from Poisson’s equation of a periodic cell in Fourier space:

Away from self-consistency there will be a deviation . This generates a change in the potential . Long cells have vectors with small G, and changes to the potential from get greatly amplified. The solution is to screen out these low Fourier components with the static dielectric function  ε(q=0). This quantity it is not given by the standard self-consistency procedure; however a model ε may be used. The full potential code lmf uses the Lindhard function, but the ASA has no plane-wave representation of the density, so constructing a model dielectric function is not so simple. A model construction is nevertheless implemented (it acts when OPTIONS_SCR=4), and it protects against these fluctuations, but at the same time it conceals true physical changes in the moments, and convergence is still slow.

It is better to calculate ε(q=0) explicitly in the ASA. lm (and lmgf) have this ability in an ASA sense; it can make a discretized form, expressed as a matrix that reflects the induced charge in channel from a disturbance in channel . Setting OPTIONS_SCR=11 turns on this mode; it causes a file psta.ext to be created. With lm, execution stops after this file is written. Once written lm will read this file to screen before mixing with . Perform these steps with the following:

lmstr asa                                                             &larr; Structure constants
touch nb.asa; rm -f {nb,fe,e}*.asa save.asa mixm.asa sv.asa log.asa   &larr; Fresh start
lm asa -vnit=-1 --pr30,20 --rs=1,0 -vscr=11 > out.psta

The last step makes the ASA polarizability, psta.asa. Note that it uses -vnit=-1 --rs=1,0 which causes lm to read rsta.asa and to create the sphere potentials. It does this because ctrl.asa contains this line:

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

Since «nit» is negative, BEGMOM evaluates to nonzero, which tells lm to start with P, Q and make potential parameters. The absolute value of «nit» tells lm how many cycles to run (observe how «nit» is used in ctrl.asa).

-vscr=11 causes token SCR to evaluate to 11, which tells lm to generate psta.asa. The 10s digit creates an on-site contribution, which is calculated from the second derivative of the total energy with respect to the moments (something like U).

With psta in hand, we can converge the density to self-consistency with ease:

mpix -n 16 lm -vmet=3 ctrl.asa -vccor=f -vnit=50 -vscr=2 --zerq~all~qout --pr40,21 --rs=0,1 > out.lm
lm asa --rsedit~rs~save@nopot@fn=rsta5~q
mv rsta5.asa rsta.asa.lm

The final two steps are not needed, but they preserve the moments in a compact form so a self-consistent density can be easily reconstructed.

To see the effect of the screening, look at the first time the screened density is constructed. You should see something like

 SCRMOM: unscreened dv=3.62e2  screened dv=1.06e0
 SCRMOM: unscreened dq=7.37e-1  screened dq=2.40e-1
 ...
 a table follows

You can see that screening strongly damps out the change in potential (dv), while the change in density is much less so.

To see how lm converges with iteration, do:

grep DQ out.lm4

and you should see

 PQMIX:  read 0 iter from file mixm.  RMS DQ=6.43e-2
 PQMIX:  read 1 iter from file mixm.  RMS DQ=4.75e-2  last it=6.43e-2
 PQMIX:  read 2 iter from file mixm.  RMS DQ=3.48e-2  last it=4.75e-2
 PQMIX:  read 3 iter from file mixm.  RMS DQ=8.38e-3  last it=3.48e-2
 PQMIX:  read 4 iter from file mixm.  RMS DQ=6.68e-3  last it=8.38e-3
 PQMIX:  read 5 iter from file mixm.  RMS DQ=5.71e-3  last it=6.68e-3
 PQMIX:  read 6 iter from file mixm.  RMS DQ=4.06e-3  last it=5.71e-3
 PQMIX:  read 0 iter from file mixm.  RMS DQ=1.96e-3  last it=4.06e-3
 PQMIX:  read 1 iter from file mixm.  RMS DQ=1.52e-3  last it=1.96e-3
 PQMIX:  read 2 iter from file mixm.  RMS DQ=1.24e-3  last it=1.52e-3
 PQMIX:  read 3 iter from file mixm.  RMS DQ=3.25e-4  last it=1.24e-3
 PQMIX:  read 4 iter from file mixm.  RMS DQ=1.84e-4  last it=3.25e-4
 PQMIX:  read 5 iter from file mixm.  RMS DQ=1.16e-4  last it=1.84e-4
 PQMIX:  read 6 iter from file mixm.  RMS DQ=8.40e-5  last it=1.16e-4
 PQMIX:  read 0 iter from file mixm.  RMS DQ=4.18e-5  last it=8.40e-5
 PQMIX:  read 1 iter from file mixm.  RMS DQ=3.12e-5  last it=4.18e-5
 PQMIX:  read 2 iter from file mixm.  RMS DQ=2.43e-5  last it=3.12e-5
 Jolly good show !  You converged to rms DQ=  0.000055

It is instructive to see what self-consistency does. Output below the line GETZV near the bottom shows the following:

  1. Inspect the Nb atoms first (last iteration)
    grep  'ATOM=Nb' out.lm | tail -12
    

    You should see

    ATOM=Nb   Z=41  Qc=36  R=3.065511  Qv=0.0007  mom=0.00258  a=0.025  nr=443
    ATOM=Nbx  Z=41  Qc=36  R=2.741494  Qv=-1.147863  mom=-0.19055  a=0.025  nr=435
    ATOM=Nb2  Z=41  Qc=36  R=3.065511  Qv=-0.090433  mom=-0.01461  a=0.025  nr=443
    ATOM=Nbx2 Z=41  Qc=36  R=2.741494  Qv=-1.396411  mom=-0.20131  a=0.025  nr=435
    ATOM=Nb3  Z=41  Qc=36  R=3.065511  Qv=0.08968  mom=0.00535  a=0.025  nr=443
    ATOM=Nbx3 Z=41  Qc=36  R=2.741494  Qv=-1.467382  mom=-0.18838  a=0.025  nr=435
    ATOM=Nb4  Z=41  Qc=36  R=3.065511  Qv=-0.108543  mom=0.00441  a=0.025  nr=443
    ATOM=Nbx4 Z=41  Qc=36  R=2.741494  Qv=-1.16373  mom=-0.1876  a=0.025  nr=435
    ATOM=Nb5  Z=41  Qc=36  R=3.065511  Qv=-0.002989  mom=0.0051  a=0.025  nr=443
    ATOM=Nbx5 Z=41  Qc=36  R=2.741494  Qv=-1.16258  mom=-0.16304  a=0.025  nr=435
    ATOM=Nb6  Z=41  Qc=36  R=3.065511  Qv=-0.125127  mom=-0.00785  a=0.025  nr=443
    ATOM=Nbx6 Z=41  Qc=36  R=2.741494  Qv=-1.221631  mom=-0.16872  a=0.025  nr=435
    

    “Frontier” Nb atom (labelled Nbx) show deficiencies of charge (Qv), indicating a charge transfer to the Fe. This is not entirely accurate; the Nb spehres are smaller and some of the charge resides in the E spheres). Also they have a small induced negative magnetic moment. “Bulk” Nb atoms — those one layer removed from the interface are already almost bulk-like, with Qv close to zero. This is a reflection of the very short screening length in a metal.

  2. Inspect the Fe atoms next (last iteration)
    grep  'ATOM=Fe' out.lm | tail -12
    
    ATOM=Fex  Z=26  Qc=18  R=2.500719  Qv=-0.364473  mom=2.21537  a=0.025  nr=391
    ATOM=Fe   Z=26  Qc=18  R=2.720097  Qv=0.278605  mom=2.35703  a=0.025  nr=399
    ATOM=Fex2 Z=26  Qc=18  R=2.500719  Qv=-0.256053  mom=1.92812  a=0.025  nr=391
    ATOM=Fe2  Z=26  Qc=18  R=2.720097  Qv=0.279149  mom=2.27022  a=0.025  nr=399
    ATOM=Fex3 Z=26  Qc=18  R=2.500719  Qv=-0.462283  mom=1.98707  a=0.025  nr=391
    ATOM=Fe3  Z=26  Qc=18  R=2.720097  Qv=0.234948  mom=2.41888  a=0.025  nr=399
    ATOM=Fex4 Z=26  Qc=18  R=2.500719  Qv=-0.34065  mom=1.63231  a=0.025  nr=391
    ATOM=Fe4  Z=26  Qc=18  R=2.720097  Qv=0.292501  mom=2.32164  a=0.025  nr=399
    ATOM=Fex5 Z=26  Qc=18  R=2.500719  Qv=-0.175108  mom=2.04966  a=0.025  nr=391
    ATOM=Fex6 Z=26  Qc=18  R=2.500719  Qv=-0.268051  mom=1.98809  a=0.025  nr=391
    ATOM=Fex7 Z=26  Qc=18  R=2.500719  Qv=-0.252931  mom=1.94343  a=0.025  nr=391
    ATOM=Fex8 Z=26  Qc=18  R=2.500719  Qv=-0.159016  mom=1.97154  a=0.025  nr=391
    

    There are only three Fe layers. The total self-consistent magnetic moment is about 23.8 μB, just shy 2 μB per Fe atom (inspect the last line of save.asa). This is a bit less than that of bulk Fe, which is 2.28 μB per Fe atom in this approximation. It shows that the Fe moment isn’t quite rigid, and gets extinguished somewhat at the Fe/Nb interface. The table shows that the moments are quenched on average by about 10 % on the frontier layers, but mostly bulk-like in the middle layer.

In summary, there is significant charge transfer at the Nb/Fe interface (as expected since work functions are different and Fermi levels align) but most of the disruption to the potential occurs on the frontier layers. Also a reduction in the local Fe moment is found there.

ASA Energy bands

Make the energy bands in a manner similar to the full LDA. As a first step, extract the a list of orbitals in the Hamiltonian associated with Nb, and ones with Fe, so as to color the bands according to the Nb or Fe character.

lm asa --pr50 --quit=ham | grep -A 70 'Orbital positions'

Inspect the the table. We can associate orbitals 1:134,243:360 with Nb, and orbitals 135:242 with Fe.

Generate the energy bands

lmchk ctrl.asa  --syml~n=41~q=0,1,0,0,0,0,1,0,-1
lm -vmet=3 ctrl.asa -vccor=f -vnit=1 --pr40,21 --quit=rho --band~col=1:134,243:360~col2=135:242~fn=syml

and create and view postscript figure for the majority and minority bands

echo -5,5,5,10 | plbnds -spin1 -fplot~sh~scl=.7~shft=-.28 -ef=0 -scl=13.6 -dat=up -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0 bnds.asa
open fplot.ps
echo -5,5,5,10 | plbnds -spin2 -fplot~sh~scl=.7~shft=.45 -ef=0 -scl=13.6 -dat=dn -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0 bnds.asa
open fplot.ps

Majority (left) and minority (right) energy bands of the Nb/Fe superlattice in the ASA approximation.   Red: Nb character; green: Fe character

ASA energy bands of the Fe/Nb interface

Majority (left) and minority (right) energy bands of the Nb/Fe superlattice from the full LDA poential.

LDA bands of the Nb/Fe superlattice

Compare to the full LDA band structure. You can see that they are very similar.

5. Crystal Green’s function calculation of the Nb/Fe superlattice

The ASA suite also has a crystal Green’s function code, lmgf. This page has a basic tutorial describing lmgf and its operation. It is of particular interest here mainly as a segue into the principal layer Green’s function code, lmpg.

This section should follow the preceding ones. To run it standalone, the following setup should be sufficient:

rm *
~/lmf/testing/test.nbfeasa --quiet 1 2
lmstr asa
cp ~/lmf/testing/nbfe/rsta.asa.lm .

Note: the ctrl file is already set up for lmgf because the --gf was included in its construction, Section 3.

Note: steps of this test have been bundled into test.nbfeasa To perform these steps automatically, do

~/lmf/testing/test.nbfeasa [--MPIK] 4

lmgf operates in a manner approximately similar to lm but it generates Green’s functions instead of wave functions. The program flow of the two are similar; only the “crystal” branch is replaced by a branch that obtains multipole moments through Green’s functions. Multipole moments provide a compact way to represent the charge density in the ASA: the ability to generate them enables lmgf to perform self-consistent calculations.

If the Hamiltonians generated by lm and lmgf were identical, they should produce identical charge densities. This is only approximately true: the third order parameterization of the potential function P are not identical; see Sec 2.13.4. in Ref. 1 for details.

As noted in the Introductory tutorial, lmgf requires some additional input, namely a GF category and an EMESH tag in the BZ category. It is strongly advised that you run through that tutorial before doing this section.

The elliptical energy contour, normally used for Brillouin zone integrations needed, e.g. to find the Fermi level and the output moments, requires information about the Fermi level. It is easiest to get an estimate from lm, since the two methods are closely linked. From the last output of lm, the Fermi level is displayed in a table with a line

 BZINTS: Fermi energy:      0.102488; 156.000000 electrons;  D(Ef):  337.606
 ...

The Fermi level appropriate to lmgf should be close to 0.102 Ry, but not identically so because lmgf uses a different integration scheme and the third order Green’s function is close to, but identical with the 3rd order ASA hamiltonian. Before running lmgf we can make a new psta file that corresponds to the dielectric function of the self-consisten density instead of the trial one

cp rsta.asa.lm rsta.asa
rm -f save.asa mixm.asa sv.asa log.asa vshft.asa
lm asa -vnit=-1 --pr30,20 --rs=1,0 -vscr=11 > out.psta2

Run lmgf to get the Fermi level satisfying charge neutrality. This is a nonlinear problem, so EF must be found iteratively (see also the Introductory tutorial). Here we choose EF to be zero (specified in the input file) but add a constant potential shift in the entire system to satisfy charge neutrality.

rm -f save.asa mixm.asa sv.asa log.asa vshft.asa
mpix -n 16 lmgf ctrl.asa -vccor=f -vnit=1 -vscr=2 --pr40,21 -vgamma=t --quit=rho > out.lmgf

–quit=rho tells lmgf to stop after EF has been found. The output should show a series of tables, the first being:

 GFASA:  integrated quantities to efermi = 0
 ...
         deviation from charge neutrality: -54.032927

This shows that EF needs to be shifted to accommodate 54 electrons. It uses a fast method (a Pade approximation to ) to find an estimate for the shift. It iterates the Pade approximation until charge neutrality is found. It settles on a converged value

 gfasa:  potential shift this iter = -0.10527.  Cumulative shift = -0.10527
 shift (-0.11) exceeded Pade tolerance (0.03) ... redo GF pass

The shift is too large for the Pade to be reliable, so lmgf remakes in the presence of this shift.

After another cycle generating the proper we find

         deviation from charge neutrality: 1.042307

so the Pade approximation captured 98% of the missing charge. This round the Pade yields

 gfasa:  potential shift this iter = 0.002996.  Cumulative shift = -0.102274

which is below the tolerance  0.03 (0.03 is a default tolerance, used if no number is explicitly specified) that would trigger yet another lmgf pass and the program stops.

Inspect file vshft.asa. lmgf required a constant potential shift,  vconst=-0.1022739, which is very close to −EF generated by lm. The close agreement should give use some confidence that the Green’s function is well described.

The RMS DQ is printed in this line:

 PQMIX:  read 3 iter from file mixm.  RMS DQ=1.38e-3

It is small enough that further self-consistency is not really necessary (self-consistent potential generated by lm and lmgf are very similar). Neverthess we can make it self-consistent:

rm -f save.asa mixm.asa sv.asa log.asa      &larr; Fresh start
mpix -n 16 lmgf ctrl.asa -vccor=f -vnit=50 -vscr=2 --pr40,21 -vgamma=t --rs=0,1 >> out.lmgf
lm asa --rsedit~rs~save@nopot@fn=rsta6~q    &larr;  preserves density for Section 6
mv rsta6.asa rsta.asa.lmgf                  &larr;  Restart file for density
head -1 vshft.asa > vshft.asa.gf            &larr;  preserves potential shift for Section 6

It should converge smoothly in about 8 iterations, with a converged potential shift  vconst=-0.1014863. Total energies and system magnetic moments in lmgf and lm are also very similar:

grep ^c out.lm out.lmgf
out.lm:c met=3 ccor=0 nit=50 scr=2 mmom=23.7892192 ehf=-122072.0118495 ehk=-122072.011848
out.lmgf:c ccor=0 nit=50 scr=2 gamma=1 mmom=23.8092855 ehf=-122072.0114363 ehk=-122072.0114346

which shows that the hamiltonians are similar.

Spectral Function from lmgf

A Green’s functions does not supply quasiparticle levels, but it does contain the spectral function . The analog of the band structure can be got by drawing as a “heat map” on a grid. Since in this case there are no alloy or many-body scattering processes, the spectral function is fully coherent with unit Z factor and sharp peaks at the QP levels. They show up as bright spots on a heat map.

lmgf will make data for a heat map in lieu of the band structure, using the --band switch. By passing the Fermi level within the --band switch, lmgf will shift the spectral function to make EF coincide with 0. This is useful for comparison to the band structures drawn previously.

lmchk ctrl.asa  --syml~n=41~q=0,1,0,0,0,0,1,0,-1
mpix -n 16 lmgf ctrl.asa -vccor=f -vnit=50 -vscr=2 --pr40,21 -vgamma=t --band:fn=syml
SpectralFunction.sh  -w -5:5
open asa_UP.png
open asa_DN.png

Majority (left) and minority (right) spectral functions bands of the Nb/Fe superlattice

NbFe SF up figure.

NbFe SF dn figure.

Spectral functions correspond well to the energy bands generated by lm, Note in particular the narrow Fe d bands are bunched together and show greater intensity (red). They correspond to the Fe d bands, seen as nearly flat bands in green near the Fermi level.

6. Nb/Fe/Nb Trilayer

The current system has two layers of Nb, followed by three Fe layers, followed by two more Nb layers. Since boundary conditions are periodic the geometry is

··· Nb ¦ Fe ¦ Nb ¦ Nb ¦ Nb ¦ Fe ¦ Nb ¦ Nb ¦ ···

lmpg works for a trilayer geometry:

··· Left ¦ Left ¦ Left ¦ Active ¦ Right ¦ Right ¦ Right ¦ ···

In our case Left and Right are both (3×1) reconstructed Nb. The density should correspond to that of bulk Nb, calculated in Section 1.

To make the transition, the following preparatory steps must be taken:

  • A PGF category must be added to the ctrl file, analogous to the GF category. This includes:
    The 3rd lattice vector out of the plane must be defined for the left end lead (PGF_PLATL), and
    the 3rd lattice vector out of the plane must be defined for the right end lead (PGF_PLATR).

  • lmpg cannot yet handle inequivalent L-cutoffs. In prior sections we used LMX=1 for the E sites; this must be converted to LMX=2.

  • Principal layers assigned (information embedded in the site file)

  • For the periodic case with a larger basis, a new polarizability file constructed the self-consistent density remade

  • A trial restart file constructed. The one generated for the periodic case can be adapted to the trilayer case.

This step is normally taken following Section 5. To run it as a standalone test do:

~/lmf/testing/test.nbfeasa --quiet 1 2
cp ~/lmf/testing/nbfe/rsta.asa.lmgf ~/lmf/testing/nbfe/vshft.asa.gf .

test.nbfeasa will perform the steps of this section automatically, with the following:

~/lmf/testing/test.nbfeasa [--MPIK] 5

It is useful as a sanity check to check that you’ve done all the steps in the proper order.

Section 6 consists of the following steps:

a. Use blm to add PGF category to the ctrl file and also set LMX=2 (output = ctrl). b. Use lmscell to assign principal layers (site4)
c. Run lm-vscr=11 to make the static ASA response function in this basis (output = psta)
d. Run lmgf-vscr=2 to converge to self-consistency in larger basis
e. Use lmpg--rsedit to create a restart fle for a trilayer, suitable for lmpg. It imports the restart file for the active region from the periodic case created by lmgf and the restart file for padding layers from bulk Nb. (output = rsta.asa)
f. Make a starting potential from the restart file.
g. Make L- and R- end leads self-consistent
h. Make the active region self-consistent (output = rsta9.asa)

a. Remake ctrl.asa with the following:

blm --pr35 --nl=3 --pgf~platl=0,0,1~platr=0,0,1 --gf --nfile --rdsite=site7 --rdspec --usefiler --nosort --shorten=no --nk=9,4,2 --mag --asa~pfree~rmaxs=7.2 --wsite@fn=null asa --omax=.18,.26,.26
cp actrl.asa ctrl.asa
sed -i.bak -e 's/LMX=1/LMX=2 IDMOD=0,0,1/' ctrl.asa
rm ctrl.asa.bak

blm is invoked in a manner similar to the construction for the crystalline case but we do not need to worry about finding E spheres or adjusting sphere radii. We can just use what we got for the periodic case. The main addition is the --pgf switch: It tells blm to add a PGF category, specifying tags PGF_PLATL and PGF_PLATR (see here). The reconstructed Nb layers are periodic normal to the plane with lattice vector (0,0,1) (third lattice vector in site3.nbasa, Section 1.

b. Define principal layers. The principal layer technique relies on the hamiltonian being short ranged; then the matrix becomes sparse. The normal axis must be partitioned into principal layers, constructed so that the hamiltonian extends only to neighboring PL, making the hamiltonian tridiagonal in the PL.

The ASA hamiltonian is indeed short-ranged; the hamiltonian can reasonably be truncated at about 7 a.u. — second neighbors for bcc Nb and Fe (see Figures 5 and 15 of Ref. 1).

lmscell has a utility to assemble PL indices automatically given a specified range of the hamiltonian (which amount to the range of structure constants in the ASA). See here for documentation of the stack editor. The following reads structural data from site5.asa, finds the PL and writes the updated information to site3.asa:

lmscell -vpgf=1 -vfile=5 ctrl.asa --stack~setpl:withinit:rangea=1.1~wsite@fn=site4

The stack editor relies on the fact that lmscell was previously invoked with initpl, to make site5.asa.

c. Remake psta.asa in the larger basis (same steps as in Section 3).

cp rsta.asa.lmgf rsta.asa
lmstr -vfile=4 asa > out.lmstr2
touch -vfile=4 nb.asa; rm -f {nb,fe,e}*.asa save.asa mixm.asa sv.asa log.asa
lm -vfile=4 asa -vnit=-1 --pr30,20 --rs=1,0 -vscr=11 > out.psta3

d. Make a new self-consistent density in the enlarged basis

rm -f save.asa mixm.asa sv.asa log.asa vshft.asa       &larr; Fresh start
cp vshft.asa.gf vshft.asa                              &larr; Not required, but saves work since vshft for large basis similar
lmgf -vfile=4 ctrl.asa -vccor=f -vnit=50 -vscr=2 --pr40,21 -vgamma=t --rs=0,1 > out.lmgf2
lm asa --rsedit~rs~save@nopot@fn=rsta7~q > out.lm5     &larr; Not required, keeps restart file for self-consistent lmgf
cp rsta7.asa rsta.asa.lmgf2                            &larr; Restart file for self-consistent lmgf density
head -1 vshft.asa > vshft.asa.gf2                      &larr; Not required, preserves shift for self-consistent lmgf

e. Create a restart file for the trilayer geometry

cp rsta.asa.lmgf2 rsta.asa
cat vshft.asa.gf2 | awk '-F[= ]' '{printf " ef=%s  vconst=%s vconst(L)=%s vconst(R)=%s\n", $3,$6,$6,$6}' > vshft.asa
lmpg -vpgf=1 -vfile=4 ctrl.asa --rsedit~rs@3d~rs@ib=61@src=1:3@fn=rsta2~rs@ib=64@src=1:3@fn=rsta2~save@nopot@fn=rsta8~q > out.lmpg0
cp rsta8.asa rsta.asa.pgf0
cp rsta.asa.pgf0 rsta.asa
rm -f save.asa mixm.asa sv.asa log.asa vshft.asa
check rsta.asa.pgf0 :  max deviation (0) within tol (2e-6) of /home/ms4/lmf/b/./testing/nbfe/rsta.asa.pgf0 ? ... yes
lmpg -vpgf=1 -vfile=4 ctrl.asa -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rs=1,0 > out.lmpg1
lmpg -vpgf=2 -vfile=4 ctrl.asa -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rs=1,0 > out.lmpg2

The restart file editor is used in a similar manner to Section 4, apart from the new feature **rs~3d$$. For restart file editor instructions, see the Appendix.

f. Create a starting potential for both end leads and the active region. Note -vnit=0 –rs=1,0 causes lmpg to read P, Q from the restart file, make potential parameters needed to make G, and exit.

rm -f save.asa mixm.asa sv.asa log.asa vshft.asa
check rsta.asa.pgf0 :  max deviation (0) within tol (2e-6) of /home/ms4/lmf/b/./testing/nbfe/rsta.asa.pgf0 ? ... yes
lmpg -vpgf=1 -vfile=4 ctrl.asa -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rs=1,0 > out.lmpg1
lmpg -vpgf=2 -vfile=4 ctrl.asa -vccor=f -vnit=0 -vscr=11 --pr40,21 -vgamma=t --rs=1,0 > out.lmpg2

g. Make the end leads self-consistent.

lmstr -vpgf=1 -vfile=4 asa > out.lmstr2
rm -f save.asa mixm.asa sv.asa log.asa vshft.asa
lmpg -vpgf=2 -vfile=4 ctrl.asa -vccor=f -vnit=20 -vscr=0 --pr35,21 -vgamma=t >> out.lmpg2

Note -vpgf=2, which assigns PGF_MODE to 2 (see the ctrl file). The function of that tag can be found in the documentation for lmpg, or use the on-line help:

lmpg --input | grep -A 10 PGF_MODE

h. Make the active region self-consistent. Note -vpgf=1 in this case.

rm -f save.asa mixm.asa sv.asa log.asa
/home/ms4/lmf/b/lmpg -vpgf=1 -vfile=4 ctrl.asa -vccor=f -vnit=50 -vscr=2 --zerq~all~qout --pr35,21 -vgamma=t --iactiv=no --rs=0,1 > out.lmpg
lmpg -vpgf=1 -vfile=4 ctrl.asa --rsedit~rs~save@nopot@fn=rsta9~q > out.lmpgf

The last step is not needed, but it makes compressed restart file, rsta9.asa which is sufficient to reconstruct the density. Thanks to the presence of psta.asa, convergence to self-consistency is smooth.

Footnotes and references

1 Dimitar Pashov, Swagata Acharya, Walter R. L. Lambrecht, Jerome Jackson, Kirill D. Belashchenko, Athanasios Chantis, Francois Jamet, Mark van Schilfgaarde, Questaal: a package of electronic structure methods based on the linear muffin-tin orbital technique, submitted to CPC. Preprint https://arxiv.org/abs/1907.06021