# Properties of the lmf basis set

This tutorial describes the **lmf** basis set, and various kinds of cutoffs that affect convergence in the basis.

### Table of Contents

- Preliminaries
- Definition of the basis set
- Tutorial
- Additional exercises

Note, the **sed** commands here are optional, you can also enter these values by hand.

```
blm init.bi2te3
cp actrl.bi2te3 ctrl.bi2te3
sed -i s/nkabc=0/nkabc=3/ ctrl.bi2te3
sed -i s/gmax=0/gmax=4.4/ ctrl.bi2te3
sed -i 's/autobas\[pnu=1 loc=1 mto=4 lmto=4\]/autobas\[pnu=1 loc=1 lmto=3 mto=1 gw=0\]/' ctrl.bi2te3
```

```
lmchk bi2te3
```

Generating basp files/Obtaining GMAX values

```
lmfa bi2te3
mv basp0.bi2te3 basp.bi2te3
lmfa bi2te3
```

```
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=2
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=3
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=4
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=5
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=6
cat save.bi2te3 | vextract h nkabc ehf | tee dat
```

Small basis self-consistent LDA/Charge density convergence

```
rm -f mixm.bi2te3
rm -f save.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3
```

```
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
lmf ctrl.bi2te3 -vgmax=6 -vnkabc=3 --quit=band
```

```
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp ctrl.bi2te3 ctrl.bak
cat ctrl.bak | sed s/LMXA=./LMXA=6/ > ctrl.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp ctrl.bak ctrl.bi2te3
```

```
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp ctrl.bi2te3 ctrl.bak
cat ctrl.bak | sed 's/LMXA=/KMXA=5 LMXA=/' > ctrl.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp ctrl.bak ctrl.bi2te3
```

```
rm -f mixm.bi2te3
rm -f save.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --optbas
cp basp.bi2te3 basp.bak
cp basp2.bi2te3 basp.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp basp.bak basp.bi2te3
```

Increasing the number of envelope functions

```
sed -i 's/autobas\[pnu=1 loc=1 lmto=3 mto=1 gw=0\]/autobas\[pnu=1 loc=1 lmto=4 mto=4 gw=0\]/' ctrl.bi2te3
lmfa bi2te3
cp basp.bi2te3 basp.bak
cp basp0.bi2te3 basp.bi2te3
sed -r -i "s/PZ=(([ ]{1,})?[0-9]+(\.[0-9]{1,})?){3}//" basp.bi2te3
lmfa bi2te3
rm -f mixm.bi2te3
rm -f rst.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3
```

Adding APWs

Note: for these commands you need to modify the *ctrl* file, please see the linked section for more information.

```
rm -rf out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=0 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=1 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=2 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=3 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=4 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=5 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=6 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=7 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
cat save.bi2te3 | vextract i pwemax ehf ehk | tee etot.bi2te3
fplot -frme 0,.8,0,.5 -frmt th=3,1,1 -tmy 2.5 -vehf=-126808.313403 -s circ -ord '(x2-ehf)*1000' etot.bi2te3
```

Modifying the input file for GW

```
blm --gw bi2te3
```

### Preliminaries

The input file structure is briefly described in this lmf tutorial for PbTe, which you may wish to go through first.

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

### Definition of the basis set

The **lmf** basis set consist of smooth Hankel functions envelope functions , which get augmented by solutions of the radial wave equation in augmentation spheres (partial waves). An additional set of local orbitals may be added to make the basis more complete in the augmentation region.

#### Envelope functions

The **lmf** basis set begins with smooth envelope functions. These functions are smooth Hankel functions, centered on an atom **R**. They are defined by smoothing radius and an energy . Smooth Hankels have the Slater-Koster form,

Here *L* is a compound index denoting the *l* and *m* angular momentum quantum numbers; the are radial functions defined by Eqns. (6-8) on this page, and the are spherical harmonics.

The are convolutions of ordinary Hankel functions of energy and Gaussian functions of smoothing radius . These two parameters define the shape of the envelope: controls the asymptotic decay for large *r* (, ) and the smoothing radius demarcates the approximate point of transition from Gaussian-like shape at small *r* () to asymptotic behavior.

Each site has its own family of . While it would be possible to have any number of on a particular site, each with unique values of and , in practice the Questaal codes allow two types () of for a particular site **R** and angular momentum *l*. For the first envelope (), you define pair through parameters **RSMH** and **EH**: Each *l* gets its own **RSMH** and **EH**. You can elect to limit the basis to one envelope function for a particular *l*, or choose a second envelope. The second set () is optional: parameters are defined through **RSMH2** and **EH2**. A single (“single kappa”) has the advantage of making for a small basis, but its accuracy is limited unless the system is fairly close-packed. At least for low *l*, i.e. for orbitals whose atomic levels are below or not far above the Fermi level, it is recommended that the “double kappa” basis (two functions per *l*) be used.

This flexibility in choosing is both a blessing and a curse. A blessing because the flexibility allows for more variational freedom, keeping the basis at low rank while maintaining high accuracy. But it is a curse because of the added burdens imposed on the user to determine the parameters. Usually you can allow the atom program **lmfa** to automatically generate parameters for basis set.

Thus we can identify the entire family of envelopes by , labeling whether the first or second envelope, the site, and the *L*. Note that for each *L*, there are 2*l*+1 basis functions.

#### Augmentation

Each smooth envelope is “augmented” in a sphere around each nucleus by *partial waves*, which are numerical solution of the Schrodinger equation in the sphere up to an augmentation radius (specified by tag **R** in the species category). For each *l* (and *m* for a given *l*) up to a cutoff **lmxa**, the is replaced (“augmented”) by a linear combination of partial waves and . These partial waves form the Hilbert space of essentially exact solutions to the Schrodinger equation, to linear order in energy around some linearization energy . This energy is normally allowed to float in the course of the self-consistency cycle, to minimize the linearization error. See also the description of partial waves on this page. Suitable linear combinations of and are taken in each sphere so as to make the augmented continuous and differentiable.

#### Local orbitals

It may be the linear method is not sufficient and that a third partial wave is required to make the basis complete over a sufficiently wide energy window. A conventional local orbital is realized by integrating the Schrodinger equation at another energy either far above or far below the energy used to make and , to obtain another partial wave . A suitable combination of and is subtracted to make the value and slope of the local orbital vanish at the augmentation radius. With such a construction it need not have an interstitial part. An “extended” local orbital (suitable for semicore states far below the Fermi level) is constructed by attaching a single smooth Hankel function, whose energy and smoothing radius are varied to match value and slope of the partial wave. You specify the local orbital through tag **PZ** in the species category. **lmfa** can automatically find deep local orbitals (one principal quantum number below the valence partial waves) that may be needed to make the basis reasonably complete. You may also select high-lying local orbitals (one principal quantum number above the valence partial waves). These are usually less relevant in DFT; however, high-lying *d* local orbitals on transition metals were necessary to obtain good agreement with other codes in the DeltaCodes validation exercise. Note that addition of local orbitals increases the rank of the Hamiltonian.

#### Hilbert space and rank of the lmf Hamiltonian

The Hilbert space of the **lmf** basis then consists of the following:

In the interstitial, smooth Hankel functions . They can be expanded in plane waves to make matrix elements of the potential.

*Note:*For a given set or you must include all*l*’s up to some basis cutoff, i.e.*l*=0,1,…**lmxb**.**lmxb**need not be the same for the set of envelopes as for the first.In augmentation sphere

**R**up to the augmentation*l*cutoff**lmxa**, the pair of partial waves (, ). The*dominant*partial wave is ; mostly it attaches on to smoothed Hankels centered at**R**. mostly attaches on to the (tails) of smooth Hankels centered at other sites and makes a*smaller contribution*. Finally there may possibly be local orbitals in some*l*channels. Local orbitals may be high-lying (far above the linearization energy) or deep, to include high-lying (semi)core levels in the valence. This extra set is specified through parameter**PZ**in the main input (*ctrl*file) file or the*basp*file.In augmentation sphere

**R**above**lmxa**, the tails of smoothed Hankels at sites other than**R**. The special manner in which augmentation is done enables the**lmf**basis set to converge very rapidly with**lmxa**— much faster than traditional all-electron basis sets. See below.

The total rank of the hamiltonian is then the number of you specify on all sites, plus any local orbitals specified. The size of the basis (number of first kappa, second kappa, and local orbitals) is printed in a table in **lmf**’s standard output.

### Tutorial

#### 1. Building the input file

This step is essentially identical to the first step in the PbTe tutorial. An abbreviated version is presented here.

Cut and paste the following into *init.bi2te3*.

```
# from http://cst-www.nrl.navy.mil/lattice/struk/c33.html
# Bi2Te3 from Wyckoff
% const a=4.3835 c=30.487 uTe=0.788 uBi=0.40
LATTICE
SPCGRP=R-3M
UNITS=A
A={a} C={c}
SITE
ATOM=Te C=0 0 0
ATOM=Te C=0 0 {uTe}
ATOM=Bi C=0 0 {uBi}
```

This tutorial explains how the input files *init.ext* and *ctrl.ext* are structured.

To create the skeleton input file invoke **blm**:

```
$ blm init.bi2te3
$ cp actrl.bi2te3 ctrl.bi2te3
```

There 2 Bi atoms and 3 Te atoms in the unit cell. Two of the Te atoms are symmetry equivalent.

This template will not work as is; three essential pieces of information which **blm** does not supply are missing, to rectify this:

- You must specify plane-wave cutoff
**GMAX**. This can be obtained by running the**lmfa**tool as seen in the PbTe tutorial. For the purpose of this tutorial we will use the value**GMAX=4.4**. A more detailed description of**GMAX**is given below. - You must specify a valid
**k**mesh for Brillouin zone integration. The correct**k**mesh size is dependent on your problem, you can find an example in the PbTe tutorial. For this tutorial we will use a**k**mesh of divisions,**nkabc=3**. A more detailed description of**nkabc**is given below. - You must define a basis set which can be done manually or automatically, as described next.

Note: if you are using the **GMAX** and **k**-mesh values listed above, you can use the following commands to automatically edit your *ctrl.bi2te3* file:

```
$ sed -i s/nkabc=0/nkabc=3/ ctrl.bi2te3
$ sed -i s/gmax=0/gmax=4.4/ ctrl.bi2te3
```

In either case, to edit these values you simply need to find the **gmax** and **nkabc** variables in the *% const* section of your *ctrl.bi2te3* file and modify them to the desired value. The **sed** command above does this for you.

**blm** creates a family of tags belonging to AUTOBAS to enable other programs to automatically find a basis set for you. We will use this tag, which sets up a standard minimal basis:

```
autobas[pnu=1 loc=1 lmto=3 mto=1 gw=0]
```

This is an alias for the tag in the **HAM** category

```
AUTOBAS[PNU=1 LOC=1 LMTO=3 MTO=1 GW=0]
```

Note that you must modify *ctrl.bi2te3* a little; the default gives [… LMTO=5 MTO=4], which makes a double kappa basis. Once again, you can use the following command to perform this replacement automatically:

```
$ sed -i 's/autobas\[pnu=1 loc=1 mto=4 lmto=4\]/autobas\[pnu=1 loc=1 lmto=3 mto=1 gw=0\]/' ctrl.bi2te3
```

**lmfa** calculates wave functions and atomic densities for free atoms. It also has a mode that automatically generates information for basis sets, using tokens in **AUTOBAS** to guide it. This information is written to a file *basp0.ext*. **AUTOBAS** specifies set of conditions that enable **lmfa** to automatically build the basis set for you, as described below. Parameters are generated, but you can modify them as you like.

Tokens in **AUTOBAS** tell **lmfa** to do the following:

- PNU=1 Calculate the logarithmic derivative parameter
*P*for the free atom._{l} - Calculated parameters are saved in file
*basp0.ext*as**P=…**. Nothing about**P**is written if**PNU=0**. - LOC=1 Look for shallow cores to be explicitly treated as valence electrons, as local orbitals.
- Shallow cores that meet specific criteria are identified and written to
*basp0.ext*as**PZ=…**. No search is made if**LOC=0** - LMTO=3 Pick a default choice about the size of basis. LMTO=3 is a standard minimal basis.
- Run
`lmfa --input`

and look for**HAM_AUTOBAS_LMTO**to see what other choices there are.**lmfa**will pick some defaults for the*l*-cutoff, e.g.*spd*or*spdf*depending on the value of**LMTO**. - MTO=1 Choose 1-kappa basis set (single orbital per
*l*channel). - Small basis for fast calculations. For better quality calculations, it is recommended you use
**MTO=2**or**MTO=4**. - GW=0 Create a setup for an LDA calculation.
- If
**GW=1**, tailor basis for a GW calculation, selecting stricter criteria for including shallow cores as valence, and the size of basis.

These tokens thus define some reasonable default basis for you. **lmfa** *writes* *basp0.ext*. This file is never read, but **lmf** will *read* *basp.ext* and use this information when assembling the basis set. The two files have the same structure and you can copy *basp0.ext* to *basp.ext*. What **lmfa** generates is not cast in concrete. You are free to adjust the parameters to your liking, e.g. add a local orbital or remove one from the basis.

The **AUTOBAS** tokens tell **lmf** what to read from *basp.ext*. It uses tokens in a manner similar, but not identical, to the way **lmfa** uses them:

- PNU=1 Read parameters P for all species present in
*basp.ext*. - If PNU=0, these parameters will not be read.
- LOC=1 tells
**lmf**to read local orbital parameters**PZ**. - Since these parameters may also be specified by the input file,

LOC=1 tells**lmf**to give precedence to parameters specified by ctrl file

LOC=2 tells**lmf**to give precedence to parameters specified by basp.

**LMTO=** is not used by **lmf**.

**MTO=1**controls what part of**RSMH**and**EH**is read from*basp.bi2te3*.- LMTO=1 or 3 tells lmf to read 1-kappa parameters specified by basp

LMTO=2 or 4 tells lmf to read 2-kappa parameters specified by basp

LMTO=1 or 2 tells lmf that parameters in the ctrl file take precedence

LMTO=2 or 4 tells lmf that parameters in the basp file take precedence - GW=0 tunes the basis for an LDA calculation
- If
**GW=1**, tune basis for a**GW**calculation. For example log derivative parameters P are floated a little differently in the self-consistency cycle. They are weighted to better represent unoccupied states, at a slight cost to their representation of occupied states.

#### 2. Checking sphere overlaps

Sphere overlaps can be checked using **lmchk**. To do this copy the template *actrl.bi2te3* to the input file and run **lmchk**:

```
$ lmchk bi2te3
```

You should see the site positions for all the atoms:

```
Site Class Rmax Hcr Position
1 1 Te 2.870279 2.009195 0.00000 0.00000 0.00000
2 3 Te2 2.870279 2.009195 -0.50000 -0.86603 1.46162
3 3 Te2 2.870279 2.009195 0.50000 0.86603 -1.46162
4 2 Bi 2.856141 1.999299 0.50000 0.86603 0.80309
5 2 Bi 2.856141 1.999299 -0.50000 -0.86603 -0.80309
```

and a table of overlaps

```
ib jb cl1 cl2 Pos(jb)-Pos(ib) Dist sumrs Ovlp % summt Ovlp %
1 4 Te Bi 2.391 -4.142 3.841 6.134 5.726 -0.41 -6.6 4.008 -2.13 -34.7
1 5 Te Bi -2.391 -4.142 -3.841 6.134 5.726 -0.41 -6.6 4.008 -2.13 -34.7
2 4 Te2 Bi -2.391 -4.142 -3.149 5.726 5.726 0.00 0.0* 4.008 -1.72 -30.0
3 5 Te2 Bi 2.391 -4.142 3.149 5.726 5.726 0.00 0.0* 4.008 -1.72 -30.0
```

By default, **blm** makes the spheres as large as possible without overlapping. In this case the Bi and Te radii are nearly the same.

The packing fraction is printed

```
Cell volume= 1141.20380 Sum of sphere volumes= 492.34441 (0.43143)
```

Generally larger packing fractions are better because the augmentation partial waves are quite accurate. **0.43** is a fairly good packing fraction.

#### 3. The atomic density and basis set parameters

If you did not do so already copy *actrl.bi2te3* to *ctrl.bi2te3* and change **[… LMTO=4 MTO=4]** → **[… LMTO=3 MTO=1]**).

Invoke **lmfa**:

```
$ lmfa bi2te3
```

The primary purpose of **lmfa** is to generate a free atom density. A secondary purpose is to supply additional information about the basis set in an automatic way. All of this information can be supplied manually in the input file, but the **blm** provides a minimum amount of information. **lmfa** generates *basp0.bi2te3* which contains

```
BASIS:
Te RSMH= 1.615 1.681 1.914 1.914 EH= -0.888 -0.288 -0.1 -0.1 P= 5.901 5.853 5.419 4.187
Bi RSMH= 1.674 1.867 1.904 1.904 EH= -0.842 -0.21 -0.1 -0.1 P= 6.896 6.817 6.267 5.199 5.089 PZ= 0 0 15.936
```

Every species gets one line. This file specifies a basis set consisting of *spdf* orbitals on Te sites, and *spdf* orbitals on Bi sites, and a local *5d* orbital on Bi. See the PbTe tutorial for further description of these parameters.

*Note:* Remember that **lmf** reads from *basp.ext*, not *basp0.ext*.

The partitioning between valence and core states is something that requires a judgement call. **lmfa** has made a default choice: from *basp0.bi2te3* you can see that **lmfa** selected the 6*s*, 6*p*, 6*d*, 5*f* states, populating them with charges 2, 3, 0, 0, making the total sphere charge Q=0. You can override the default, e.g. choose the 5*d* over the 6*d* with SPEC_ATOM_P; override the *l* channel charges with SPEC_ATOM_Q.

**lmfa** does the following to find basis set parameters:

- Automatically finds deep states to include as valence electrons.
- Selects sphere boundary condition for partial waves
- Finds envelope function parameters

The process is essentially the same as described in the PbTe tutorial; it is described in some detail there and in the annotated **lmfa** output. The main difference is that in this case a smaller, single-kappa basis was specified (**LMTO=3**); **lmfa** makes (**RMSH,EH**) instead of the double kappa (**RMSH,EH**; **RMSH2,EH2**). Later we will improve on the basis by adding APW’s.

With your text editor insert lines from *basp0.bi2te3* in the **SPEC** category of the ctrl file, viz

```
ATOM=Te Z= 52 R= 2.870279 LMX=3 LMXA=3
RSMH= 1.615 1.681 1.914 1.914 EH= -0.888 -0.288 -0.1 -0.1 P= 5.901 5.853 5.419 4.187
ATOM=Bi Z= 83 R= 2.856141 LMX=3 LMXA=4
RSMH= 1.674 1.867 1.904 1.904 EH= -0.842 -0.21 -0.1 -0.1 P= 6.896 6.817 6.267 5.199 5.089 PZ= 0 0 15.936
```

Alternatively make **lmf** read these parameters from *basp.bi2te3*. Copy *basp0.bi2te3* to *basp.bi2te3*, and modify it as you like. *basp.ext* is read after the main input file is read, if it exists. According to which of following tokens is present, their corresponding parameters will be be read from the basp file, superseding prior values for these contents:

AUTOBAS | lmfa writes, lmf reads |

MTO | RSMH,EH (and RSMH2,EH2 if double kappa basis) |

P | P |

LOC | PZ |

If this information is supplied in both the ctrl file and the basp file, values of **MTO** and **LOC** tell **lmf** which to use. In this tutorial we will work just with the basp file.

```
$ cp basp0.bi2te3 basp.bi2te3
```

The atm file was created by **lmfa** without prior knowledge that the 5*d* local orbital is to be included as a valence state (via a local orbital). Thus it incorrectly partitioned the core and valence charge. You must do one of the following:

Remove

**PZ=0 0 15.936**from*basp.bi2te3*. It will no longer be treated as a valence state. Removing it means the remaining envelope functions are much smoother, which allows you to get away with a coarser mesh density. Whether you need it or not depends on the context: with*GW*, for example, this state is a bit too shallow to be treated with Fock exchange only (which is how cores are handled in*GW*).Copy

*basp0.bi2te3*to*basp.bi2te3*and run**lmfa**over again:

```
$ cp basp0.bi2te3 basp.bi2te3
$ lmfa bi2te3
```

With the latter choice **lmfa** operates a little differently from the first pass. Initially the Bi 5*d* was part of the core (similar to the Pb 5*d* in the Pbte tutorial; now it is included in the valence.

#### 4. GMAX and NKABC

##### 4.1 Setting GMAX

**blm** makes no attempt to automatically assign a value to the plane wave cutoff for the interstitial density. This value determines the mesh spacing for the charge density. **lmf** reads this information through **HAM_GMAX** or **EXPRESS_gmax**, or **HAM_FTMESH**. It is a required input; but **blm** does not automatically pick a value because its proper choice depends on the smoothness of the basis. In general the mesh density must be fine as the most strongly peaked orbital in the basis.

**lmfa** makes **RSMH** and **EH** and can determine an appropriate value for **HAM_GMAX** based on them. In the present instance, when the usual 6*s*, 6*p*, 6*d*, 5*f* states are included **lmfa** recommends **GMAX=4.4** as can be seen by inspecting the first **lmfa** run. If you keep the 5*d* in the valence and run **lmfa** second time you will see this output

```
FREEAT: estimate HAM_GMAX from RSMH: GMAX=4.4 (valence) 8.1 (local orbitals)
```

The valence states still require only **GMAX=4.4** but the 5*d* state is strongly peaked but **GMAX=8.1** for the local orbitals. The 5*d* state is strongly peaked at around the atom, and requires more plane waves to represent reasonably, even a smoothed version of it, than the other states. The difference between 8.1 and 4.4 is substantial, and it reflects the additional computational cost of including deep core-like states in the valence. This is the all-electron analog of the “hardness” of the pseudopotential in pseudopotential schemes. If you want high-accuracy calculations (especially in the *GW* context), you will need to include these states as valence. Including the Bi 5*d* is a bit of overkill for LDA calculations however. If you eliminate the Bi 5*d* local orbital you can set GMAX=4.4 and significantly speed up the execution time. It is what this tutorial does.

##### 4.2 Setting NKABC

**blm** assigns the initial *k*-point mesh to zero. Note the following lines in *ctrl.bi2te3*:

```
% const nkabc=0 gmax=0
...
# Brillouin zone
nkabc= {nkabc} # 1 to 3 values
```

*Note:* **nkabc** is simultaneously a *variable* and a *tag* here. This can be somewhat confusing. The expression **{nkabc}** gets parsed by the preprocessor and is turned into the value of *variable* **nkabc** (see how **nit** gets turned into **10** in the PbTe tutorial), so after preprocessing, the argument following *tag* is a number.

**EXPRESS_nkabc** (alias **BZ_NKABC**) governs the mesh of *k*-points. An appropriate choice will depend strongly on the context of the calculation and the sytem of interest; the density-of-states at the Fermi level; whether Fermi surface properties are important; whether you want optical properties as well as total energies well described; the precision you need; the integration method, and so on. Any automatic formula can be dangerous, so **blm** will not choose an operational default.

The most reliable way determine the appropriate mesh is to vary **nkabc** and monitor the convergence. We don’t need a self-consistent calculation for that: the k-convergence will not depend strongly whether the potential is converged or assembled from free atoms.

Do the following (assuming no 5*d* local orbital)

```
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=2
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=3
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=4
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=5
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=6
```

The meaning of `--quit=band`

is described here.

Total energies are written to the save file *save.bi2te3*. It should read:

```
h gmax=4.4 nkabc=2 ehf=-126808.2361717 ehk=-126808.1583957
h gmax=4.4 nkabc=3 ehf=-126808.3137885 ehk=-126808.2492178
h gmax=4.4 nkabc=4 ehf=-126808.3168406 ehk=-126808.2505156
h gmax=4.4 nkabc=5 ehf=-126808.3165536 ehk=-126808.2497121
h gmax=4.4 nkabc=6 ehf=-126808.3164058 ehk=-126808.2494041
```

You can use the **vextract** tool to conveniently extract a table of the Harris-Foulkes total energy as a function of nkabc

```
$ cat save.bi2te3 | vextract h nkabc ehf | tee dat
```

You can plot the data, or just see by inspection that the energy is converged to less than a mRy with 4×4×4 *k* mesh and about 0.1 mRy with a 5×5×5 *k* mesh. A detailed analysis of *k* point convergence is given here.

#### 5. Self-consistent LDA calculation, small basis

In the following we will use **nkabc=3**, though after convergence is reached we might consider increasing it a little. Before doing the calculation you can use your text editor to set **nkabc=3** and **gmax=4.4**, or continue to set assign values from the command line.

The density can be made self-consistent with

```
$ rm -f mixm.bi2te3 save.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3
```

##### 5.1 Convergence in the charge density

You should **lmf** reach self-consistency in 9 iterations.

The Harris-Foulkes and Kohn-Sham total energies are **ehf=-126808.3137885** and **ehk=-126808.2492178** from the Mattheis construction. At self-consistency they come together: **ehf=-126808.2950696** and **ehk=-126808.2950608** (the large integer part comes almost entirely from the core states.)

The self-consistent value is 18 mRy *less binding* than the Harris-Foulkes energy of the Mattheis construction, and 46 mRy *more binding* than the corresponding Kohn-Sham energy. That the two initial functionals bracket the self-consistent result, and that the HF is generally closer to the final result than the HK functional, is typical behavior. (The HF functional is generally more stable.)

##### 5.2 Convergence in G cutoff

The *G* cutoff, **4.4u** Ry^{1/2} yields precisions in plane-wave expansions of envelope functions as shown this table:

```
sugcut: make orbital-dependent reciprocal vector cutoffs for tol=1.0e-6
spec l rsm eh gmax last term cutoff
Te 0 1.61 -0.89 4.603 3.31E-06 1635*
Te 1 1.68 -0.29 4.662 5.08E-06 1635*
Te 2* 1.91 -0.10 4.273 1.15E-06 1539
Te 3 1.91 -0.10 4.471 1.71E-06 1635*
Bi 0 1.67 -0.84 4.441 1.29E-06 1635*
Bi 1 1.87 -0.21 4.183 1.08E-06 1411
Bi 2* 1.90 -0.10 4.297 1.37E-06 1539
Bi 3 1.90 -0.10 4.497 2.06E-06 1635*
```

**gmax=4.4** isn’t quite large enough to push the error below tolerance (**1.0e-6**), but it is pretty close. You can check how well the total energy is converged by increasing **gmax**:

```
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
$ lmf ctrl.bi2te3 -vgmax=6 -vnkabc=3 --quit=band
```

and comparing **ehf** in the last two lines of *save.bi2te3*. You should find that the energy is converged to about 0.1 mRy.

##### 5.3 Convergence in the forces

The exact forces are the change in total energy with respect to displacement of the nucleus. As they are calculated in the Questaal package, some correction terms are added to accelerate convergence with respect to deviations *n*^{in}−*n* in the guessed density from the self-consistent one. See the annotation of **lmf** output.

*Note:* The forces are not exact derivatives of the total energy. This is because the change in shape of the augmented partial waves and is not taken into account as a nucleus displaces. The effect is usually very small. Forces should be exactly consistent with the energy if the shape of the partial waves is frozen, however, which you can do with **HAM_FRZWF**.

To what extent the basis set affects the forces is taken up in the Additional Exercises.

##### 5.4 Convergence in LMXA

In an augmented wave method, envelope functions centered a different site is must be expanded around a local site, as one-center expansions in spherical harmonics *Y _{lm}*.

**LMXA**is a cutoff that truncates the one-center expansion to a finite

*l*=

**LMXA**. The input file reads:

```
SPEC
ATOM=Te Z= 52 R= 2.870279 LMX=3 LMXA=3
ATOM=Bi Z= 83 R= 2.856141 LMX=3 LMXA=4
```

**LMXA=3** and **LMXA=4** are very low *l* cutoffs for an augmented wave method. They are about half of what is needed for a reasonably converged calculation with traditional augmentation. However, the Questaal suite has a unique augmentation procedure that converges very rapidly with *l*. Indeed, it can be seen as a justification for pseudopotential methods that limit the pseudopotential to very low *l*, e.g. *l*=2.

It is a useful exercise to see how well converged the total energy is using the default values **LMXA=3** and **LMXA=4**. First, run **lmf** with the unaltered ctrl file:

```
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
```

Set both LMXA to 6: and run **lmf** again

```
$ cp ctrl.bi2te3 ctrl.bak
$ cat ctrl.bak | sed s/LMXA=./LMXA=6/ > ctrl.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
```

When the restart file is read you should see

```
site 1, species Te : augmentation lmax changed from 3 to 6
site 1, species Te : inflate local density from nlm= 16 to 49
```

Compare the last two lines of *save.bi2te3*. You should be able to confirm that the energy change is 0.25 mRy, By playing around with the two **LMXA** you can confirm that increasing the Te **LMXA** to 4, nearly all the error is eliminated.

Before continuing, be sure to restore the original ctrl file

```
$ cp ctrl.bak ctrl.bi2te3
```

##### 5.5 Convergence in KMXA

Ordinary Hankel functions can be expanded in Bessel functions around a remote site. This follows from the fact that both are solutions of the same second order differential equation, one being regular at the origin and the other being regular at ∞. Smooth Hankel functions are more complicated: the one-center expansion has no such elementary function, but it can be written as a linear combination of Laguerre polynomials *P _{kl}* of half integer order; see Eq. (12.17) in this J. Math. Phys. paper.

The polynomials are cut off at a fixed order given by **KMXA**. **blm** doesn’t write **KMXA** to the template file; instead a default value is used, which is written to this table

```
species data: augmentation density
spec rmt rsma lmxa kmxa lmxl rg rsmv kmxv foca rfoca
Te 2.870 1.148 3 3 3 0.718 1.435 15 1 1.148
Bi 2.856 1.142 4 3 4 0.714 1.428 15 1 1.142
```

You can check the convergence by in **KMXA** by supplying an input for it. First run **lmf** with the unmodified file:

```
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
```

Set **KMXA** to 5 and run **lmf** again:

```
$ cp ctrl.bi2te3 ctrl.bak
$ cat ctrl.bak | sed 's/LMXA=/KMXA=5 LMXA=/' > ctrl.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
```

When the restart file is read you should see an indication that **KMXA** has been increased:

```
site 1, species Te : augmentation kmax changed from 3 to 5
site 2, species Te : augmentation kmax changed from 3 to 5
```

Compare the last two lines of *save.bi2te3*. You should be able to confirm that the energy change is 0.15 mRy.

*Note:* Before continuing, be sure to restore the original ctrl file.

```
$ cp ctrl.bak ctrl.bi2te3
```

#### 6. Optimizing the basis set

##### 6.1 Tuning the envelope function parameters

The envelope function shape generated by **lmfa** is tailored to the atom. They are very good basis sets for that case, but less so for the crystal.

Ideally a basis of the size could be constructed that yields energies converged almost as well as we accomlish for the free atom. This is difficult to do in practice, though it is the aim of the Jigsaw Puzzle Orbitals basis.

The best we can with the existing **lmf** basis is to optimize the energy by tuning **RSMH** and **EH**. This cannot be done analytically, but switch **--optbas** causes **lmf** to perform this optimization by brute-force minimization of the Harris-Foulkes energy.

First we need to regenerate the self-consistent density we generated previously:

```
$ rm -f mixm.bi2te3
$ rm -f save.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3
```

Then we can use the **--optbas** switch:

```
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --optbas > out
```

This will loop through **RSMH** in each species, minimizing the total energy for each one by one. It does not vary **EH** because the total energy is less sensitive to it. **lmf** will print out a table of the parameters it will optimize:

```
LMFOPB: optimizing energy wrt 8 parameters, etol=5e-4:
spec l type start
1:Te 0 rsm 1.615
1:Te 1 rsm 1.681
1:Te 2 rsm 1.914
1:Te 3 rsm 1.914
2:Bi 0 rsm 1.674
2:Bi 1 rsm 1.867
2:Bi 2 rsm 1.904
2:Bi 3 rsm 1.904
```

Each cycle carries out non self-consistent calculations to get the Harris-Foulkes total energy. After completing **lmf** prints out

```
LMFOPB: before optimization ehf=-126808.2951. After optimization ehf=-126808.3092
Exit 0 Optimization complete. See file basp2
```

By tuning the parameters the energy decreased by 0.0141 Ry. To see which orbitals contributed the most, do:

```
$ grep LMFOPB out | grep optimized
```

You should see

```
LMFOPB: var #1 optimized to 1.602 (7 iter): ehf=-126808.295097, delta ehf=-2.38e-5
LMFOPB: var #2 optimized to 1.617 (5 iter): ehf=-126808.298989, delta ehf=-0.00392
LMFOPB: var #3 optimized to 1.778 (6 iter): ehf=-126808.299166, delta ehf=-1.77e-4
LMFOPB: var #5 optimized to 1.605 (5 iter): ehf=-126808.301, delta ehf=-9.12e-4
LMFOPB: var #6 optimized to 1.658 (7 iter): ehf=-126808.307945, delta ehf=-0.00695
LMFOPB: var #7 optimized to 2.245 (5 iter): ehf=-126808.308856, delta ehf=-9.11e-4
LMFOPB: var #8 optimized to 1.679 (6 iter): ehf=-126808.309198, delta ehf=-3.42e-4
```

The orbitals that mad the most difference were the Te 5*p* (**1.681** → **1.617**) and Bi 5*p* (**1.867** → **1.658**).

Optimized parameters where the total energy decreased by more than **etol** (**5e-4**) were modified; the optimized parameters are written to *basp2.bi2te3*

Next, we are going to modify the basp file. You should make a backup of the current file first:

```
$ cp basp.bi2te3 basp.bak
```

You can copy *basp2.bi2te3* to *basp.bi2te3*, but there is little point in optimizing any but the Te 5*p* and Bi 6*p*. Instead we use a text editor and modify *basp.bi2te3* in those channels:

```
BASIS:
Te RSMH= 1.615 1.617 1.914 1.914 EH= -0.888 -0.288 -0.1 -0.1 P= 5.901 5.853 5.419 4.187
Bi RSMH= 1.674 1.658 1.904 1.904 EH= -0.842 -0.21 -0.1 -0.1 P= 6.896 6.817 6.267 5.199 5.089
```

Run **lmf** again

```
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
```

and confirm that **ehf** comes out close (-126808.3086945) to the optimized value (-126808.309198).

You can at this point make the density self-consistent again.

*Note:* Before continuing, be sure to restore the original basp file:

```
$ cp basp.bak basp.bi2te3
```

In sum,

- without optimizing the basis, you should achieve total energy
**ehf=-126808.2950696**at self-consistency. - With the
**--optbas**switch it should become**ehf=-126808.306962**.

##### 6.2 Increasing the number of envelope functions

**lmfa** uses the **lmto** tags to make default choices for the size of basis. The tutorial used a relative minimal basis, **lmto=3**.

```
autobas[pnu=1 loc=1 lmto=3 mto=1 gw=0]
```

With your editor, modify this line to read

```
autobas[pnu=1 loc=1 lmto=4 mto=4 gw=0]
```

This will increase the basis, most notably include two envelope functions per *l* channel (double-kappa)

Run **lmfa** again and note how the basp file changed:

```
$ lmfa bi2te3
$ diff basp0.bi2te3 basp.bi2te3
```

Remove the **PZ=…** as before and copy *basp0.bi2te3* to *basp.bi2te3*. Save the original file in a backup.

```
$ nano basp0.bi2te3
$ cp basp.bi2te3 basp.bak
$ cp basp0.bi2te3 basp.bi2te3
```

Since we modified the *basp* file we need to re-run **lmfa**:

```
$ lmfa bi2te3
```

Run **lmf** to self-consistency

```
$ rm -f mixm.bi2te3
$ rm -f rst.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3
```

Note that size of the basis has increased: With **lmto=3** there are 80 orbitals

```
Makidx: hamiltonian dimensions Low, Int, High, Negl: 80 0 18 27
suham : 98 augmentation channels, 98 local potential channels Maximum lmxa=4
```

which increases to 125 orbitals with **lmto=5**:

```
Makidx: hamiltonian dimensions Low, Int, High, Negl: 125 0 71 54
kappa Low Int High L+I L+I+H Neglected
1 80 0 18 80 98 27
2 45 0 53 45 98 27
all 125 0 71 125 196 54
suham : 98 augmentation channels, 98 local potential channels Maximum lmxa=4
```

At self-consistency you should find that the total energy converges to **ehf=-126808.313403** Ry.

##### 6.3 Adding APWs : the PMT method

You can also use augmented plane waves (APWs) to improve on the basis set. The combination of smooth Hankel functions with APW’s is called the PMT basis set. Adding APW’s provides a systematic way of making the basis of envelope functions complete.

These tags in the **HAM** category control how many APWs are added:

```
PWMODE={pwmode} PWEMIN=0 PWEMAX={pwemax} # For APW addition to basis
```

To use APW’s some tolerances in the ctrl file should be tightened. If you invoke **blm** with **--pmt** these extra tolerances are included automatically. That is the easiest way to make the changes but here we just add include tolerances **HAM_TOL** and **EWALD_TOL** by hand. Insert two lines at the end of the **HAM** category:

```
HAM
...
FORCES={so==0} ELIND=-0.7
TOL=1d-16
EWALD TOL=1d-16
```

One problem with the **PMT** method is that the smoooth-Hankel and APW basis set span nearly the same hilbert space. If you add too many plane waves the overlap matrix does not remain positive definite. Tightening the tolerances above helps, and so does increasing the fineness of the density mesh, as these points are also used for the PW expansion of the basis.

*Note:* as a last resort, you can stabilize the overlap with the **HAM_OVEPS** switch, but it is better to reduce the rank of the LMTO or APW basis sets.

**PWMODE=11** adds APWs in a standard way: envelope functions of the form *e*^{i(k+G)·r} are added to the smooth Hankel basis. You must specify the PW cutoffs. Here you can specify both the lower and upper bounds since the smooth Hankels will efficiently cover most of the lower part of the space.

With these extra considerations, run **lmf** as follows

```
rm -f out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=0 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=1 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=2 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=3 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=4 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=5 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=6 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=7 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
```

and observe how the total energy decreases with **pwemax**. Use the **vextract** tool:

```
cat save.bi2te3 | vextract i pwemax ehf ehk | tee etot.bi2te3
```

Draw a picture of the total energy relative to the double-kappa value. The **fplot** command shown in the box below will generate a postscript file; the figure actually shown has been elaborated a little. The red circle shows the self-consistent double-kappa result of Sec. 6.2. The light grey line follows the PMT procedure as above, but taking for the LMTO part a smaller, optimized *spd* single kappa basis (see Additional exercises). Numbers in parenthesis are the number of LMTO’s. The red circle is the 2-kappa basis without APWs; the purple is the LAPW basis without LMTOs.

```
$ fplot -frme 0,.8,0,.5 -frmt th=3,1,1 -tmy 2.5 -vehf=-126808.313403 -s circ -ord '(x2-ehf)*1000' etot.bi2te3
```

The figure shows that the double-kappa basis (additional 45 orbitals) stabilizes the single-kappa *spdf* basis by about 18 mRy, and that the same can be accomplished with APWs with a plane-wave cutoff of 2 Ry (additional 48-62 APWs). By further increasing the APW cutoff, another 5 mRy can be gained. For most purposes this extra gain is not important. Note that the APW basis with *E*_{max}=7 Ry is quite large: 343-370 APWs.

Note that the APW basis is generally less efficient than the LMTO basis. To reach a precision comparable to the 2-kappa basis (125 orbitals) starting from the 1-kappa *spd* basis, about 160 APWs are needed, or about 200 orbitals all told. The power of the PMT method can be compared against a straight LAPW basis. About 300 APWs are needed to achieve the convergence of the double kappa basis. See purple line in the Figure.

To see how many orbitals the APW basis adds, do:

```
$ grep ndham out.lmf
```

#### Modifying the input file for GW

*GW* calculations demand more of the basis set because unoccupied states are important. To set up a job in preparation for a *GW* calculation, invoke **blm** as :

```
$ blm --gw bi2te3
```

Compare *actrl.bi2te3* generated with the **--gw** switch to one without. One important difference will be that the default basis parameters are modified because **AUTOBAS** becomes:

```
AUTOBAS[PNU=1 LOC=1 LMTO=5 MTO=4 GW=1]
```

The basis is similar to **LMTO=4** but **EH** has been set a little deeper. This helps the QS*GW* implementation interpolate between *k*-points. The larger basis makes a minor difference to the valence bands; but the conduction bands change, especially the higher in energy you go.

Note also that the LMTO f orbitals are more efficient in converging the total energy per extra orbital added than plane waves are.

### Additional exercises

Reduce the smooth Hankel basis set of Section 5 by eliminating the

*f*orbitals. This reduces the basis set to 9 orbitals/atom. You should find that the total energy is**ehf= -126808.271025**without**--optbas**, and**ehf=-126808.280117**with it.Add APWs to this basis, and observe that the total energy converges to the same value. Note also that the LMTO

*f*orbitals are more efficient in converging the total energy than plane waves are. You should get something similar to the grey line in the Figure of Section 6.3.Try an all LAPW basis : use

**pwmode=12**. You should get something similar to the purple line in the Figure of Section 6.3. Results start to get reasonable around**pwemax=4.5**.At self-consistency, the forces should reach convergence for given basis set, as described in the annotated

**lmf**output. Forces should be derivatives of the total energy wrt nuclear displacements. As the basis set becomes complete, the forces should stop changing. Confirm that this is the case by varying the basis set. Bi_{2}Te_{3}has 5 atoms; two of the Te atoms have one force equal and opposite, and the two Bi another force, also equal and opposite. You should find something similar to the following:Basis F(Te) F(Bi) 1-kappa (optimized) 10.6 15.7 2-kappa 7.5 13.1 1-kappa(opt)+ APW(pwemax=2) 8.0 13.3 1-kappa(opt)+ APW(pwemax=6) 8.0 13.4