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

- Create the input file

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

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

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

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

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 *n*_{1}, *n*_{2}, *n*_{3}. *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 *n*_{1}, *n*_{2}, *n*_{3}. If you supply only one number (as is the case here), the code will set *n*_{2} and *n*_{3} to *n*_{1}.

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*_{1}. *n*_{2} and *n*_{3} are set to *n*_{1}, 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 , , for . (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 . The simplest way is to let **lm** choose defaults. In the absence of another specification, **lm** reads preset values for the and from a lookup table and sets . 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 , , (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 , , with a specification of the boundary condition through the “continuous principal quantum number” , and generates a unique potential consistent with those moments. For a given , 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 . Alternatively, the can be specified implicitly though the boundary condition at the augmentation radius. Usually the ASA codes specify implicitly from continuous principal quantum number , which fixes the logarithmic derivative , which fixes . In sum, complete information about the sphere is given from the and the . 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

refers to site, and is a compound index of and quantum numbers corresponding to an envelope function centered at site . In the ASA there is one basis function per partial wave, so there is a one-to-one correspondence between the envelope at and the partial wave at the same site. The rank of is the number of 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 .

Self-consistency proceeds by alternating between the crystal branch which takes potential parameters as input and generates as output; and the sphere branch, which generates the potential parameters from a given and .

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

### 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 increaases the rank of the Hmailtonian. You may also need to freeze or constrain the (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)
```