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

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

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

1. Create the input file
blm --express=0 --asa --wsitex --findes init.pbte
cp actrl.pbte ctrl.pbte

1. Checks
lmchk ctrl.pbte


The following alternate route to build empty spheres from section 2.1 is optional

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

1. Self-consistency
sed -i s/nkabc=0/nkabc=4/ ctrl.pbte
lmstr ctrl.pbte
lm -vnit=0 ctrl.pbte
lm -vnit=20 ctrl.pbte
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 of the init file 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 position with POS=, which specifies positions in Cartesian coordinates, in units of the lattice constant.

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  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). The switches blm knows about are described on the command line switches page. This particular group has the following meaning:

The --express switch
--express=0 controls the brevity of the input file (smaller numbers make more verbose files with more information). It is worth experimenting with this switch to find which style of control file you are most confortable with. See the discussion in Additional Exercises.
The --asa switch
tailors the ctrl file for the ASA. To see how it affects the ctrl file, try running blm without --asa.
The --wsitex switch
causes blm to write site positions in crystal coordinate.
The --findes switch
causes blm to find empty spheres to fill the unit cell, or to place floating orbitals.

The last switch is necessary when using ASA in open systems, because the sum of sphere volumes must equal the cell volume. The ASA only works well when sites are closely packed. Close packing of open systems can be artificially accomplished by adding “empty” sites with Z=0.

blm has a number of other command-line switches; to see what they are, type:

blm --h


or refer to the command line switches page.

Now the only thing left to do is to rename 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.

Note

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

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}


to

  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 tool

lmstr ctrl.pbte


As yet the ctrl file is incomplete, because we have not specified the k mesh. Change the nkabc variable within the control file to nkabc=4. This variable controls the k-mesh density. Use your text editor to change:

% const nkabc=0


to

% const nkabc=4


or simply make the substitution with the sed editor:

sed -i s/nkabc=0/nkabc=4/ ctrl.pbte


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 n1, n2, n3. ctrl.pbte has this line:

BZ    METAL={met}  NKABC={nkabc}


Note: if you have run blm without --express=0, the information appears in the EXPRESS category as

  nkabc=  {nkabc}                  # 1 to 3 values


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

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 n1. n2 and n3 are set to n1, 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.

The penultimate step is to generate 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_{0l}$, $Q_{1l}$, $Q_{2l}$ for $l=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 l’s are required solve the Schrodinger equation to a given precision.)

The are several ways to specify initial estimates $Q_{0{\ldots}2l}$. The simplest way is to let lm choose defaults. In the absence of another specification, lm reads preset values for the $Q_{0l}$ and $P_l$ from a lookup table and sets $Q_{1l}{=}Q_{2l}{=}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 the self-consistent values.

Invoke lm executable with zero number of iterations as:

lm -vnit=0 ctrl.pbte


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

lm has two main branches: an “sphere” part that loops over each site independently, and a “crystal” part that generates hamiltonians and constructs eigenfunctions and derived quantities. The setup step (lm -vnit=0) executed just the “sphere” branch. It takes as input multipole moments $Q_{0l}$, $Q_{1l}$, $Q_{2l}$ with a specification of the boundary condition through the “continuous principal quantum number” $P_l$, and generates a unique potential $v(r)$ consistent with those moments. 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_\mu$. Alternatively, the $\varepsilon_\mu$ can be specified implicitly though the boundary condition at the augmentation radius. Usually the ASA codes specify $\varepsilon_\mu$ implicitly from continuous principal quantum number $P_l$, which fixes the logarithmic derivative $D_l$, which fixes $\varepsilon_\mu$. In sum, complete information about the sphere is given from the $P_l$ and the $Q_{0{\ldots}2l}$. One important output derived from the partial waves are the “potential parameters”, the most important of which are the band “center of gravity” C and “bandwidth” Δ.

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

$\mathbf{R}$ refers to site, and $L$ is a compound index of $l$ 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}2l}$.

Self-consistency proceeds by alternating between the crystal branch which takes potential parameters as input and generates $Q_{0{\ldots}2l}$ as output; and the sphere branch, which generates the potential parameters from a given $Q_{0{\ldots}2l}$ and $P_l$.

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

lm -vnit=20 ctrl.pbte


You should see “Jolly good show” at the end of the standard output. This indicates that self-consistency has been achieved.

Note: -vnit=20 does not directly set the number of iterations. This number is determined by token NIT in this line:

ITER  MIX=B2,b=.3,k=7  NIT={nit}  CONVC=1e-5


The preprocessor sees the curly brackets ({nit}), treats the contents as an expression (a trivial one in this case) and substitutes the result (20) for {nit}.

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}2l}$ from the atom files, clears the log file log.pbte, and overwrites it with linearization parameters P and the $Q_{1{\ldots}2l}$ , in a form suitable for the ctrl file. Append log.pbte to it with, e.g.

cat log.pbte >> ctrl.pbte


In this way the ctrl file retains the essential information for the self-consistent density.

1. 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 increaases the rank of the Hmailtonian. You may also need to freeze or constrain the $P_l$ (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).
  conv=   1e-5                     # Convergence tolerance (energy)