Nb/Fe Superlattice

The electronic structure of a Nb/Fe superlattice is set up, using Questaal’s superlattice editor lmscell.

Nb and Fe are both bcc but are strongly lattice mismatched, which considerably complicates matters. An interface is constructed by rotating both Fe and Nb off the cubic axes, so that the two lattices normal to the interface approximately align. The Fe is slightly sheared normal to the interface to make the two align exactly.

The tutorial begins by building reconstructed elemental Nb and Fe unit cells separately and determinng the electronic structure. Then the Nb/Fe superlattice is constructed, and the atomic positions relaxed. The relaxed lattice is needed to set up a follow-on tutorial for Landauer-Buttiker transport through a Nb/Fe/Nb trilayer.

A separate tutorial also explains how to construct superlattices with lmscell. The superlattice in that case is simpler; it is written for QSGW calculations.

Table of Contents

• Preliminaries
• Introduction
  • Basic strategy
• Part 1. Elemental Nb and elemental Fe in reconstructed unit cells
  • Structural information for elemental Nb and Fe
  • Reconstructed Nb, self-consistent LDA
    • Energy bands of reconstructed Nb
  • Reconstructed Fe, self-consistent LDA
    • Energy bands of reconstructed Fe
• Part 2. The Nb(2)/Fe(3)/Nb(2) superlattice
  • Create a simple Nb(2)/Fe(3)/Nb(2) superlattice
  • Strategy to make practical Nb/Fe/Nb superlattices
  • 1. Ideal Nb(2)/Fe(3)/Nb(2) superlattice
  • 2. Crude estimate for relaxation of the Nb(2)/Fe(3)/Nb(2) superlattice
  • 3 Relaxed Nb(2)/Fe(3)/Nb(2) superlattice
    • 3a. Self-consistency
    • 3b. Constrained relaxation
    • 3c. Optimizing the Nb/Fe layer spacing
    • 3d. LDA Energy Bands and layer resolved DOS of the superlattice
• Part 3. Longer period superlattices
• Additional exercises
• Footnotes and references

touch init.nb; rm *.nb *.fe;
cp $HOME/lmf/fp/test/nbfe/{init.nb,init.fe} .

rm -f out.blm out.lmscell out.lmfa out.lmf out.lmchk
blm --nfile --ctrl=ctrl --gw nb
echo p 1 0 0  0 0 1  1 1 0 | lmscell nb --wsite~short~fn=site2 >: out.lmscell
echo m 1 0 0 0 3 0 0 0 1 | lmscell -vfile=2 nb --wsite~short~fn=site3 '--rot=(0,1,0)pi/4' >> out.lmscell
blm --autobas~loc=9 --nfile --gmax=5.8 --nk=12,12,12 --ctrl=ctrl --gw nb >> out.blm
lmfa -vfile=3 ctrl.nb >> out.lmfa
cat basp0.nb > basp.nb
grep Nb basp0.nb | sed s/Nb/Nbx/ >> basp.nb
lmfa -vfile=3 ctrl.nb >> out.lmfa
$lmf -vnit=20 -vfile=3 ctrl.nb >> out.lmf
lmchk ctrl.nb -vfile=3 '--syml~n=11~q=0,1,0,0,0,0,1/sqrt(2),0,0' >> out.lchk
$lmf -vfile=3 ctrl.nb --band~col=1,17,seq=31,61~col2=2:4,18:20,seq=32,62~col3=5:9,21:25,26:30,seq=35,65~fn=syml >> out.lmf
echo -6,6,5,10 | plbnds -fplot~sh~scl=.5 -ef=0 -scl=13.6 -dat=green -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 bnds.nb
open fplot.ps

blm --nfile --ctrl=ctrl --mag --gw fe >> out.blm
echo p 1 0 0  0 0 1  1 1 0 | lmscell fe --wsite~short~fn=site2 >> out.lmscell
echo m 1 1 0 -3 1 0 0 0 1 | lmscell '--scala=sqrt(4/3)' '--rot=(1,0,1)-acos(1/sqrt(3)),y:pi/4' -vfile=2 fe --wsite~short~fn=site3 >> out.lmscell
lmscell -vfile=3 ctrl.fe -vz=6.24000171/6.226 '--stack~scale=1/z~stretch=z^3~sort@h3~wsite@short@fn=site4' >> out.lmscell
blm --nit=30 --autobas~loc=9 --nfile --nk=10,5,13 --ctrl=ctrl --gw --gmax=8.4 --mag fe >> out.blm
lmfa -vfile=4 ctrl.fe >> out.lmfa
cat basp0.fe > basp.fe
$lmf -vfile=4 ctrl.fe >> out.lmf
cp syml.nb syml.fe
$lmf -vfile=4 ctrl.fe  --band~col=1,17,seq=26,51,76~col2=2:4,18:20,seq=27,52,77~col3=5:9,21:30,seq=35,60,95~fn=syml >> out.lmf
echo -6,6,5,10 | plbnds -spin2 -fplot~sh~scl=.5~shft=.5  -ef=0 -scl=13.6 -dat=dn -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 bnds.fe
open fplot.ps

mcx -cmpf~fn1=ctrl.nb~fn2=$HOME/lmf/fp/test/nbfe/ctrl.nb
mcx -cmpf~fn1=basp.nb~fn2=$HOME/lmf/fp/test/nbfe/basp.nb
mcx -cmpf~fn1=site3.nb~fn2=$HOME/lmf/fp/test/nbfe/site3.nb
mcx -cmpf~fn1=ctrl.fe~fn2=$HOME/lmf/fp/test/nbfe/ctrl.fe
mcx -cmpf~fn1=site4.fe~fn2=$HOME/lmf/fp/test/nbfe/site4.fe
mcx -cmpf~fn1=basp.fe~fn2=$HOME/lmf/fp/test/nbfe/basp.fe

rm -f *.nb *.nbfe out.{blm,lmchk,lmscell}
cp $HOME/lmf/fp/test/nbfe/{ctrl.nb,site3.nb,site4.fe,basp.nb,basp.fe} .

lmscell -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp1=0 -vdp2=0 -vdp3=0 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=3@dpos=dp1/2,dp2/2,-dz+dp3/2~~file@site3.nb@dup=2@dpos=dp1,dp2,dp3~show~basp@fn=basp.nbfe@basp.nb@basp.fe~shorp3~pad=dp1,dp2,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~' >> out.lmscell
cp pos.nb pos0.nbfe
cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=.16 --nosort --shorten=no --wpos=pos0 --ctrl=ctrl nbfe >> out.blm
lmchk nbfe --mino~dxmx=.001~xtol=.0001~nkill=10~maxit=200~4:10,15:21 --wpos=pos  >> out.lmchk
mcx -f3f12.6 pos.nbfe pos0.nbfe -- > dpos30.nbfe
lmscell -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp3=0 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=3@dpos=0,0,-dz+dp3/2~~file@site3.nb@dup=2@dpos=0,0,dp3~setlbl@spid=Nbx@targ=4:6,nbas-5:nbas-3~basp@fn=basp.nbfe@basp.nb@basp.fe~shorp3~pad=0,0,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~' >> out.lmscell
cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=8 --mag --frzwf --gw nbfe '--wsite~ib=\!4:10,15:21~rlx=000~long~fn=sitex' >> out.blm
cp sitex.nbfe site0.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=9 --mag --frzwf --gw nbfe --rdpos=dpos30 --rsham~el=-.1,-0.9~rscut=1.6 --dhftol=45 --molstat~xtol=0~gtol=0.01 >> out.blm
cat actrl.nbfe | sed s:'\(ATOM=Fe.*\):\1 MMOM=0,0,2:' > ctrl.nbfe
lmscell -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp3=-.10 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=3@dpos=0,0,-dz+dp3/2~~file@site3.nb@dup=2@dpos=0,0,dp3~setlbl@spid=Nbx@targ=4:6,nbas-5:nbas-3~shorp3~pad=0,0,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~' >> out.lmscell
cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=8 --mag --frzwf nbfe '--wsite~ib=\!4:10,15:21~rlx=000~long~fn=sitex' >> out.blm
cp sitex.nbfe site10.nbfe
lmscell -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp3=-.20 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=3@dpos=0,0,-dz+dp3/2~~file@site3.nb@dup=2@dpos=0,0,dp3~setlbl@spid=Nbx@targ=4:6,nbas-5:nbas-3~shorp3~pad=0,0,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~' >> out.lmscell
cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=8 --mag --frzwf nbfe '--wsite~ib=\!4:10,15:21~rlx=000~long~fn=sitex' >> out.blm
cp sitex.nbfe site20.nbfe

mcx -cmpf~fn1=ctrl.nbfe~fn2=$HOME/lmf/fp/test/nbfe/ctrl.nbfe
mcx -cmpf~fn1=site0.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site0.nbfe
mcx -cmpf~fn1=dpos30.nbfe~fn2=$HOME/lmf/fp/test/nbfe/dpos30.nbfe
mcx -cmpf~fn1=basp.nbfe~fn2=$HOME/lmf/fp/test/nbfe/basp.nbfe
mcx -cmpf~fn1=site10.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site10.nbfe
mcx -cmpf~fn1=site20.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site20.nbfe

rm -f *.nbfe spec.ni spec.nb site4.ni out.0 out.10 out.20
cp $HOME/lmf/fp/test/nbfe/{ctrl.nbfe,site0.nbfe,site10.nbfe,site20.nbfe,basp.nbfe,dpos30.nbfe} .

rm -f mixm.nbfe hssn.nbfe log.nbfe save.nbfe rst.nbfe
lmfa -vfile=0 ctrl.nbfe >> out.lmfa
mpirun -n 14 lmf ctrl.nbfe -vfile=0 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos --rdpos=dpos30 --shorten=2,2,1 --wpos=posx >> out.0
mpirun -n 14 lmf ctrl.nbfe -vfile=0 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos -vfrzwf=t -vminx=10 --wpos~fn=posx~shorten@0,0,1 >> out.0
cp posx.nbfe pos0.nbfe

mcx -f3f14.8 -vdp3=-.1 pos0.nbfe -av:nr,3 p3max -e3 0 0 'dp3*x3/p3max' pos0.nbfe -+ > posnow.nbfe
rm -f mixm.nbfe hssn.nbfe log.nbfe save.nbfe rst.nbfe
mpirun -n 14 lmf ctrl.nbfe -vfile=10 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos --rpos=posnow --shorten=2,2,1 > out.10
mpirun -n 14 lmf ctrl.nbfe -vfile=10 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos -vfrzwf=t -vminx=10 --rpos=posnow --shorten=no --rs=11,1,1 --wpos~fn=posx~shorten@0,0,1 >> out.10
cp posx.nbfe pos10.nbfe

mcx -f3f14.8 -vdp3=-.2 pos0.nbfe -av:nr,3 p3max -e3 0 0 'dp3*x3/p3max' pos0.nbfe -+ > posnow.nbfe
rm -f mixm.nbfe hssn.nbfe log.nbfe save.nbfe rst.nbfe
mpirun -n 14 lmf ctrl.nbfe -vfile=20 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos --rpos=posnow --shorten=2,2,1 > out.20
mpirun -n 14 lmf ctrl.nbfe -vfile=20 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos -vfrzwf=t -vminx=10 --rpos=posnow --shorten=no --rs=11,1,1 --wpos~fn=posx~shorten@0,0,1 >> out.20
cp posx.nbfe pos20.nbfe

mcx -cmpf~fn1=pos0.nbfe~fn2=$HOME/lmf/fp/test/nbfe/pos0.nbfe
mcx -cmpf~fn1=pos10.nbfe~fn2=$HOME/lmf/fp/test/nbfe/pos10.nbfe
mcx -cmpf~fn1=pos20.nbfe~fn2=$HOME/lmf/fp/test/nbfe/pos20.nbfe

lmchk ctrl.nbfe -vfile=20 '--syml~n=21~q=0,1,0,0,0,0,1/sqrt(2),0,0'
mpix -n 16 lmf ctrl.nbfe -vfile=20 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 --band~col=1:180,541:720~col2=181:540~fn=syml
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.nbfe
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.nbfe
open fplot.ps

touch ctrl.nbfe; rm -f *.nbfe out.lmchk out.lmscell out.lmfa out.lmf
cp $HOME/lmf/fp/test/nbfe/{ctrl.nbfe,site0.nbfe,site10.nbfe,site20.nbfe,basp.nbfe,dpos30.nbfe} .

lmchk ctrl.nbfe -vfile=20 --wpos=posx >> out.lmchk
mcx -f3f14.8 pos20.nbfe posx.nbfe  -- > dpos20.nbfe

foreach n (4 5 6 7 8 9)
echo 
echo "... make site$n and dpos$n"
lmscell -vn=$n -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp3=-.20 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=n@dpos=0,0,-dz+dp3/2~~file@site3.nb@dup=2@dpos=0,0,dp3~setlbl@spid=Nbx@targ=4:6,nbas-5:nbas-3~shorp3~pad=0,0,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~' >> out.lmscell
cp sitein.nb sitein.nbfe
blm -vn=$n --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=8 --mag --frzwf nbfe '--wsite~ib=\!4:10,11+4*1:11+4*n+6~rlx=000~long~fn=sitex' >> out.blm
cp sitex.nbfe site$n.nbfe
mcx -f3f14.8 -vn=$n -vm=n-2 [ 0:'4*m'-1 '-array[1,3]' 0,0,0 '?i>1' -rcat ] -a ins dpos20.nbfe  -split x 1,11,15,nr+1 1,nc+1 x11 ins x31 -rcat -rcat > dposx.nbfe
cp dposx.nbfe dpos$n.nbfe
end

mcx -cmpf~fn1=site4.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site4.nbfe
mcx -cmpf~fn1=site7.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site7.nbfe
mcx -cmpf~fn1=dpos7.nbfe~fn2=$HOME/lmf/fp/test/nbfe/dpos7.nbfe

Preliminaries

This tutorial makes extensive use of many of the Questaal packages. including lmscell, lmplan lmf, and miscellaneous utilities such as vextract. They are assumed to be in your path, or in the directory where you cloned the Questaal package (assumed here to be $HOME/lmf).

It 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 fewer, or run in serial mode.

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

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

Introduction

This tutorial builds a complex Nb/Fe superlattice, using the LDA code lmf to relax the structure of a periodic Nb / Fe superlattice. Even though Fe and Nb are bcc lattices, they are severely lattice mismatched. However, this tutorial shows how Nb and Fe can be reconstructed so Nb and Fe cells can be stacked on top of one another with good lattice matching. The disruption at the interface is large, however, and there is a large attendant lattice relaxation, which must be dealt with.

The tutorial then switches to the ASA approximation. The ASA-LDA code lm code is then applied to the relaxed superlattice to confirm the reliability of the ASA.

Next the crystal Green’s function package lmgf as setup to the layer Green’s function package lmpg. This latter code is used to perform transmission calculations through the interface. It requires a non-periodic trilayer geometry (… Nb / Fe / Nb …), which is related to but distinct from the Nb/Fe superlattice.

Basic strategy

To stack unit cells of Nb(110) and Fe(111) on top of one another, the lattice vectors in the stacking plane must match. First the following elemental cells are constructed:

• a special reconstruction of bcc Nb with three atoms in the (110) plane. The crystal rotated by π/4 around y to render this plane perpendicular to z.
• a different reconstruction of bcc Fe with four atoms in a (111) plane. In this case the z axis is rotated to put the (111) plane normal to z, and the Fe planes are slightly stretched to exactly match lattice constants of the reconstructed Nb and Fe cells (coincident site lattices in the plane normal to the new z direction in each case).

Next the superlattice is made.

• a rigid superlattice of 7 planes, Nb(2)/Fe(3)/Nb(2), is assembled assuming that the frontier Nb/Fe planes are spaced at the average of the Nb/Nb and Fe/Fe planes.
• A crude estimate for how the atoms a the interface relax is made by minimizing some empirical function of overlaps between spheres drawn around each site.
• The lattice is relaxed subject to constraints described below. Constraints are imposed so the superlattice is built in a lego-like fashion, layer by layer. This greatly simplifies extension to larger period Nb(n)/Fe(m)/Nb(n) lattices.
• The assumption of Nb/Fe spacing is removed by minimizing the total energy, repeating the relaxation at several values of the Fe/Nb spacing.

Not described here, but in a follow-on tutorial, Landauer-Buttiker transport through a Nb/Fe/Nb trilayer is described.

• Repeat the self-consistent calculation of the interface code with the ASA-LDA code lm. lm is more approximate than the full LDA code lmf, but it is very efficient, and provides a stepping stone to the transport code, lmpg.
  The interface has an additional complication in that open spaces appear at the interface, which the ASA must deal with by packing with additional empty spheres.
• Use lmpg to perform transmission calculations.

The prescription is in outline:

• Use blm to build input files for elemental Nb and Fe in simple bcc structures.
• Use the supercell maker lmscell to generate different reconstructions and rotations of the simple bcc Nb and Fe lattices that can be together with the Nb(110) and Fe(111) planes.
• The reconstructed elemental systems will be made self-consistent and their bands in the plane of the interface examined.
• Use the superlattice maker lmscell --stack to stack reconstructed Nb and Fe lattices into superlattices.
• Use the LDA code lmf to relax the interface, since the interface has a significant disruption of atomic environments relative to bcc. Relaxation proceeds in several steps; see Steps 2, 3

Part 1. Elemental Nb and elemental Fe in reconstructed unit cells

This section requires files init.nb and init.fe. Cut and paste these files from boxes below, or copy them the Questaal distribution:

cp $HOME/lmf/fp/test/nbfe/{init.nb,init.fe} .

Note: The steps in this section have been bundled into a test script. To confirm what you do is working as it should, you can check your steps against those of the script

$HOME/lmf/fp/test/test.nbfe --MPIK=16 1

Three files are required to use lmf to make the self-consistent electronic structure for a specified lattice. All are autogenerated once structural information is supplied.

• A ctrl file, which is the main input file. We will autogenerate this file using the blm utility.
• A site file with structural data. (The ctrl file can also hold structural information, but we keep it separate so only the site changes as the structure changes.) blm also makes this file.
• A basp file, which defines the basis set. The free-atom code lmfa makes this file.

We need to assemble files in several contexts:

• For the end leads, Nb
• For the central region, Fe
• Fe/Nb superlattices, built so that the interface can be relaxed
• A non-periodic trilayer (… Nb Nb / Fe / Nb Nb …) suitable for study of transport of an active region between semi-infinite leads

Structural information for elemental Nb and Fe

Copy the boxes below into files init.nb and init.fe. They contain the requisite structural information for elemental bcc Nb and Fe.

# Init file for Nb.  Lattice constant at 0K
LATTICE
    ALAT=6.226 PLAT= -.5 .5 .5  .5 -.5 .5  .5 .5 -.5
SPEC
    ATOM=Nb   PZ=0,0,0
SITE
    ATOM=Nb   X=0 0 0

# Init file for Fe.  Lattice constant at 0K
LATTICE
    ALAT=5.404 PLAT= -.5 .5 .5  .5 -.5 .5  .5 .5 -.5
SPEC
    ATOM=Fe  MMOM=0,0,2.2
SITE
    ATOM=Fe   X=0 0 0

Notes:

• Both crystals are in the bcc lattice structure; however, the volume differential is large: much too large to grow Fe on Nb epitaxially without reconstructions of both lattices and a rotation of one structure against the other.
• Fe is ferromagnetic. We must give the starting point some estimate of the moment to allow it to find the ferromagnetic solution.

Also, it is not necessary to do this, but if you want the steps to follow exactly consistent with the present tutorial, create the following file site4.fe.

% site-data vn=3.0 fast io=15 nbas=4 alat=6.226 plat= -0.70710678 0.5 0.0 2.12132034 1.5 0.0 -0.35594409 -0.50338096 0.61651325
#                        pos
 Fe        0.0000000   0.0000000   0.0000000
 Fe        0.0000000   0.7500000   0.0000000
 Fe        0.7071068   1.0000000   0.0000000
 Fe        1.4142136   1.2500000   0.0000000

The lmscell instruction below should make this file, but it may possible change the site order.

Reconstructed Nb, self-consistent LDA

To get started touch init.nb; rm -f *.nb ; cp /home/ms4/lmf/b/./fp/test/nbfe/{init.nb} .

Build a template ctrl and site file for bcc Nb. This file is incomplete but sufficient for steps needed for the reconstruction.

blm --nfile --ctrl=ctrl --gw nb

--nfile tells blm to insert the following lines into ctrl.nb :

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

This allows Questaal codes to read structural data from siten.ext, when preprocessor variable file is assigned to n. These lines have no effect if file is not assigned, in which case Questaal codes will read from the usual file site.ext.

It turns out that the under a (3×1) reconstruction the [001] planes of Nb can be nearly epitaxially joined with a reconstruction of Fe with four atoms/plane. Note that the ratio of areas of equivalent planes in the two systems is 6.226²/5.404² = 1.327, or nearly 4/3. Anticipating this in advance, the tutorial builds directly a Nb superlattice that will join onto Fe, with three atoms in a plane of Nb and 4 in a plane of Fe. The area lattice mismatch will be only 4/3/1.327 = 1.0045.

echo p 1 0 0  0 0 1  1 1 0 | lmscell nb --wsite~short~fn=site2

Inspect site file site2.nb. The third lattice vector of points along the z axis. site2.nb. still has only one atom/cell; this is an equivalent, though less symmetric way of writing the primitive unit cell.

Next, make a (3×1) reconstruction of this lattice

echo m 1 0 0 0 3 0 0 0 1 | lmscell -vfile=2 nb --wsite~short~fn=site3 '--rot=(0,1,0)pi/4'

Inspect site3.nb, and note that it has 3 atoms/cell. Use lmplan to examine some of its characteristics:

lmplan -vfile=3 nb

 SPLCLS:  1 species split into 2 classes
 Species  Class      Sites...
 Nb       1:Nb       1
          2:Nb2      2 3

 Class        Qtot       Qbak       Vmad     Vh(Rmax)    V(Rmax)
 Nb         0.000000   0.000000   0.000000   0.000000  -0.700000
 Nb2        0.000000   0.000000   0.000000   0.000000  -0.700000
 Sum Q=0.000000  Emad=0.000000(0.000000)  Vmtz=-0.700000 (up-dn=-0.700000)

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

    1       1:Nb          1    0      0.00000   0.00000   0.00000   0.00000   0.70711
    2       2:Nb2         1    0      0.70711   0.50000   0.00000   0.00000   0.00000
    3       2:Nb2         1    0      1.41421   1.00000   0.00000   0.00000   0.00000

Notes:

• The planes defined by the first two vectors P₁×P₂ stack lie normal to the [0,0,1] direction (‘normal’ vector in the output)). If --rot=(0,1,0)pi/4 had been omitted, the normal direction would be [1,0,1]. The z axis is in fact the [1,0,1] direction of the bcc lattice.
• All three atoms lie in a single (001) plane.
• Spacing between planes is 0.70711, which is 1/√2.
• The 48 symmetry operations get reduced to 8, and the three Nb atoms are no longer symmetry-equivalent. This is artificial since this superlattice is still elemental Nb, but it is unimportant for our purposes.

Complete the input file for an lmf calculation and make the density self-consistent. To complete the input file the G cutoff parameter is needed to specify the mesh density. Use lmfa to generate a good estimate:

lmfa -vfile=3 ctrl.nb

Read off the suggested value 5.8 from the end of the output. (Note the value 8.4 would be needed if we include local orbitals. We will suppress them,cre as explained shortly.)

We also need a k mesh, which is chosen conservatively since the calculation is fast anyway. The steps below make a complete ctrl.nb and rerun lmfa since ctrl.nb has changed (In this case it isn’t necessary, but in it’s a safe practice in general)

blm --nfile --gmax=5.8 --nk=12,12,12 --ctrl=ctrl --gw nb
lmchk ctrl.nb --wspec~fn=fpspec
lmfa -vfile=3 ctrl.nb

lmchk writes file fpspec.nb, to be used in a later section.

lmfa creates a basis file, and writes it to basp0.ext. It automatically found local orbitals for the Nb 4p state and Nb 5d state. For high accuracy GW calculations they can matter, but these extensions are an unnecessary complication the present purposes. From lmfa’s output you can see that the Nb p state sits at -2.54 Ry, a moderately shallow core state. Copy this autogenerated file to basp.nb, stripping the PZ tag for the local orbitals:

cat basp0.nb | sed 's/ *PZ.*//' > basp.nb

Before proceeding to self-consistency it is useful to check that the basis is well constructed. To see how to do this and things to look at, see the PbTe tutorial.

Make the density self-consistent

lmfa -vfile=3 ctrl.nb
mpix -n 16 lmf -vnit=20 -vfile=3 ctrl.nb

The forces should all be identically zero, but there are slight deviations from zero, because e.g. because the k mesh doesn’t possess the full cubic symmetry. In fact, the forces are very small.

One way to confirm the density is nearly self-consistent is to compare the Harris-Foulkes and Kohn-Sham total energies. These functionals are different but they should approach the same value at self-consistency. You should find convergence in 10 iterations with

c nit=20 file=3 ehf=-22895.0474411 ehk=-22895.0474405

Energy bands of reconstructed Nb

We are interested in interfaces and energy bands in the plane normal to the interface. Draw bands one lines in a plane normal to (1,0,1). This plane is defined by the lines (0,1,0) and (1,0,-1). Vectors (0,1,0) and (1/2,0,-1/2) are the smallest vectors that touch a boundary point on the Brillouin zone (half-integer multiples of reciprocal lattice vectors). To confirm this, use the cart2lat feature of the superlattice stack maker

lmscell ctrl.nb -vfile=3 '--stack~cart2lat@q@x=0,1,0~cart2lat@q@x=1/sqrt(2),0,0~q'

Create syml.nb with

lmchk ctrl.nb -vfile=3 '--syml~n=31~q=0,1,0,0,0,0,1/sqrt(2),0,0'

syml.nb should be created in symmetry lines format.

Further, it is very useful to know the orbital characteristic of the states, particularly at the Fermi surface where transport takes place. For this purpose, we use a Mulliken decomposition of states into the s, p, d orbital character, and draw bands with color weights, e.g. red for s, green for p, and blue for d.

To make this decomposition we need to where the respective orbitals reside in the LMTO hamiltonian. lmf prints this table out:

lmf -vfile=3 ctrl.nb --pr55 --quit=ham | grep -A 5 ' Orbital positions in hamiltonian, resolved by l'

You should see the following table:

 Site  Spec  Total    By l ...
   1   Nb    1:33   1:1(s)   2:4(p)   5:9(d)   10:16(f) 17:17(s) 18:20(p) 21:25(d) 26:28(p) 29:33(d)
   2   Nb   34:66   34:34(s) 35:37(p) 38:42(d) 43:49(f) 50:50(s) 51:53(p) 54:58(d) 59:61(p) 62:66(d)
   3   Nb   67:99   67:67(s) 68:70(p) 71:75(d) 76:82(f) 83:83(s) 84:86(p) 87:91(d) 92:94(p) 95:99(d)

The s, p, d orbitals are conveniently grouped by using the Questaal standard syntax for Integer Lists. The following creates and makes a postscript figure for the energy bands:

mpix -n 16 lmf -vfile=3 ctrl.nb --band~col=1,17,seq=34,67~col2=2:4,18:20,seq=35,68~col3=5:9,21:25,29:33,seq=38,71~fn=syml
echo -6,6,5,10 | plbnds -fplot~sh~scl=.5 -ef=0 -scl=13.6 -dat=green -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 bnds.nb
open fplot.ps

Color weights group the 6 s, 18 p, and 30 d orbitals into first, second, and third weights, respectively. The second line makes a postscript file fplot.ps of the bands, using the plbnds utility and fplot utility.

Note that the bands at the Fermi level are of mixed p and d character.

To compare to the ASA, see this page.

Reconstructed Fe, self-consistent LDA

Build a ctrl and site file for bcc Fe, in the same manner as for Nb, remembering that Fe is spin polarized

blm --nfile --ctrl=ctrl --mag --gw fe
echo p 1 0 0  0 0 1  1 1 0 | lmscell fe --wsite~short~fn=site2

Reconstruction of the Fe lattice is more complicated. First, the reconstruction involves four atoms in a plane (making the ratio of areas in the Nb and Fe planes nearly the same) Also the z axis must be rotated around the (1,0,1) line to the (1,1,1).

echo m 1 1 0 -3 1 0 0 0 1 | lmscell '--scala=sqrt(4/3)' '--rot=(1,0,1)-acos(1/sqrt(3)),y:pi/4' -vfile=2 fe --wsite~short~fn=site3

Inspect site3.fe and see how the atoms stack up:

lmplan -vfile=3 fe

Notes:

• The planes stack this cell normal to the [0,0,1] direction, as in the Nb case. But there are four atoms instead of 3, and they lie in a single plane.
• The first two lattice vectors are identical to those of site3.nb, so the Nb and Fe cells can be stacked on top of one another.
• There is slight mismatch in the lattice constant, 6.24/6.226 = 1.00225 which must be taken into account before Fe can be stacked onto Nb.

site3.nb and site4.fe are conveniently designed to stack into superlattices. To see this in detail, look at the neighbor table from site4.fe

lmchk -vfile=3 fe --shell:tab:r=1.0

There are 8 neighbors in the first shell, 6 in the second as usual for the bcc lattice. The lattice constant has been enlarged by √(4/3), so the NN and 2NN distances are now 0.75 and √(3/4) in units of the new alat, instead of the usual √(3/4) and 1. Now use lmscell to make a 3-layer stack

lmscell -vfile=3 ctrl.fe --stack~file@site3.fe@dup=2~wsite@fn=site12

Run lmplan and lmchk on this stack and confirm (1) it has three layers, and (2) every site has the usual bcc coordination.

lmplan -vfile=12 fe
lmchk -vfile=12 fe --shell:tab:r=1.0

You can remove site12.fe, as it is a temporary file made only to demonstrate the stacker will work as needed

To account for the remaining Nb/Fe lattice mismatch, we will stretch it Fe and keep Nb fixed, because Fe is much thinner. Make a new file site4.fe taking this into account.

lmscell -vfile=3 ctrl.fe -vz=6.24000171/6.226 '--stack~scale=1/z~stretch=z^3~sort@h3~wsite@short@fn=site4'~

Note: You can alternatively copy site4.fe from the box at the top of this tutorial.

The shear reduces the point group operations from 16 to 4; also the 8 NN split into two groups with slightly different bond lengths:

lmplan -vfile=4 fe
lmchk -vfile=12 fe --shell:tab:r=1.0

Construct file basp.fe and remake ctrl.fe. Because of the odd orientation of the supercell, we should make the k mesh different along each axis. To get an estimate of an appropriate partitioning, do the following. lmf can accept a negative argument for tag BZ_NKABC: its absolute value is an estimate for the total number of k points. lmf will then look for an optimal division between nk₁, nk₂, nk₃ that approximates the total number of k points and makes the spacing along the three axes as uniform as it can.

lmchk ctrl.fe --wspec~fn=fpspec
lmfa -vfile=4 ctrl.fe
cat basp0.fe | sed 's/ *PZ.*//' > basp.fe
blm --nit=20 --nfile --nk=-600 --ctrl=ctrl --gw --gmax=8.4 --mag fe
lmf -vfile=4 ctrl.fe --quit=ham

You should see this line

 BZMESH: 198 irreducible QP from 650 ( 10 5 13 )  shift= F F F

This gives us an optimal mesh to select.

lmfa supplies the value 8.4 for gmax.

Remake the ctrl file with this information

blm --nit=30 --nfile --nk=10,5,13 --ctrl=ctrl --gw --gmax=8.4 --mag fe

Strip off the autogenerated specification of a 4d local orbital and make the density self-consistent

mpix -n 16 lmf -vfile=4 ctrl.fe

The self-consistent magnetic moment should come out to about 8.68, or about 2.17 a.u./Fe atom. This is close to the observed magnetic moment, and slightly below a carefully converged LDA calculation. In this case the HF and HK total energies are:

c file=4 mmom=8.6771739 ehf=-10164.1622129 ehk=-10164.1622143

Energy bands of reconstructed Fe

For the Mulliken decomposition make the table or orbital positions

lmf -vfile=4 ctrl.fe --pr55 --quit=ham | grep -A 5 ' Orbital positions in hamiltonian, resolved by l'

and make the bands. Use the same symmetry lines files as for Nb (endpoints don’t touch zone boundaries in this case):

cp syml.nb syml.fe
mpix -n 16 lmf -vfile=4 ctrl.fe  --band~col=1,17,seq=26,51,76~col2=2:4,18:20,seq=27,52,77~col3=5:9,21:25,seq=30,55,80~fn=syml
echo -6,6,5,10 | plbnds -spin1 -fplot~sh~scl=.5~shft=-.1 -ef=0 -scl=13.6 -dat=up -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 bnds.fe
cp fplot.ps feup.ps
open feup.ps
echo -6,6,5,10 | plbnds -spin2 -fplot~sh~scl=.5~shft=.5  -ef=0 -scl=13.6 -dat=dn -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 bnds.fe
cp fplot.ps fedn.ps
open feup.ps

feup.ps and fedn.ps contain majority and minority spin bands of Fe, respectively.

Bands near E_F are dominated by states of d character; there is some p character as well. Note that within a given spin channel the d states cluster into two groups, yielding the classical two-peak structure in the density-of-states.

To compare against the ASA, see this tutorial.

Sanity checks: the following files should be nearly identical to those in the distribution, $HOME/lmf/fp/test/nbfe :

mcx -cmpf~fn1=ctrl.nb~fn2=$HOME/lmf/fp/test/nbfe/ctrl.nb
mcx -cmpf~fn1=site3.nb~fn2=$HOME/lmf/fp/test/nbfe/site3.nb
mcx -cmpf~fn1=ctrl.fe~fn2=$HOME/lmf/fp/test/nbfe/ctrl.fe
mcx -cmpf~fn1=site4.fe~fn2=$HOME/lmf/fp/test/nbfe/site4.fe

Part 2. The Nb(2)/Fe(3)/Nb(2) superlattice

This part constructs files for the Nb(2)/Fe(3)/Nb(2) superlattice, and relaxes the structure. First input files are built; then the structure is relaxed.

The steps in this Part have been bundled into a test script. To confirm what you do is working as it should, you can check your steps against those of the script. Job 2 makes the input files; job 3 performs the self-consistency and lattice relaxation.

$HOME/lmf/fp/test/test.nbfe 2
$HOME/lmf/fp/test/test.nbfe --MPIK=14 3

Part 2 requires the files generated in Part 1. Run Part 1 before doing Part 2, or copy these files from the Questaal distribution for the setup step

rm -f *.nb *.nbfe out.{blm,lmchk,lmscell}
cp $HOME/lmf/fp/test/nbfe/{ctrl.nb,site3.nb,site4.fe,basp.nb,basp.fe} .

lmscell –stack can assemble superlattices from individual site files with equivalent lattice constant and first two lattice vectors. The supercell is constructed by stacking the new atoms above the end of the prior lattice and the third lattice vector to the supercell lattice vectors. Note that site3.nb and site4.fe were specially constructed to have this property. There is no requirement that the third lattice vector from different constituents point in the same direction.

Create a simple Nb(2)/Fe(3)/Nb(2) superlattice

Use the superlattice stack maker to make the following superlattice:

lmscell -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=3~file@site3.nb@dup=2~show~basp@fn=basp.nbfe@basp.nb@basp.fe~show@planes'

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

    1       1:Nb          1    0      0.00000   0.00000   0.00000   0.00000   0.70711  ← layer Nb-L1
    2       1:Nb          1    0      0.70711   0.50000   0.00000   0.00000   0.00000
    3       1:Nb          1    0      1.41421   1.00000   0.00000   0.00000   0.00000
    4       1:Nb          2    0     -0.70711   0.00000   0.00000   0.70711   0.70711  ← layer Nb-L2
    5       1:Nb          2    0      0.00000   0.50000   0.00000   0.70711   0.00000
    6       1:Nb          2    0      0.70711   1.00000   0.00000   0.70711   0.00000
    7       2:Fe          3    0     -1.41421   0.00000   0.00000   1.41421   0.70711  ← layer Fe-1
    8       2:Fe          3    0     -1.41421   0.75000   0.00000   1.41421   0.00000
    9       2:Fe          3    0     -0.70711   1.00000   0.00000   1.41421   0.00000
   10       2:Fe          3    0      0.00000   1.25000   0.00000   1.41421   0.00000
   11       2:Fe          4    0     -1.77016  -0.50338   0.00000   2.03073   0.61651  ← layer Fe-2
   12       2:Fe          4    0     -1.77016   0.24662   0.00000   2.03073   0.00000
   13       2:Fe          4    0     -1.06305   0.49662   0.00000   2.03073   0.00000
   14       2:Fe          4    0     -0.35594   0.74662   0.00000   2.03073   0.00000
   15       2:Fe          5    0     -2.12610  -1.00676   0.00000   2.64724   0.61651  ← layer Fe-3
   16       2:Fe          5    0     -2.12610  -0.25676   0.00000   2.64724   0.00000
   17       2:Fe          5    0     -1.41899  -0.00676   0.00000   2.64724   0.00000
   18       2:Fe          5    0     -0.71189   0.24324   0.00000   2.64724   0.00000
   19       1:Nb          6    0     -2.48205  -1.51014   0.00000   3.26375   0.61651  ← layer Nb-R2
   20       1:Nb          6    0     -1.77494  -1.01014   0.00000   3.26375   0.00000
   21       1:Nb          6    0     -1.06783  -0.51014   0.00000   3.26375   0.00000
   22       1:Nb          7    0     -3.18915  -1.51014   0.00000   3.97086   0.70711  ← layer Nb-R1
   23       1:Nb          7    0     -2.48205  -1.01014   0.00000   3.97086   0.00000
   24       1:Nb          7    0     -1.77494  -0.51014   0.00000   3.97086   0.00000

This command performs the following sequence (a breakdown of the --stack instruction is given in the next section

• Starts with an initial lattice given by site3.nb.
• Stacks the Nb lattice (site3.nb) on top of the original so that the current cell has 6 atoms in 2 planes. The bottom plane will form the “anchor” and become the end lead in the transport. The top will become the frontier plane at the first Nb/Fe interface.
• Stacks three layers (@dup=3) of Fe (site4.fe) on top of the existing structure. At this stage the superlattice has 3×2+3×4 = 18 sites.
• Stacks the Nb lattice site3.nb twice again on top, the first forming the second Fe/Nb interface, and the second layer forming the upper anchor layer. This completes the superlattice with 3×2+3×4+3×2 = 24 sites.
• Assembles file basp.nbfe from basp.nb and basp.fe. Thus the superlattice will use the same basis as the elemental systems.
• Shows the planes in the supercell stack.

lmscell -vfile=3 ctrl.nb --stack
file site3.nb
show
file site4.fe dup=3 dpos=-.01,0,-.01
show
file site3.nb dup=2
show
basp@fn=basp.nbfe@basp.nb@basp.fe
wsite short fn=sitein

lmscell -vfile=3 ctrl.nb '--stack~file@site3.nb@dup=2~show~wsite@short@fn=site4'

lmchk --shell -vfile=1 nb
lmchk --shell -vfile=4 nb

This lattice has 12 Nb and 12 Fe atoms, arranged in the following planes.

                                    Table 1
  Layer    Nb-L1     Nb-L2     Fe-1       Fe-2     Fe-3      Nb-R2     Nb-R1
  sites    1-3       4-6       7-10       11-14    15-18     19-21     22-24
  force    small     large     large      small    large     large     small

Strategy to make practical Nb/Fe/Nb superlattices

The lattice just generated is too simple for the following reasons:

1. Spacings at the two interfaces are not symmetrically distributed.
2. Frontier atoms (at the interface) will have large forces, as Table 1 shows.
3. The superlattice was constructed so that the volume should be 12×120.669456+12×78.907090 (12× sum of elemental Fe + Nb volumes). You can confirm from the output that this is the volume of the superlattice. But there is no reason why this must be so; in fact all three components of the third lattice vector P₃ can change.

The rigid lattice must be relaxed, which the following section accomplish within some constraints. Constraints are imposed so the superlattice is built in a lego-like fashion, layer by layer. That way construction can be automatic and easily generated for any number of layers. It will be particularly important for Landauer-Buttiker transport, where there the end layers become semi-infinite and must have the structure of the bulk Nb crystal.

This is greatly facilitated by restricting atomic relaxation to atoms on frontier planes. It is an approximation, but it is a reasonable one; and it can be checked by looking at the forces ramaining after the constrained relaxation. 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.

1. Ideal Nb(2)/Fe(3)/Nb(2) superlattice

If you do not already have these files in your working directory, you can run this step and the next one independently of prior ones by copying these from the Questaal distribution:

rm -f *.nb *.nbfe out.{blm,lmchk,lmscell}
cp $HOME/lmf/fp/test/nbfe/{ctrl.nb,site3.nb,site4.fe,basp.nb,basp.fe} .

The output of the previous lmscell instruction shows that the spacing between Nb-L2 and Fe-1 is Δh=0.70711=1/√2, while the spacing between Fe-3 and Nb-R2 is significantly smaller, Δh=0.61651. This spacing should be distributed equally across the interfaces. We do this by modifying the lmscell instruction to shift the Fe planes by half of this difference. Since the planes are stacked on [001], we can do it as follows:

lmscell -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp1=0 -vdp2=0 -vdp3=0 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=3@dpos=dp1/2,dp2/2,-dz+dp3/2~~file@site3.nb@dup=2@dpos=dp1,dp2,dp3~show~basp@fn=basp.nbfe@basp.nb@basp.fe~shorp3~pad=dp1,dp2,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~'

lmscell makes sitein.nb and basp.nbfe. We could use lmfa to autogenerate the latter, but for consistency we will stick with the Fe and Nb basis sets calculated for their respective elements.

2. Crude estimate for relaxation of the Nb(2)/Fe(3)/Nb(2) superlattice

You need to perform step 1 before doing this one.

In the preceding section lmscell provided for a general xyz translation of the relative positions of Nb and Fe blocks, (dp1,dp2,dp3), though all were set to 0 there. In the steps that follow we take dp1=dp2=0, but construct three site files for dp3=0, −0.1, −0.2, in preparation for determining the optimal value of dp3, accomplished in section 3c. Files site0, site10, and site20 are made corresponding to dp3=0, −0.1, −0.2 without relaxation, as well as a crude estimate for the relaxation for dp3=0 (dpos30).

Table 1’s second line shows the relative magnitude forces expected for the ideal, rigid lattice, as would apply if the internuclear forces were dominate by nearest neighbors. The rigid lattice must be relaxed, which this section accomplishes within some constraints. The idea is structure the stacks so the superlattice is built in a lego-like fashion, layer by layer. That way construction can be automatic and easily generated for any number of layers. This is greatly facilitated by restricting atomic relaxation to atoms on frontier planes. This is an approximation, but it is a reasonable one; and it can be checked by looking at the forces ramaining after the constrained relaxation. 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 lattice can be straightforwardly relaxed from the rigid structure generated in the previous section We can improve on this by shifting frontier atoms around to push apart pairs with small bond lengths. This can be accomplished with the sphere overlap minimizer utility in lmchk, which penalizes small bond lengths. This should distribute the forces (reduce the largest ones, while increasing some of the smallest ones).

To do this, make a temporary input file from sitein.nb generated in section 1, where the spheres are allowed to overlap by as much as 16%. Then use the --mino utility in lmchk and allow a subset of atoms to move, to minimize a function of the sphere overlaps.

cp pos.nb pos0.nbfe
cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=.16 --nosort --shorten=no --wpos=pos0 --ctrl=ctrl nbfe
lmchk nbfe --mino~dxmx=.001~xtol=.0001~nkill=10~maxit=200~4:10,15:21 --wpos=pos

A rather larger sphere overlap of 0.16 allows atomic sphere radii to overlap by as much as (16%), which supplies plenty of scope for the overlap minimizer to reduce the penalty function. Relaxed positions are written to file pos.nbfe.

As for arguments to --mino, dxmx, xtol, nkill and maxit are explained on the command line arguments web page. They control details of the (very nonlinear) minimization process.\ 4:10,15:21 is a site list specifying which sites are allowed to move. Correspondence between sites and layers are shown in the table above.

It is convenient for later steps, to work with shifts in position relative to the rigid lattice. Make this file:

mcx -f3f12.6 pos.nbfe pos0.nbfe -- > dpos30.nbfe

The ‘3’ in the name refers to 3 Fe layers; the ‘0’ refers to the crude estimate for relaxation. These extra labels are there to disambiguate other position files.

Now we are in a position to make the ctrl and site files. The next steps build the latter. For bookeeping reasons we will create the site file for the ideal, rigidly stacked lattice, and when needed, read in the estimates for the relaxed positions from pos,nbfe. Create site0.nbfe :

cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=8 --mag --frzwf --gw nbfe '--wsite~ib=\!4:10,15:21~rlx=000~long~fn=sitex' >> out.blm
cp sitex.nbfe site0.nbfe

blm reads structural data from sitein.nbfe, and generates actrl.nbfe and site0.nbfe. (We designate the name site0 because this corresponds to 0 shift in the Nb/Fe interplanar spacing relative to the average. Soon site10 and site20 will be made, corresponding respectively to shifts of −-0.1 and −-0.2.)

%ifdef file | file==0
  file=   site{file}
%endif

We could use actrl.nbfe blm generated just now for the input file, but the spheres would be very small owing to very close spacing between some Fe and Nb. Instead we construct the ctrl file again, reading the approximate relaxed positions from dpos30.nbfe.

blm --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=9 --mag --frzwf --gw nbfe --rdpos=dpos30 --rsham~el=-.1,-0.9~rscut=1.6 --dhftol=45 --molstat~xtol=0~gtol=0.01
cat actrl.nbfe | sed s:'\(ATOM=Fe.*\):\1 MMOM=0,0,2:' > ctrl.nbfe

Arguments following nbfe are different from before; they are explained in the drop-down box above.

blm doesn’t know about magnetic moments so we have to modify the autogenerated actrl.nbfe to give it a reasonable starting guess for the moment. Bulk Fe has a moment of 2.2 μ_B; the second line adds a tag to the species category supplying a guess of 2.0 μ_B for all the Fe atoms.

Repeat the preceding commands to make ideal site files with Nb/Fe interplanar spacing reduced by 0.1 and 0.2, respectively. These commands make site10.nbfe and site20.nbfe.

lmscell -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp3=-.10 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=3@dpos=0,0,-dz+dp3/2~~file@site3.nb@dup=2@dpos=0,0,dp3~setlbl@spid=Nbx@targ=4:6,nbas-5:nbas-3~shorp3~pad=0,0,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~' >> out.lmscell
cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=8 --mag --frzwf nbfe '--wsite~ib=\!4:10,15:21~rlx=000~long~fn=sitex' >> out.blm
cp sitex.nbfe site10.nbfe
lmscell -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp3=-.20 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=3@dpos=0,0,-dz+dp3/2~~file@site3.nb@dup=2@dpos=0,0,dp3~setlbl@spid=Nbx@targ=4:6,nbas-5:nbas-3~shorp3~pad=0,0,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~' >> out.lmscell
cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=8 --mag --frzwf nbfe '--wsite~ib=\!4:10,15:21~rlx=000~long~fn=sitex' >> out.blm
cp sitex.nbfe site20.nbfe

Sanity checks: the following files should be nearly identical to those in the distribution, $HOME/lmf/fp/test/nbfe :

mcx -cmpf~fn1=ctrl.nbfe~fn2=$HOME/lmf/fp/test/nbfe/ctrl.nbfe
mcx -cmpf~fn1=site0.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site0.nbfe
mcx -cmpf~fn1=dpos30.nbfe~fn2=$HOME/lmf/fp/test/nbfe/dpos30.nbfe
mcx -cmpf~fn1=basp.nbfe~fn2=$HOME/lmf/fp/test/nbfe/basp.nbfe
mcx -cmpf~fn1=site10.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site10.nbfe
mcx -cmpf~fn1=site20.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site20.nbfe

3 Relaxed Nb(2)/Fe(3)/Nb(2) superlattice

This step does the following:

1. Makes the density self-consistent for crudely relaxed positions generated in test 2, Nb/Fe interplanar spacing taken as average of Nb and Fe (site0.nbfe)
2. Relaxes the lattice to minimize forces
3. Does the same for the Nb/Fe interplanar spacing reduced −0.1 (site10.nbfe)
4. Does the same for the Nb/Fe interplanar spacing reduced −0.2, which is near the minimum-energy point (site20.nbfe)

Setup: run steps 1 and 2 of Part 2, or copy necessary files from the Questaal distribution

rm -f *.nbfe spec.ni spec.nb site4.ni out.0 out.10 out.20
cp $HOME/lmf/fp/test/nbfe/{ctrl.nbfe,site0.nbfe,site10.nbfe,site20.nbfe,basp.nbfe,dpos30.nbfe} .

3a. Self-consistency

Before relaxing the lattice it is useful to obtain a self-consistent density to get a sense of idea of what is transpiring at the interface and the forces there.
Note this step takes about 15 minutes to execute. The number of MPI processes cannot exceed 14.

rm -f mixm.nbfe hssn.nbfe log.nbfe save.nbfe
lmfa -vfile=0 ctrl.nbfe
mpirun -n 14 lmf ctrl.nbfe -vfile=0 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos --rdpos=dpos30 --shorten=2,2,1 --wpos=posx >> out.0
cp posx.nbfe pos0.nbfe

Inspect the forces in the last iteration (out.0):

                                           Table 2
  ib           estatic                  eigval                    total
   1      3.39     9.66    34.57     -6.42   -14.66   -60.34     -3.03    -4.99   -25.78
   2     -3.07    -6.81    34.65      5.34     8.88   -61.64      2.27     2.07   -26.99   ← layer Nb-L1
   3     -0.75    -1.01    27.46      0.00     2.36   -35.75     -0.67     1.35    -8.28

   4    -29.70   -20.99   -78.37     33.25    15.17   108.43      3.55    -5.82    30.06
   5      0.10     0.48  -182.65      0.00     0.12   240.87      0.00     0.60    58.22   ← layer Nb-L2
   6     30.71    22.91   -82.57    -34.05   -17.19   114.25     -3.35     5.71    31.68

   7      1.64    -1.35     3.05     -1.56     1.07   -16.26      0.00    -0.28   -13.21
   8     -5.62   -31.16    39.33      1.55    35.42   -60.10     -4.07     4.26   -20.78   ← layer Fe-1
   9     -1.30    -1.20    80.03      1.78     0.67   -90.98      0.47    -0.53   -10.95
  10      4.23    29.18    38.91     -0.48   -34.46   -53.70      3.74    -5.28   -14.79

  11      3.38     3.21   -22.98      0.83    -0.34    34.78      4.21     2.87    11.80
  12      2.93     1.58     1.12     -3.14    -1.93     0.00     -0.22    -0.34     1.17   ← layer Fe-2
  13     -0.99    -1.30   -20.08     -4.31    -1.36    30.92     -5.30    -2.67    10.84
  14     -1.06    -2.40   -21.56     -1.21     2.74    42.84     -2.27     0.35    21.28

  15     17.44    -1.90    20.65    -24.60     3.73   -21.01     -7.16     1.84    -0.36
  16      0.00     0.41    46.69     -1.66    -0.42   -65.50     -1.66     0.00   -18.81   ← layer Fe-3
  17    -24.83     3.59    29.84     39.46    -6.28   -36.62     14.63    -2.69    -6.78
  18     -1.70    -2.98    47.76      3.10     4.28   -56.41      1.40     1.30    -8.64

  19      1.05    -2.00     2.53     -2.76     4.02   -14.12     -1.71     2.02   -11.59
  20      6.47    14.87    23.61    -12.21   -33.38   -42.69     -5.74   -18.50   -19.08   ← layer Nb-R2
  21     -1.08   -11.54    10.66      2.62    29.99   -22.95      1.54    18.45   -12.29

  22     -1.20     2.24    -9.56      0.86    -0.52    15.88     -0.34     1.72     6.33
  23      2.43     1.08   -11.73      0.40    -6.66    23.74      2.83    -5.59    12.01   ← layer Nb-R1
  24     -2.53    -4.57   -11.38      3.22     8.77    26.29      0.69     4.19    14.91

The largest force (site 5) is about 58 mRy/a.u. We can expect the frontier atoms would have much larger forces had the calculation been done for the rigid lattice instead. (To confirm this redo the calculation omit the --rdpos=dpos30 switch. You should find that the largest force is ~80 mRy/a.u. (site 9), while the end layers Nb-L1 and Nb-R1 would have smaller forces.)

3b. Constrained relaxation

(Complete step 3a before doing this one.)

Use lmf to perform the (constrained) lattice relaxation. Note that the site files were constructed so that only atoms in the frontier layers are allowed to relax.
Note this step takes about 2 hours to execute.

mpirun -n 14 lmf ctrl.nbfe -vfile=0 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos -vfrzwf=t -vminx=10 --wpos~fn=posx~shorten@0,0,1 >> out.0
cp posx.nbfe pos0.nbfe

It uses a Broyden algorithm to relax the structure for ten iterations. It writes the positions file for what would be iteration six into posx.nbfe. It is copied to pos0.nbfe to make the naming convention consistent with site0.nbfe.

Try grep RELAX out.0 to see how the gradient evolved in successive relaxation steps:

 RELAX: (warning) failed to read Hessian; set to unity
 RELAX:  completed Broyden iter 1;  max shift=0.01000  |g|=0.0907
 RELAX:  completed Broyden iter 2;  max shift=0.01000  |g|=0.0824
 RELAX:  completed Broyden iter 3;  max shift=0.01000  |g|=0.0716
 RELAX:  completed Broyden iter 4;  max shift=0.01000  |g|=0.0616
 RELAX:  completed Broyden iter 5;  max shift=0.01000  |g|=0.049
 RELAX:  completed Broyden iter 6;  max shift=0.01000  |g|=0.0374
 RELAX:  completed Broyden iter 7;  max shift=0.01000  |g|=0.0232
 RELAX:  completed Broyden iter 8;  max shift=0.01000  |g|=0.0123
 RELAX:  converged to tolerance;  max shift=0.00292  |g|=0.00405

|g| is the length of the 42 component force “vector”; it drops smoothly as relaxation proceeds. Evolution of total energy also shows a steady drop, tapering off as iterations proceed (egrep -E '^[Ccx]' out.0)

c file=0 conv=1e-4 convc=.001 beta=.1 mmom=23.6913094 ehf=-122072.4593395 ehk=-122072.455455
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=23.7494935 ehf=-122072.4593708 ehk=-122072.4581496
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=23.7626436 ehf=-122072.4685081 ehk=-122072.4680317
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=23.7189084 ehf=-122072.4750251 ehk=-122072.4749764
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=23.611728 ehf=-122072.4803413 ehk=-122072.4802643
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=23.5291273 ehf=-122072.4847483 ehk=-122072.4846855
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=23.4060968 ehf=-122072.488403 ehk=-122072.488258
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=23.2702541 ehf=-122072.4908124 ehk=-122072.4907438
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=23.1521162 ehf=-122072.4923798 ehk=-122072.4922336
c file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=22.9615162 ehf=-122072.4927374 ehk=-122072.4926169
C file=0 conv=1e-4 convc=.001 beta=.1 frzwf=1 minx=10 mmom=22.9615162 ehf=-122072.4927374 ehk=-122072.4926169

Maximum Harris force = 47.4 mRy/au (site 5)  Max eval correction = 8.2   ← iter 1 rlx 0
Maximum Harris force = 47.9 mRy/au (site 5)  Max eval correction = 7.3   ← iter 2 rlx 0
Maximum Harris force = 117 mRy/au (site 5)  Max eval correction = 748.4  ← iter 1 rlx 1
...
Maximum Harris force = 45.1 mRy/au (site 5)  Max eval correction = 21    ← last iter rlx 1
...
Maximum Harris force = 23.4 mRy/au (site 14)  Max eval correction = 40.5 ← last iter rlx 5
...

The maximum force on an individual atom can be seen from grep 'Maximum Harris force' out.0

 Maximum Harris force = 56.7 mRy/au (site 5)  Max eval correction = 6.9  ← before relaxation
...
 Maximum Harris force = 25.3 mRy/au (site 14)  Max eval correction = 3.3 ← after relaxation, iteration 5

The force on site 14 is a bit on the large side; it will not change with further relaxation because it is in the constrained middle layer. We do not relax this layer in keeping with the “lego” strategy.

Final site positions were written to file posx.nbfe and copied to pos0.nbfe. These positions do not correspond to the last iteration were forces were calculated, but the ersatz next iteration. It is not necessary that you do this, but you can run lmf at these new positions and confirm that the relaxation was nearly complete.

rm -f mixm.nbfe hssn.nbfe log.nbfe save.nbfe
mpix -n 14 lmf ctrl.nbfe -vfile=0 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos -vfrzwf=t -vminx=0 --rpos=pos0 --rs=11,1,1 >> out.0 &

For the meaning of --rs=11,1,1, run lmf --h or see command line arguments documentation.

It is worthwhile to check that sphere overlaps after relaxation. Do lmchk ctrl.nbfe -vfile=0 --rpos=pos0 --pr20 and you should see this line

 OVMIN, 468 pairs:  fovl = 3.52798e-13   <ovlp> = 0.8%   max ovlp = 0.8%

There is a small overlap. Thanks to Questaal’s three-component augmentation spheres can overlap up to about 5% with very small error, and sometimes a bit larger (see Sec 3.6 of the Questaal methods paper).

3c. Optimizing the Nb/Fe layer spacing

(Complete step 3b before doing this one.)

There is no good justification to assume that the Nb-Nb spacing is 0.61651/2+0.70711/2. Instead, this spacing should be determined by minimizing the total energy. In fact all three components of the third lattice vector P₃ should be varied; however components 1 and 2 lie normal to the propagation direction, do not affect the volume substantially and are less important than the z component.

The lmscell instruction that assembled the site file of the superlattice is set up modify the z component of P₃ through the variable p3. To find p3 that minimizes the total energy we can repeat steps 3a and 3b with changes for different values of p3 until the minimum energy is found.

Also since we have found one relaxed structure from DFT (p3=0), a better choice for the relaxation can be made using that information. A reasoanble guess is to reduce the third component of each basis vector the dp3=−0.1, scaled in proportion to its height. This exactly preserves the relation between the top layer and the bottom layer, as it should since periodic boundary conditions are used. Do this with the mcx calculator

         mcx -f3f14.8 -vdp3=-.1 pos0.nbfe -av:nr,3 p3max -e3 0 0 'dp3*x3/p3max' pos0.nbfe -+ > posnow.nbfe

-av:nr,3 p3max assigns scalar p3max to the element (nr,3) read from pos0.nbfe. It creates a shift array (0,0,dp3*x3/p3max), adds it to pos0.nbfe, and writes the result to posnow.nbfe.

Repeat the steps in 3a, 3c, substituting site10.nbfe for site0.nbfe and read site positions from posnow.nbfe

rm -f mixm.nbfe hssn.nbfe log.nbfe save.nbfe rst.nbfe
mpirun -n 14 lmf ctrl.nbfe -vfile=10 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos --rpos=posnow --shorten=2,2,1 > out.10
mpirun -n 14 lmf ctrl.nbfe -vfile=10 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos -vfrzwf=t -vminx=10 --rpos=posnow --shorten=no --rs=11,1,1 --wpos~fn=posx~shorten@0,0,1 >> out.10
cp posx.nbfe pos10.nbfe

You can confirm that the ansatz encoded in posnow.nbfe was pretty reasonable : egrep -E '^[Ccx]' out.10 shows the energy gained by relaxation was only 3 mRy.

Do the steps once again for dp3=−0.2

mcx -f3f14.8 -vdp3=-.2 pos0.nbfe -av:nr,3 p3max -e3 0 0 'dp3*x3/p3max' pos0.nbfe -+ > posnow.nbfe
rm -f mixm.nbfe hssn.nbfe log.nbfe save.nbfe rst.nbfe
mpirun -n 14 lmf ctrl.nbfe -vfile=20 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos --rpos=posnow --shorten=2,2,1 > out.20
mpirun -n 14 lmf ctrl.nbfe -vfile=20 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos -vfrzwf=t -vminx=10 --rpos=posnow --shorten=no --rs=11,1,1 --wpos~fn=posx~shorten@0,0,1 >> out.20
cp posx.nbfe pos20.nbfe

To find the optimal p3 we need to make a table of the relaxed total energies. Do this with First look at the converged magnetic moment and total energy for each of the p3.

egrep -E '^[Ccx]' out.0 out.10 out.20 | grep :C

If the Fe local magnetic M retained its bulk value, we would expect M=12×2.2 μ_B = 26.4 μ_B for the superlattice, since it has 12 Fe atoms. The actual moment is somewhat less than that. To see a table of individual moments, do

lmf ctrl.nbfe -vfile=20 --rsedit~rs~a

                                              Table 3
   site          z     rmt     nr   a  lmxl   qtrue     mtrue     q2        m2
   1 Nb        41.0  2.40242  425 0.025  4   2.98558   0.02844   1.88649   0.01674
   2 Nb        41.0  2.40242  425 0.025  4   2.98285   0.00444   1.88481   0.00296   ← layer Nb-L1
   3 Nb        41.0  2.40242  425 0.025  4   3.01698   0.02203   1.91298   0.01299					 
																				                       
   4 Nbx       41.0  2.31392  421 0.025  4   2.82051  -0.12580   1.77776  -0.07388					 
   5 Nbx       41.0  2.31392  421 0.025  4   2.80527  -0.07955   1.76174  -0.04604   ← layer Nb-L2
   6 Nbx       41.0  2.31392  421 0.025  4   2.81867  -0.12256   1.77689  -0.07225					 
																				                       
   7 Fe        26.0  2.08458  377 0.025  4   6.37478   1.96629   2.29314   0.41430					 
   8 Fe        26.0  2.08458  377 0.025  4   6.44111   1.69123   2.35841   0.35192   ← layer Fe-1 
   9 Fe        26.0  2.08458  377 0.025  4   6.56111   1.39020   2.47715   0.28296					 
  10 Fe        26.0  2.08458  377 0.025  4   6.43785   1.77174   2.35584   0.37183					 
																				                       
  11 Fe        26.0  2.08458  377 0.025  4   6.43294   2.05272   2.35206   0.42267					 
  12 Fe        26.0  2.08458  377 0.025  4   6.39846   2.24441   2.31373   0.46674   ← layer Fe-2 
  13 Fe        26.0  2.08458  377 0.025  4   6.41570   2.23589   2.33349   0.46471					 
  14 Fe        26.0  2.08458  377 0.025  4   6.46188   2.11634   2.38500   0.43881					 
																				                       
  15 Fe        26.0  2.08458  377 0.025  4   6.44852   1.74429   2.37155   0.36311					 
  16 Fe        26.0  2.08458  377 0.025  4   6.43273   1.72522   2.35703   0.35984   ← layer Fe-3 
  17 Fe        26.0  2.08458  377 0.025  4   6.43686   1.77967   2.35954   0.37103					 
  18 Fe        26.0  2.08458  377 0.025  4   6.45338   1.75519   2.38119   0.36365					 
																				                       
  19 Nbx       41.0  2.31392  421 0.025  4   2.77547  -0.11428   1.74272  -0.06709					 
  20 Nbx       41.0  2.31392  421 0.025  4   2.76843  -0.10233   1.73722  -0.06041   ← layer Nb-R2
  21 Nbx       41.0  2.31392  421 0.025  4   2.76949  -0.08910   1.73796  -0.05284					 
																				                       
  22 Nb        41.0  2.40242  425 0.025  4   2.99457   0.00973   1.89460   0.00616					 
  23 Nb        41.0  2.40242  425 0.025  4   2.99594   0.00507   1.89518   0.00359   ← layer Nb-R1
  24 Nb        41.0  2.40242  425 0.025  4   2.96189  -0.01609   1.86601  -0.00846

Two features stand out:

1. qtrue (charge inside the augmentation sphere) shows a slight transfer of charge from the frontier Nb layer to the frontier Fe layer. This is the dipole that must be set up to align the bulk Fe and Nb Fermi levels (chemical potentials). Fe has a deeper work function than Nb, so we expect the charge to flow from Nb to Fe.
2. The system magnetisation drops from about 23 for p3=0, dropping to 21.7 for p3=−0.2. This is to be expected, since as Fe planes are pushed against Nb planes, the coupling increases, thus increasing the kinetic energy which competes with the magnetic exchange interactions. The latter become less effective and M drops. It is apparent from column mtrue that the drops on layers Fe-1 and Fe-3 are substantial, while the middle (Fe-2) layer is bulk-like.

Extract a table total energy as a function of p3

egrep -E '^[Ccx]' out.0 out.10 out.20 | grep :C | vextract . file ehf > dat

You should see

  0 -122072.4927374
  10 -122072.5281779
  20 -122072.5391974

The first argument corresponds to x1=−100×p3; the second to the total energy. Fit p3 that minimizes the energy (in a parabolic approximation)

echo q | pfit 2 dat

The minimum energy falls at x1=19.5, showing that dp3=−0.2 is close to the minimum point. We will take this structure as a good approximation to the fully relaxed one.

3d. LDA Energy Bands and layer resolved DOS of the superlattice

The large unit cell ensures the bands will be spaghetti-like, so it is useful to know whether a band is mostly Fe-like or Nb-like. When making bands we can assign color weights to distinguish orbital character. For this a list of orbitals in the superlattice is needed:

lmf ctrl.nbfe -vfile=20 -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos --rpos=posnow --shorten=2,2,1  --quit=ham --pr50 | grep -A 50 'Orbital positions'

Orbitals 1:180 and 541:720 belong to Nb while 181:540 belong to Fe. Make bands with color weights

lmchk ctrl.nbfe -vfile=20 '--syml~n=21~q=0,1,0,0,0,0,1/sqrt(2),0,0'
mpix -n 16 lmf ctrl.nbfe -vfile=20 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 --band~col=1:180,541:720~col2=181:540~fn=syml

Draw the 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.nbfe
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.nbfe
open fplot.ps

The Fe states cluster into two groups bands, with a pseudogap separating them, as occurs with elemental Fe.

Rather remarkably, there is little Fe character at E_F: both groups both fall E_F for the majority spin, and for the minority spin E_F falls in the pseudogap.

The Nb-derived (red) bands clearly show families of free-electron-like (approximately parabolic) bands. These bands are likely to carry most of the transport. There are multiple parabolic bands closely spaced in k: this spread Δk is a consequence of the reconstruction. In reality the reconstruction is likely somewhat disordered. However, disorder on top of the reconstruction already there is not likely to increase Δk significantly. Thus Δk seen in the figure approximately represents the broadening in k because of nonidealities at the interface. Note that Δk is relatively small compared to the entire Brillouin zone, particularly near E_F. Thus, the often-used approximation that k is completely smeared out over the Fermi surface is pretty far removed from the actual situation.

Another way to see this is to inspect the density-of-states resolved by the 6 layers of the superlattice. To do this, use lmf with the --pdos switch as shown below . Note that a finer k mesh is used to render the DOS smooth. lmdos combines the raw output of lmf into data readable by pldos, which makes input for fplot to draw a postscript file.

mpix -n 16 lmf ctrl.nbfe -vnk1=18 -vnk2=8 -vfile=20 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 '--pdos~mode=0~group=1:3;4:6;7:10;11:14;15:18;19:21;22:24' --quit=rho
lmdos ctrl.nbfe -vnk1=18 -vnk2=8 -vfile=20 --shorten=no '--pdos~mode=0~group=1:3;4:6;7:10;11:14;15:18;19:21;22:24' --dos:npts=501:window=-.8,.5
echo 5 5.99 -4 4 | pldos -ef=0 -esclxy=13.6 -fplot~open~tmy=2~dmin=3~xl=E -lst='1;3;5;7;9;11;13' -lst2 dos.nbfe
fplot -f plot.dos
open fplot.ps

Sanity checks: the following files should be nearly identical to those in the Questaal distribution, $HOME/lmf/fp/test/nbfe :

mcx -cmpf~fn1=pos0.nbfe~fn2=$HOME/lmf/fp/test/nbfe/pos0.nbfe
mcx -cmpf~fn1=pos10.nbfe~fn2=$HOME/lmf/fp/test/nbfe/pos10.nbfe
mcx -cmpf~fn1=pos20.nbfe~fn2=$HOME/lmf/fp/test/nbfe/pos20.nbfe

Part 3. Longer period superlattices

This section uses relaxation from Part 2 to make Nb(2)/Fe(n)/Nb(2) superlattices. n can be any integer larger than 3.

Setup: run steps in Section 3 of Part 2, or copy the following files from the Questaal distribution

touch ctrl.nbfe; rm -f *.nbfe out.lmchk out.lmscell out.lmfa out.lmf
cp $HOME/lmf/fp/test/nbfe/{ctrl.nbfe,site0.nbfe,site10.nbfe,site20.nbfe,basp.nbfe,dpos30.nbfe} .

Create file dpos20.nbfe that holds the relaxation of the basis relative to the unrelaxed Nb(2)/Fe(3)/Nb(2) structure

lmchk ctrl.nbfe -vfile=20 --wpos=posx >> out.lmchk
mcx -f3f14.8 pos20.nbfe posx.nbfe  -- > dpos20.nbfe

The “lego” strategy is to create a site file for the Nb(2)/Fe(n)/Nb(2) unrelaxed structure as was done for the Nb(2)/Fe(3)/Nb(2) structure in Part 2 The relaxation was confined to frontier planes, and the change in positions is contained in dpos20.nbfe. By inserting rows of zeros corresponding to the insertion of extra planes relative to the Nb(2)/Fe(3)/Nb(2) reference, the analog of dpos20.nbfe is easily made.

In the steps below, you can set n to any number. For concreteness, select n=5 for the Nb(2)/Fe(5)/Nb(2) structure and make site5.nbfe and pos5.nbfe:

set n=5     ← Use this in tcsh. The corresponding bash instruction is :  n=5
lmscell -vn=$n -vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2 -vdp3=-.20 -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dup=n@dpos=0,0,-dz+dp3/2~~file@site3.nb@dup=2@dpos=0,0,dp3~setlbl@spid=Nbx@targ=4:6,nbas-5:nbas-3~shorp3~pad=0,0,dp3~shorten@mode=2,2,1~wpos@fn=pos~wsite@fn=sitein~' >> out.lmscell
cp sitein.nb sitein.nbfe
blm -vn=$n --nfile --ewald=1d-12 --rdsite --omax=-.07 --nosort --shorten=no --nit=20 --nk=9,3,1 --gmax=8 --mag --frzwf nbfe '--wsite~ib=!4:10,11+4*1:11+4*n+6~rlx=000~long~fn=sitex' >> out.blm
cp sitex.nbfe site5.nbfe
mcx -f3f14.8 -vn=$n -vm=n-2 [ 0:'4*m'-1 '-array[1,3]' 0,0,0 '?i>1' -rcat ] -a ins dpos20.nbfe  -split x 1,11,15,nr+1 1,nc+1 x11 ins x31 -rcat -rcat > dposx.nbfe
cp dposx.nbfe dpos5.nbfe

Sanity checks: the following files should be nearly identical to those in the Questaal distribution, $HOME/lmf/fp/test/nbfe :

mcx -cmpf~fn1=site5.nbfe~fn2=$HOME/lmf/fp/test/nbfe/site5.nbfe
mcx -cmpf~fn1=dpos5.nbfe~fn2=$HOME/lmf/fp/test/nbfe/dpos5.nbfe

Additional exercises

You can make the densities of larger period superlattices self-consistent, and inspect the evolution with layer spacing

Choose some n and do

lmfa -vfile=$n ctrl.nbfe
mpix -n 14 lmf ctrl.nbfe -vfile=$n -vconv=1e-4 -vconvc=1e-3 -vbeta=.1 --rhopos --rdpos=dpos$n --shorten=2,2,1 > out.$n

You should find

1. the forces for sites allowed to relax at the frontier layer remain small (<10 mRy/a.u.); the constrained atoms in the constrained layers adjacent to them have somewhat larger forces than they did in the Nb(2)/Fe(3)/Nb(2) case, and forces in layers removed at least once from frontier layers have small forces.

2. the magnetic moment is approximately linear in n, with a change in M between 8 and 9 as a layer is added. The value of M at the frontier changes a little as n increases, but M in atoms removed from the frontier are consistently between 2.15 and 2.25.

3. the total energy is very nearly linear in n.

Footnotes and references

A detailed description of Questaal and lmf can be found in the following reference:
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, Comp. Phys. Comm. 249, 107065 (2020).