# 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 QS*GW* calculations.

### Table of Contents

- Preliminaries
- Introduction
- Part 1. Elemental Nb and elemental Fe in reconstructed unit cells
- Part 2. The Nb(2)/Fe(3)/Nb(2) superlattice
- Part 3. Longer period superlattices
- Additional exercises
- Footnotes and references

**Part 1. Elemental compounds**

Requires files *init.nb* and *init.fe*. Copy and paste them to your working directory, or do:

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

**Bulk Nb**

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

**Bulk Fe**

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

Pass checks

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

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

*Sections 1 and 2 make input files* for unrelaxed Nb(2)/Fe(3)/Nb(2) superlattice: *ctrl.nbfe*, *site0.nbfe*, *dpos30.nbfe*, *basp.nbfe*, *site10.nbfe*, *site20.nbfe*.

Setup:

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

Make site files for **dp3=0, -0.1, -0.2**, and crude estimate for relaxation (*dpos30.nbfe*)

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

Pass checks

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

*Section 3: Self-consistency and relaxation*

Setup:

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

Relax **dp3=0**

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

Relax **dp3=−0.1**

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

Relax **dp3=−0.2**

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

Pass checks

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

Make and draw the energy bands

```
lmchk ctrl.nbfe -vfile=20 '--syml~n=21~q=0,1,0,0,0,0,1/sqrt(2),0,0'
mpix -n 16 lmf ctrl.nbfe -vfile=20 --shorten=no -vbeta=.1 --rhopos -vfrzwf=t -vnit=1 --band~col=1:180,541:720~col2=181:540~fn=syml
echo -5,5,5,10 | plbnds -spin1 -fplot~sh~scl=.7~shft=-.28 -ef=0 -scl=13.6 -dat=up -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0 bnds.nbfe
open fplot.ps
echo -5,5,5,10 | plbnds -spin2 -fplot~sh~scl=.7~shft=.45 -ef=0 -scl=13.6 -dat=dn -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0 bnds.nbfe
open fplot.ps
```

**Part 3. Longer period superlattices**

Setup

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

Create file *dpos20.nbfe*

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

Make *siten.nbfe* and *posn.nbfe* for *n*=4,5,6,7,8,9 (*tcsh* style)

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

Selected checks

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

### Preliminaries

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

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

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

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

### Introduction

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

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

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

#### Basic strategy

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

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

Next the superlattice is made.

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

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

- Repeat the self-consistent calculation of the interface code with the ASA-LDA code
*lm*.*lm*is more approximate than the full LDA code*lmf*, but it is very efficient, and provides a stepping stone to the transport code,*lmpg*.

The interface has an additional complication in that open spaces appear at the interface, which the ASA must deal with by packing with additional empty spheres. - Use
*lmpg*to perform transmission calculations.

The prescription is in outline:

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

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

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

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

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

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

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

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

We need to assemble files in several contexts:

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

#### Structural information for elemental Nb and Fe

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

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

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

*Notes:*

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

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

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

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

#### Reconstructed Nb, self-consistent LDA

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

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

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

`--nfile`

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

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

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

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

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

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

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

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

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

```
lmplan -vfile=3 nb
```

```
SPLCLS: 1 species split into 2 classes
Species Class Sites...
Nb 1:Nb 1
2:Nb2 2 3
Class Qtot Qbak Vmad Vh(Rmax) V(Rmax)
Nb 0.000000 0.000000 0.000000 0.000000 -0.700000
Nb2 0.000000 0.000000 0.000000 0.000000 -0.700000
Sum Q=0.000000 Emad=0.000000(0.000000) Vmtz=-0.700000 (up-dn=-0.700000)
Gtplan: 1 different planes normal to ( 0.000000 0.000000 1.000000)
Projection of normal onto P1, P2, P3: 0.000000 0.000000 0.707107
ib Class Plane PL x-x.h y-y.h z-z.h h del h
1 1:Nb 1 0 0.00000 0.00000 0.00000 0.00000 0.70711
2 2:Nb2 1 0 0.70711 0.50000 0.00000 0.00000 0.00000
3 2:Nb2 1 0 1.41421 1.00000 0.00000 0.00000 0.00000
```

*Notes:*

- The planes defined by the first two vectors
**P**_{1}×**P**_{2}stack lie normal to the [0,0,1] direction (‘normal’ vector in the output)). If`--rot=(0,1,0)pi/4`

had been omitted, the normal direction would be [1,0,1]. The*z*axis is in fact the [1,0,1] direction of the bcc lattice. - All three atoms lie in a single (001) plane.
- Spacing between planes is 0.70711, which is 1/√2.
- The 48 symmetry operations get reduced to 8, and the three Nb atoms are no longer symmetry-equivalent. This is artificial since this superlattice is still elemental Nb, but it is unimportant for our purposes.

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

```
lmfa -vfile=3 ctrl.nb
```

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

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

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

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

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

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

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

Make the density self-consistent

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

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

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

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

##### Energy bands of reconstructed Nb

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

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

Create *syml.nb* with

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

*syml.nb* should be created in symmetry lines format.

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

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

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

You should see the following table:

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

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

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

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

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)),y:pi/4' -vfile=2 fe --wsite~short~fn=site3
```

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

```
lmplan -vfile=3 fe
```

*Notes:*

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

You should see this line

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

This gives us an optimal mesh to select.

*lmfa* supplies the value **8.4** for **gmax**.

Remake the *ctrl* file with this information

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

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

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

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

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

##### Energy bands of reconstructed Fe

For the Mulliken decomposition make the table or orbital positions

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

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

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

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

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

Bands near *E _{F}* are dominated by states of

*d*character; there is some

*p*character as well. Note that within a given spin channel the

*d*states cluster into two groups, yielding the classical two-peak structure in the density-of-states.

To compare against the ASA, see this tutorial.

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

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

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

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

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

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

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

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

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

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

Use the superlattice stack maker to make the following superlattice:

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

```
Gtplan: 7 different planes normal to ( 0.000000 0.000000 1.000000)
Projection of normal onto P1, P2, P3: 0.000000 0.000000 4.677967
ib Class Plane PL x-x.h y-y.h z-z.h h del h
1 1:Nb 1 0 0.00000 0.00000 0.00000 0.00000 0.70711 ← layer Nb-L1
2 1:Nb 1 0 0.70711 0.50000 0.00000 0.00000 0.00000
3 1:Nb 1 0 1.41421 1.00000 0.00000 0.00000 0.00000
4 1:Nb 2 0 -0.70711 0.00000 0.00000 0.70711 0.70711 ← layer Nb-L2
5 1:Nb 2 0 0.00000 0.50000 0.00000 0.70711 0.00000
6 1:Nb 2 0 0.70711 1.00000 0.00000 0.70711 0.00000
7 2:Fe 3 0 -1.41421 0.00000 0.00000 1.41421 0.70711 ← layer Fe-1
8 2:Fe 3 0 -1.41421 0.75000 0.00000 1.41421 0.00000
9 2:Fe 3 0 -0.70711 1.00000 0.00000 1.41421 0.00000
10 2:Fe 3 0 0.00000 1.25000 0.00000 1.41421 0.00000
11 2:Fe 4 0 -1.77016 -0.50338 0.00000 2.03073 0.61651 ← layer Fe-2
12 2:Fe 4 0 -1.77016 0.24662 0.00000 2.03073 0.00000
13 2:Fe 4 0 -1.06305 0.49662 0.00000 2.03073 0.00000
14 2:Fe 4 0 -0.35594 0.74662 0.00000 2.03073 0.00000
15 2:Fe 5 0 -2.12610 -1.00676 0.00000 2.64724 0.61651 ← layer Fe-3
16 2:Fe 5 0 -2.12610 -0.25676 0.00000 2.64724 0.00000
17 2:Fe 5 0 -1.41899 -0.00676 0.00000 2.64724 0.00000
18 2:Fe 5 0 -0.71189 0.24324 0.00000 2.64724 0.00000
19 1:Nb 6 0 -2.48205 -1.51014 0.00000 3.26375 0.61651 ← layer Nb-R2
20 1:Nb 6 0 -1.77494 -1.01014 0.00000 3.26375 0.00000
21 1:Nb 6 0 -1.06783 -0.51014 0.00000 3.26375 0.00000
22 1:Nb 7 0 -3.18915 -1.51014 0.00000 3.97086 0.70711 ← layer Nb-R1
23 1:Nb 7 0 -2.48205 -1.01014 0.00000 3.97086 0.00000
24 1:Nb 7 0 -1.77494 -0.51014 0.00000 3.97086 0.00000
```

This command performs the following sequence (a breakdown of the `--stack`

instruction is given in the next section

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

Try running *lmscell* in interactive mode. You can watch how it builds the stack one block at a time. Instructions **file**, **show**, **basp**, **wsite** are documented in the manual

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

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

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

#### Strategy to make practical Nb/Fe/Nb superlattices

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

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

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

This is greatly facilitated by restricting atomic relaxation to atoms on frontier planes. It is an approximation, but it is a reasonable one; and it can be checked by looking at the forces ramaining after the constrained relaxation. In any case the real material will be highly non-ideal, which is very difficult to model, so there is little point in trying to allow more general relaxations.

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

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

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

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

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

`-vdz0='1/sqrt(2)' -vdz=dz0/2-.61651/2`

assigns to variable**dz**half the difference in Nb/Nb and Fe/Fe spacing.

**dp1**,**dp2**,**dp3**aren’t used now, but later it will be used to ajust the translation of the top layer, as it is not determined by symmetry.`-vfile=3 ctrl.nb`

tells*lmscell*to start from the reconstructed 3-atom Nb cell (a single plane) as a base.`--stack`

invokes the superlattice editor. It is invoked in batch mode with**~**separating commands and does the following:`file@site3.nb`

stacks a second the Nb layer (*site3.nb*) on top of the original`file@site4.fe@dup=3@dpos=dp1,dp2,-dz+dp3/2`

stacks 3 reconstructed Fe layer on top of Nb. The shift causes the Nb-Fe spacing to take the average of the Nb-Nb and Fe-Fe spacings, plus a shift**p3/2**. There is also a lateral translation given by**dp1**,**dp2**.`file@site3.nb@dup=2`

stacks two Nb layers (*site3.nb*) on top. This completes the stack.`basp@fn=basp.nbfe@basp.nb@basp.fe`

tells*lmscell*to assemble a*basp*file from files of the constituent elements.`pad=dp1,dp2,dp3`

adds**dp1,dp2,dp3**to the third lattice vector. This will be used in Step 3.`shorp3`

adds linear combinations of the first two lattice vectors to the third one, to shorten it.`shorten@mode=2,2,1`

adds lattice vectors to site positions. In the**P**_{1}and**P**_{2}directions, they are shortened to the minimum length. In the**P**_{3}direction lattice vectors are added to the position is between 0 and 1 in units of**P**_{3}.`~wpos@fn=pos`

writes site positions to file*pos.nb*. It isn’t used in this section but will be in the a later one.`wsite@fn=sitein`

writes a*site*file.

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

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

You need to perform step 1 before doing this one.

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

Table 1’s second line shows the relative magnitude forces expected for the ideal, rigid lattice, as would apply if the internuclear forces were dominate by nearest neighbors. The rigid lattice must be relaxed, which this section accomplishes within some constraints. The idea is structure the stacks so the superlattice is built in a lego-like fashion, layer by layer. That way construction can be automatic and easily generated for any number of layers. This is greatly facilitated by restricting atomic relaxation to atoms on frontier planes. This is an approximation, but it is a reasonable one; and it can be checked by looking at the forces ramaining after the constrained relaxation. In any case the real material will be highly non-ideal, which is very difficult to model, so there is little point in trying to allow more general relaxations.

The main thing is that the frontier atoms are not pressed too close together or kept too far apart. The fact that the displacements are limited only to the interface greatly simplifies the adoption the lego-like, layer-by-layer construction, since each layer is essentially independent of the others.

The lattice can be straightforwardly relaxed from the rigid structure generated in the previous section We can improve on this by shifting frontier atoms around to push apart pairs with small bond lengths. This can be accomplished with the sphere overlap minimizer utility in *lmchk*, which penalizes small bond lengths. This should distribute the forces (reduce the largest ones, while increasing some of the smallest ones).

To do this, make a temporary input file from *sitein.nb* generated in section 1, where the spheres are allowed to overlap by as much as 16%. Then use the `--mino`

utility in *lmchk* and allow a subset of atoms to move, to minimize a function of the sphere overlaps.

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

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

As for arguments to `--mino`

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

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

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

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

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

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

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

Command-line switches do the following:

`--nfile`

tells*blm*to insert the following lines into the*ctrl*file. This makes it possible to select different structures from the command line, through the choice of the*site*file.

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

`--ewald=1d-12`

: is probably not necessary here, but setting a tight tolerance on Ewald summations for envelope functions keeps the overlap matrix positive definite, (a problem which may arise if cells become long).`--rdsite`

: tells*blm*to read structural data from*sitein.nbfe*`--omax=-.07`

: but because of the messy interface, lattice relaxations will be large. Keeping the augmentation spheres a little bit apart will be needed so that when the lattice is fully relaxed the sphere overlaps are not too severe.`--nosort`

: suppresses re-ordering of site positions`--shorten=no`

: suppresses automatic shortening of site positions by adding multiples of lattice translation vectors.`--nit=20`

: variable controlling the number of iterations stop at if convergence is not reached. If modifies the input file via tag**EXPRESS_nit**, an alias for ITER_NIT`--nk=10,3,2`

:*k*mesh for BZ_NKABC Brillouin zone integration. This choice spaces*k*along the three axes approximately equally.`--gmax=8`

: specify mesh density cutoff HAM_GMAX.*lmfa*can find it for you (which is how the value was arrived at).`--mag`

: tells*blm*to create an input file for a magnetic system.`--frzwf`

: Add a variable**frzwf**which if set, suppresses updating partial waves through the tag HAM_FRZWF`--wsite~ib=\!4:10,15:21~rlx=000~long~fn=sitex`

: writes a file*sitex.nbfe*sets the rlx flag to 000 for all sites that are*not*frontier layers.**4:10,15:21**is a list of all sites in the frontier layers, as noted above. The preceding shriek tells the list maker to include all sites*except*those listed.*Note*`--wsite~ib=\!4...`

is specific to*tcsh*. For*bash*, use`--wsite~ib=!4...`

.

In the *blm* instruction that follows these switches are also used:

`--rdpos=dpos30`

: the contents of*dpos30.nbfe*, to the rigid lattice`--wsite:none`

: tells*blm*not to write a*site*file`--rsham~el=-.1,-0.9~rscut=1.6`

: supplies information for the basis. This is mainly a placeholder for the future JPO form.`--dhftol=45`

: Adds convergence criteria (ITER_DHFTOL) for self-consistency. This allows charge convergence tolerance to be set loosely (see below)`--molstat~xtol=0~gtol=0.01`

: Adds DYN category for molecular statics, using gradients only as a check for convergence

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

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

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

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

*μ*for all the Fe atoms.

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

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

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

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

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

This step does the following:

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

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

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

##### 3a. Self-consistency

Before relaxing the lattice it is useful to obtain a self-consistent density to get a sense of idea of what is transpiring at the interface and the forces there.

*Note* this step takes about 15 minutes to execute. The number of MPI processes cannot exceed 14.

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

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

Table 2 ib estatic eigval total 1 3.39 9.66 34.57 -6.42 -14.66 -60.34 -3.03 -4.99 -25.78 2 -3.07 -6.81 34.65 5.34 8.88 -61.64 2.27 2.07 -26.99 ← layer Nb-L1 3 -0.75 -1.01 27.46 0.00 2.36 -35.75 -0.67 1.35 -8.28 4 -29.70 -20.99 -78.37 33.25 15.17 108.43 3.55 -5.82 30.06 5 0.10 0.48 -182.65 0.00 0.12 240.87 0.00 0.60 58.22 ← layer Nb-L2 6 30.71 22.91 -82.57 -34.05 -17.19 114.25 -3.35 5.71 31.68 7 1.64 -1.35 3.05 -1.56 1.07 -16.26 0.00 -0.28 -13.21 8 -5.62 -31.16 39.33 1.55 35.42 -60.10 -4.07 4.26 -20.78 ← layer Fe-1 9 -1.30 -1.20 80.03 1.78 0.67 -90.98 0.47 -0.53 -10.95 10 4.23 29.18 38.91 -0.48 -34.46 -53.70 3.74 -5.28 -14.79 11 3.38 3.21 -22.98 0.83 -0.34 34.78 4.21 2.87 11.80 12 2.93 1.58 1.12 -3.14 -1.93 0.00 -0.22 -0.34 1.17 ← layer Fe-2 13 -0.99 -1.30 -20.08 -4.31 -1.36 30.92 -5.30 -2.67 10.84 14 -1.06 -2.40 -21.56 -1.21 2.74 42.84 -2.27 0.35 21.28 15 17.44 -1.90 20.65 -24.60 3.73 -21.01 -7.16 1.84 -0.36 16 0.00 0.41 46.69 -1.66 -0.42 -65.50 -1.66 0.00 -18.81 ← layer Fe-3 17 -24.83 3.59 29.84 39.46 -6.28 -36.62 14.63 -2.69 -6.78 18 -1.70 -2.98 47.76 3.10 4.28 -56.41 1.40 1.30 -8.64 19 1.05 -2.00 2.53 -2.76 4.02 -14.12 -1.71 2.02 -11.59 20 6.47 14.87 23.61 -12.21 -33.38 -42.69 -5.74 -18.50 -19.08 ← layer Nb-R2 21 -1.08 -11.54 10.66 2.62 29.99 -22.95 1.54 18.45 -12.29 22 -1.20 2.24 -9.56 0.86 -0.52 15.88 -0.34 1.72 6.33 23 2.43 1.08 -11.73 0.40 -6.66 23.74 2.83 -5.59 12.01 ← layer Nb-R1 24 -2.53 -4.57 -11.38 3.22 8.77 26.29 0.69 4.19 14.91

The largest force (site 5) is about 58 mRy/a.u. We can expect the frontier atoms would have much larger forces had the calculation been done for the rigid lattice instead. (To confirm this redo the calculation omit the `--rdpos=dpos30`

switch. You should find that the largest force is ~80 mRy/a.u. (site 9), while the end layers Nb-L1 and Nb-R1 would have smaller forces.)

##### 3b. Constrained relaxation

(Complete step 3a before doing this one.)

Use *lmf* to perform the (constrained) lattice relaxation. Note that the *site* files were constructed so that only atoms in the frontier layers are allowed to relax.

*Note* this step takes about 2 hours to execute.

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

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

`-vconv=1e-4`

: Set a tolerance rougher than the default for convergence in the energy before shifting atoms`-vconvc=1e-3`

: Set a tolerance rougher than the default for convergence in the charge density before shifting atoms. This is further explained below.`-vbeta=.1`

: specify the charge mixing parameter**beta**.**beta=0.1**is probably conservative but it is safe.`-vfrzwf=t`

: suppresses updating of the shape of the partial waves. The constraint is not essential, but it stabilizes convergence and most of the change is already realized in the initial self-consistent cycle which did allow the potential to float. (Keeping it frozen to the shpae of the free atom is essentially what a pseudopotential does.) Since the atomic relaxation will only slightly modify the this potential, it is safer and more stable with little loss of accuracy. Also without this switch the forces are not quite exact. Finally the error it can always be checked after the relaxation by letting the potential float once more.`-vminx=10`

: Perform at most 10 relaxation steps. It is sufficient to fully relax the lattice, as the discussion shows.`-wpos=posx`

: Writes the position file after every move, including the last one (the positions for the ersatz last+1 iteration).

Try `grep RELAX out.0`

to see how the gradient evolved in successive relaxation steps:

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

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

)

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

Parameters in **DYN_MSTAT** govern *lmf*’s relaxation algorithm. Tokens are explained in the input file; you can also get a quick summary with `lmf --input | grep -A 3 DYN_MSTAT`

.

**DYN_MSTAT_MODE=6**invokes a Broyden method.- Convergence is checked only on the forces (gradients) since
**XTOL=0**. **STEP=.010**limits the steps size to a relatively conservative displacement. It could be made larger, but the disruption to the potential becomes greater and it is safer not keep displacements fairly small.- Because of the constraints (see discussion of
**rlx**above) only frontier atoms are allowed to move. There are 10 and 14 Nb and Ni atoms per layer; thus there are 48×3 = 144 degrees of freedom in the relaxation.

*lmf* has two nested loops: the inner loop converged the charge at the given atomic positions, the outer loop relaxes the nuclei. This disrupts the potential, so that it has to be converged again at the updated positions. You can see this by doing `grep 'Maximum Harris' out.0`

.

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

Lines show maximum force on any site and an estimate for the correction to this force.

- When that correction is on the order of the force itself, usually the density is converged enough that there is little advantage to driving self-consistency much farther. Correction is small in first two lines above because the density was already nearly self-consistent from the initial run of
*lmf*. - In the third line the atoms had been moved; it is the first iteration of the inner loop. The correction is huge and the forces not very meaningful.
- The next line [last iter rlx 1] is the last step in the inner loop. In this case
*lmf*has no difficulty converging. If self-consistency has not been reached by the last iteration the*save*file’s line for that relaxation step would begin with**x**, rather than**c**as it does in this case. - The energy and charge convergence tolerances (
**conv=1e-4 convc=.001**) are set fairly loosely. This is safe because an extra condition is imposed (note**dhftol**in the*ctrl*file). The correction term to the Harris-Foulkes forces must fall below this value. 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 refined is all we need at this stage. As it happens in this case the relaxation proceeds smoothly anyway.

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

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

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

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

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

For the meaning of `--rs=11,1,1`

, run `lmf --h`

or see command line arguments documentation.

It is worthwhile to check that sphere overlaps after relaxation. Do `lmchk ctrl.nbfe -vfile=0 --rpos=pos0 --pr20`

and you should see this line

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

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

##### 3c. Optimizing the Nb/Fe layer spacing

(Complete step 3b before doing this one.)

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

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

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

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

`-av:nr,3 p3max`

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

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

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

You can confirm that the ansatz encoded in *posnow.nbfe* was pretty reasonable : `egrep -E '^[Ccx]' out.10`

shows the energy gained by relaxation was only 3 mRy.

Do the steps once again for **dp3=−0.2**

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

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

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

If the Fe local magnetic *M* retained its bulk value, we would expect *M*=12×2.2 *μ _{B}* = 26.4

*μ*for the superlattice, since it has 12 Fe atoms. The actual moment is somewhat less than that. To see a table of individual moments, do

_{B}```
lmf ctrl.nbfe -vfile=20 --rsedit~rs~a
```

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

Two features stand out:

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

Extract a table total energy as a function of **p3**

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

You should see

```
0 -122072.4927374
10 -122072.5281779
20 -122072.5391974
```

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

```
echo q | pfit 2 dat
```

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

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

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

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

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

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

Draw the bands:

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

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

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

Rather remarkably, there is little Fe character at *E _{F}*: both groups both fall

*E*for the majority spin, and for the minority spin

_{F}*E*falls in the pseudogap.

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

*k*is completely smeared out over the Fermi surface is pretty far removed from the actual situation.

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

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

Majority and minority DOS of NbFe resolved by layer.

Each of the 7 panels displays majority DOS (indicated by *D*(*E*)>0) and a minority DOS (indicated by *D*(*E*)<0).

The seven panels correspond to the site groupings of Table 3

layer Nb-L1

layer Nb-L2

layer Fe-1

layer Fe-2

layer Fe-3

layer Nb-R2

layer Nb-R1\

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

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

### Part 3. Longer period superlattices

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

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

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

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

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

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

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

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

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

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

### Additional exercises

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

Choose some *n* and do

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

You should find

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

the magnetic moment is approximately linear in

*n*, with a change in*M*between 8 and 9 as a layer is added. The value of*M*at the frontier changes a little as*n*increases, but*M*in atoms removed from the frontier are consistently between 2.15 and 2.25.the total energy is very nearly linear in

*n*.

### Footnotes and references

A detailed description of Questaal and *lmf* can be found in the following reference:

Dimitar Pashov, Swagata Acharya, Walter R. L. Lambrecht, Jerome Jackson, Kirill D. Belashchenko, Athanasios Chantis, Francois Jamet, Mark van Schilfgaarde, *Questaal: a package of electronic structure methods based on the linear muffin-tin orbital technique*, Comp. Phys. Comm. **249**, 107065 (2020).