# Detailed lmf tutorial: self-consistent LDA calculation for PbTe

This tutorial carries out a self-consistent density-functional calculation for PbTe using the **lmf** code. It has a purpose similar to the basic tutorial on Si but provides much more detail. See also the Fe tutorial, an LDA+QSGW calculation for a ferromagnetic metal.

It synchronizes with an ASA tutorial on the same system, enabling a comparison of the ASA and full potential methods, and forms a starting point for other tutorials, e.g. on optics.

Make an input and site files:

```
nano init.pbte
blm --nk=6 --gmax=8.2 --ctrl=ctrl init.pbte
```

Free atomic density and basis parameters

```
lmfa ctrl.pbte #use lmfa to make basp file, atm file and to get gmax
cp basp0.pbte basp.pbte #copy basp0 to recognised basp prefix
lmfa ctrl.pbte #remake atomic density with updated valence-core partitioning
```

Alternatively perform the same function with one step

```
lmfa ctrl.pbte --usebasp #use lmfa to make basp file, atm file and to get gmax
```

```
lmf ctrl.pbte -vnkabc=6 -vgmax=8.2
```

### Table of Contents

- 1. Building the input file
- 2. How the ctrl file is organized
- 3. The
**EXPRESS**category - 4. Determining what input an executable seeks
- 5. Initial setup: free atomic density and parameters for basis
- 6. Remaining Inputs
- 7. Sanity Checks
- 8. Self consistency
- Other Resources
- FAQ
- Additional exercises

### Preliminaries

Some of the basics are covered in the basic lmf tutorial for Si. It is easier to read but less detailed. See also the the tutorial on building input files.

The standard outputs from running this tutorial are explained in the annotation of *lmfa* output and annotation of *lmf* output. Many details omitted here are given in the annotated outputs.

Executables *blm*, *lmchk*, *lmfa*, and **lmf** are required and are assumed to be in your path.

### 1. Building the input file

(See also the tutorial on building input files).

PbTe crystallizes in the rocksalt structure with lattice constant *a* = 6.428 Å. You need the structural information in the box below to construct the main input file, *ctrl.pbte*. Start in a fresh working directory and cut and paste the box’s contents to *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
```

The primitive lattice vectors are in row format (the first row contains the *x*, *y* and *z* components of the first lattice vector and so forth). In the **SITE** section, the a tom type and coordinates are shown. **X=** specifies the site coordinates. They are specified as fractional multiples of lattice vectors **PLAT** (sometimes called “direct” representation). You can also use Cartesian coordinates; instead of **X=** you would use **POS=** (see additional exercises below). Positions in Cartesian coordinates are in units of **ALAT**, like the lattice vectors.

*Note :* You can also import structural data from other sources. See this page for options.

Use the *blm* tool as in the box below to create the input file (*ctrl.pbte*) and the site file (*site.pbte*):

```
blm init.pbte
cp actrl.pbte ctrl.pbte
```

Or just simply

```
blm init.pbte --ctrl=ctrl
```

*Note :* If you are preparing for a later QS*GW* calculation, use `blm --gw init.pbte`

. For quick reference to the command-line switches *blm* recognizes, try `blm --h`

. This page offers more complete documentation.

### 2. How the ctrl file is organized

In this tutorial, *blm* use “standard” mode to create the input file *ctrl.pbte*. (The basic tutorial creates a simpler form using `blm --express init.si`

). Standard mode makes limited use of the preprocessing capabilities of the Questaal input system : it uses algebraic variables which can be modified on the command line. Thus `lmf -vnit=10 ...`

sets variable **nit** to 10 before doing anything else. Generally lines in the input file:

- which begin with ‘
**#**’ are comment lines and are ignored. (More generally, text following a `**#**’ in any line is ignored). - beginning with ‘
**%**’ are directives to the preprocessor.

The beginning of the ctrl file generated by *blm* should look like the following:

```
# Variables entering into expressions parsed by input
% const nit=10
% const met=5
% const so=0 nsp=so?2:1
% const lxcf=2 lxcf1=0 lxcf2=0 # for PBE use: lxcf=0 lxcf1=101 lxcf2=130
% const pwmode=0 pwemax=3 # Use pwmode=1 or 11 to add APWs
% const nkabc=0 gmax=0
```

**% const** tells the preprocessor that it is declaring one or more variables. **nit**, **met**, etc, used in expressions later on. The parser interprets the contents of brackets **{…}** as algebraic expressions: The contents of **{…}** is evaluated and the numerical result is substituted for it. Expression substitution works for input lines proper, and also in the directives.

For example this line

```
metal= {met} # Management of k-point integration weights in metals
```

becomes

```
metal= 5
```

because **met** is a numerical expression (admittedly a trivial one). It evaluates to 5 because **met** is declared as an algebraic variable and assigned value 5 near the top of the ctrl file. The advantage is that you can do algebra in the input file, and you can also re-assign values to variables from the command line, as we will see shortly.

Lines corresponding to actual input are divided into **categories** and *tokens* within the categories.

A category begins when a character (other than **%** or **#**) occurs in the first column. Each token belongs to a category; for example in box below **IO** contains three tokens, *SHOWMEM*, *IACTIV* and *VERBOS* :

```
IO SHOWMEM=f
IACTIV=f VERBOS=35,35
```

(Internally, a complete identifier (aka *tag*) would be **IO_IACTIV=**, though it does not appears in that form in the ctrl file.)

This link explains the structure of the input file in more detail.

### 3. The **EXPRESS** category

*blm* normally includes an **EXPRESS** category in *ctrl.pbte*.

```
EXPRESS
# Lattice vectors and site positions
file= site
# Basis set
gmax= {gmax} # PW cutoff for charge density
autobas[pnu=1 loc=1 lmto=5 mto=4 gw=0]
```

Tags in the **EXPRESS** category are effectively aliases for tags in other categories, e.g. **EXPRESS_gmax** corresponds to the same input as **HAM_GMAX**. If you put a tag into **EXPRESS**, it will be read there and ignored in its usual location; thus in this instance adding **GMAX** to the **HAM** category would have no effect.

The purpose of **EXPRESS** is to simplify the input file, collecting the most commonly used tags in one place.

### 4. Determining what input an executable seeks

Executables accept input from two primary streams : tags in the ctrl file and additional information through command-line switches. Each executable reads its own particular set, though most executables share many tags in common.

Usually an input file contains only a small subset of the tags an executable will try to read; defaults are used for the vast majority of tags.

There are four special modes designed to facilitate managing input files. For definiteness consider the executable *lmfa*.

```
$ lmfa --input
$ lmfa --help
$ lmfa --showp
$ lmfa --show | lmfa --show=2
```

`--input`

puts *lmfa* in a special mode. It doesn’t attempt to read anything; instead, it prints out a (large) table of all the tags it would try to read, including a brief description of the tag, and then exits.

See here for further description.

`--help`

performs a similar function for the command line arguments: it prints out a brief summary of arguments effective in the executable you are using.

See annotated lmfa output for further description.

`--showp`

reads the input through the preprocessor, prints out the preprocessed file, and exits.

See the annotated lmf output for a comparison of the pre- and post-processed forms of the input file in this tutorial.

`--show`

tells *lmfa* to print out tags as it reads them (or the defaults it uses).

It is explained in the annotated lmf output.

### 5. Initial setup: free atomic density and parameters for basis

- The basp file
- In order for
*lmf*carry out a self-consistent calculation for the crystal, it requires information about the basis set, and an initial trial density.*lmfa*computes densities for free atoms, which are overlapped to form an initial trial density (called the Mattheis construction).*lmfa*will also automatically generate a basis set, and save it in a file called*basp.ext*. What is stored in this file depends on the contents of**HAM_AUTOBAS**, but typically three pieces of information are kept:- Parameters
**PZ**specifying local orbitals - Parameters
**(RSMH,EH)**defining the shape of the smooth Hankel envelope functions - Parameters
**P**defining the boundary conditions of partial waves on augmentation boundaries

- Parameters

As will be described in more detail in the pages below, *lmfa* prepares the following.

Make a self-consistent atomic density for each species.

Fit the density outside the augmentation radius.

*lmf*needs this information to overlap atomic densities for an initial trial density.

Information about the augmented and interstitial parts of the density are written to file*atm.pbte*.Provide a reasonable estimate for the Gaussian smoothing radius

*r*and hankel energy_{s}*ε*) that fix the shape of the smooth Hankel envelope functions for*l*=0, 1,…. The*l*cutoff is determined internally, depending on the setting of**HAM_AUTOBAS_LMTO**.

These parameters are written to file*basp0.pbte*as**RSMH**and**EH**.Provide a reasonable estimate for boundary conditions that fix linearization energies, parameterized by the logarithmic derivative parameter

*P*, aka the “continuous principal quantum number.”_{l}

These parameters are written to*basp0.pbte*as**P**.Decide on which shallow cores should be included as local orbitals.

Local orbitals are written*basp0.pbte*as nonzero values of**PZ**.Supply an estimate for the interstitial density plane wave cutoff

**GMAX**.

*lmfa* will provide all of this information automatically. It will write atomic density information to *atm.pbte* and basis set information to template *basp0.pbte*. The Questaal suite reads from *basp.pbte*, but *lmfa* writes to *basp0* to avoid overwriting a file you may want to preserve. You can edit *basp.pbte* and customize the basis set.

As a first step, do:

```
lmfa ctrl.pbte #use lmfa to make basp file, atm file and to get gmax
cp basp0.pbte basp.pbte #copy basp0 to recognised basp prefix
```

The output is annotated in some detail here. It begins with a header:

```
LMFA: nbas = 2 nspec = 2 verb 35
pot: XC:BH
autogen: mto basis(4), pz(1), pnu(1) Autoread: pz(1)
```

The **pot** line says that *lmfa* the potential will be made from the Barth-Hedin functional. To use a GGA, see here.

The **autogen** line says that *lmfa* will make the basis set information (points 3-5 outlined above).

The next few sections amplify on these three points. Point 6 is discussed here.

#### 5.1 Local orbitals

Part of *lmfa*’s function is to identify *local orbitals* that extend the linear method. Linear methods are reliable only over a limited energy window; certain elements may require an extension to the linear approximation for accurate calculations. This is accomplished with local orbitals. *lmfa* will automatically look for atomic levels which, if certain criteria are satisfied it designates as a local orbital, and includes this information in the basp0 file. The annotated lmfa output explains how *lmfa* analyzes core states for local orbitals.

Inspect *basp.pbte*. Note in particular this text connected with the Pb atom:

```
PZ= 0 0 15.9336
```

(The same information can be supplied in the input file, through **SPEC_ATOM_PZ**.)

*lmfa* is suggesting that the Pb 5*d* state is shallow enough that it be included in the valence. Since this state is far removed from the Fermi level, we would badly cover the Hilbert space spanned by Pb 6*d* state were we to use Pb 5*d* as the valence partial wave. (In a linear method you are allowed to choose a single energy to construct the partial wave; it is usually the “valence” state, which is near the Fermi level.)

This problem is resolved with local orbitals : these are partials wave at an energy far removed from the Fermi level. The three numbers following **PZ** correspond to specifications for local orbitals in the *s*, *p*, and *d* channels. Zero indicates “no local orbital;” there is only a *d* orbital in this case.

**15.9336** is actually a compound of **10** and the “continuous principal quantum number” **5.9336**. The 10’s digit tells *lmf* to use an “enhanced” local orbital as opposed to the usual variety found in most density-functional codes. Enhanced orbitals append a tail so that the density from the orbital spills into the interstitial. You can specify a “traditional” local orbital by omitting the 10, but the extended kind is more accurate, and there is no advantage to doing so.

The continuous principal quantum number (**5.9336**) specifies the number of nodes and boundary condition. The large fractional part of *P* is large for core states, typically around 0.93 for shallow cores. *lmfa* determines the proper value for the atomic potential. In the self-consistency cycle the potential will change and *lmf* will update this value.

*lmfa* automatically selects the valence-core partitioning; the information is given in *basp.pbte*. You can set the partitioning manually by editing this file.

*Note:* high-lying states can also be included as local orbitals; they improve on the Hilbert space far above the Fermi level. In the LDA they are rarely needed and *lmfa* will not add them to *basp.pbte*. But they can sometimes be important in *GW* calculations, since in contrast to the LDA, unoccupied states also contribute to the potential.

After *basp.pbte* has been modified, you must run *lmfa* a second time:

```
lmfa ctrl.pbte
```

This is necessary whenever the valence-core partitioning changes through the addition or removal of a local orbital. Even though *lmfa* writes the atomic to *atm.pbte*, this file will change when partitioning between core and valence will change with the introduction of local orbitals, as described next. This is because core and valence densities are kept separately.

Alternatively, you can cause *lmfa* to remake the valence-core partitioning after the LO have been found, and write directly to file *basp.pbte*.

```
lmfa ctrl.pbte --usebasp
```

It saves effort since you don’t need to copy *basp0.pbte* to *basp.pbte* or run *lmfa* twice. On the other hand it gives you less flexibility since you no longer control which orbitals are to be treated as LO.

To see how to control and monitor lmfa’s search for local orbitals, see the annotated lmfa output.

##### Relativistic core levels

Normally *lmfa* determines the core levels and core density from the scalar Dirac equation. However there is an option to compute the core levels from the full Dirac equation.

Tag **HAM_REL** controls how the Questaal package manages different levels of relativistic treatment. Run `lmfa --input`

and look for **HAM_REL**. You should see:

```
HAM_REL opt i4 1, 1 default = 1
0 for nonrelativistic Schrödinger equation
1 for scalar relativistic Schrödinger equation
2 for Dirac equation (ASA only for now)
10s digit 1: compute core density with full Dirac equation
10s digit 2: Like 1, but neglect coupling (1,2) pairs in 4-vector
```

Set **HAM_REL=11** to make *lmfa* calculate the core levels and core density with the full Dirac equation.

You might want to see the core level eigenvalues; they can shift significantly relative to the scalar Dirac solution. Also, *l* is no longer a good quantum number so there can be multiple eigenvalues connected with the scalar Dirac *l*. To see these levels, invoke *lmfa* with a sufficiently high verbosity. In the present instance insert **HAM REL=11** into *ctrl.pbte* and do

```
lmfa --pr41 ctrl.pbte
```

You should see the following table:

```
Dirac core levels:
nl chg <ecore(S)> <ecore(D)> <Tcore(S)> <Tcore(D)> nre
1s 2 -6461.412521 -6461.420614 9160.575645 9160.568216 439
ec(mu) -6461.420614 -6461.420614
2s 2 -1154.772794 -1154.777392 2201.484620 2201.485036 473
ec(mu) -1154.777392 -1154.777392
3s 2 -277.137428 -277.136313 700.148783 700.160432 501
ec(mu) -277.136313 -277.136313
4s 2 -62.683976 -62.678557 231.671152 231.686270 531
ec(mu) -62.678557 -62.678557
5s 2 -10.589828 -10.580503 60.826909 60.833608 567
ec(mu) -10.580503 -10.580503
2p 6 -990.094400 -1001.984462 1702.510726 1772.365432 475
ec(mu) -948.389636 -1109.174115 -948.389636 -1109.174115 -948.389636 -948.389636
3p 6 -229.993746 -232.623198 568.649082 585.156080 505
ec(mu) -220.667558 -256.534478 -220.667558 -256.534478 -220.667558 -220.667558
4p 6 -47.246014 -47.902771 184.751871 189.523363 537
ec(mu) -44.969950 -53.768412 -44.969950 -53.768412 -44.969950 -44.969950
5p 6 -6.300710 -6.422904 43.507054 44.670581 577
ec(mu) -5.869706 -7.529300 -5.869706 -7.529300 -5.869706 -5.869706
3d 10 -182.032939 -182.146340 501.452676 502.171493 509
ec(mu) -179.091564 -186.728504 -179.091564 -186.728504 -179.091564 -186.728504 -179.091564 -186.728504 -179.091564 -179.091564
4d 10 -29.432703 -29.453418 150.979227 151.198976 545
ec(mu) -28.796634 -30.438595 -28.796634 -30.438595 -28.796634 -30.438595 -28.796634 -30.438595 -28.796634 -28.796634
5d 10 -1.566638 -1.562069 23.907636 23.945913 605
ec(mu) -1.485638 -1.676716 -1.485638 -1.676716 -1.485638 -1.676716 -1.485638 -1.676716 -1.485638 -1.485638
4f 14 -9.755569 -9.751307 117.412788 117.457023 569
ec(mu) -9.592725 -9.962749 -9.592725 -9.962749 -9.592725 -9.962749 -9.592725 -9.962749 -9.592725 -9.962749 -9.592725 -9.962749 -9.592725 -9.592725
qcore(SR) 78.000000 qcore(FR) 78.000000 rho(rmax) 0.00000
sum ec : -25841.9031 (SR) -25934.9233 (FR) diff -93.0203
sum tc : 48113.1010 (SR) 48677.3220 (FR) diff 564.2210
```

The scalar Dirac Pb 5*d* eigenvalue (**-1.566638 Ry**) gets split into 6 levels with energy **-1.485638 Ry** and four with **-1.676716 Ry**. The mean (**-1.56207 Ry**) is close to the scalar Dirac value. In the absence of a magnetic field a particular *l* will split into two distinct levels with degeneracies 2*l* and 2*l*+2, respectively.

The bottom part of the table shows how much the free atom’s total energy changes as a consequence of the fully relativistic Dirac treatment.

#### 5.3 Automatic determination of envelope function parameters

Some details of the basis set (envelope functions, augmentation, local orbitals) are explained in this tutorial.

*lmfa* loops over each species, generating a self-consistent density.

Given a density and corresponding potential, *lmfa* will construct some estimates for the basis set, namely the generation of envelope function parameters **RSMH** and **EH** (and possibly **RSMH2** and **EH2**, depending on the setting of **HAM_AUTOBAS_MTO**), analyzing which cores should be promoted to local orbitals, and reasonable estimates for the boundary condition of the partial wave.

- Envelope functions
- The envelope functions (smoothed Hankel functions) are characterized by
**RSMH**and**EH**.**RSMH**is the Gaussian “smoothing radius” and approximately demarcates the transition between short-range behavior, where the envelope varies as , and asymptotic behavior where it decays exponentially with decay length , where is one of the**EH**.*lmfa*finds an estimate for**RSMH**and**EH**by fitting them to the “interstitial” part of the atomic wave functions (the region outside the augmentation radius). - Fitting the smooth Hankel function to the numerically tabulated exact function is usually quite accurate. For Pb, the error in the energy (estimated from the single particle sum) is 0.00116 Ry — very small on the scale of other errors. The fitting process is described in more detail in the annotated lmfa output.
*lmf*requires**RSMH**and**EH**. Those generated by*lmfa*are reasonable, but unfortunately not optimal choices for the crystal, as explained in the annotated lmfa output. You can change them by hand, or optimize them with*lmf*’s optimizing function,`--opt`

. To make an accurate basis, a second envelope function is added through**RSMH2**and**EH2**. (*lmfa*automatically does this, depending on the setting of**HAM_AUTOBAS_MTO**). Alternatively you can add APW’s to the basis. For a detailed discussion on how to select the basis set, see this tutorial.

*Note 1:*The envelope functions [ for the (**RSMH,EH**) group, for the (**RSMH2,EH2**) group], are augmented by partial waves in augmentation spheres. Thus the**lmf**basis set consists of augmented smooth Hankel functions.*Note: 2:*The new Jigsaw Puzzle Orbital basis is expected significantly improve on the accuracy of the existing Questaal basis. High quality envelope functions are automatically constructed that continuously extrapolate the accurate augmented partial waves smoothly into the interstitial; the kinetic energy of the envelope functions are continuous across the augmentation boundary.

- Local orbitals
*lmfa*searches for core states which are shallow enough to be treated as local orbitals, using the core energy and charge spill-out of the augmentation radius (**rmt**) as criteria; see annotated lmfa output.- When it was run for the first time,
*lmfa*singled out the Pb 5*d*state, using information from the table below taken from*lmfa*’s standard output. Once local orbitals are specified*lmfa*is able to appropriately partition the valence and core densities. This is essential because the two densities are treated differently in the crystal code. Refer to the annotated lmfa output for more details.

```
Find local orbitals which satisfy E > -2 Ry or q(r>rmt) > 5e-3
l=2 eval=-1.569 Q(r>rmt)=0.0078 PZ=5.934 Use: PZ=15.934
l=3 eval=-9.796 Q(r>rmt)=3e-8 PZ=4.971 Use: PZ=0.000
```

- Boundary conditions
- The free atomic wave function satisfies the boundary condition that the wave function decay as
*r*→∞. Thus, the value and slope of this function at**rmt**are determined by the asymptotic boundary condition. This boundary condition is needed for fixing the linearization energy of the partial waves in the crystal code.*lmfa*generates an estimate for this energy and encapsulates it into the “continuous principal quantum number”, saved as**P**in*basp0.pbte*(normally**P**will updated in the self-consistency cycle).

Refer to the annotated lmfa output for more details.

#### 5.4 Fitting the interstitial density

*lmfa* fits valence and core densities to a linear combination of smooth Hankel functions. This information will be used to overlap free-atomic densities to obtain a trial starting density. This is explained in the annotated lmfa output.

#### 5.5 Estimate for GMAX

After looping over all species *lmfa* writes basis information to *basp0.pbte*, atomic charge density data to file *atm.pbte*, and exits with the following printout:

```
FREEAT: estimate HAM_GMAX from RSMH: GMAX=4.6 (valence) 8.2 (local orbitals)
```

This is the *G* cutoff **EXPRESS_gmax** or **HAM_GMAX** needed to determine the interstitial mesh spacing. Two values are printed, one determined from the shape of valence envelope functions (**4.3**) and, if local orbitals are present the largest value found from their shape, as explained in the annotated lmfa output.

### 6. Remaining Inputs

*k* mesh

We are almost ready to carry out a self-consistent calculation. Try the following:

```
lmf ctrl.pbte
```

*lmf* stops with this message:

```
Exit -1 bzmesh: illegal or missing k-mesh
```

We haven’t yet specified a *k* mesh. You must supply it yourself since there are too many contexts to supply a sensible default value. Information is supplied through **BZ_NKABC**, which value is governed by variable **nkabc** in this input file. See this page for more details on how the Brillouin mesh is constructed, and how **nkabc** affects **BZ_NKABC**.

In this case a *k*-mesh of divisions is more than adequate. With your text editor or the sed command-line editor change **nkabc=0** in the ctrl file to **nkabc=6**, e.g.

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

or alternatively run *blm* again with the switch `--nk=6`

; or alternatively assign variable **nkabc** on the command line using **-vnkabc=6** (which is what this tutorial will do).

#### Charge density mesh

We also haven’t specified the *G* cutoff for the density mesh. *blm* does not determine this parameter automatically because it is sensitive to the selection of basis parameters, which local orbitals are included. *lmfa* conveniently supplies that information for us, based in the shape of envelope functions it found. In this case the valence *G* cutoff is quite small (**4.3**), but the Pb 5*d* local orbital is a much sharper function, and requires a larger cutoff (**8.2**). You must use the larger of the two.

*Note:* if you change the shape of the envelope functions you must take care that **gmax** is large enough. This is described in the lmf output below.

Change variable **gmax=0** in the ctrl file, or alternatively add a variable to the command line (**-vgmax=8.2**), as we do in the next section. Or, you can run *blm* again, with command-line argument **-gmax=8.2**.

### 7. Sanity Checks

This step shows how perform some sanity checks on the input. It is optional, but always advisable when you are getting started with a new system. Run *lmf* with a high verbosity, quitting after hamiltonian parameters are generated.

```
lmf ctrl.pbte -vnkabc=6 -vgmax=8.2 --pr55 --quit=ham
```

The output should confirm :

- the
*k*mesh is 6×6×6 divisions, which generates 16 points on the irreducible wedgeBZMESH: 16 irreducible QP from 216 ( 6 6 6 ) shift= F F F

**gmax=8.2**yields a charge density mesh of 20×20×20 divisionsMSHSIZ: mesh has 20 x 20 x 20 divisions; length 0.429, 0.429, 0.429

and this mesh is accurate enough to converge the plane-wave expansion of the smooth Hankel basis to about 10

^{−6}(estimated error is printed in table**sugcut**)- The basis consists of 30 Pb orbitals (
*spdfspd*+LO) and 25 Te orbitals (*spdfspd*)

Also notice the table `Orbital positions ...` . It locates the position of each orbital in the Hamiltonian, which is very useful when bands are to be plotted with color weights.

### 8. Self consistency

Carry out a self-consistent calculation as follows:

```
lmf ctrl.pbte -vnkabc=6 -vgmax=8.2
```

*lmf* will iterate up to 10 iterations. The cycle is capped to 10 iterations because of the following lines in *ctrl.pbte*, which before and after preprocessing read:

```
before preprocessing after preprocessing
% const nit=10 |
... |
EXPRESS | EXPRESS
... |
nit= {nit} | nit= 10
```

#### Initialization steps

*lmf* begins with some initialization steps. Each step is explained in more detail in the annotated lmf output.

- Read basis set parameters from
*basp.pbte*. This information can also be given via the ctrl file, depending on the settings of**EXPRESS_autobas**. - Informational printout about computing conditions, lattice structure, and atomic parameters such as augmentation radii and
*l*-cutoffs - Automatic determination of crystal symmetry
- Setup for the Brillouin zone integration
- Construction of the mesh for interstitial density and potential
- Assemble and display information about the size and constituents of the basis set
- Read or assemble an input density
- A QS
*GW*potential Σ^{0}may be read in.

#### Self-consistent cycle

Each iteration of the self-consistency cycle begins with

```
--- BNDFP: begin iteration 1 of 10 ---
...
--- BNDFP: begin iteration 2 of 10 ---
...
```

One iteration consists of the following steps. The standard output is annotated in some detail here.

- Construct the potential and matrix elements.
- Interstitial and local parts of the potential are made.
- Partial waves and are integrated from the potential subject to the boundary conditions.
- Matrix elements of the partial waves (kinetic energy, potential energy, overlap) are assembled for the Kohn-Sham hamiltonian.
- Matrix elements of the interstitial potential . for envelope functions .

This is sufficient to make the Kohn-Sham hamiltonian. Other matrix elements may be made depending on circumstances, matrix elements for optics or for spin-orbit coupling (

**HAM_SO=t**).

- Makes an initial pass through the irreducible
*k*points in the Brillouin zone to obtain the Fermi level and obtain integration weights for each band and*k*point into a binary file*wkp.pbte*. In general until the Fermi level is known, the weights assigned to each eigenfunction are not known, so the charge density cannot be assembled. How labor is divided between the first and second pass depends on**BZ_METAL**. See here for further discussion. - Makes a second pass to accumulate the output mesh and local densities. For the latter essential information is retained as coefficients of the local density matrix (a compact form).
- Assembles the local densities from the density matrix.
- Symmetrizes the density.
- Finds new logarithmic derivative parameters
**pnu**by floating them to the band center-of-gravity - Computes the Harris-Foulkes and Kohn-Sham total energies and forces.
Mixes the input and output densities to form a new trial density. A segment of the output reads:

`mixrho: sought 2 iter from file mixm; read 0. RMS DQ=2.16e-2`

**DQ=2.17e-2**is a measure of the root mean square deviation*n*^{out}−*n*^{in}. At self-consistency this number should be small.- Checks for convergence.

*lmf* should converge to self-consistency in 10 iterations.

At the end of the self-consistency cycle the density is written to *rst.pbte*

```
iors : write restart file (binary, mesh density)
```

and a check is made for convergence. No check is made in the first iteration because there is no prior iteration to compare the change in total energy. The second iteration reads:

```
diffe(q)= 0.004805 (0.018562) tol= 0.000010 (0.000030) more=T
i nkabc=6 gmax=8.2 ehf=-55318.1657749 ehk=-55318.1568595
```

Two checks are made: against the change (**0.004805**) in total energy and the RMS DQ (**0.018562**). When both checks fall below tolerances self-consistency is reached. In this case it occurs in iteration 10, where the convergence check reads:

```
diffe(q)= 0.000000 (0.000005) tol= 0.000010 (0.000030) more=F
c nkabc=6 gmax=8.2 ehf=-55318.1620975 ehk=-55318.1620958
```

The first line prints out the change in Harris-Foulkes energy relative to the prior iteration and some norm of RMS change in the charge density *n*^{out}−*n*^{in}, followed by the tolerances required for self-consistency.

The last line prints out variables specified on the command line, and total energies from the Harris-Foulkes and Kohn-Sham functionals. Theses are different functionals but they should approach the same value at self-consistency. The **c** at the beginning of the line indicates that this iteration achieved self-consistency with the tolerances specified. See the annotated output for more details.

### Other Resources

Click here to see annotated standard output from

*lmfa*, and here to see annotated standard output from*lmf*.An input file’s structure, and features of the programming language capability, is explained in some detail here. The full syntax of categories and tokens can be found in the input file manual.

This tutorial more fully describes some important tags the

*lmf*reads, and this one presents alternative ways to build input files from various sources such as the VASP*POSCAR*file.This tutorial more fully explains the

*lmf*basis set. There is a corresponding tutorial on the basics of a self-consistent ASA calculation for PbTe. A tutorial on optics can be gone through after you have finished the present one.This document gives an overview of some of

*lmf*’s unique features and capabilities.The theoretical formalism behind the

*lmf*is described in detail in this book chapter: M. Methfessel, M. van Schilfgaarde, and R. A. Casali, ``A full-potential LMTO method based on smooth Hankel functions,’’ in*Electronic Structure and Physical Properties of Solids: The Uses of the LMTO Method*, Lecture Notes in Physics,**535**, 114-147. H. Dreysse, ed. (Springer-Verlag, Berlin) 2000.

### FAQ

How does

*lmf*iterate to self-consistency?It mixes the input density

*n*^{in}with output density*n*^{out}generated by*lmf*, to construct a new input density*n*^{in}. This process is repeated until*n*^{out}=*n*^{in}(within a specified tolerance). The actual mixing algorithm can be quite involved; see this page.The gap is small and Pb is a heavy element. Doesn’t spin-orbit coupling affect the band structure?

Yes, it does. The bandgap will change significantly when spin-orbit coupling is added.

The LDA is supposed to underestimate bandgaps. But the PbTe bandgap looks pretty good. Why is that?

This turns out to be largely an accident. If spin orbit coupling is included, the bandgap appears to be pretty good, but in fact levels L

_{6}^{+}and L_{6}^{−}that form the valence and conduction band edges are inverted in the LDA. See Table I of this paper. As the paper notes, they are well described in QS*GW*.How do you know where the band edges are?

PbTe is has a quite simple band structure with high symmetry. It’s a good bet that the band edges are on high-symmetry lines. But in general the position of band edges can be quite complex. A slightly more complicated case is Si. See this tutorial.

Is there an easy way to calculate effective masses?

Yes, once you know where the band edge is. See this tutorial.

The augmentation

*l*cutoff in the ctrl file is 4 for Pb and 3 for Te. LAPW methods typically require much higher*l*cutoffs (6 or 8) to be well converged. Why is this? Is it a worry?Questaal uses a different kind of augmentation. It a difference in the basis set (smoothed Hankels vs LAPWs; Questaal in fact has an LAPW basis that may be combined with the augmented smooth Hankels), but because of the different way augmentation is carried out. Both kinds of augmentation converge to the same answer in the limit of large

*l*, but Questaal’s augmentation is much more rapidly convergent. It can be incorporated into APW basis sets as well. See this tutorial.

### Additional exercises

Try self-consistent calculations with the Pb 5

*d*in the valence as a local orbital. Repeat the calculation but remove the**PZ**part from*basp.pbte*.Specify symops manually.

Turn on spin orbit coupling and observe how the band structure changes.

Try rotations

k-convergence. Try BZ_BZJOB