Questaal Home
Navigation

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

Preliminaries

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.

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 are 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 noneriodic 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 nonperiodic 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 elmeental 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.

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

Notes:

  • 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 anway. 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 unecessary 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
open fplot.ps

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

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 downfolding, 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
open fplot.ps

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

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

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

  -122072.2671369
  -122072.2768195
  -122072.2819049

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

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

FP and ASA energy bands of Fe compared

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.

ASA band structure of the superlattice

rm -f mixm.asa log.asa

Footnotes and references


Edit This Page