Self-Consistent ASA calculation for PbTe

This tutorial covers the basics for running a self-consistent LDA calculation in the Atomic Spheres Approximation (ASA), starting with the creation of an input file. There are two main sections:

Building a suitable input file for an ASA-LDA calculation
Running a self-consistent calculation

This page presents an overview of the ASA, and Other Resources provide references for many details on the ASA.

There is also an introductory tutorial for the full-potental program lmf as well as a more detailed tutorial for PbTe. See also a follow-on tutorial to plot the energy bands for PbTe in the ASA.

Create the input file

blm --express=0 --asa --wsitex --findes --nk=4 init.pbte
cp actrl.pbte ctrl.pbte

Checks

lmchk ctrl.pbte

Self-consistency

lmstr ctrl.pbte
lm -vnit=0 ctrl.pbte
lm ctrl.pbte > out.asasc
lm ctrl.pbte >> out.asasc
lmctl ctrl.pbte
cat log.pbte >> ctrl.pbte

1. Building the input file

Under normal atmospheric conditions PbTe crystallises in the rocksalt structure with lattice constant a=6.428 Å. To build an input file, the first step is to construct file init.ext (ext is replaced by a name of your choosing, usually related to the material being studied, pbte in this case). init.pbte will contain the structural information needed for the calculations demonstrated here. For PbTe it will look similar to the following. Copy the content of the box below into init.pbte:

LATTICE
    ALAT=6.427916  UNITS=A
    PLAT=    0.0000000    0.5000000    0.5000000
             0.5000000    0.0000000    0.5000000
             0.5000000    0.5000000    0.0000000
SITE
    ATOM=Pb   X=     0.0000000    0.0000000    0.0000000
    ATOM=Te   X=     0.5000000    0.5000000    0.5000000

File init.pbte shown above has two sections: LATTICE and SITE. The lattice section includes information regarding the lattice structure such as lattice constant (ALAT=) and its units (UNITS=, Å in this case) plus the primitive lattice translation vectors (PLAT=). The site section includes the basis information. In the case of PbTe there are two atoms, one Pb and one Te at the position indicated by “X=”. X= indicates positions in crystal coordinates, which is a compact name for position vectors expressed as linear combinations (i.e. fractional multiples of) of the three primitive lattice vectors.

Alternatively, you can supply the positions with POS=, which specifies positions in Cartesian coordinates, in units of the lattice constant ALAT

Now that init.pbte has been created, blm will create the ctrl file, which is the primary input file for the Questaal suite. Invoke it this way:

blm --express=0 --asa --wsitex --findes --nk=4 init.pbte

This command will generate two new files, namely site.pbte, which contains structural information in a form the Questaal package can read, and the main input file ctrl.pbte, which contains all other information needed to carry out a self-consistent calculation (LDA-ASA calculation in this case). blm has a many command-line switches; to see what they are, type:

blm --h

or refer to the command line switches page. This group used here have the following meaning:

--express controls the brevity of the input file (smaller numbers make more verbose files with more information, so --express=0 is the most verbose). It is worth experimenting with this switch to find which style you are most confortable with. See the discussion in Additional Exercises.
--asa tailors the ctrl file for the ASA. To see how it affects the ctrl file, try running blm without --asa.
--wsitex causes blm to write site positions in crystal coordinates.
--nk=4 specify how many divisions of k points to use for numerical integrations over the Brillouin zone. 4 divisions is a small number; you should check convergence by trying larger meshes; see below.
--findes causes blm to find empty spheres to fill the unit cell, or to place floating orbitals.

--findes is necessary when using ASA in open systems, because the sum of sphere volumes must equal the volume of the unit cell. This is a necessary constraint to conserve charge (very important!) since the ASA has no interstitial; but it means limits the ASA to cases when sites are reasonably well packed. Good packing of open systems can be artificially accomplished by adding “empty” sites with Z=0. --findes will locate them for you automatically.

Ideally the sphere overlaps should be less than 15%. In complex geometries this is not realistic, and experience shows 18% still yields good results. lmchk will check sphere overlaps for you (see Section 2 below). The trick of adding empty spheres is also essential at complex interfaces where misalignment of adjacent planes results in large voids. This tutorial illustrates such a case.

Copy actrl.pbte to ctrl.pbte, which is the name of the main input file.

cp actrl.pbte ctrl.pbte

blm writes to actrl, rather than ctrl, to avoid overwriting a file you may wish to keep.

Notes about the ctrl file

Lines which begin with # are comment lines and are ignored. (More generally, text following # in any line is ignored).
Lines beginning with % directives to the preprocessor. Directives can perform various functions similar to a normal programming language, such as assigning variables, evaluating expressions, conditionally readings some lines, and repeated loops over sections of input.

2. Further setup and checks

The ASA requires that the sum-of-sphere volumes equal the unit cell volume. To check neighbor tables, sphere overlaps and packing fraction, do:

lmchk ctrl.pbte

If you used blm as above to find the empty spheres (ES) you can skip section 2.1 and jump directly to section 3. Section 2.1 is included to explain how to add empty spheres to an already existing ctrl file.

2.1 Adding Empty Spheres to an already-existing ctrl file

Begin by running blm with a slightly different set of switches:

blm --express=0 --asa --wsitex --addes init.pbte
cp actrl.pbte ctrl.pbte

--addes differs from --findes in that in the former case, blm prepares the input file for adding empty spheres, without actually adding them, while in the latter case blm addes them automatically.

To check ES and volume packing invoke:

lmchk ctrl.pbte

The full output can be viewed by clicking here (waiting for file storing system), however the important informations regarding the volume packing and overlap are towards the end of the standard output (stdout) and in this particular case it can be seen (in the case of no ES) in one line

Cell volume= 448.07190   Sum of sphere volumes= 301.06511 (0.67191)

Here the cell volume, sum of all sphere volumes and their ratio are presented, the latter of which has to be equal to 1 for an ASA calculation. Another important value is the overlap percentage, which in this case is given by

  OVMIN, 38 pairs:  fovl = 4.24366e-7   <ovlp> = 8.7%   max ovlp = 8.7%

This line tells us about the average and maximum sphere overlaps. Generally, for ASA the overlaps should be kept below 16% where possible. For full potential calculations generally 5% and below is safe. For GW calculations it should be below 2%. As the sum of sphere volumes is less than the cell volume, we have to add artificial atoms with Z=0 (“empty spheres”) to meet this requirement.

The appropiate space filling spheres can be found using lmchk by invoking:

lmchk --findes --wsitex ctrl.pbte

We have added two command-line switches: the first is to invoke the procedure to find the empty spheres and the second is to write the information into a new site file called essite.pbte. At the end of the stdout generated by lmchk you will see the following snippet:

 ... Final sphere sizes and packing fraction (2 new spheres)

 SCLWSR:  mode = 30  vol = 448.072 a.u.   Initial sphere packing = 81.3%  scaled to 100%
 constr omax1=  16.0  18.0  20.0 %    omax2=  40.0 100.0 100.0 %
 actual omax1=   8.7  12.1   0.0 %    omax2=  16.0  24.6   0.0 %

 spec  name        old rmax    new rmax     ratio
    1   Pb          3.300000    3.300000    1.000000
    2   Te          3.300000    3.300000    1.000000
    3   E           1.959808    2.598601    1.325947

which indicates two new empty spheres have been found, and the new sphere packing is 100%. The control file has to be changed to reflect the new basis. First, change:

  NBAS=2+{les?0:0}  NL=4  NSPEC=2+{les?0:0}

  NBAS=2+{les?2:0}  NL=4  NSPEC=2+{les?1:0}

Quantities in brackets {…} are algebraic expressions expr. In particular, {les?2:0} uses C-like syntax consisting of three expressions separated by ‘?’ and ‘:’. The first (les in this case) is evaluated. If les is nonzero, the result of the entire bracket is the result of the second expression (2 in this case) and otherwise the result of the third (0). In the line above, the contents of tag NBAS= can be interpreted as “if les>0 then NBAS=4 else NBAS=2”. Similarly the contents of NSPEC can be interpreted as “if les>0 NSPEC=3 else NSPEC=2”. “les” is a variable that is defined in the control file through the following line

  % const nit=10 les=1

Here we have also defined nit with value of 10. Finally, we have to pass the information about the empty sphere sites to the control file. We do this by commenting both instances of the line beginning FILE=site and uncommenting the line FILE=essite, as essite.pbte has the new appropiate information.

Note: there are two instances of FILE=site : the first occurs in the STRUC category, where lattice information is read. The second occurs in the SITE category, where basis information is read. Be sure to comment/uncomment both instances.

The last step is to copy the new species information from file poses.pbte file to the SPEC category within the ctrl file, including the new empty spheres.

3. Self-consistency

Before a self consistent calculation can be performed the real-space structure constants have to be generated. They are made once, for a given structure, with a separate utility

lmstr ctrl.pbte

Note that we specified a 4×4×4 k mesh. This should be checked.

Quantities such as the total energy and charge density require an integration over the Brillouin zone. How to perform this integration depends on whether the system is a metal or insulator. (The metal case is more complicated because it must take the discontinuous change in occupation at the Fermi level. The metal case doesn’t concern us here because PbTe is an insulator.)

Whether metal or insulator, standard practice is to take a uniform mesh of points in the first Brillouin zone. Thus the mesh is specified by the number of divisions along each the reciprocal lattice vectors; call them n₁, n₂, n₃. ctrl.pbte has this line:

BZ    METAL={met}  NKABC={nkabc}

met and nkabc are preprocessor variables that are translated into numbers by the preprocessor before being parsed by the input file reader. Default values are set at the top of ctrl.pbte in this line:

% const met=(asa?1:5) nkabc=4

The mesh is specified through NKABC, which is a tag followed by three integers n₁, n₂, n₃. If you supply only one number (as is the case here), the code will set n₂ and n₃ to n₁.

The tag actually reads NKABC={nkabc}. The curly brackets around nkabc tell the file preprocessor to parse the contents of {\hellip;} as an algebraic expression; in this case the result is just the value of the preprocessor variable nkabc. If variables nkabc and met have been assigned to 4 and 0 respectively, file preprocessor will transform this line into

BZ    METAL=0  NKABC=0

The input parser detects only a single expression, which becomes n₁. n₂ and n₃ are set to n₁, as noted before.

Note: there is a additional freedom in choosing the k mesh (not used in this tutorial). Aecond triplet of parameters can specify whether the mesh passes through the origin along a particular G vector, or is staggered. See this page for more details.

Before performaing a self-consistent calculation, we must guess a starting potential. In the ASA, information needed to determine the potential is contained in a compact form. It is completely determined by multipole moments $Q_{0;\ell}$ , $Q_{1;\ell}$ , $Q_{2;\ell}$ for $\ell=0{\ldots}3$ . (Inspect the input file. Note that blm selected LMX=3 for Pb and Te. This is because both are rather large atoms: the larger the radius of the augmentation sphere, the more $\ell$ ’s are required solve the Schrodinger equation to a given precision.)

The are several ways to specify initial estimates $Q_{0{\ldots}2;\ell}$ . The simplest way is to let lm choose defaults. In the absence of another specification, lm reads preset values for the $Q_{0;\ell}$ and $P_\ell$ from a lookup table and sets $Q_{1;\ell}{=}Q_{2;\ell}{=}0$ . This will generate a crude initial estimate for the density, but it is usually sufficient for a starting guess, which can be iterated to self-consistency. (Sometimes the estimate is too crude and alternatives must be chosen.)

Invoke lm executable with zero number of iterations as:

lm -vnit=0 ctrl.pbte

This command starts from $Q_{0;\ell}$ , $Q_{1;\ell}$ , $Q_{2;\ell}$ (default values are taken since none are specified explicitly) and makes a trial potential from it.

Note: -vnit=0 does not directly set the number of iterations. It assigns a value 0 to preprocessor variable, translates this line of the input

      NIT=    {abs(nit)}

into

      NIT=    0

It also translates this line

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

into

START CNTROL=1 BEGMOM=1

BEGMOM governs program flow in lm; specifically BEGMOM=1 tells lm to start with moments $Q_{0{\ldots}2;\ell}$ to make potential parameters which define the ASA one-particle hamiltonian.

lm has two main branches:

a “sphere” branch that loops over each site independently, takes as input multipole moments $Q_{0;\ell}$ , $Q_{1\ell}$ , $Q_{2;\ell}$ with a specification of the boundary condition through the “continuous principal quantum number” $P_\ell$ , and generates a unique potential $v(r)$ consistent with those moments. One important output derived from the partial waves are the “potential parameters”, which define the ASA hamiltonian. The most important of these are the band “center of gravity” C and “bandwidth” Δ.
a “crystal” branch that generates hamiltonians and constructs eigenfunctions and derived quantities (in particular the $Q_{0;\ell}$ , $Q_{1\ell}$ , $Q_{2;\ell}$ ) from potential parameters.

The setup step (lm -vnit=0) executes just the “sphere” branch because the number of iterations is zero, and BEGMOM=1.

For a given $v(r)$ , partial waves can be numerically integrated for any energy. The LMTO method is a linear one, so the full energy dependence of the partial wave is approximated by the partial wave and its energy derivative at some linearization energy $\varepsilon_\nu$ . Alternatively, the $\varepsilon_\nu$ can be specified implicitly though the boundary condition at the augmentation radius. Usually the ASA codes specify $\varepsilon_\nu$ implicitly from continuous principal quantum number $P_\ell$ , which fixes the logarithmic derivative $D_l$ , which fixes $\varepsilon_\nu$ . In sum, complete information about the sphere is given from the $P_\ell$ and the $Q_{0{\ldots}2;\ell}$ .

At the end of this branch, other quantities may by calculated depending on what is sought. Examples are matrix elements of the gradient operator for optics calculations, matrix elements of the spin-orbit hamiltonian, and matrix elements of an external field.

The ASA Hamiltonian depends only on potential parameters and structure constants S. The latter depend only on nuclear positions; they are calculated beforehand for a given structure with the utility lmstr. Indeed the simplest form of the ASA hamiltonian is $H_{\mathbf{R}L,\mathbf{R}'L'} = C_{\mathbf{R}L}\delta_{\mathbf{R}L,\mathbf{R}'L'} + \sqrt{\Delta}_{\mathbf{R}L} S_{\mathbf{R}L,\mathbf{R}'L'} \sqrt{\Delta}_{\mathbf{R}L,\mathbf{R}'L'}$

$\mathbf{R}$ refers to site, and $L$ is a compound index of $\ell$ and $m$ quantum numbers corresponding to an envelope function centered at site $\mathbf{R}$ . In the ASA there is one basis function per partial wave, so there is a one-to-one correspondence between the envelope at $\mathbf{R}$ and the partial wave $\phi_l$ at the same site. The rank of $H$ is the number of $\mathbf{R}L$ channels.

For a known set of potential parameters, lm can run the crystal branch to generate wave functions, from which it can make quantities of physical interest, e.g. the dielectric function. Unless it is run in a special mode (e.g. to generate energy bands for plotting), the crystal branch performs an integration over the Brillouin zone to generate multipole moments $Q_{0{\ldots}2;\ell}$ .

Self-consistency proceeds by alternating between the crystal branch which takes potential parameters as input and generates $Q_{0{\ldots}2;\ell}$ as output; and the sphere branch, which generates the potential parameters from a given $Q_{0{\ldots}2;\ell}$ and $P_\ell$ . $PPAR\stackrel{crystal}{\longrightarrow} (P_\ell, Q_{0{\ldots}2;\ell})\stackrel{sphere}{\longrightarrow} PPAR \stackrel{crystal}{\longrightarrow} (P_\ell, Q_{0{\ldots}2;\ell})\stackrel{sphere}{\longrightarrow}$

Between the crystal and sphere branches, some other steps are taken. One is to float (adjust) the $P_\ell$ to the band center-of-gravity: by changing the boundary conditions the linearization energy changes, and thus $Q_{1{\ldots}2;\ell}$ change. The band center is defined as the $P_\ell$ for which $Q_{1;\ell}{=}0$ . This is generally desirable; however, for partial waves far removed from the Fermi level, floating $P_\ell$ to the point $Q_{1;\ell}{=}0$ can be dangerous, and cause “ghost bands”. For such partial waves it is necessary to limit the shifts in $P_\ell$ . Token SPEC_ATOM_IDMOD will freeze $P_\ell$ for any l you specify. You can also set HAM_PMIN to −1: this will prevent and $P_\ell$ from falling below its free electron value.

Another is to mix input and output moments for an estimate of the self-consistent set. How this mixing is carried out is governed by parameters in the ITER_MIX tag (or EXPRESS_mix if present).

Note that lm can start directly with the crystal branch if potential parameters are given. Whether it starts in the crystal branch or the sphere branch is governed by the token START_BEGMOM.

For a self consistent LDA-ASA calculation invoke lm

lm ctrl.pbte > out.asasc

If we have reached self-consistency, the RMS change in (output-input) moments $Q_{0{\ldots}2;\ell}$ (DQ) would be very small. We can check its progress with iteration with

grep DQ out.asasc

You should see the following

 PQMIX:  read 0 iter from file mixm.  RMS DQ=1.67e-1
 PQMIX:  read 1 iter from file mixm.  RMS DQ=1.17e-1  last it=1.67e-1
 PQMIX:  read 2 iter from file mixm.  RMS DQ=8.37e-2  last it=1.17e-1
 PQMIX:  read 3 iter from file mixm.  RMS DQ=5.13e-3  last it=8.37e-2
 PQMIX:  read 4 iter from file mixm.  RMS DQ=2.74e-3  last it=5.13e-3
 PQMIX:  read 5 iter from file mixm.  RMS DQ=1.64e-3  last it=2.74e-3
 PQMIX:  read 6 iter from file mixm.  RMS DQ=1.10e-3  last it=1.64e-3
 PQMIX:  read 0 iter from file mixm.  RMS DQ=5.19e-4  last it=1.10e-3
 PQMIX:  read 1 iter from file mixm.  RMS DQ=3.45e-4  last it=5.19e-4
 PQMIX:  read 2 iter from file mixm.  RMS DQ=2.35e-4  last it=3.45e-4

Indeed the (rms) DQ has become small, but it is on the cusp of whether it fully self-consistent or not. Run lm again:

lm ctrl.pbte >> out.asasc

The end of the standard output should look like this:

 SV:   2  1.853e-5 1.0000  2.745e-5  -55318.14249975 0.000000 B 81.4
 it 2 of 10        ehf=  -55318.1424998   ehk=  -55318.1424997
 From last iter    ehf=  -55318.1425057   ehk=  -55318.1425014
 diffe(q)=  0.0000059 (0.0000185)    tol= 0.0000100 (0.0000300)   more=F
c ehf=-55318.1424998 ehk=-55318.1424997
 cpudel    ...   Time this iter:  time(s):  0.0307   total:  0.104s

 Jolly good show !  You converged to rms DQ=  0.000027

In two additional iterations we obtained good convergence. “Jolly good show” indicates that self-consistency has been achieved (within the tolerances specified in the ctrl file).

Multipole moments, potential, potential parameters, and derived quantities such as matrix elments of the spin orbit coupling or the momentum operator are stored in atom files. There is one atom file per class; in this case there are files pb.pbte, te.pbte, e.pbte. You can inspect these files.

As a final step, you can collect the self-consistent moments generated by lm and add them to the ctrl file. Type

lmctl ctrl.pbte

lmctl reads the $Q_{1{\ldots}2;\ell}$ from the atom files, clears the log file log.pbte, and overwrites it with parameters $P\ell$ _ and $Q_{1{\ldots}2;\ell}$ , in a form suitable for the ctrl file. Append log.pbte to it with, e.g.

cat log.pbte >> ctrl.pbte

Now the ctrl file retains essentially complete information for to generate the self-consistent density.

As a next step, you might want to draw the energy bands.

Additional Exercises

With ctrl.pbte generated by blm, lm downfolds the Pb f and Te f states: they are removed from the basis in a particular way. Explore what happens if you turn off the downfolding (set 2 in tag IDXDN to 1). This increases the rank of the Hamiltonian. You may also need to freeze or constrain the $P_\ell$ (use tag ATOM_SPEC_IDMOD or possibly HAM_PMIN) to prevent ghost bands from appearing. Also compare how the bands change if the f states are removed all together (set LMX=2).

Try running blm without the --express=0 switch. actrl.pbte will appear somewhat different, but the two input files should generate the practically the same results. They should be identical except in this latter case the calculation run for one additional iteration, because an extra convergence criterion is included in the ctrl file, namely

  conv=   1e-5                     # Convergence tolerance (energy)