# Nb/Ni Superlattice

The electronic structure of a Nb(110)/Ni(111) superlattice is created and the interface relaxed, using Questaal’s superlattice editor *lmscell*.

JMRAM devices typically consist films of Nb, oriented on the (110) direction, followed by a Cu spacer layer, followed by Ni, with another Cu spacer layer and finally capped with Nb. Cu and Ni are oriented along (111). Ni and Cu are both fcc with room temperature lattice constants 3.524 Å and 3.615 Å, respectively. Nb on the other hand has a lattice constant of 3.300 Å. To make an A/B superlattice unit cells of the A and B regions in the plane normal to the interface must be coincident. This is a challenging requirement, in part because the lattices are not matched, but also the (110) and (111) are not compatible.

The tutorial begins by rotating Nb [110] axis and the Ni [111] axis respectively to [001]. Then superlattice vectors are identified that are nearly coincident. By applying a small shear to the Ni, the lattice vectors are made exactly coincident.

At this point superlattices can be constructed. Some are made and the site positions at the frontier layers relaxed.

A separate tutorial takes the structural information made in this tutorial to investigate Landauer-Buttiker transport through a Nb/Ni/Nb trilayer.

*Note: at present that tutorial is written for Nb[110]/Ni(110)/Nb([110] trilayer*

### Table of Contents

- Nb/Ni(3
*n*)/Nb superlattices - Preliminaries
- Introduction
- 1. Reconstructed unit cells for elemental Nb and Ni
- 2. Lattice structure of the Nb(
*m*)/Ni(*n*)/Nb(*m*) superlattice - 3. Relaxing the Nb(2)/Ni(3)/Nb(2) superlattice
- Footnotes and references

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

** Section 1. Coincident site Nb and Ni site lattices**

Creates files *site20.nb* and *site42.ni*, and *spec.nb* and *spec.ni*, and also *sitea.nb*, *siteb.nb* and *sitea.ni*, *siteb.ni* and *sitec.ni*.

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

```
touch init.nb; rm *.nb *.ni;
cp ../input.fp.Nb110.Ni111/init.nb ../input.fp.Nb110.Ni111/init.ni .
```

Make *site20.nb*,*spec.nb*, *sitea.nb* and *siteb.nb*

```
blm --nfile --ctrl=ctrl --omax=.14 nb
lmscell --plx~m~-1,-4,0,2,-2,0,1,1,2 '--rot=(1,-1,0)-pi/2,z:pi/4,z:acos(5/sqrt(43))' nb --wsitex~fn=site2 --sort:"x3 x2"
sed -i.bak 's:\(0.00000001\):0.0:g' site2.nb
cp site2.nb site20.nb
lmchk ctrl.nb --wspec > out.lmchk.nb
lmscell -vfile=20 ctrl.nb '--stack~rmsite@targ=11:20~pad@0,0,-1/sqrt(2)~wsite@fn=sitea'
lmscell -vfile=20 ctrl.nb '--stack~rmsite@targ=1:10~addpos@dx=0,0,-1/sqrt(2)~pad@0,0,-1/sqrt(2)~wsite@fn=siteb'
```

Make *site42.ni*,*spec.ni*, *sitea.ni*, *siteb.ni* and *sitec.ni*

```
blm --nfile --ctrl=ctrl --omax=.14 ni
lmscell '--rot=z:pi/4,x:acos(sqrt(1/3)),z:acos(4/sqrt(19))' --plx~m~3,0,2,-1,0,4,-1,3,-1 ni --wsitex~fn=site2 --sort:"x3 x2"
echo 3.27871926 0.0 0.0 1.82998284 2.15665546 0.0 0.0 0.0 1.41421356 | mcx -f9f12.6 . -s6.226 | sed 's/rows 1 cols 9 real/rows 3 cols 3 real/' > pnb
echo 3.082207 0.0 0.0 1.62221421 1.9668302 0.0 0.0 0.0 1.73205081 | mcx -f9f12.6 . -s6.643 | sed 's/rows 1 cols 9 real/rows 3 cols 3 real/' > pni
cp ctrl.ni ctrl.ni.bak
echo STRUC >> ctrl.ni
echo ' DEFGRD= ' `mcx -f9f12.8 pni -i pnb -x -t -a:nr=1 s s -w:nohead . | sed 's:[0-9.]*$:1:'` >> ctrl.ni
lmchk -vfile=2 ni --wsitex~fn=site3
cp ctrl.ni.bak ctrl.ni
lmchk -vfile=2 ni | grep vol; lmchk -vfile=3 ni | grep vol
sed -i.bak 's:\(plat.*[1-9.]*\)$:\1/1.00983:' site3.ni
lmchk -vfile=2 ni | grep vol; lmchk -vfile=3 ni | grep vol
lmchk -vfile=3 ni --wsitex~scala=6.226/6.643~fn=site42
sed -i.bak 's:\(plat.*1.83\).*$:\1:' site42.ni
lmchk ctrl.ni --wspec > out.lmchk.ni
lmscell -vfile=42 ctrl.ni '--stack~rmsite@targ=15:42~pad@0,0,-2*.61~wsite@fn=sitea'
lmscell -vfile=42 ctrl.ni '--stack~rmsite@targ=1:14,29:42~addpos@dx=0,0,-.61~pad@0,0,-2*.61~wsite@fn=siteb'
lmscell -vfile=42 ctrl.ni '--stack~rmsite@targ=1:28~addpos@dx=0,0,-2*.61~pad@0,0,-2*.61~wsite@fn=sitec'
```

** Section 2. Initial estimate for Nb(110)/Ni(111)**

Creates files *site20.nb* and *site42.ni*.

#### Nb/Ni(3*n*)/Nb superlattices

This section requires the files *ctrl.nb*, *site20.nb*, *spec.nb*, *site42.ni*, *spec.ni*.

Set initial conditions with

```
touch spec.nb; rm *.nb *.ni
cp ../input.fp.Nb110.Ni111/{ctrl.nb,site20.nb,spec.nb,site42.ni,spec.ni} .
```

Confirm Nb/Ni and Ni/Nb have equal spacing

```
lmscell -vfile=20 ctrl.nb '--stack~file@sitea.ni@dpos=0,0,-1/2/sqrt(2)+.61/2~wsite@fn=sitex'
cp sitex.nb sitein.dat
blm --nfile --rdsite=sitein --rdspec~fn=spec.nb~fn=spec.ni --usefiler --nosort --shorten=no dat
cp actrl.dat ctrl.dat
lmplan dat
```

### Preliminaries

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

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

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

### Introduction

This tutorial builds a Nb(110)/Ni(111) superlattice, using Questaal tools *blm*, *lmscell* and *lmchk*. Since Nb[110] and Ni[111] are poorly lattice matched, unit cells in the plane normal to [110] and [111] respectively must be found that are nearly coincident. Within the coincident lattice, bonds across the interface are very irregular, and there is a large attendant lattice relaxation which must be dealt with. This is done in two steps: a trial relaxation that minimizes some measure of overlaps between atomic density, and subsequently using the LDA code *lmf* to relax the structure by minimizing the forces.

#### Basic strategy

To stack unit cells of Nb and Ni on top of one another, the two lattice vectors normal the stacking plane must match (they must be coincident). This entails:

- a reconstruction of bcc Nb with 10 atoms in a plane normal to [110]. The Nb lattice is also rotated orienting the [110] axis along
*z*. - a reconstruction of fcc Ni with 14 atoms in a plane normal to [111]. The Ni lattice is also rotated orienting [111] axis along
*z*. - The Ni lattice is sheared slightly to make Nb and Ni exactly coincident in the plane normal to
*z*. It is also strained along*z*to partially accommodate the volume change from the in-plane shear, according to Poisson’s ratio.

The prescription is in outline:

- Use
*blm*to build input files for elemental Nb and Ni in reconstructed lattices. - Use the supercell maker
*lmscell*in combination with*blm*to generate different reconstructions of the elemental lattices that can be joined together. - Use the superlattice editor
`lmscell --stack`

to stack reconstructed Nb and Ni lattices into superlattices. - Use
*lmchk*to make a rough estimate of the atomic relaxation of frontier atoms at the Nb/Ni interfaces, followed by the LDA code*lmf*to perform a constrained relaxation the interface, allowing atoms in the frontier planes to relax until forces on them are small.

### 1. Reconstructed unit cells for elemental Nb and Ni

This tuturial requires only two *init* files corresponding to element Nb and Ni to complete.

This tutorial concludes with a self-consistent *finish*

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, Ni
- Ni/Nb superlattices, built so that the interface can be relaxed in a constrained manner. The constraint is necessary for later use, when superlattices are converged into trilayer (… Nb Nb / Ni / Nb Nb …) suitable for study of transport of an active region. The semi-infinite left and right leads are required to have a periodic lattice.

#### 1a. Structural information for elemental Nb and Ni

Copy the boxes below into files *init.nb* and *init.ni*. They contain the requisite structural information for elemental bcc Nb and fcc Ni. Here we use the 0K lattice constants.

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

and

```
# Init file for Ni. Lattice constant at 0K
LATTICE
ALAT=6.643 PLAT= .5 .5 .0 .5 .0 .5 .0 .5 .5
SPEC
ATOM=Ni MMOM=0,0,0.6 PZ=0,0,0
SITE
ATOM=Ni X=0 0 0
```

#### 1b. Reconstructed Nb, and Ni, coincident site lattices

##### Initial reconstruction of Nb

Nb is bcc, and it will be rotated so the [110] direction points along *z* which will be axis of the superlattice.

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

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

`--nfile`

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

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

The other switches supply information *blm* cannot determine on its own: the *k* mesh and the plane wave cutoff for the interstitial density.

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

Next we rotate the crystal with the following (admittedly odd) construction. The total rotation consists of a sequence of three separate ones: first, a rotation of the [110] axis to the [001], then a 45 degree rotation the new *z* axis, finally another rotation around the same axis. `--rot=(1,-1,0)-pi/2`

rotates [110] to *z* which will be the superlattice axis. The latter rotations will be made clear when we align Ni and Nb.

```
lmscell --noplx '--rot=(1,-1,0)-pi/2,z:pi/4,z:acos(5/sqrt(43))' nb --wsitex~fn=site1
```

Next, make a superlattice of this rotated cell. It reads site information from *site1.nb* which *lmscell* generated in the previous instruction.

```
lmscell -vfile=1 --plx~m~-1,-4,0,2,-2,0,1,1,2 nb --wsitex~fn=site2 --sort:"x3 x2"
```

This superlattice is artfully constructed from a linear combination of integer multiples bcc lattice vectors. These particular multiples yield the following lattice vectors:

```
plx plx/plat qlx
3.27872 0.00000 0.00000 -1.00000 -4.00000 0.00000 0.30500 -0.25880 0.00000 ← PLAT(1)
1.82998 2.15666 0.00000 2.00000 -2.00000 0.00000 0.00000 0.46368 0.00000 ← PLAT(2)
0.00000 0.00000 1.41421 1.00000 1.00000 2.00000 0.00000 0.00000 0.70711 ← PLAT(3)
```

This superlattice consists of two planes 10 atoms each. The first plane passes through the origin, and the second through , which is the spacing between (110) planes in the bcc lattice. To confirm this, do

```
lmplan -vfile=2 nb --pledit~q
```

bcc planes are stacked in an ABAB… pattern. A and B planes are equivalent but B is translated in the lateral direction. The third plane at a height , coincides with the first without a lateral translation. Thus this rotation/superlattice combination puts **PLAT(3)** along [001] with a length .

##### Initial reconstruction of Ni

Ni is fcc, and it will be rotated so the [110] direction points along *z*.

Make a template *ctrl* and *site* file for fcc Ni:

```
blm --nfile --ctrl=ctrl ni
```

In the Ni case we will rotate the [111] to *z* using `--rot=z:pi/4,x:acos(sqrt(1/3))`

. An additional rotation `z:acos(4/sqrt(19))`

is made for reasons that will soon become clear. We also need a superlattice of this rotated cell, with different superlattice vectors than the Nb case. We can accomplish both the rotation and the superlattice in one step as follows:

```
lmscell '--rot=z:pi/4,x:acos(sqrt(1/3)),z:acos(4/sqrt(19))' --plx~m~3,0,2,-1,0,4,1,1,1 ni --wsitex~fn=site1 --sort:"x3 x2"
```

Compare the first lines of *site1.ni* and *site2.nb*. The lattice vectors **PLAT(1)** and **PLAT(2)** (when scaled by their respective lattice constants) are very similar.

This particular superlattice consists of 14 atoms all in a single plane normal to the new *z*, which is [111] in the unrotated frame. To confirm this, do

```
lmplan -vfile=1 ni --pledit~q
```

Ni is fcc whose planes are stacked ABC. We will use as a basic unit a superlattice with three planes (42 atoms) that include all the ABC planes. As needed we will split the ABC layers into single planes.

```
lmscell '--rot=z:pi/4,x:acos(sqrt(1/3)),z:acos(4/sqrt(19))' --plx~m~3,0,2,-1,0,4,-1,3,-1 ni --wsitex~fn=site2 --sort:"x3 x2"
```

You can confirm that **PLAT(1)** and **PLAT(2)** in files *site1.ni* and *site2.ni* are the same; only the third axis now points along *z*.

*lmplan* confirms that the new cell of 42 atoms consists of 3 planes, normal to [111] in the unrotated frame.

```
lmplan -vfile=2 ni --pledit~q
```

You can confirm that sites 15:28 are uniformly translated relative to sites 1:14, and sites 29:42 are translated by twice that amount. Each plane is spaced by in units of the lattice constant.

##### Nb and Ni, coincident site lattices

Note that *lmscell* created the superlattices using `--wsitex`

instead of `--wsite`

. This saves the basis vectors in the form of fractional multiples of the lattice vectors rather than Cartesian coordinates. This will be important for the shear which we need to do next, because to stack Nb and Ni into a superlattice **PLAT(1)** and **PLAT(2)** have to be exactly coincident.

From each of the files *site2.nb* and *site2.ni* , extract the lattice vectors and scale by their respective lattice constants.

```
echo 3.27871926 0.0 0.0 1.82998284 2.15665546 0.0 0.0 0.0 1.41421356 | mcx -f9f12.6 . -s6.226 | sed 's/rows 1 cols 9 real/rows 3 cols 3 real/' > pnb
echo 3.082207 0.0 0.0 1.62221421 1.9668302 0.0 0.0 0.0 1.73205081 | mcx -f9f12.6 . -s6.643 | sed 's/rows 1 cols 9 real/rows 3 cols 3 real/' > pni
```

The `sed`

instruction is appended to render *pnb* and *pni* into the 3×3 matrices they actually are. Inspect them and you should see:

```
% rows 3 cols 3 real
20.413306 0.000000 0.000000 11.393473 13.427337 0.000000 0.000000 0.000000 8.804894
% rows 3 cols 3 real
20.475101 0.000000 0.000000 10.776369 13.065653 0.000000 0.000000 0.000000 11.506014
```

This shows that the superlattice vectors are nearly coincident, with the first vector in each case pointing along *x*. This explains why these particular superlattice vectors and rotations were chosen.

We must shear one or the other lattice to render the first two vectors **PLAT(1)** and **PLAT(2)** exactly coincident. (**PLAT(3)** has no such constraint and we will deal with it later.) We elect to shear the Ni lattice, since real Ni layers are disordered anyway.

Do the following:

```
mcx -f3f12.8 pni -i pnb -x
```

You should get

```
0.99698194 0.00000000 0.00000000
0.04972026 1.02768205 0.00000000
0.00000000 0.00000000 0.76524277
```

This matrix tells us how much we have to shear Ni to match Nb. For now we are interested only in the 2×2 subblock since the third axis is not dictated by matching planes of the two lattices. For now let’s restrict the deformation to the (**PLAT(1)**, **PLAT(2)**) in the plane normal to *z*.

Temporarily add these lines to *ctrl.ni*:

```
STRUC
DEFGRD= 0.99698194 0.00000000 0.00000000 0.04972026 1.02768205 0.00000000 0.00000000 0.00000000 1
```

and create a new file *site3.ni* with the transformed structural data.

```
lmchk -vfile=2 ni --wsitex~fn=site3
```

Remove or comment out **DEFGRD** before continuing.

Token **DEFGRD** supplies a general linear transformation of the lattice: it scales both lattice and basis vectors by a 3×3 matrix. Compare lattice vectors of *site2.nb* and *site3.ni*: you should see that (**PLAT(1)**, **PLAT(2)**) are coincident (remember that **alat** is different in the two cases)

```
head -1 site2.nb; head -1 site3.ni
```

Finally we must decide how much to scale **PLAT(3)**. Were Ni incompressible we would know how much we should strain Ni along *z*: any volume change generated by **DEFGRD** would be compensated by a tensile strain along *z*. To determine the volume change induced by **DEFGRD** do the following (be sure you have removed or commented out the **DEFGRD** line first!):

```
lmchk -vfile=2 ni | grep vol; lmchk -vfile=3 ni | grep vol
```

You should find that the volume changed by this ratio:

```
dval 3153.75620/3078.09526 = 1.02458
```

As a sanity check, copy *site3.ni* to *site3.ni.bak* and scale the 9^{th} element of **plat** by 1/1.02458, or simply use *sed*:

```
sed -i.bak 's:\(plat.*[1-9.]*\)$:\1/1.02458:' site3.ni
```

Run *lmchk* again and confirm the unit cell volumes of *site2.ni* and *site3.ni* are the same. This is the required strain for the incompressible case, which corresponds to a Poisson ratio of 0.5. The opposite limit, Poisson ratio of 0, means that change in lateral area has no effect on tensile strain along *z* (cork is approximately like that). You can look up Poisson ratios for elements here : https://periodictable.com/Properties/A/PoissonRatio.v.html; it is 0.31 for Ni, close to 3/5.

Since Ni isn’t incompressible, this strain is tool large. This volume increase should be reduced by 2/5, which corrresponds to a strain `dval 1+.02458*2/5`

= 1.00983 . Replace the 9^{th} element of **plat** with `dval 1.73205081/1.00983`

= 1.71519. To do it on the command-line, invoke

```
cp site3.ni.bak site3.ni
sed -i.bak 's:\(plat.*[1-9.]*\)$:\1/1.00983:' site3.ni
```

Run *lmchk* as before and confirm the ratios of cell volumes are `dval 3123.05566/3078.09526`

= 1.01461 . Thus the ~2.5% strain is reduced to a ~1.5% strain.

To make the final reconstructed Ni *site* file, do the following

```
lmchk -vfile=3 ni --wsitex~scala=6.226/6.643~fn=site42
```

Command line argument `scala`

causes *lmchk* to scale the lattice constant **ALAT** by a factor and the lattice and basis vectors by its reciprocal. The structure is unchanged, only the distribution between **ALAT** and **PLAT** is changed. In this context **ALAT** becomes the same for Ni and Nb. The name *site42.ni* was used simply to mark how many atoms are in the unit cell.

Note **plat(3,3)** is 1.83006861, which is very close to 1.83, which is a multiple of 0.61. For convenience replace **plat(3,3)** with 1.83, since we will need to divide this cell into three subcells corresponding to each of the three planes. Use your text editor or use *sed*:

```
sed -i.bak 's:\(plat.*1.83\).*$:\1:' site42.ni
```

### 2. Lattice structure of the Nb(*m*)/Ni(*n*)/Nb(*m*) superlattice

The Nb(*m*)/Ni(*n*)/Nb(*m*) superlattice (*m* and *n* refer to numbers of planes in each region) has an Nb/Ni and an Ni/Nb interface; also periodic boundary conditions apply and there is a Nb/Nb interface as well. Provided the right layer ends in a B plane and the left begins with an A plane, the Nb/Nb interface will just be a continuum of 4 bcc Nb planes along *z* ([110] in the unrotated lattice). Periodic boundary conditions requires that cannot be translated with respect to each other, so we impose the constraint that the left Nb lead begins with an untranslated A plane and the right B lead ends with an untranslated B plane.

The Nb/Ni and Ni/Nb interfaces are more complicated. The main difficulty is that stacking ideal Nb and Ni planes will bring some Ni-Nb bonds much too close together while leaving void regions. Atoms can be relaxed, but there are some constraints we want to follow. We would like to 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. 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.

When stacking Ni on Nb, we are immediately confronted with these issues:

- What is the spacing between Nb and Ni planes?

As a first guess we will take the average of the Nb-Nb and Ni-Ni spacing. - Should the Ni lattice be translated relative to the Nb?

In principle this question would be resolved by lattice relaxations. However, the interfaces are complicated and there maybe many metastable configurations where forces vanish.

To handle any number of Ni layers we need to consider three distinct configurations:

- the
**Nb/**(**A’B’C’**)_{n}**/Nb**family with 3*n*Ni planes - the
**Nb/**(**A’B’C’**)_{n}**A’/Nb**family with 3*n*+1 Ni planes - the
**Nb/**(**A’B’C’**)_{n}**A’B’/Nb**family with 3*n*+2 Ni planes

Within any family the interfaces will be the same, but between families they will not, and each family must be dealt with independently of the others.

#### Nb/Ni(3*n*)/Nb superlattices

This section makes a *site* file for the **Nb(AB)/**Ni(**A’B’C’**)**/Nb(AB)** superlattice, the simplest lattice with 3*n* Ni planes. It requires the files *site20.nb*, *spec.nb*, *site42.ni*, *spec.ni*.

If the files below are not already present in your working directory create initial conditions with

```
cp ../input.fp.Nb110.Ni111/{ctrl.nb,site20.nb,spec.nb,site42.ni,spec.ni} .
```

We will distinguish the three kinds of lattices by denoting them as 3*n*, 3*n*+1 and 3*n*+2 superlattices, respectively. This section builds the simplest superlattice in the 3*n* family. Within this family even without any translations or reconstructions, there are four possibilities: on the left side the frontier Nb layer can be either **A** or **B**; similarly for the right side.

We take the most obvious choice: **AB|A’B’C’|AB**. It turns out to be the best choice among the four possibilities. Build the stack with the superlattice editor *lmscell* starting with the reconstructed Nb lattice as a base.

```
lmscell -vfile=20 ctrl.nb -vdz='1/2/sqrt(2)-.61/2' --stack~file@site42.ni@dpos=0,0,-dz~file@site20.nb~wsite@fn=sitex
```

This makes the **AB|A’B’C’|AB** superlattice, using the average of Nb-Nb and Ni-Ni spacings as the spacing between Nb and Ni layers. The superlattice structure is written to file *sitex.nb*.

`-vfile=20 ctrl.nb`

tells*lmscell*to build the stack starting from the recontructed 20-atom Nb cell as a base.The Nb-Nb spacing is and the Ni-Ni spacing is . is half this difference.

`--stack`

invokes the superlattice editor. It is invoked in batch mode with**~**separating commands and does the following:`~file@site42.ni@dpos=0,0,-dz`

stacks the 42-atom reconstructed Ni layer on top of Nb. Absent the modifier**dpos**, the first Ni plane would be placed at the boundary of the 20-atom reconstructed Nb layer, which lies at a height above the topmost Nb layer. Rigidly shifting the 42 Ni atoms by causes the Nb-Ni spacing to take the average of the Nb-Nb and Ni-Ni spacings.`~file@site20.nb`

stacks the 20atom recontructed Nb layer on top of the stack. Since the Ni was shifted by , the second Ni-Nb is automatically at the average of Ni-Ni and Nb-Nb spacings.`~wsite@fn=sitex`

writes the file*sitex.nb*, which is a*site*file of the superlattice.

Use *blm* to construct a temporary *ctrl* and *site* file with extension *.dat*. At this stage we only want a *ctrl* file with enough information to examine the lattice structure. *blm* can read structural information from a *site* file, a feature which we use here:

```
cp sitex.nb sitein.dat
blm --nfile --rdsite=sitein --rdspec~fn=spec.nb~fn=spec.ni --usefiler --nosort --shorten=no --ctrl=ctrl dat
```

Command-line arguments to *blm* are documented on this page.

`--nfile`

was explained above`--rdsite=sitein`

tells*blm*to read structural data from*sitein.dat*`--rdspec~fn=spec.nb~fn=spec.ni`

supplies*blm*with species-specific information gleaned from*spec.nb*and*spec.ni*, notably in this case sphere radii.`--usefiler`

tells*blm*to accept without change the muffin-tin radius (**R**) it read from*spec.nb*and*spec.ni*. It is useful in this context because will obtain a crude estimate for reconsructiong by minimizing overlaps between spheres on neighboring sites.`--nosort`

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

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

writes the template input file directly to*ctrl.dat*(instead of the default*actrl.dat*)

This lattice has 20+42+20 = 82 sites. Use *lmscell* to see the structure of stack

```
echo q | lmplan dat
```

Inspect the table below this line

```
ib Class Plane PL x-x.h y-y.h z-z.h h del h
```

**delh** indicates the change in projection along z compared to the preceding site. **delh** is always zero except for the first atom in a new plane:

```
1 1:Nb 1 0 0.00000 0.00000 0.00000 0.00000 0.70711
...
11 21:Nb11 2 0 2.05873 0.10783 0.00000 0.70711 0.70711
...
21 2:Ni 3 0 0.00000 0.00000 0.00000 1.36566 0.65855
...
35 30:Ni15 4 0 0.43389 0.05135 0.00000 1.97566 0.61000
...
49 58:Ni29 5 0 0.86779 0.10270 0.00000 2.58566 0.61000
...
63 41:Nb21 6 0 0.00000 0.00000 0.00000 3.24421 0.65855
...
73 61:Nb31 7 0 2.05873 0.10783 0.00000 3.95132 0.70711
```

This confirms that planes are stacked (Nb, Nb, Ni, Ni, Ni, Nb, Nb) with Nb-Nb spacings **0.70711**, Ni-Ni spacings **0.61** and Nb-Ni spacings **0.65855** (the average of **0.70711** and **0.61**). Note also the very first site has a spacing **0.70711**. This is the spacing between the last plane, sites 73:82, and the first plane, sites 1:10. It is , as it should be because of periodic boundary conditions.

Next we need to examine Ni-Nb bonds, and relax frontier sites to space bonds more evenly. The simplest way to do this is to minimize sphere overlaps. *blm* automatically determined sphere radii in the elemental compounds by enlarging the sphere radius to match the Wigner Seitz cell volume, and thus making them overlap: 14% in the bcc case, 11% in the close-packed fcc case. (*blm* imported these radii into *ctrl.dat* via `--rdspec`

and `--usefiler`

.) Such geometry violation is always present in ASA calculations, but in this context we want to use the sphere overlaps as a measure of how bonds are spaced.

Do the following

```
lmchk dat --wpos=pos0
```

Near the end you should see this line

```
OVMIN, 1426 pairs: fovl = 3.78727e-5 <ovlp> = 18.3% max ovlp = 36%
```

The maximum overlap is large, but slightly less than what it would be if we had chosen interfaces other than **AB|A’B’C’|AB**, e.g. **ABA|A’B’C’|AB** or **AB|A’B’C’|BAB**.

*lmchk* has a facility to move atoms that minimize sphere overlaps. It does this by minimizing an objective function (e.g. the sixth power of the sum of all positive overlaps). We will use this facility to get an initial estimate of the reconstruction. As noted above, we do not want to move all sites, only the frontier ones so we can build any 3*n* lattice lego-style. Do the following

```
lmchk dat --mino~dxmx=.001~xtol=.0001~nkill=10~maxit=200~11:34,49:72 --wsite~short~fn=site82 --wpos=pos
```

Note that only sites 11:34,49:72 move (compare to table above).

Command-line arguments to *lmchk* are documented on this page. Note that **~** separates arguments.

`--mino~dxmx=.001~xtol=.0001~nkill=10~maxit=200~11:34,49:72`

tells*lmchk*to move a selected list of sites to minimize the objective function.

The last argument selects the list (see here for syntax of integer lists).

**11:34,49:72**encompasses all the sites in the frontier layers (see table above).

**dxmx**,**xtol**,**nkill**and**maxit**are explained on the command line arguments web page.`--wsite~short~fn=site82`

tells*lmchk*to write*site*file*site82.dat*.`--wpos=pos`

tells*lmchk*to write site positions to file*pos.dat*.

Near the end you should find the following:

```
minimized: fovl = 6.96742e-6 <ovlp> = 13.8% max ovlp = 19.6%
```

The maximum overlap is reduced to 19.6% from an initial 36%. Since there is a heavy energy penalty for atoms being too close together, this is a much better starting point than the unrelaxed lattice.

*lmchk* wrote the unrelaxed and relaxed positions to files *pos0.dat* and *pos.dat* respectively To see what atoms moved do

```
mcx pos.dat pos0.dat --
```

Displacements are significant, on the order of 10% of the bond length. Note that if *all* the sites were included in the list, the maximum overlap would fall to 17.3% — not much lower than the constrained case. This provides a rough justification for constraining the relaxation to frontier atoms.

### 3. Relaxing the Nb(2)/Ni(3)/Nb(2) superlattice

The previous section supplied a suitable structure to perform a density-functional calculation of the Nb(2)/Ni(3)/Nb(2) superlattice.

Next we need to construct an input file suitable for *lmf*. *blm* will accomplish this reading structural information from *site82.dat* generated in the previous section.

```
cp site82.dat sitein.nbni
blm --autobas~loc=5 --rsham~el=-.3,-1.2~rscut=1.6 --dhftol=45 --molstat~xtol=0~gtol=0.01 --nfile --rdsite --omax=-.01 --nosort --shorten=no --nit=20 --nk=4,5,2 --gmax=10 --mixb=.1 --frzwf --wsite~ib=1:10,35:48,73:82~rlx=000~long~fn=site82 nbni
cp actrl.nbni ctrl.nbni
```

Command-line arguments to *blm* are documented on this page and more details are provided there.

`--autobas~loc=5`

Modifies AUTOBAS_LOC to designate low-lying local orbitals as of the conventional type.`--rsham~el=-.3,-1.2`

Adds category STR for real-space strux. This will play a role when the JPO basis is available.`--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`--nfile`

was explained above`--rdsite`

tells*blm*to read structural data from*sitein.nbni*`--omax=-.01`

select sphere radii so that the maximum sphere overlap is −1%`--nosort`

suppress re-ordering of site positions`--shorten=no`

suppress 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=4,5,2`

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

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

Specify the charge mixing parameter**beta**.**beta=0.2**was found empirically to be a more efficient choice than the default**beta=0.3**.`--frzwf`

Add a variable**frzwf**which if set, suppresses updating partial waves through the tag HAM_FRZWF`--wsite~ib=1:10,35:48,73:82~rlx=000~long~fn=site82`

writes a file*site82.nbni*sets the rlx flag to 000 for all sites that are not frontier layers. This restricts relaxation to frontier atoms.

With an input file (*ctrl.nbni*) and structure file (*site82.nbni*), we are in a position to perform self-consistent LDA calculations. We will have to repeat them in the ASA context, but the ASA cannot be used to relax the lattice. Thus this step is used mainly to obtain lattice displacments needed by the ASA. It also can serve as a benchmark to validate the quality of the ASA.

- This invocation of
*blm*makes a*ctrl.nbni*to be used in self-consistent LDA calculations. Some care now has to be taken concerning chemical attributes in addition to structural ones, e.g. choice of augmentation sphere radii.*blm*found good choices for the sphere radii automatically: they were made as large as possible subject to the constraint that overlap between any pair not exceed −1%.

To confirm this, try`lmchk ctrl.nbni -vfile=82`

These lines should appear in the output:`Cell volume= 7949.71676 Sum of sphere volumes= 4464.56181 (0.5616) ... OVMIN, 1404 pairs: fovl = 0 <ovlp> = 0% max ovlp = -1%`

See this page for general considerations in how to select sphere radii. Here we would like to make them as large as we can while avoiding inaccuracies from geometry violation when spheres overlap. As that page notes

*lmf*can tolerate a limited amount of overlap. By starting with and initial −1%, we may expect some positive overlaps to appear in the relaxation process. As we will see the maximum overlap will be about 5.5% when the system is relaxed, which is acceptable. Also the average sphere packing is 56%, which is quite healthy, so the basis should be accurate. **--molstat**adds provision in the input file for molecular statics, to relax the lattice to equilibrium positions. There is a tutorial for a simple case, which you might want to consult before continuing with this one. This is a restricted relaxation for reasons explained above Only the interfacial atoms will be free to move (information is fed to*lmf*via*site82.ni*, the**rlx**column).*blm*sets**rlx**to zero for all sites*not*on a frontier layer.- The
**k**mesh,**4,5,2**specifies the division of the reciprocal lattice into microcells. The requisite mesh density cannot be reliably predicted in advance: some trial and error is needed. This mesh was found to accurate be fairly accurate, certainly enough that the forces will be reliable. We are using tetrahedron integration and as the calculation prints out the “Bloechl correction,” which gives a measure of the mesh quality. (It is about 2 mRy/au in this case, which for 82 atoms is quite good.) The cell has 82 atoms: it is a bit expensive so we want keep the mesh relatively light.**4,5,2**was selected because it was found empirically to be dense enough to be accurate, and also space lengths of the edges of parallelpipeds as evenly as possible. You can do determine this by hand, but*lmf*can do it for you. Temporarily insert a line`nkabc= -50`

near the beginning of the EXPRESS category, at any rate before the existing tag

**nkabc**. (Always the first tag is read; later ones are ignored.) Then do`lmf ctrl.nbni -vfile=82 --quit=ham`

and look for this line:`BZMESH: 22 irreducible QP from 40 ( 4 5 2 ) shift= F F F`

If you did this check, be sure to remove the temporary line

`nkabc= -50`

. `--rsham`

was included anticipating that we will use Pashov’s redesign of*lmf*.*lmf*will use alternative set of critical blocks that are significantly faster than the ones they replace, and moreover*lmf*can be run on multicore machines with both OMP and BLAS threading as well as MPI parallelization over**k**. This instruction is put there in anticipation of the Jigsaw Puzzle Orbital basis. Even though it does not work as of this writing, when it does work no change to the input file will be needed. In particular these lines supply information for the JPO basis. At present they supply Hankel energies EL. (The JPO basis does not permit an*l*dependent EH.)`% ifdef rsham>1 STR RMAXA={rscut} ENV[MODE=1 NEL=2 EL={el1} {el2}] % endif`

We use

**el1=−0.3**and**el2=−-1.2**.

We need a starting density and basis set both supplied by *lmfa*:

```
lmfa -vfile=82 ctrl.nbni --tbeh -vrsham=2
```

`-vfile=82`

tells *lmfa* to read site file *site82.nbni*, as can seen by inspecting *ctrl.nbni*. `--tbeh -vrsham=2`

turn on Pashov’s fast *lmf* and also tell *lmf* to use a real-space transformation of the basis set. *lmfa* also reads them but only to ensure that parameters **EH** and **EH2** of the traditional *lmf* basis are assigned values −0.3 and −-1.2. *lmfa* writes the basis file to *basp0.nbni* and information about the atomic density to *atm.nbni*.

We are ready to run *lmf*. Initially we want to make the system self-consistent without relaxing any atoms. Then we can see how large the forces are with our simple but rude estimate for lattice displacements, and to monitor the precision of the **k** point integration.

With 82 atoms, running in serial mode is not a viable option. Exactly how you run *lmf* will depend on the hardware available to you. *lmf* will parallelize over **k** points, but we need to know how many. Extract that information this way:

```
lmf ctrl.nbni -vfile=82 --quit=ham | grep BZMESH
```

There are 22 irreducible QP. If you have a machine with at least 22 cores you can parallelize over all of them. Suppose you have a single node machine with 16 cores, and *lmf* is compiled with Intel’s MKL library (this tutorial was written using such a machine). You could do either one of these:

```
mpix -n 11 lmf -vfile=82 --v8 --tbeh -vrsham=2 ctrl.nbni > out
env OMP_NUM_THREADS=2 MPI_NUM_THREADS=2 | mpix -n 11 lmf -vfile=82 --v8 --tbeh -vrsham=2 ctrl.nbni > out
```

The latter uses 22 threads which is more than 16, but the extra cost of swapping threads is made up for my making better use of the 16 available cores.

*Caution:* When this tutorial was written the `--v8`

option did not compute the forces correctly, so conventional algoriths are used:

```
mpix -n 11 lmf -vfile=82 --v8 --tbeh -vrsham=2 ctrl.nbni > out
```

Run *lmf* again using the (nearly) self-consistent *rst.nbni* file, this time allowing frontier atoms to relax.

```
rm -f mixm.nbni
mpix -n 11 lmf -vfile=82 -vfrzwf=t -vconv=1e-4 -vconvc=1e-3 -vminx=5 --wpos=pos6 ctrl.nbni > out1
```

This uses a Broyden algorithm to relax the structure for five iterations, sufficient reduce the maximum force from an initial 63 mRy/au to about 25 mRy/au.

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

: Perform at most 5 relaxation steps. It gets pretty close to the minimum energy point; a second pass will be taken to refine the positions.`-wpos=pos6`

: Writes the position file after every move, including the last one.

Try `grep RELAX out1`

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.0100 |g|=0.238
RELAX: completed Broyden iter 2; max shift=0.01000 |g|=0.152
RELAX: completed Broyden iter 3; max shift=-0.0100 |g|=0.114
RELAX: completed Broyden iter 4; max shift=-0.0100 |g|=0.0534
RELAX: completed Broyden iter 5; max shift=0.00466 |g|=0.0192
```

**|g|** is the length of the 144 component force “vector”; it should drop from an initial value **0.238** to **0.0192** in the 5^{th} iteration.

The total energy also drops significantly as `egrep -E '^[cx]' out1`

shows:

```
c file=82 frzwf=1 conv=1e-4 convc=.001 minx=5 beta=.2 ehf=-432804.2616707 ehk=-432804.2598271
x file=82 frzwf=1 conv=1e-4 convc=.001 minx=5 beta=.2 ehf=-432804.3137953 ehk=-432804.3135852
x file=82 frzwf=1 conv=1e-4 convc=.001 minx=5 beta=.2 ehf=-432804.3451069 ehk=-432804.3452563
x file=82 frzwf=1 conv=1e-4 convc=.001 minx=5 beta=.2 ehf=-432804.3711576 ehk=-432804.371441
x file=82 frzwf=1 conv=1e-4 convc=.001 minx=5 beta=.2 ehf=-432804.3830901 ehk=-432804.3833868
```

The final gradient is reasonably small, but it needs a little refinement: the change in energy betwen steps 4 and 5 is still significant.

- Concerning the relaxation algorithm
- 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.

- Iterations to convergence
*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`egrep -E 'DQ|Maximum Harris|more=F' out5`

.

This extracts lines with RMS change in charge**DQ**, the maximum value of the Harris force, and summary lines terminating The RMS change in the charge**DQ**drops steadily until the last iteration at given positions (**more=F**). the inner cycle. Atoms are then moved, the potential is disrupted and it the density has to be converged again.`Maximum Harris force = 63.2 mRy/au (site 64) Max eval correction = 0.9 ← iter 1 rlx 0 Maximum Harris force = 63.2 mRy/au (site 64) Max eval correction = 0.8 ← iter 2 rlx 0 Maximum Harris force = 233 mRy/au (site 64) Max eval correction = 355.5 ← iter 1 rlx 1 ... Maximum Harris force = 36.8 mRy/au (site 64) Max eval correction = 111.8 ← iter 20 rlx 1 ...`

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 [iter 20 rlx 1] is the last step in the inner loop (iterations are limited to 20, even if self-consistency has not been reached). The density is nevertheless reasonably converged because the Harris-Foulkes and Kohn-Sham forces are nearly equal. But when the convergence criteria were not met, the summary line e.g.
`x file=82 frzwf=1 conv=1e-4 convc=.001 minx=5 beta=.2 ehf=-432804.3830901 ehk=-432804.3833868`

has an

**x**in the first column, rather than**c**. - 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.

- 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

Allow the system to relax for five more iterations. Note that in the previous run site positions were written to file *pos6*, and the contents of this file are positions of iteration 6 had the cycle continued. Forces for these positions are not calculated, so they are different from the sites positions in the *rst* file which retains positions corresponding to the minimum value of known forces. Nevertheless the contents of *pos6* are likely closer to equilbrium. To take advantage of it, we tell *lmf* to read site positions from *pos6* and include `--rs=11,1,1`

to override the default setting, which would read site positions from *rst.nbni*.

```
rm -f mixm.nbni hssn.nbni
mpix -n 11 lmf -vfile=82 -vfrzwf=t -vconv=1e-4 -vconvc=1e-3 -vminx=5 --rpos=pos6 --wpos=pos --rs=11,1,1 ctrl.nbni > out2
```

Now the forces converge to tolerance: `grep RELAX out6`

```
RELAX: completed Broyden iter 1; max shift=-0.0083 |g|=0.0274
RELAX: completed Broyden iter 2; max shift=-0.0052 |g|=0.0574
RELAX: completed Broyden iter 3; max shift=-0.0028 |g|=0.0181
RELAX: completed Broyden iter 4; max shift=0.00133 |g|=0.0102
RELAX: converged to tolerance; max shift=0.00194 |g|=0.00956
```

and the energy gained by the relaxation is only a few mRy.

Frontier atoms continue to have significant forces, e.g. atom 48 (a frontier Ni atom), as much as 25 mRy/au. If need be we could reduce this force by including more layers in the relaxation. But, this complicates our ability to assemble stacks “lego” style, and since the real interface is probably rather non-ideal and almost certainly different from this one, it isn’t worth the trouble to reduce the residual forces.

### Footnotes and references

^{1} 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).