Questaal Home
Navigation

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

Bulk Nb

Copy and paste init.nb file to your working directory.

blm --nfile --gmax=5.5 --nk=12,12,12 --ctrl=ctrl --gw nb
echo p 1 0 0 0 0 1 1 1 0 | lmscell nb --wsite~short~fn=site2
echo m 1 0 0 0 3 0 0 0 1 | lmscell -vfile=2 nb --wsite~short~fn=site3
sed -i 's/ *PZ.*//' ctrl.nb
lmfa -vfile=3 ctrl.nb
cat basp0.nb | sed 's/ *PZ.*//' > basp.nb
mpix -n 16 lmf -vfile=3 ctrl.nb

Energy bands

lmchk ctrl.nb -vfile=3 --syml~n=31~q=0,1,0,0,0,0,1/2,0,-1/2
mpix -n 16 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,seq=35,65~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

Bulk Fe

Copy and paste init.fe file to your working directory.

blm --nfile --ctrl=ctrl --mag --gw fe
echo p 1 0 0  0 0 1  1 1 0 | lmscell fe --wsite~short~fn=site2
echo m 1 1 0 -3 1 0 0 0 1 | lmscell '--scala=sqrt(4/3)' '--rot=(1,0,1)-acos(1/sqrt(3))' -vfile=2 fe --wsite~short~fn=site3
# lmscell -vfile=3 ctrl.fe -vz=6.24000171/6.226 '--stack~scale=1/z~stretch=z^3~sort@h3~wsite@short@fn=site4'
lmfa -vfile=4 ctrl.fe
cat basp0.fe | sed 's/ *PZ.*//' > basp.fe
blm --nit=30 --nfile --nk=10,5,13 --ctrl=ctrl --gw --gmax=8.4 --mag fe
sed -i 's/ *PZ.*//' ctrl.fe
mpix -n 16 lmf -vfile=4 ctrl.fe

Energy bands

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

Nb/Fe superlattice

Make sitein.nbfe, basp.nbfe, ctrl.nbfe, site.nbfe

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~wsite@fn=sitein~'
cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.01 --nosort --shorten=no --nit=20 --nk=9,4,2 --gmax=8 --mag --frzwf nbfe
cat actrl.nbfe | sed s:'\(ATOM=Fe.*\):\1 MMOM=0,0,2:' > ctrl.nbfe

Make self-consistent for dpos=-.02,0,-.02

lmscell -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dpos=-.02,0,-.02@dup=3~file@site3.nb@dup=2~show~basp@fn=basp.nbfe@basp.nb@basp.fe~wsite@fn=sitein~'
cp sitein.nb site2.nbfe
rm -f mixm.nbfe log.nbfe rst.nbfe
lmfa ctrl.nbfe
rm -f rst.nbfe mixm.nbfe log.nbfe
mpix -n 13 lmf ctrl.nbfe -vbeta=.1 --rhopos -vfrzwf=t -vnit=10 -vfile=2 > out.m02 &
cp rst.nbfe rst.nbfe.bk

Fully relax the lattice

blm --dhftol=45 --conv=1e-4 --convc=1e-3 --molstat --nfile --ewald=1d-12 --rdsite --omax=-.01 --nosort --shorten=no --nit=20 --nk=9,4,2 --gmax=8 --mag --frzwf --wsite@fn=null nbfe
cat actrl.nbfe | sed s:'\(ATOM=Fe.*\):\1 MMOM=0,0,2:' | sed s:'XTOL=.001 GTOL=0:XTOL=0 GTOL=0.01:'  > ctrl.nbfe
rm -f save.nbfe mixm.nbfe hssn.nbfe
mpix -n 13 lmf ctrl.nbfe -vfile=2 -vbeta=.1 --rhopos -vfrzwf=t -vnit=20 -vminx=5 --wpos=pos > out.dorelax &
cp rst.nbfe rst.nbfe.bk2

Adjust positions in the relaxed lattice to remove translation vectors, confirm self-consistent

lmf ctrl.nbfe -vfile=2 -vbeta=.1 --rhopos -vfrzwf=t -vnit=0 --rs=1,0 --wpos=pos
lmscell ctrl.nbfe -vfile=2 --stack~cmppos@shorten@wdx=posx@fn=pos~q
lmscell ctrl.nbfe -vfile=2 --stack~addpos@fn=posx.nbfe~wpos@fn=posfr~q
blm --dhftol=20 --molstat --nfile --ewald=1d-12 --rdsite --omax=-.01 --nosort --shorten=no --nit=20 --nk=9,4,2 --gmax=8 --mag --frzwf --wsite@fn=null nbfe
cat actrl.nbfe | sed s:'\(ATOM=Fe.*\):\1 MMOM=0,0,2:' | sed s:'XTOL=.001 GTOL=0:XTOL=0 GTOL=0.01:'  > ctrl.nbfe
rm -f save.nbfe mixm.nbfe hssn.nbfe
mpix -n 13 lmf --shorten=no ctrl.nbfe -vfile=2 -vbeta=.1 --rhopos -vfrzwf=t -vnit=20 --rpos=posfr --rs=11,1,1 > out.fullrelax &
cp rst.nbfe rst.nbfe.bk3

Relax the lattice constraining the Nb at the end

lmscell ctrl.nbfe -vfile=2 --stack~addpos@targ=4:21@src=4:21@fn=posx.nbfe~cmppos@fn=pos@shorten~wpos@fn=poscr0~q
lmscell ctrl.nbfe -vfile=2 --stack~setrlx@ijk=000@targ=1:3,22:24,2~wsite@fn=site3~q
cp rst.nbfe.bk2 rst.nbfe
rm -f mixm.nbfe log.nbfe mixm.nbfe hssn.nbfe save.nbfe
mpix -n 13 lmf ctrl.nbfe -vfile=3 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=20 --rs=11,1,1 --rpos=poscr0 -vminx=10 > out.doconstrainedrelax &

lmf ctrl.nbfe -vfile=2 -vbeta=.1 --rhopos -vfrzwf=t -vnit=0 --rs=1,0 --wpos=poscr
lmscell ctrl.nbfe -vfile=3 --rpos=poscr --stack~wsite@short@fn=site4~q
rm -f mixm.nbfe log.nbfe mixm.nbfe hssn.nbfe save.nbfe
mpix -n 13 lmf ctrl.nbfe -vfile=4 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 --rs=11,0 > out.constrainedrelax
cp rst.nbfe rst.nbfe.bk4
lmchk ctrl.nbfe -vfile=4 --syml~n=41~q=0,1,0,0,0,0,1,0,-1
mpix -n 13 lmf ctrl.nbfe -vfile=4 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 --band~col=1:150,367:516~col2=151:366~fn=syml

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 ~/lm).

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

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 and Fe on top of one another, the lattice vectors in the stacking plane must match. This entails:

  • a (3×1) reconstruction of bcc Nb(/tutorial/application/nbfe_superlattice/#reconstructed-nb-self-consistent-lda) with three atoms in the plane
  • a special reconstruction of bcc Fe with four atoms in a plane
  • a rotation of the z axis in Fe around the (1,0,1) axis to the (1,1,1)
  • a slight stretching of the Fe planes to exactly match lattice constants of the reconstructed Nb and Fe cells

The prescription is in outline:

  • Use blm to build input files for elemental Nb and Fe in reconstructed bcc structures.
  • Use the supercell maker lmscell to generate different reconstructions of the simple bcc Nb and Fe lattices that can be joined together.
    These bulk 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.
  • 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.

Elemental Nb and elemental Fe in reconstructed unit cells

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

  • 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
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 to be 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.5 0.5 0.5 1.5 1.5 -1.5 0.18425022 -0.50338096 0.68763117
#                        pos
 Fe        0.0000000   0.7500000   0.0000000
 Fe        0.5000000   1.0000000  -0.5000000
 Fe        1.0000000   1.2500000  -1.0000000
 Fe        0.0000000   0.0000000   0.0000000

The only difference should be in the site order (lmscell has changed since this tutorial was written)

Reconstructed Nb, self-consistent LDA

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.2262/5.4042 = 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

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.707107    0.000000    0.707107)
 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.50000   0.50000  -0.50000   0.00000   0.00000
    3       2:Nb2         1    0      1.00000   1.00000  -1.00000   0.00000   0.00000

Notes:

  • The planes defined by the first two vectors P1×P2 stack lie normal to the [1,0,1] direction (h in the table). All three atoms lie in a single (101) 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 in elemental Nb, but 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.5 --nk=12,12,12 --ctrl=ctrl --gw nb
lmfa -vfile=3 ctrl.nb

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

blm also automatically put a high-lying LO on the Nb d. It also isn’t needed here, so strip it from the ctrl file:

sed -i 's/ *PZ.*//' ctrl.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 -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. The two energies for iteration 10 should read:

   it 10  of 10    ehf=  -22895.047312   ehk=  -22895.047298
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/2,0,-1/2

Create syml.nb with

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

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:25   1:1(s)   2:4(p)   5:9(d)   10:16(f) 17:17(s) 18:20(p) 21:25(d)
   2   Nb   26:50   26:26(s) 27:29(p) 30:34(d) 35:41(f) 42:42(s) 43:45(p) 46:50(d)
   3   Nb   51:75   51:51(s) 52:54(p) 55:59(d) 60:66(f) 67:67(s) 68:70(p) 71:75(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=26,51~col2=2:4,18:20,seq=27,52~col3=5:9,21:25,seq=30,55~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.

Energy bands of Nb. Red, green, blue correspond to s, p, and d orbital character.

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))' -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 (1,0,1) direction: the first two lattice vectors are the same as in the Nb case. Now there four atoms, and they lie in a single plane
  • The first two lattice vectors are identical to those of site3.nb, so the 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. We will assume that Fe is much thinner than Nb and stretch it, keeping the Nb fixed.

Make a new file site4.fe taking this into account. Note: skip this step if you want to use file site4.fe that exactly reproduces this tutorial.

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

Note that the shear reduces the point group operations from 16 to 4:

lmplan -vfile=4 fe

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 nk1, nk2, nk3 that approximates the total number of k points and makes the spacing along the three axes as uniform as it can.

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

sed -i 's/ *PZ.*//' ctrl.fe
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:

  it 27  of 30    ehf=  -10164.162235   ehk=  -10164.162234
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.

Majority (left) and minority (right) energy bands of Fe. Red, green, blue correspond to s, p, and d orbital character.

energy bands of Fe in the plane normal to [1,0,1]

Bands near EF 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.

The Nb/Fe superlattice

lmscell –stack can assemble superlattices from individual site files with identical lattice constant and the first two lattice vectors. The supercell is constructed by adding along the third lattice vector. 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.

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~wsite@fn=sitein~'

This command does the following

  • 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 next 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. Now 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.
  • Writes the superlattice site file to sitein.nb.

Try running lmscell in interactive mode. You can watch how it builds the stack one block at a time:

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

Also try making just three layers of Nb.

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

This makes a superlattice of two planes of Nb, three atoms/plane, and writes the structure to site4.nb. In the last table lmscell prints to stdout, h is the projection along the (1,0,1) line. The atoms are ordered into three planes. This lattice is still merely a superlattice of 9 Nb atoms. Compare the two lattices:

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

The original cell has a volume of 120.669456 a.u., and all the spheres touch (0.00% overlap), with a standard neighbor list of (8,6,12) pairs in the first, second, third shells.

The supercell has volume 1086.025100 a.u. = 9×120.669456, and all the spheres continue to exactly touch, and the neighbor list of each Nb atom is that of bcc.

Make ctrl.nbfe from this structural data. Note that blm doesn’t know about magnetic moments so we have to modify the autogenerated actrl.nbfe to give it a reasonable starting point. The steps below make the density roughly self-consistent. This will provide idea of what is transpiring at the interface.

Rigid stacking of Nb and Fe planes

Setup:

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~wsite@fn=sitein~'

lmscell makes sitein.nbfe. We could have use the autogenerated basp0.nbfe for the basis, but for consistency we will stick with the Fe and Nb basis sets calculated for their respective elements. lmscell makes basp.nbfe from basp.nb and basp.fe.

cp sitein.nb sitein.nbfe
blm --nfile --ewald=1d-12 --rdsite --omax=-.01 --nosort --shorten=no --nit=20 --nk=9,4,2 --gmax=8 --mag --frzwf nbfe
cat actrl.nbfe | sed s:'\(ATOM=Fe.*\):\1 MMOM=0,0,2:' > ctrl.nbfe
rm -f rst.nbfe mixm.nbfe log.nbfe
lmfa ctrl.nbfe

blm reads structural data from sitein.nbfe, and generates actrl.nbfe. Command-line switches modify the ctrl file in the following ways:

  • --ewald=1d-12: long cells put greater strains keeping the overlap matrix positive definite. Setting a tight tolerance on Ewald summations for envelope functions keeps the overlap matrix positive definite.
  • --rdsite : tells blm to read structural data from sitein.nbfe
  • --nk=9,4,2 : is a mesh with approximately equal spacing between k points
  • --omax=-.01 : is not necessary, but because of the messy interface, keeping the augmentation spheres a little bit apart is a safety measure that simplifies and stabilizes convergence.
  • --frzwf : is not necessary, but it keeps the basis set frozen (at a slight loss in accuracy), which simplifies and stabilizes convergence.

For the remaining switches see blm’s command line documentation.

Make approximately self-consistent:

mpix -n 13 lmf ctrl.nbfe -vbeta=.1  --rhopos -vfrzwf=t -vnit=10 > out.00  &

After self-consistency is reached inspect the forces at the end of out.00:

  ib           estatic                  eigval                    total
   1   -8.07  -10.32   -2.57     4.93    9.88   -0.22    -3.13   -0.44   -2.80
   2   -3.96   10.14   -8.77     2.33  -10.36    6.95    -1.63   -0.22   -1.82
   3   -3.51   -0.37   -4.55     6.89    0.33    8.65     3.38   -0.03    4.11

   4  -49.43   13.19  -32.23    57.48  -29.30   39.60     8.05  -16.11    7.37
   5  -39.59    0.63  -38.68    43.03   -0.25   41.43     3.43    0.38    2.75
   6  -32.91  -12.32  -49.96    39.84   29.35   58.20     6.93   17.03    8.24

   7   -5.88  -36.08    0.86    -0.80   52.63   -7.02    -6.68   16.54   -6.15
   8    9.17   -2.57   11.61     8.04    1.87    6.91    17.21   -0.70   18.51
   9    2.67   39.79   -4.62    -5.10  -56.80    2.38    -2.42  -17.02   -2.24
  10    3.81   -1.21    0.55   -20.73    0.76  -17.78   -16.92   -0.44  -17.23

  11    1.52   -0.53   -5.08    -2.22   -1.61    5.53    -0.69   -2.14    0.45
  12   -2.32   -1.53    1.22     0.66   -0.37   -0.79    -1.66   -1.90    0.44
  13   -1.12    2.62    5.65     1.13   -1.20   -7.30     0.01    1.43   -1.66
  14    6.39   -0.24   -2.96    -6.30    2.25    1.03     0.09    2.01   -1.93

  15   60.21   -2.02   59.86   -91.44    2.01  -91.02   -31.23   -0.01  -31.16
  16   38.35    7.02   60.68   -46.23   -8.16  -96.06    -7.88   -1.14  -35.38
  17   55.02   -3.03   56.32   -69.10    3.94  -71.52   -14.08    0.91  -15.20
  18   55.21   -3.17   37.61   -89.96    5.58  -44.25   -34.75    2.41   -6.65

  19  -25.96    0.22  -26.36    41.22    0.30   42.10    15.26    0.52   15.75
  20  -45.68   18.16  -34.06    64.61  -40.46   59.30    18.93  -22.30   25.24
  21  -37.07  -18.19  -46.60    61.85   39.35   63.89    24.78   21.16   17.30

  22    7.15    8.58    4.51     0.71   -6.90    4.68     7.86    1.67    9.19
  23    6.36   -8.97    7.83     2.27    6.80   -2.01     8.63   -2.17    5.82
  24    9.64    0.20    9.73    -3.12    0.35   -2.65     6.51    0.55    7.08

Forces in the middle layers (sites 1-3 for Nb and 22-24 for Nb, 11-14 for Fe) are rather small, as might be hoped and expected. However, forces on the frontier atoms are much larger. Consider the projection of the forces along the [101] line, averaged by layer:

  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
  force    small     small    medium(+)   small   large(-)  large(+)   small

The largest planar average forces occur between Fe 3 and Nb 3. They push in opposite senses, wanting to separate. Thus, if the Fe block is rigidly shifted to a smaller projection h along [101] the energy should go down. (There is a counterbalance: forces on the Fe 1 layer will increase). The reason is easily seen: do

lmscell ctrl.nbfe --stack~show@planes~q

The spacing between Nb 2 and Fe 1 is Δh=0.70711, while the spacing between Fe 3 and Nb 3 is significantly smaller, Δh=0.61651.

We can improve the unrelaxed stacking by rigidly shifting the entire Fe block a variable amount and finding the point of minimum energy. This step is not necessary but it makes the full relaxation easier, and moreover it is a natural step in anticipation of the fact that we will later constrain positions of the anchor Nb layers. This will be necessary when it is used as part of an infinitely repeating lead, as it must have the structure of bulk Nb.

Repeat construction of the superlattice, displacing the Fe block by (-0.01,0,-0.01) and again by (-0.02,0,-0.02). Note that files site.ext, site1.ext, site2.ext hold structural data for the three displacements.

lmscell -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dpos=-.01,0,-.01@dup=3~file@site3.nb@dup=2~show~basp@fn=basp.nbfe@basp.nb@basp.fe~wsite@fn=sitein~'
cp sitein.nb site1.nbfe
rm mixm.nbfe log.nbfe rst.nbfe
mpix -n 13 lmf ctrl.nbfe -vbeta=.1 --rhopos -vfrzwf=t -vnit=10 -vfile=1 > out.m01 &

lmscell -vfile=3 ctrl.nb '--stack~file@site3.nb~file@site4.fe@dpos=-.02,0,-.02@dup=3~file@site3.nb@dup=2~show~basp@fn=basp.nbfe@basp.nb@basp.fe~wsite@fn=sitein~'
cp sitein.nb site2.nbfe
rm mixm.nbfe log.nbfe rst.nbfe
mpix -n 13 lmf ctrl.nbfe -vbeta=.1 --rhopos -vfrzwf=t -vnit=10 -vfile=2 > out.m02 &

This command command extracts the energy of the last iteration for each of the three displacements and should yield the following:

grep ^x save.nbfe | vextract . ehf
  -122072.2671338
  -122072.2768167
  -122072.2819031

To see the maximum force on any atom at the last iteration for each of the three cases do the following. It should yield

find -maxdepth 1 -type f  -name 'out.*0*'  -exec  sh -c 'echo "$(grep 'Maximum' $1 | tail -n 1),$1"' _ {} \;
 Maximum Harris force = 44.1 mRy/au (site 15)  Max eval correction = 18.2,./out.00
 Maximum Harris force = 38.2 mRy/au (site 8)  Max eval correction = 14.7,./out.m02
 Maximum Harris force = 35 mRy/au (site 15)  Max eval correction = 16.9,./out.m01

The global maximum force drops for the initial displacement (dpos=-.01,0,-.01), but begins to rise again at dpos=-.02,0,-.02. Compare the forces at the last iteration in out.m02 to out.00. It reveals that indeed forces on layers Fe 3 and Nb 3 are significantly smaller, while forces on the Fe 1 layer have increased.

Note that with this restricted relaxation some Nb/Fe atoms have moved closer together. There is now a 3% overlap between atoms 5 and 8, as can be seen from the output of

lmchk ctrl.nbfe
lmchk -vfile=2 ctrl.nbfe

A parabolic fit of the total shows that the energy would be minimized near  disp=(-0.026,0,-0.026). It’s close enough to the structure with  disp=(-0.02,0,-0.02), so we fully relax the lattice starting from the last structure.

Relaxing the superlattice

To prosecute the full relaxation, use blm to add some features to ctrl.nbfe

blm --dhftol=45 --conv=1e-4 --convc=1e-3 --molstat --nfile --ewald=1d-12 --rdsite --omax=-.01 --nosort --shorten=no --nit=20 --nk=9,4,2 --gmax=8 --mag --frzwf --wsite@fn=null nbfe
cat actrl.nbfe | sed s:'\(ATOM=Fe.*\):\1 MMOM=0,0,2:' | sed s:'XTOL=.001 GTOL=0:XTOL=0 GTOL=0.01:'  > ctrl.nbfe

Note: saving rst.nbfe isn’t necessary, but by hanging onto it we can restart at this point to carry out restricted relaxations, which we will carry out after the full one

The blm switches modify the ctrl file in the following ways:

  • --conv=1e-4  : Set a tolerance rougher than the default for convergence in the energy before shifting atoms
  • --convc=1e-3 : Set a tolerance rougher than the default for convergence in the charge density before shifting atoms
  • --dhftol=45  : Set an additional requirement for convergence : an upper limit on the correction term to the Harris-Foulkes forces lmf estimates a correction to the raw forces from the Harris functional; this correction must fall below 45 for the convergence criterion to be satisfied.
  • --molstat : Add lines that turn on molecular statics (lattice relaxation) if preprocessor variableminx  is set.
  • --ewald=1d-12: long cells put greater strains keeping the overlap matrix positive definite. Setting a tight tolerance on Ewald summations for envelope functions keeps the overlap matrix positive definite.
  • --frzwf: is not necessary, but it keeps the basis set frozen (at a slight loss in accuracy), which makes relaxations more stable.
  • --wsite@fn=null : is there only to prevent blm from overwriting site.nbfe.

For the remaining switches see the previous invocation of blm or the command line documentation.

Setting soft limits for --conv and --convc while imposing an extra condition --dhftol=45 is suitable for molecular statics. We don’t care so much whether the energy is well converged, but that the forces are well described. The forces are good enough for a rough relaxation, which since it will be repeated with some constraints is all we need at this stage.

cp rst.nbfe rst.nbfe.bk
rm -f save.nbfe mixm.nbfe hssn.nbfe
mpix -n 13 lmf ctrl.nbfe -vfile=2 -vbeta=.1 --rhopos -vfrzwf=t -vnit=20 -vminx=5 --wpos=pos > out.dorelax

This uses a Broyden algorithm to relax the structure for five iterations. This is sufficient to that the maximum force is 6 mRy/au, as can be seen by inspecting the final forces table in out.dorelax.

The mean force drops with each iteration, and the total energy goes down. Using the command below you should see the following:

egrep ^c\|RELAX out.dorelax
grep ^c out.dorelax  | vextract . mmom ehf | awk '{printf " 'mom' %10.6f  'ehf' %10.6f\n", $1, $2}'

From the first command |g| is the length of the 72 component force “vector”; it should drop from an initial value 0.0893 to 0.0142 in the 5th iteration, and ehf should drop each iteration.

The next steps do the following:

  1. dump the (relaxed) positions contained in rst.nbfe into file pos.nbfe
  2. calculate the shift relative to starting point (removing lattice translation vectors) and write it to file posx.nbfe
  3. construct a “fully relaxed” positions file posfr.nbfe
lmf ctrl.nbfe -vfile=2 -vbeta=.1 --rhopos -vfrzwf=t -vnit=0 --rs=1,0 --wpos=pos
lmscell ctrl.nbfe -vfile=2 --stack~cmppos@shorten@wdx=posx@fn=pos~q
lmscell ctrl.nbfe -vfile=2 --stack~addpos@fn=posx.nbfe~wpos@fn=posfr~q

Make the system self-consistent at the final position the relaxation algorithm found, reverting to tighter tolerances in the convergence; this is what blm does. Note also that shortening of positions by lattice translations is suppressed. This is not necessary but it simplifies synchronization.

cp rst.nbfe rst.nbfe.bk2
blm --dhftol=20 --molstat --nfile --ewald=1d-12 --rdsite --omax=-.01 --nosort --shorten=no --nit=20 --nk=9,4,2 --gmax=8 --mag --frzwf --wsite@fn=null nbfe
cat actrl.nbfe | sed s:'\(ATOM=Fe.*\):\1 MMOM=0,0,2:' | sed s:'XTOL=.001 GTOL=0:XTOL=0 GTOL=0.01:'  > ctrl.nbfe
rm -f save.nbfe mixm.nbfe hssn.nbfe
mpix -n 13 lmf --shorten=no ctrl.nbfe -vfile=2 -vbeta=.1 --rhopos -vfrzwf=t -vnit=20 --rpos=posfr --rs=11,1,1 > out.fullrelax &
cp rst.nbfe rst.nbfe.bk3

Constrained relaxation of the superlattice

We repeat the relaxation, this time constraining the “anchor” planes Nb 1 and Nb4 to their ideal bcc positions. We do this because our objective is to study transport through a non-periodic trilayer ( Nb / Fe / Nb ), where Nb extends to infinity on the left and right and should have the structure of perfect bcc.

To accomplish this we must (1) restore positions of the end layers to their ideal positions and (2) redo the relaxation under the constraint that atoms in these layers don’t move. For the former, do

lmscell ctrl.nbfe -vfile=2 --stack~addpos@targ=4:21@src=4:21@fn=posx.nbfe~cmppos@fn=pos@shorten~wpos@fn=poscr0~q

This makes a positions file poscr0.nbfe. Compare it to the fully relaxed positions file posfr.nbfe to confirm that sites 1-3 and 21-24 have been restored to their ideal positions:

mcx poscr0.nbfe posfr.nbfe  --

For the latter, you can tell lmf not to move specific atoms when performing a relaxation. Note the  rlx column in site2.nbfe. It is a sequence of three integers (each 0 or 1) that perform the same function as SITE_ATOM_RELAX, namely flag which of the (x,y,z) components of position should be allowed to move when relaxations are carried out. It can be read via the site file rather than the ctrl file. Use your text editor, or let lmfscell make this change for you:

lmscell ctrl.nbfe -vfile=2 --stack~setrlx@ijk=000@targ=1:3,22:24,2~wsite@fn=site3~q

Repeat the relaxation, this time with constraints

cp rst.nbfe.bk2 rst.nbfe
rm -f mixm.nbfe log.nbfe mixm.nbfe hssn.nbfe save.nbfe
mpix -n 13 lmf ctrl.nbfe -vfile=3 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=20 --rs=11,1,1 --rpos=poscr0 -vminx=10 > out.doconstrainedrelax &

Inspect the final forces table and confirm there remains a residual strain at the end layers. This is to be expected, and the strain is small enough that we can neglect it.

Extract the relaxed positions from rst.nbfe, and create a file site3.nbfe. Then confirm that the system is self-consistent

lmf ctrl.nbfe -vfile=2 -vbeta=.1 --rhopos -vfrzwf=t -vnit=0 --rs=1,0 --wpos=poscr
lmscell ctrl.nbfe -vfile=3 --rpos=poscr --stack~wsite@fn=site4~q
rm -f mixm.nbfe log.nbfe mixm.nbfe hssn.nbfe save.nbfe
mpix -n 13 lmf ctrl.nbfe -vfile=4 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 --rs=11,0 > out.constrainedrelax
LDA Energy Bands and layer resolved DOS of the superlattice

To distinguish Fe and Nb character, read the orbitals table:

We need a pair of symmetry lines in the plane perpendicular to (1,0,1). Usually bands and on faces of the Wigner Seitz cell, that is at half-integer multiples of reciprocal lattice vectors. To confirm that points (0,1,0) and (1,0,−1) (perpendicular to (1,0,1)) are also the smallest half-integer multiples of reciprocal lattice vectors, use the superlattice editor and do

lmscell ctrl.nbfe -vfile=4 --stack~cart2lat@q@x=0,1,0~cart2lat@q@x=1,0,-1~q

This shows that q=(0,1,0) and q=(1,0,−1) are indeed near half-integer multiples of reciprocal lattice vectors, namely (0.5,1.5,-1.5) and (-1.0,3.0,-5.5). Make the symmetry lines file with:

lmchk ctrl.nbfe -vfile=4 --syml~n=41~q=0,1,0,0,0,0,1,0,-1

It is useful to know whether a band is mostly Fe-like or Nb-like. Get a table of orbitals in the superlattice:

lmf ctrl.nbfe -vfile=4 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 --pr55 --quit=ham | grep -A 50 'Orbital positions'

Orbitals 1:150 and 367:516 belong to Nb while 151:366 belong to Fe. Make bands with color weights

mpix -n 13 lmf ctrl.nbfe -vfile=4 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 --band~col=1:150,367:516~col2=151:366~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

Majority (left) and minority (right) energy bands of Fe.   Red: Nb character; green: Fe character

LDA bands of the Nb/Fe superlattice

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 EF: both groups both fall EF for the majority spin, and for the minority spin EF 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 EF. 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=4 --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=4 --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' -lst2 dos.nbfe
fplot -f plot.dos

Majority and minority DOS of NbFe resolved by layer

Footnotes and references