Questaal Home

Nb/Fe/Nb Metallic multilayer, Electronic Structure and Landauer-Buttiker Transport

The electronic structure of a Nb/Fe/Nb multilayer is set up, using Questaal’s superlattice editor, and transmission computed through a Nb/Fe/Nb trilayer.

Nb and Fe are both bcc but are strongly lattice mismatched, which considerably complicates matters.

The tutorial begins by building a reconstructed Nb/Fe superlattice, and relaxing it in the LDA using the lmf code.

The tutorial switches to the ASA, where transport through a Nb/Fe/Nb trilayer can be computed using the layer code lmpg.

A separate tutorial also shows how to build superlattices. Its purpose is to explain how to calculate the potential of superlattices in a QuasiParticle Self-Consistent GW framework

Table of Contents

Bulk Nb

Bulk Fe

An Nb/Fe superlattice

to be completed


This tutorial makes extensive use of many of the Questaal packages. including lmscell, lmplan lmf, lm, lmgf, lmpg, 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 to be ~/lm here).

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


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 much match. This entails:

  • a (3×1) reconstruction of bcc Nb 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
    ALAT=6.226 PLAT= -.5 .5 .5  .5 -.5 .5  .5 .5 -.5
    ATOM=Nb   X=0 0 0
# Init file for Fe.  Lattice constant at 0K
    ALAT=5.404 PLAT= -.5 .5 .5  .5 -.5 .5  .5 .5 -.5
    ATOM=Fe  MMOM=0,0,2.2
    ATOM=Fe   X=0 0 0


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

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}

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


  • The planes defined by the first two vectors P1×P2 stack lie normal to the [1,0,1] direction. All three atoms lie in a single (101) plane.
  • 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. 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

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.

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). Copy the following to syml.nb

31 0 1 0      0 0 0
31 0 0 0     .5 0 -.5

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:30   1:1(s)   2:4(p)   5:9(d)   10:16(f) 17:17(s) 18:20(p) 21:25(d) 26:30(d)
   2   Nb   31:60   31:31(s) 32:34(p) 35:39(d) 40:46(f) 47:47(s) 48:50(p) 51:55(d) 56:60(d)
   3   Nb   61:90   61:61(s) 62:64(p) 65:69(d) 70:76(f) 77:77(s) 78:80(p) 81:85(d) 86:90(d)

The s, p, d orbitals are conveniently grouped by using extensions to the Questaal standard syntax for Integer Lists. The following band calculation

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

makes bands and groups the 6 s, 18 p, and 30 d orbitals into first, second, and third weights, respectively. The second line makes a postscript file 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.

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

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

Reconstructed Nb, self-consistent ASA-LDA

It is convenient here to repeat the calculation in the ASA approximation, as we will need the ASA for transport and need to confirm its reliability for bulk Nb.

cp init.nb init.nbasa
blm --nfile --nk=12,12,12 --asa nbasa
cat actrl.nbasa | sed 's/IDXDN.*/GRP2=1/'  | awk '{ if ($1 == "OPTIONS") {print "HAM   PMIN=-1 QASA=1"}; print}' > ctrl.nbasa

blm makes an acceptable ctrl file, but we make a few tweaks in anticipation of the fact that ultimately we will use lmpg for transport. This code doesn’t allow down-folding, and it uses KKR style accumulation of weights. Compare the changes and refer to the input file documentation for the meaning of each new tag.

Check the sphere overlaps

cp site3.nb site3.nbasa
lmchk -vfile=3 ctrl.nbasa

The MT sphere has been enlarged so that the sphere volume matches the unit cell volume, as the ASA requires.

Make the ASA density self-consistent in the reconstructed structure

lmstr -vfile=3 ctrl.nbasa
lm -vnit=0 -vfile=3 ctrl.nbasa
lm -vnit=30 -vfile=3 ctrl.nbasa
cp syml.nb syml.nbasa
lm -vnit=30 -vfile=3 ctrl.nbasa --band~fn=syml

Draw the ASA bands:

echo -6,6,5,10 | plbnds -fplot~sh~scl=.5 -ef=0 -scl=13.6 --nocol -dat=red  red

The ASA energy bands are quite similar to the full LDA ones, and sufficiently close for the present purposes. Note the bands cross the Fermi surface at similar points with similar slope.

Blue: LDA bands generated by lmf. Red: ASA-LDA generated by lm.

FP and ASA energy bands of Nb compared

Reconstructed Fe, self-consistent LDA

Build a ctrl and site file for bcc Fe, in the same manner as for Nb

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


  • 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.
lmscell -vfile=3 ctrl.fe -vz=6.24000171/6.226 '--stack~scale=1/z~stretch=z^3~sort@h3~wsite@short@fn=site4'

Construct file basp.fe and remake ctrl.fe:

lmfa -vfile=4 --basfile=basp ctrl.fe
blm --nit=20 --nfile --nk=10,5,12 --ctrl=ctrl --gw --gmax=7.9 --mag fe

lmfa supplies the value 7.9 for gmax. The k mesh was selected to make spacing between points on the three axes approximately equal.

Make the density self-consistent

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

The self-consistent magnetic moment should come out to about 8.67, or about 2.17 a.u./Fe atom. This is close to the observed magnetic moment.

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=31,61,91~col2=2:4,18:20,seq=32,62,92~col3=5:9,21:25,seq=35,65,95~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
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 and 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 classical two-peak structure in the density-of-states.

Reconstructed Fe, self-consistent ASA-LDA

It is convenient here to repeat the calculation in the ASA approximation, as we will need the ASA for transport and need to confirm its reliability for bulk Nb.

cp init.fe init.feasa
blm --nfile --nk=12,12,12 --mag --asa feasa
cat actrl.feasa | sed 's/\(MMOM.*\)/\1 GRP2=1/'  | awk '{ if ($1 == "OPTIONS") {print "HAM   PMIN=-1 QASA=1"}; print}' > ctrl.feasa

Check the sphere overlaps

cp site4.fe site4.feasa
lmchk -vfile=4 ctrl.feasa

The MT sphere has been enlarged so that the sphere volume matches the unit cell volume, as the ASA requires.

Make the ASA density self-consistent in the reconstructed structure

lmstr -vfile=4 ctrl.feasa
lm -vnit=0 -vfile=4 ctrl.feasa
lm -vnit=30 -vfile=4 ctrl.feasa
cp syml.fe syml.feasa
lm -vnit=30 -vfile=4 ctrl.feasa --band~fn=syml

Draw the ASA bands:

echo -6,6,5,10 | plbnds -fplot~sh~scl=.5 -ef=0 -scl=13.6 --nocol -dat=red  red

As with Nb, the ASA energy bands are quite similar to the full LDA ones, and indeed significantly closer to each other than to bands computed from a higher level theory such as GW.

Majority (left) and minority (right) energy bands of Fe. blue: LDA generated by lmf; red: ASA-LDA generated by lm.

FP and ASA energy bands of Fe compared

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
file site4.fe dup=3 dpos=-.01,0,-.01
file site3.nb dup=2
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 the density roughly self-consistent. This will provide idea of what is transpiring at the interface.

Rigid stacking of Nb and Fe planes

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
mpix -n 13 lmf ctrl.nbfe -vbeta=.1  --rhopos -vfrzwf=t -vnit=10 > out.00  &

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.

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 &

grep ^x save.nbfe | vextract . ehf

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


Inspect the forces in out.m02. 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 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
cp null.nbfe site.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 component of 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:

egrep ^c\|RELAX out.dorelax

Here |g| is the length of the 72 component force “vector”.

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. (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 latter, you can tell lmf not to move specific atoms when performing a relaxation. Edit site2.nbfe and note the last column:

#                  pos                                  vel                                   eula                   vshft  PL rlx
 Nb  0.0000000   0.0000000   0.0000000   0.0000000    0.0000000    0.0000000   0.0000000    0.0000000    0.0000000 0.000000  0 000

The three digits in the rlx column 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.

For the first three sites and also the last three, use your text editor to set rlx to 000. Then 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  --

After you have edited the rlx column in site2.nbfe, repeat the relaxation

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=2 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=20 --rs=11,1,1 --rpos=poscr0 -vminx=5 > 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=2 --rpos=poscr --stack~wsite@short@fn=site3~q
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=1 --rs=11,0,1 > out.constrainedrelax
LDA Energy Bands of the superlattice

To distinguish Fe and Nb character, get 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=3 --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). Cut and past the following box into file syml.nbfe

51 0 1 0    0 0 0
35 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=3 --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=3 --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
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

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.

Self-consistent ASA Calculation of the Nb/Fe interface

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

Create a directory asa below your working directory to keep the files distinct.

A fist cut at a ctrl file can be obtained with

blm --pr35 --gf --nfile --rdsite=site4 --nosort --shorten=no --nk=9,4,2 --mag --asa~pfree~rmaxs=7 --wsite@fn=null nbfe -vfile=4 --findes~rmin=.9~shorten=0,0,1~omax=0~sclwsr=20~spec~nspmx=4~nesmx=36 --omax=.18,.26,.26

However, there are some subtleties in how to deal with location of empty spheres we defer for nnow about how ctrl.asa and site.asa. Copy the boxes below into ctrl.asa and site.asa.)


# Autogenerated from init.nbfe using:
# blm --pr35 --gf --nfile --rdsite=site4 --nosort --shorten=no --nk=9,4,2 --mag --asa~pfree~rmaxs=7 --wsite@fn=null nbfe -vfile=4 --findes~rmin=.9~shorten=0,0,1~omax=0~sclwsr=20~spec~nspmx=4~nesmx=36 --omax=.18,.26,.26

# Variables entering into expressions parsed by input
% const nit=10 asa=t les=1
% const met=(asa?1:5)
% const nsp=2 so=0
% const lxcf=2 lxcf1=0 lxcf2=0     # for PBE use: lxcf=0 lxcf1=101 lxcf2=130
% const ccor=T sx=0 gamma=sx scr=4 # ASA-specific
% const rmaxs=7                    # Range for structure constants
% const gfmode=1 nz=16 ef=0 c3=t   # lmgf-specific variables
% const nk1=9 nk2=4 nk3=2 beta=.3

VERS  LM:7 ASA:7 # FP:7
# Lattice vectors and site positions
%ifdef file>=1
  file=   site{file}
  file=   site
# file= essite

# Self-consistency
  nit=    {abs(nit)}               # Maximum number of iterations
#  mix=    A2,b={beta},n=5,k=5,0
#   mix=    A7,b={beta},n=3,k=3;A2,b={beta/10},n=2,k=2,elind=0.001
#   mix=    A7,b={beta},n=7,k=7;A2,b={beta/10},n=2,k=2,elind=0
  mix=    B2,b={beta},k=7          # Charge density mixing parameters
  conv=   1e-5                     # Convergence tolerance (energy)
  convc=  3e-5/3                   # tolerance in RMS (output-input) density

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

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

SYMGRP  find
HAM   PMIN=-1                      # Constrain Pnu > P_free

      ASA[ CCOR={ccor} TWOC=F ADNF=F GAMMA={gamma}] SCR={scr} SX={sx}
STR   RMAXS={rmaxs}                # Range of strux

GF    MODE={gfmode} GFOPTS={?~c3~p3;~p2;}
BZ    EMESH={nz},10,-1,{ef},.5,.3  # nz-pts;contour mode;emin;emax;ecc;bunching

# SCLWSR=11 WSRMAX=3.3 OMAX1=0.18 0.26 0.26
  ATOM=Nb       Z= 41  R= 3.065511  LMX=2  IDXDN=0,0,0,2 # IDMOD=0,1
  ATOM=Nbx      Z= 41  R= 2.730129  LMX=2  IDXDN=0,0,0,2 # IDMOD=0,1
  ATOM=Fex      Z= 26  R= 2.496030  LMX=2
  ATOM=Fe       Z= 26  R= 2.719829  LMX=2
  ATOM=E        Z=  0  R= 1.715533  LMX=0
  ATOM=E1       Z=  0  R= 1.710466  LMX=0
  ATOM=E2       Z=  0  R= 1.651325  LMX=0
  ATOM=E3       Z=  0  R= 1.202364  LMX=0

# SCLWSR=11 WSRMAX=3.3 OMAX1=0.18 0.26 0.26
  ATOM=Nb       Z= 41  R= 3.065511  LMX=2  IDXDN=0,0,0,2  IDMOD=0,1
  ATOM=Nbx      Z= 41  R= 2.683855  LMX=2  IDXDN=0,0,0,2  IDMOD=0,1
  ATOM=Fex      Z= 26  R= 2.453724  LMX=2  MMOM=0,0,2
  ATOM=Fe       Z= 26  R= 2.673730  LMX=2  MMOM=0,0,2
  ATOM=E        Z=  0  R= 1.715533  LMX=0
  ATOM=E1       Z=  0  R= 1.796268  LMX=0
  ATOM=E2       Z=  0  R= 1.737788  LMX=0
  ATOM=E3       Z=  0  R= 1.202364  LMX=0

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


% site-data vn=3.0 fast io=15 nbas=60 alat=6.226 plat= -0.5 0.5 0.5 1.5 1.5 -1.5 0.55275066 -1.51014288 6.06289351
#                        pos
 Nb        0.0000000   0.0000000   0.0000000
 Nb        0.5000000   0.5000000  -0.5000000
 Nb        1.0000000   1.0000000  -1.0000000
 E3       -0.4737501  -0.7394575   1.0371136
 E3        0.5357464   0.7459879   0.0428640
 E3        0.0444594   0.2488781   0.5344094
 E3       -0.4454942  -0.2743963   1.0326326
 E3        0.5547866   0.2622121   0.0406214
 E3        0.0405665  -0.5050013   0.7628725
 E3        0.5391316  -0.0233702   0.2705096
 E3       -0.4670669  -0.9928333   1.2831687
 E3        0.2968868  -0.0115608   0.5300687
 E3        0.8001271   0.5039323   0.0363901
 Nbx       0.0038386  -0.0166966   0.9910505
 Nbx       0.9928133   1.0176382   0.0048461
 Nbx       0.5161225   0.4990571   0.5150855
 E3       -0.0164014  -0.4956354   1.2589020
 E3       -0.3739799  -0.9880550   1.6531767
 E3        0.7909603   0.0066238   0.5321685
 E3        0.4664294  -0.3020557   0.9890124
 E3       -0.0259534  -0.7485791   1.4813952
 E3        0.5611756   0.1357335   0.9604225
 E3        0.0229530  -0.2263954   1.4986452
 E2       -0.5404845  -0.6328016   2.0620827
 E1       -0.4216957  -1.1627362   1.9898101
 Fex      -0.0399873  -0.0010873   1.9603949
 Fex      -0.0259839   0.7809822   1.9780095
 Fex       0.9792297   1.2189900   0.9795957
 Fex       0.5160204   0.9996464   1.5211879
 Fe        0.1509317   0.2442758   2.6493636
 Fe        0.1626577  -0.5019115   2.6567494
 Fe        0.6537761   0.4947430   2.1670347
 Fe        1.1725712   0.7480391   1.6834229
 Fex       0.8421375  -0.0133975   2.8192889
 Fex       0.3308379  -0.2533275   3.3370495
 Fex       1.3371460   0.2434396   2.3440138
 Fex       0.3278900  -0.9987992   3.3586058
 E3       -0.3619496  -1.1309320   4.2693144
 E         0.8837140  -0.5118108   3.3420287
 E3       -0.0787140  -1.4964656   4.3085914
 E3        0.6610712  -0.7050789   3.5688063
 E3        0.2767879  -1.0210851   4.0233808
 E3        0.5795804  -1.0270209   3.8087107
 Nbx       1.5400241  -0.4828374   3.0627345
 Nbx       0.5490728  -1.5111615   4.0608585
 Nbx       1.0623082  -1.0379833   3.5524052
 E3        1.0240526  -0.5127517   3.7389527
 E3        0.5165514  -1.0225847   4.2903860
 E3        0.0116892  -1.4962605   4.7978325
 E3        0.3010275  -1.4799238   4.5283927
 E3        1.2996213  -0.5585176   3.5297990
 E3        1.0168605  -0.2808990   4.0027595
 E3       -0.0005736  -1.2373961   5.0284631
 E3       -0.0037828  -1.7497071   5.0358071
 E3        0.5024477  -0.7478516   4.5295767
 E3        0.5170890  -1.2513967   4.5213959
 E3        1.0304754  -0.7795029   4.0211892
 Nb        0.5527507  -1.5101429   5.0628935
 Nb        1.0527507  -1.0101429   4.5628935
 Nb        1.5527507  -0.5101429   4.0628935

Remove any atom files that might be present and Make the system self-consistent:

touch nb.asa
rm {nb,nbx,fe,fex,e}*.asa
lmstr asa
lm asa -vnit=0 --pr30,20
rm -f save.asa mixm.asa sv.asa log.asa
mpix -n 16 lm -vmet=3 ctrl.asa -vccor=f -vnit=400 --iactiv -vscr=4 --zerq~all~qout --pr30,21 > out  &

Convergence to self-consistency is slow, as a consequence of the screening algorithm. This is necessary to protect against low-energy, long-wave charge oscillations along the long axis. That this is an important issue can readily be seen from Poisson’s equation of a periodic cell in Fourier space:

Away from self-consistency there will be a deviation . This generates a change in the potential . Long cells have are vectors of low G, and changes to the potential from get greatly amplified. The solution is to screen out these low Fourier components using a model dielectric function. But the ASA has no plane-wave representation of the density, so constructing a model dielectric function is not so simple.

A model construction is nevertheless implemented, and it protects against these fluctuations, but at the same time it conceals true physical changes in them. This limits the rate of convergence. It can be dealt with by some extra tricks, but for the present we simply iterate a very large number of times to circumvent the problem.

The self-consistent magnetic moment is close to 24 μB, or about 2 μB per Fe atom (inspect the last line of save.asa). This is close to the full LDA result, 23.6 μB, and is a bit less than the moment of 12 Fe atoms in bulk, 12×2.2 = 26.4 μB. It shows that the Fe moment isn’t quite rigid, and gets extinguished somewhat at the Fe/Nb interface.

As a last step run lm one final iteration, to save in file rsta.asa essential information that can reconstruct the self-consistent density:

mpix -n 16 lm -vmet=3 ctrl.asa -vccor=f -vnit=1 --iactiv -vscr=4 --pr40,21 --rs=0,1 --quit=rho > out.asa

That rsta.asa contains sufficient information to reconstruct the density, remove the atom files containing potentials and do

rm {nb,fe,e}*.asa
lm asa -vnit=0 --pr30,20 --rs=1,0
mpix -n 16 lm -vmet=3 ctrl.asa -vccor=f -vnit=1 --iactiv -vscr=4 --pr40,21 --quit=rho > out

Files out and out.asa are essentially identical.

ASA Energy bands

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

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

Inspect the the table. We can associate orbitals 1:64 and 188:244 with Nb, and orbitals 74:172 with Fe.

Generate the energy bands

cp ../syml.nbfe syml.asa
lm -vmet=3 ctrl.asa -vccor=f -vnit=1 --iactiv -vscr=4 --pr40,21 --quit=rho --band~col=1:64,188:244~col2=74:172~fn=syml

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

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

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

ASA energy bands of the Fe/Nb interface

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

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

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

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

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

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

 BZINTS: Fermi energy:      0.154488; 156.000000 electrons;  D(Ef):  323.085

The Fermi level appropriate to lmgf should be close to 0.154 Ry, but not identically so because lmgf uses a different integration scheme and the third order Green’s function is close to, but identical with the 3rd order ASA hamiltonian.

Following the Introductory tutorial run lmgf to get the Fermi level satisfying charge neutrality

lmgf ctrl.asa -vccor=f -vnit=1 -vscr=4 --pr40,21 --quit=rho -vgamma=t -vef=.154

Inspect file vshft.asa. It found a constant potential shift, vconst=-0.0004575 must be subtracted from the estimate for the Fermi level. The Fermi level then should be ef-vconst=0.1544 Ry, perhaps fortuitously close to the Fermi level generated by lm. The close agreement should give use some confidence that the Green’s function is well described.

Spectral Function from lmgf

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

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

mpix -n 16 lmgf -vmet=3 ctrl.asa -vccor=f -vnit=1 -vscr=4 --pr40,21 --quit=rho -vgamma=t  -vef=.1544 --band~ef=.1544~fn=syml2

Draw a figure with

~/lm/startup/  -w -5:5
open asa_UP.png
open asa_DN.png

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

NbFe SF up figure.

NbFe SF dn figure.

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

Nb/Fe/Nb Trilayer: lmpg

… to be completed

Footnotes and references

If this page has any errors, there is something you think is missing or unclear, or for any other issues, you can create a post here letting us know.