Questaal Home

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

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, and spec.nb and, and also sitea.nb, siteb.nb and, and

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

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

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

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
echo STRUC >>
echo ' DEFGRD= ' `mcx -f9f12.8 pni -i pnb -x -t -a:nr=1 s s -w:nohead . | sed 's:[0-9.]*$:1:'` >>
lmchk -vfile=2 ni --wsitex~fn=site3
lmchk -vfile=2 ni | grep vol; lmchk -vfile=3 ni | grep vol
sed -i.bak 's:\(plat.*[1-9.]*\)$:\1/1.00983:'
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:'
lmchk --wspec >
lmscell -vfile=42 '--stack~rmsite@targ=15:42~pad@0,0,-2*.61~wsite@fn=sitea'
lmscell -vfile=42 '--stack~rmsite@targ=1:14,29:42~addpos@dx=0,0,-.61~pad@0,0,-2*.61~wsite@fn=siteb'
lmscell -vfile=42 '--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

Nb/Ni(3n)/Nb superlattices

This section requires the files ctrl.nb, site20.nb, spec.nb,,

Set initial conditions with

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

Confirm Nb/Ni and Ni/Nb have equal spacing

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


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.


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 editorlmscell --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 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
    ALAT=6.226 PLAT= -.5 .5 .5  .5 -.5 .5  .5 .5 -.5
    ATOM=Nb   X=0 0 0


# Init file for Ni.  Lattice constant at 0K
    ALAT=6.643 PLAT= .5 .5 .0  .5 .0 .5  .0 .5 .5
    ATOM=Ni  MMOM=0,0,0.6  PZ=0,0,0
    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}

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 siten.ext, when preprocessor variable file is assigned to n. These lines have no effect if file is not assigned, in which case Questaal codes will read from the usual file site.ext.

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

 DEFGRD=  0.99698194  0.00000000  0.00000000  0.04972026  1.02768205  0.00000000  0.00000000  0.00000000  1

and create a new file 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 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

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 to and scale the 9th element of plat by 1/1.02458, or simply use sed:

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

Run lmchk again and confirm the unit cell volumes of and 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 :; 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 9th element of plat with dval 1.73205081/1.00983 = 1.71519. To do it on the command-line, invoke

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

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

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:

  1. the Nb/(A’B’C’)n/Nb family with 3n Ni planes
  2. the Nb/(A’B’C’)nA’/Nb family with 3n+1 Ni planes
  3. the Nb/(A’B’C’)nA’B’/Nb family with 3n+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(3n)/Nb superlattices

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

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

We will distinguish the three kinds of lattices by denoting them as 3n, 3n+1 and 3n+2 superlattices, respectively. This section builds the simplest superlattice in the 3n 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',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:

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

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

    3. ~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 --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
  •   supplies blm with species-specific information gleaned from spec.nb and, notably in this case sphere radii.
  • --usefiler   tells blm to accept without change the muffin-tin radius (R) it read from spec.nb and 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 3n 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,2k 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.

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

  2. --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, the rlx column). blm sets rlx to zero for all sites not on a frontier layer.

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

  4. --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 5th 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  &larr; iter 1 rlx 0
Maximum Harris force = 63.2 mRy/au (site 64)  Max eval correction = 0.8  &larr; iter 2 rlx 0

Maximum Harris force = 233 mRy/au (site 64)  Max eval correction = 355.5  &larr; iter 1 rlx 1
Maximum Harris force = 36.8 mRy/au (site 64)  Max eval correction = 111.8  &larr; 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.

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