# Introductory lmf tutorial

This tutorial provides a light introduction to running a self-consistent DFT calculation for silicon, using the main density-functional code in the Questaal suite, *lmf*. The goal is to introduce you to the different file types and the basics of running the code, as propaedeutic for other tutorials such as this more detailed tutorial for PbTe. See also a corresponding introductory tutorial for Questaal’s ASA program, *lm*.

Follow-on tutorials draw energy bands in Si, and find band edges and effective masses.

```
mkdir si && cd si
```

Create *init.si*

```
cat << END > init.si
LATTICE
# ALAT is a length scale, in atomic units
ALAT=10.26
# PLAT refers to primitive lattice vectors, in units of alat
PLAT= 0.00 0.50 0.50
0.50 0.00 0.50
0.50 0.50 0.00
# POS refers to cartesian coordinates, in units of alat
SITE
ATOM=Si POS= 0.00 0.00 0.00
ATOM=Si POS= 0.25 0.25 0.25
END
```

Make the *ctrl* and *atm* files

```
blm init.si --simple --wsite --nk,ins --ctrl=ctrl
lmfa ctrl.si
```

Make a self-consistent density

```
lmf ctrl.si > out.lmfsc
```

Alternatively, once you have created *init.si*, you can run the following one-liner to get the same result.

```
blm init.si --simple --wsite --nk,ins --ctrl=ctrl && lmfa ctrl.si && lmf ctrl.si > out.lmfsc
```

### Preliminaries

This tutorial assumes you have downloaded the repository (see this page for instructions on getting access to and downloading the Questaal package), that you have built (compiled) the executable codes and installed them in a some directory in your path, e.g. *~/bin*.

The source code for all Questaal executables can be found here.

### Tutorial

Create a new working directory and move into it, e.g. “**si**” as shown below.

```
mkdir si && cd si
```

Most Questaal codes read from the main input file, *ctrl.ext* (*ctrl.si* in this case). The *ctrl* file is roughly akin to the *POSCAR* file in VASP. *ctrl* files are most easily built from an *init* file which contains mostly structural information. *blm*’s main purpose is to create *ctrl* files; usually it reads structural data from *init* files, which take a format shown in the box below, and generate input (*ctrl*) files that many other Questaal executables read as their primary input stream. In this tutorial we create *init.si* by hand. *init* files can be made in other ways, e.g. importing structural data from *POSCAR* files using the utility *poscar2init*.

Create a file named *init.si* and copy the following lines to it:

```
LATTICE
# ALAT is a length scale, in atomic units
ALAT=10.26
# PLAT refers to primitive lattice vectors, in units of alat
PLAT= 0.00 0.50 0.50
0.50 0.00 0.50
0.50 0.50 0.00
# POS refers to cartesian coordinates, in units of alat
SITE
ATOM=Si POS= 0.00 0.00 0.00
ATOM=Si POS= 0.25 0.25 0.25
```

*Note*: be sure the **LATTICE** and **SITE** tags are in the leftmost column of *init.si*, while the other lines are indented. (Tags that appear in the first column are called categories; they are there to organize the input by categories or groups of tags that belong together.

**ALAT** is the lattice constant for Si, which is in atomic units by default (Questaal uses atomic Rydberg units). The primitive lattice vectors are listed after **PLAT** and are in row format (i.e. the first row contains the *x*, *y* and *z* components of the first lattice vector and so forth). The **SITE** section specifies the atom types and their positions in cartesian coordinates (indicated by **POS**) in units of the lattice constant in **ALAT**. More generally **ALAT** is a length scale that scales **PLAT** and **POS**, and sometimes other length-bearing quantities (the augmentation radii can be specified in units of **ALAT**, for example).

More information can be found in the Additional Exercises.

The DFT program *lmf* requires more information than the structural information supplied in *init.si*. The *blm* utility reads the *init* creates a template *actrl.si* (and typically an auxiliary file *site.si* with purely structural information). The Questaal package recognises certain names as file types such as

*ctrl.ext*the main input file, which*blm*will generate*site.ext*, an auxiliary input file with purely structural information, which in this tutorial,*blm*will also generate.*basp.ext*, an auxiliary input file defining parameters for basis, which in this tutorial,*blm*will also generate.*atm.ext*, a file which contains self-consistent densities of isolated free atoms. It is generated by*lmfa*, as shown shortly

and there are many more, depending on the context. Extension *si* labels the family of files associated with this materials system; it is something you choose. Most files read or generated by Questaal programs will append this extension.

Run *blm* as shown below. The **--simple** switch tells *blm* to make a particularly simple input file; we will see more complicated examples in later tutorials.

```
blm init.si --simple --wsite --nk,ins --ctrl=ctrl
lmfa ctrl.si
```

*blm*’s command line switches are documented here. You can quickly get an abbreviated version of the documentation by running `blm --h`

.

**--simple**: tells*blm*to create a simple input file, with minimal information**--wsite**: tells*blm*to park the structural (lattice vectors and site positions) into an independent file*site.si*. Without this switch, this information is embedded directly in*ctrl.si*.**--nk,ins**: tells*blm*to generate a default*k*mesh assuming an insulator, written to tag**NKABC**.*Warning*:*blm*’s default values are only estimates. Proper values are materials- and property- specific. The user should take care that the default values are reasonably converged.**--ctrl=ctrl**: tells*blm*to write the input file directly to*ctrl.si*. By default it will write the input to a template file,*actrl.si*, so as not to overwrite an already existing input file.

With the command as written, *blm* makes three input files (*ctrl.si*, *basp.si*, *site.si*). If you leave out **--wsite** and add **--genbas,nofile**, *blm* will make a single, standalone *ctrl.si* file.

When *blm* is invoked with switches as given, followed by *lmfa*, the two codes construct a complete set of input files: the *ctrl* file, the *basp* file, the *site* file, and the *atm* file, ready for Questaal’s band code, *lmf*.

Now look at the input file generated by *blm*. *ctrl.si* should look like the following. Note that text following a **#** on any line is treated as a comment.

```
# Autogenerated from init.si using:
# blm init.si --simple --wsite --nk,ins --ctrl=ctrl
#symgrp i*r3(-1,-1,1)::(1/4,1/4,1/4)
#symgrp spglib
ham
gmax= 0
autobas[pnu=1 loc=1 mto=4]
forces= 0
nspin= 1
so= 0
xcfun="lm_lda_bh"
iter
nit= 99
mix= b2,b=0.3,k=7
conv= 1e-5
convc= 3e-5
bz
nkabc= 6
struc
file=site
spec
atom=Si z= 14 r= 2.221355 lmx=3 lmxa=5
```

Data parsed from the *ctrl* file is identified by particular tag, which consists of a token (e.g. **GMAX**) and the category it belongs to (**HAM**). The full tag is then **HAM_GMAX**. Any text beginning in the first column marks the start of a new category, with the text naming the category. All data within a category must be indented at least one space. Input organized by *categories* and *tokens* is a convenient way to group information that belongs together.

*Note:* some tags actually have three parts, e.g. **HAM_AUTOBAS_PNU**, **HAM_AUTOBAS_LOC** (and also **SPEC_ATOM_Z**, **SPEC_ATOM_R**, **SPEC_ATOM_LMX**, **SPEC_ATOM_LMXA**). The tags belonging to **HAM_AUTOBAS** are contained in brackets, **PNU=1 LOC=1 MTO=4** in this case.

The input file can contain many more categories and tags than is shown here. All the input codes read from the *ctrl* file is documented on the ctrl file page. The categories shown above have the following meaning:

**VERS**is a category for version control. It checks that the input file is consistent with the code you are running.**HAM**contains specifications for the construction of the hamiltonian.**BZ**contains specifications for Brillouin zone integration (also energy integrations in the Green’s function context).**STRUC**contains structural data that is not specific to a particular nucleus or site. In this case it reads the information from file*site.si*.**SITE**contains structural data that is specific to a particular site. In this case it reads the information from file*site.si*.**SPEC**contains chemical information for each species (atom type) in the crystal.

The tags in *ctrl.si* are as follows:

Tag | Explanation |
---|---|

VERS_LM | Major version number for the entire Questaal package |

VERS_FP | Major version number for the full potential code lmf |

HAM_GMAX | G-vector cutoff for the plane wave representation of smoothed Hankel envelope functions and the smooth density |

HAM_AUTOBAS_PNU | PNU is an -dependent continuous principal quantum number that determines the boundary condition for the partial wave in the augmentation sphere. AUTOBAS_PNU=1 tells the input file reader to look for this information in file basp. blm created basp.si together with ctrl.si. |

HAM_AUTOBAS_LOC | PZ is an -dependent continuous principal quantum number that determines the boundary condition for low-lying core-like (and also high-lying states far above the Fermi level) states that are added to the usual linear method via the method of local orbitals. AUTOBAS_LOC=1 tells the input file reader to look for this information in file basp. In this instance there are no local orbitals. |

HAM_AUTOBAS_MTO | Tells the input file reader to look for parameters RSMH, EH, RSMH2, EH2 that define the shape of the envelope functions, in file basp. The basis consists of envelope functions made of smooth Hankel functions, augmented by numerical solutions of partial waves in the augmentation spheres. A smooth Hankel is defined by smoothing radius RSMH and an energy EH. There is one function for each , and AUTOBAS_MTO=4 tells the input file reader to seek for as many as two envelope functions on a site for each (RSMH, EH, and RSMH2, EH2) |

BZ_NKABC | The number of divisions on each of the three directions of the reciprocal lattice vectors. k points are generated on a uniform mesh along each of these axes. The parser will attempt to read three integers. If only one number is read, the missing second and third entries assume the value of the first. |

STRUC_FILE | Tells the parser to read structural information from an independent file, site.si in this case. It reads ALAT, and PLAT from the file and also determines how many sites there are (NBAS) and the number of species, i.e. atom types (NSPEC). |

SPEC_ATOM | A label identifying the species (aka atom type). These labels are used ubiquitously in the code and they must synchronize with the labels in the SITE category. |

SPEC_ATOM_Z | The atomic number associated with SPEC_ATOM |

SPEC_ATOM_R | The augmentation radius associated with SPEC_ATOM |

SPEC_ATOM_LMX | The cutoff for the smooth-Hankel basis. Here , which corresponds to a basis of spd orbitals (recall that there are two envelope functions for each ) |

SPEC_ATOM_LMXA | The cutoff for the augmentation of tails of envelope functions from sites other than the central one. Strictly SPEC_ATOM_LMXA should be twice the largest relevant in the basis. Si’s d orbitals () mostly lie above the Fermi energy, so the contribution to the density from d-d products is small. Thus LMXA=3 is sufficient. This is largely true for any element that is not a transition metal or an f shell metal. |

SITE_FILE | Tells the parser to read site information from an independent file, site.si in athis case. In this instance site.si contains only site positions. |

Note that all the **HAM_AUTOBAS** tags point the parser to read data from the *basp* file. *blm* was invoked in a way that automatically makes this file. You can tell *blm* not to create this file, but to put the information directly in the *ctrl* instead (use `--genbas~nofile`

). Similarly *blm* can put structural data directly into the *ctrl* file (use `--wsite~nofile`

).

Recall that *blm* was invoked with the switch **--nk,ins**. It tells *blm* to make a default **k** mesh on the reciprocal lattice, assuming the material is an insulator. The mesh is uniformly spaced along the three directions and discretizes the Brillouin zone into contiguous parallelpipeds; it is needed for numerical integrations over the zone. There is no good general way for *blm* to select a mesh in general, as it depends on the band structure and details of the Fermi surface, and also the precision needed for the physical property of interest; thus it is something you should check. If you omit this switch, *blm* will assign **NKABC= NULL**, and you will have to edit *ctrl.si* before runing *lmf*.

The start of the *blm* output shows some structural and symmetry information. It does this first by identifying the crystal system, which sets the point group and an upper bound to the allowed symmetry operations. The site positions may reduce the allowed symmetry operations; it seeks all group operations within the known crystal system consistent with the site positions. This determines the space group.

Next it must find the augmentation sphere radii. It does this by finding the radius corresponding to the maximum value of the potential around each sphere, and later scaling all the radii by a single factor to reach a specified overlap. In this case it used the default overlap of 0, i.e. touching spheres. You can trace what *blm* is doing in the output following **makrm0**.

Next it generates information for the **k** mesh, because *blm* was invoked with **--nk** (in this case **--nk,ins**), which it uses to estimate the total number of *k* points in the Brillouin zone and then tries to divide the reciprocal lattice on the three directions to make the spacing as equal a possible. *blm* prints what it finds, which in this case reads

```
Set k mesh for supplied nkabc = 6 ...
BZMESH: 16 irreducible QP from 216 ( 6 6 6 ) shift= F F F
```

Thus the mesh is one with 4 divisions on each axis, and with the symmetry operations it found (48 in all), the 4×4×4 = 64 **k** points reduce to 8 inequivalent points.

Next, it generates information about the basis set. What *blm* does in this respect depends on the style of input file you select (**--brief[= n], --simple**, or

**--express[=**) and whether and how

*n*]**--genbas[~options]**is invoked. In the present context it generates all the information needed :

**RSMH, EH**, continuous principal quantum numbers for the valence partial waves

**P**and local orbitals

**PZ**, if any (none in this case). All of this is saved in

*basp.si*. Further, it finds a good estimate for the plane wave cutoff

**GMAX**. After

*ctrl.si*has been created, it overwrites the token for

**GMAX**with the value it finds.

Thus *blm* and *lmfa* supply all the necessary input for a density-functional calculation, without any special knowledge required by the user, other than structural information.

A great deal more could be said about the code’s basis set and other features, but we will leave further details to later tutorials and the input file documentation.

The Questaal methods paper, Ref. 1 provides a thorough formal development of Questaal’s implementation of density-functional theory: its basis set, construction of the hamiltonian, and representation of the charge density. The theory of smooth Hankel functions is developed in Ref. 3. Questaal is an all-electron linear augmented-wave method, and so it does not suffer from approximations inherent in pseudopotential methods (such methods are also linear methods, but approximate formulations of them; see Ref. 5). See Ref. 4 for Andersen’s classic paper where the theory of linear methods was first formulated. The reliability of Questaal’s standard set approaches the “Gold Standard” Linear Augmented Plane Wave method, but because the basis is minimal the uncontrolled approximations can lead to errors, especially in very open systems where the sphere packing is low. Workarounds in the present version (**LM:7**) are to add plane waves to the basis (the **PMT** method, Ref. 6; see also this tutorial), or to add interstitial sites or points with envelope functions but no augmentation spheres (floating orbitals). Questaal’s next generation (version 8) will include Jigsaw Puzzle Orbitals, which should enable Questaal’s standard basis to be similar LAPW fidelity, but minimal at the same time.

#### Checking sphere overlaps

It is advisable before starting a calculation, as a sanity check, to see whether the augmentation spheres are properly sized. Do this with *lmchk* utility.

```
lmchk ctrl.si
```

You should see that the autogenerated muffin-tin radius (**2.221355**) causes the spheres to touch but not overlap (maximum overlap = 0%). Thanks to of its special form of augmentation, you can allow the spheres to overlap little (up to about 10%) with little adverse consquence. Making the spheres larger increases the volume of the highly accurate partial waves, but there it is compensated by an error from the geometry violation.

#### Self-consistency

We now have everything we need to run an all-electron, full-potential DFT calculation. This is done using the *lmf* program. A lot of text is produced so it will be easier to redirect the output to a file, here we call it *out.lmfsc* (appending sc to indicate a self-consistent cycle).

```
lmf ctrl.si > out.lmfsc
```

Before looking at the output, do the following:

```
grep DQ out.lmfsc
```

**DQ** is the RMS change in the charge density as the program cycles through iterations. You can see it becomes small as iterations proceed; thus the density is self-consistent.

Now look at the output, *out.lmfsc*. It begins with lattice structure, symmetry, and the **k** point mesh, and is similar to the blm output.

Then follows a presentation of the basis set:

```
Makidx: hamiltonian dimensions Low, Int, High, Negl: 50 0 94 0
...
```

The table following this line shows each Si has 13 orbitals (16 *spdf* orbitals for the first envelope functions **RSMH,EH**, 9 *spd* orbitals for the second envelope functions **RSMH2,EH2**)

Look for the following lines:

```
lmfp : no rst file ... try to overlap atomic densities
rdovfa: read and overlap free-atom densities (mesh density) ...
```

The first line above tells us about what input density is used. *lmf* first looks for a restart file *rst.si* and if it’s not found it then looks for the free atom density file *atm.si*. *lmf* then overlaps the free atom densities to form a trial density (Mattheis construction) and this is used as the input density. Next *lmf* begins the first iteration of a self-consistent cycle: calculate the potential from the input density, use this potential to solve the Kohn-Sham equations and then perform an integration over the Brillouin zone to get the output density. The input and output densities are mixed and the cycle is repeated until **DQ** becomes small. Towards the end of each cycle, the Kohn-Sham total energy is reported along with the Harris-Foulkes total energy. These two energies will be the same (or very close) at self-consistency.

Now move to the end of the file. The ** c ** in front of the Harris-Foulkes **ehf** and Kohn-Sham **ehk** energies indicates that convergence was reached (note how similar the **ehf** and **ehk** energies are). A few lines up you can see that it took 7 iterations to converge: look for “**it 7 of 10**”. At the end of each iteration the **ehf** and **ehk** total energies are printed and a check is made for self-consistency. Default values are used for the tolerances, since the tags that read them are not present in *ctrl.si*. To see what the tolerancees are, do the following:

```
lmf ctrl.si --input
```

and look for **ITER_CONV** and **ITER_CONVC**. The first is the tolerance in the change in **ehf** from the prior iteration (default = **1e-5**); the second is the root mean square change in the density **DQ** of the current iteration (default = **3e-5**). Note that by default both tolerances have to be met. To use a single tolerance you simply set the one that you don’t want to zero. (You must add a category **ITER** to *ctrl.si* and include one or both tokens **CONV**, **CONVC**.)

Further down the Fermi energy and band gap, and other key bits of information are reported in the Brillouin zone integration section. You should find something similar to the output snippet below.

```
... only filled or empty bands encountered: ev=0.177908 ec=0.222268
VBmax = 0.177908 CBmin = 0.222268 gap = 0.044360 Ry = 0.60355 eV
BZINTS: Fermi energy: 0.177908; 8.000000 electrons; D(Ef): 0.000
Sum occ. bands: -1.5637863 incl. Bloechl correction: 0.0000000
```

The bandgap, 0.604 eV is about 0.06 eV larger than the exact LDA value with the Barth-Hedin functional (computed with a larger basis set), and smaller than the experimental one, as is typical in density functional theory.

To see how the density, bandgap, and energy changes between iterations, try grepping for “**DQ**” and “**ehk=-**”. The **RMS DQ** measures the RMS change in the density between iterations; **ehk** is the Hohenberg-Kohn Sham total energy (also just called Kohn-Sham).

```
egrep 'DQ|gap|ehk=' out.lmfsc
```

As iterations proceed, the RMS **DQ**→0 and **ehf**→**ehk**. The bandgap from the Mattheis construction is larger than the self-consistent one, as is typical.

To get the cohesive energy, you need to subtract the total energy of the free atom from **ehf** *lmfa* prints this out: it is **etot=-577.640888** Ry for Si. The (negative of the) cohesive energy is than −1156.1532−2×−577.6408 = −0.87161 Ry for the unit cell. It is overbound relative to experiment. Part of the overbinding comes because the free atom must be computed spin polarized; the rest is because limits to the accuracy of the local density approximation.

And that’s it! You now have a self-consistent density and have calculated some basic properties such as the band gap and total energy. Other tutorials show how to generate energy band structures, and density-of-states, or calculate a mechanical property such as the optical mode frequency.

### Other Resources and References

This tutorial on PbTe also covers the basic self-consistency cycle, in a more detail. It has a companion tutorial for the ASA, allowing you to compare the FP and ASA methods. There is a more detailed tutorial with some description of important tags the *lmf* reads, and of the *lmf* basis set.

The output of *lmf*, taken from the PbTe tutorial, is annotated in detailed on this page. See this page for a corresponding annotated output of *lmfa*.

For a tutorial showing a self-consistent calculation of ferromagnetic metal, see this tutorial. This tutorial shows how to calculate optical properties for PbTe. See this tutorial for the calculation of the optical mode frequency in Si. This document gives an overview of the *lmf* implementation; see also References, below.

#### References

[1] Questaal’s main methods paper:

Dimitar Pashov, Swagata Acharya, Walter R. L. Lambrecht, Jerome Jackson, Kirill D. Belashchenko, Athanasios Chantis, Francois Jamet, Mark van Schilfgaarde, *Questaal: a package of electronic structure methods based on the linear muffin-tin orbital technique*, Comp. Phys. Comm. **249**, 107065 (2020). This is Questaal’s main methods paper.

[2] This paper describes the method of the ful-potential code *lmf* as it was first formulated; it is the antecedent to Questaal.

M. Methfessel, Mark 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. H. Dreysse, ed. (Springer-Verlag, Berlin)* 2000.

[3] This paper presents the formal development of smooth Hankel functions, invented by M. Methfessel.

E. Bott, M. Methfessel, W. Krabs, and P. C. Schmidt, *Nonsingular Hankel functions as a new basis for electronic structure calculations*, Journal of Mathematical Physics 39, 3393 (1998).

[4] This classic paper established the framework for linear methods in band theory:

O. K. Andersen, *Linear methods in band theory* Phys. Rev. B12, 3060 (1975)

[5] This paper presents an excellent general presentation of electronic structure within density functional theory, with a heavy focus on plane wave, pseudopotential methods.

R. M. Martin, *Electronic Structure*, Cambridge University Press (2004).

[6] These papers describe the PMT method: an enhancement of the Questaals standard basis with plane waves.

T. Kotani and M. van Schilfgaarde, *A fusion of the LAPW and the LMTO methods: the augmented plane wave plus muffin-tin orbital (PMT) method*, Phys. Rev. B81, 125117 (2010).

See also T. Kotani, H. Kino, H. Akai, *Formulation of the Augmented Plane-Wave and Muffin-Tin Orbital Method*, J. Phys. Soc. Jpn. 84, 034702 (2015).

### FAQ

Below is a list of frequently asked questions. Please get in contact if you have other questions.

- 1) What are the real and reciprocal lattices?
- A crystal is defined as an object that is unchanged by integer multiple of primitive lattice translation vectors , with integers.

are primitive lattice vectors (three members of which make the smallest nonzero volume). See Wikipedia for more details.

In Questaal’s notation are used for (dimensionless) primitive lattice vectors, with , is the lattice constant read as**ALAT**in the input file (atomic units).

The reciprocal lattice are the family of vectors defined from three primitive vectors, . is related by .

Questaal uses for reciprocal primitive lattice vectors, with .

The simplest way to relate and is to construct 3×3 matrices and .

Then . Most Questaal programs print this information near the beginning of the standard output. - 2) What is the relation between Cartesian coordinates and coordinates expressed as fractions of lattice vectors (aka crystal coordinates)?
- The connection between positions
**POS**expressed in Cartesian coordinates and**X**expressed as fractional multiples of the primitive lattice vectors is:

POS_{i}= X_{1}P_{i,1}+ X_{2}P_{i,2}+ X_{3}P_{i,3}where P_{i,j}is Cartesian component*i*of lattice vector*j*.The expression above equally applies in reciprocal space if the primitive reciprocal lattice vectors

**Q**are substituted for**P**.

Sometimes called “crystal coordinates” (VASP uses the term “direct coordinates”), they are used when site positions or**k**points are expressed in fractional multiples of the lattice vectors (**ALAT**in real space,**2π/ALAT**in**k**space).

This example converts a**k**point given in Cartesian coordinates to one in multiples of the reciprocal lattice vectors, using Questaal conventions. - 3) What is the
**Harris-Foulkes**energy? - It is a functional of the input density and input potential, rather than the output density and input potential. At self-consistency it should be the same as the standard Kohn-Sham functional, and like the Kohn-Sham functional, it is stationary at the self-consistent density. The original motivation was to provide a framework for DFT analogs of tight-binding theory, since if a good trial input density could be found, self-consistency can be avoided. In practice the Harris-Foulkes functional is valuable because it tends to be the most stable in the approach to self-consistency. But it is not necessarily a minimum point, even if it is a stationary one. See this paper by M. Foulkes and R. Haydock.
It can be useful to have both KS and HF functionals: the form is a minimum while the latter tends to be a maximum (though the latter is not guaranteed). In such a case the two functionals bounds the self-consistent total energy, and it is an interesting exercise to see how the two approach the same value. It is not trivial that the two functionals yield the same result; indeed if tolerances are not tight enough (e.g. the mesh density typically set by parameter

**GMAX**), they may differ somewhat. - 4) What is the Brillouin zone?
- For Questaal’s purpose the Brillouin zone is the Wigner Seitz cell of the reciprocal lattice, sometimes called the first Brillouin zone. Bloch’s theorem establishes that wave number
**k**is a good quantum number for eigenfunctions of the crystal, and that its values of**k**are restricted to the (first) Brillouin zone. (Eigenfunctions of any state with quantum number**k**+**G**are equivalent to those of**k**.) Typically in solid state calculations integrations must be performed over**k**, e.g. for the total energy or charge density as was done in this tutorial. In practice this is done by discretizing**k**, dividing the Brillouin zone into microcells: parallelopipeds whose edges have size . then define the discretization of**k**; they are read from token**nkabc**described in this tutorial. (**nkabc**is an alias for BZ_NKABC.) Integration over**k**is approximated by a sum over the discrete points: each microcell gets equal weight. The denser the discretized**k**mesh the more accurate the numerical integration, though the rate of convergence depends strongly on density-of-states at the Fermi level (best convergence is for insulators which have no density of states there). From a convergence point of view, the spacing between points on each direction should usually be as similar as possible, convergence is only as good as the worst case. - 5) How does
*blm*determine the augmentation sphere radii? *blm*overlaps free atom densities and adjusts the radii to make the potential as similar as possible on the surface of each sphere. You can constrain the radius to a value of your choosing with species token**CSTRMX**.- 6) What is the
*log*file? - File
*log.si*keeps a compact record of key outputs in the current directory. In successive runs, data is appended to this file.

### Additional Exercises

**1) Short version of tutorial with one-liner command**

The entire tutorial can be run with the following one-liner:

```
blm init.si --gmax=7 --brief --wsite --nk,ins --ctrl=ctrl && lmfa ctrl.si && lmf ctrl.si > out.lmfsc
```

Create a new directory and try running the above command, starting from the *init* file. N.B. You should be careful about running calculations in the same directory as certain files (such as the *rst*, *mix* and *moms* files) will be read and used. Note that here we use the more standard **--brief** rather than **--simple**. Compare this new *ctrl* file to the original one and you will see that the file is different, and provides more options and makes use of Questaal’s powerful preprocessor. with variables that can be set on the command line. Note also that **GMAX** was supplied by hand using a command-line argument. The output will be slightly different, but if you run *blm* with the auto-generate value (**--gmax=6.957**) it should be essentially identical. The changes in output give a measure of how well converged the default value (**6.957**) is.

**2) Accurate bandgap**

The bandgap printed out by the code is not the actual LDA gap, but the smallest separation between the highest occupied and lowest unoccupied state it found on a discrete 6×6×6 k mesh. The actual minimum occurs near k=(1,0,0), commonly referred to as the X point. The X point is on the 6×6×6 k-mesh, but the conduction band minimum itself is not quite at X. Rerun the calculation with a very fine **k** mesh (not self-consistently this time) and observe that the bandgap is slightly smaller. Use your text editor to set **nkabc** to 12 or 16 and do:

```
lmf ctrl.si --quit=rho
```

You should find that the gap gets smaller by about 0.13 eV, which is close to the experimentally observed splitting between the conduction minimum and the conduction band at X. Instead of running a very fine **k** mesh, you can also do a gradients minimization as described in the Extremal points and effective mass tutorial.

**3) Structural data in other formats ( POSCAR or CIF)**

If you have structural data in another format, you can simply convert to an *init* file and repeat the same steps. The *cif2init* and *poscar2init* tools convert *CIF* and *POSCAR* files, respectively, to *init* files. For example, create a new directory and copy the following lines to a file called *POSCAR*:

```
system Si
5.430
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
Si
2
Direct
0.00 0.00 0.00
0.25 0.25 0.25
```

(Take care that “Direct” starts in the first column.)

Convert to an *init* file with the following command:

```
poscar2init POSCAR > init.si
```

The output has been redirected to the file *init.si* (it is sent to standard output by default). Take a look at the converted *init* file. **UNITS=A** on the third line specifies that the lattice constant is in Angstrom units, which is the default for the *POSCAR* format. The Questaal code accepts other units in the *init* file but will convert to atomic units. Check this by running *blm* and you will see that the lattice constant in the site file has been converted to Bohr. Last thing to note is the **X=** after the atom species. The *init* file from the tutorial denotes positions by **POS=**, indicating cartesian coordinates but here **X=** indicates “direct” or “fractional” coordinates, in which vectors are defined as (noninteger) multiples of the three lattice vectors (see the following additional exercises for more). The *poscar2init* utility only accepts *POSCAR* files with coordinates specified in direct format. Also note that the *POSCAR* file must specify the atom types and number of atoms per atomic species (lines 6 and 7 in the POSCAR example above). More information, including an example importing from a **CIF** file the detailed input file tutorial.

**4) Converting between fractional and Cartesian coordinates**

There is a simple relation between “fractional” coordinates **X=** and Cartesian coordinates **POS=**, as described in the FAQ.

To compare, try running the command `blm init.si --brief --wsitex`

and you will see that **xpos** has been added to the first line of *site.si*; this indicates that the coordinates are now in fractional form. Note that in this example for silicon the Cartesian and fractional coordinates happen to be the same.

**5) Specifying the structure based on the space group ( init file)**

Run *blm* again but this time use the following *init* file:

```
# Init file for Silicon
LATTICE
SPCGRP=227
A=5.43 UNITS=A
SITE
ATOM=Si X=0 0 0
```

Si forms in the diamond cubic structure, space group 227 as you can readily determine from, e.g. this web page. Space group 227 is Fd-3m which has an underlying fcc Bravais lattice. The *init* file equally allows you to supply the space group as **SPCGRP=Fd-3m**, or supply the lattice vectors directly, with a token **PLAT**, *viz*

```
PLAT= 0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
```

**UNITS=A** tells *blm* that **A** is in Å;. If you do not specify the units, **A** is given in atomic units. (Questaal uses Atomic Rydberg units.)

Lattice vectors, lattice constant and the position of the basis vectors place in **SITE** is all that is needed to fix the structural information. Atomic numbers for each atom is inferred from the symbol (Si). Token **X=** specifies that the coordinates are in “direct” representation, that is, as fractional multiples of lattice vectors. It doesn’t matter in this case, but you can use Cartesian coordinates : use **POS=** instead of **X** (see additional exercise 4). In such a case positions are given in Cartesian coordinates in units of **ALAT**. So are the lattice vectors, if you supply them instead of the space group number.

**6) Draw the energy band structure**

To draw the energy bands, you must select particular lines in the Brillouin zone and put them in a file designed for this purpose, i.e. symmetry-lines mode.

*lmchk* has special mode designed to make files in this format; it is activated by the --syml switch. This page contains common symmetry lines and labels for any crystal structure, Information there was extracted from this paper. In this section you can find the following instruction for the face centered cubic lattice (Si has the fcc lattice, as noted in Exercise 5).

```
lmchk ctrl.si --syml~n=21~mq~wmq~lblq:G=0,0,0,L=1/2,1/2,1/2,X=0,1/2,1/2,W=1/4,1/2,3/4,K=3/8,3/8,3/4~lbl=LGXWGK
```

Labels **L**, **G**, **X**, and so on are widely used names for particular **k** points (**G** is the origin, and is a substitute for the usual label **Γ**).

**--syml** expects definitions of two or more labels and the **k** points they are associated with. These are supplied in the subtag `~lblq:...~`

. **G**, **L**, **X**, **W**, **K** are widely used labels for fcc; the **k** points indicated in the line above are in “direct” or “fractional” coordinates (see FAQ #2 above). `~mq`

tells *lmchk* that the **k** points are in direct coordinates (without this switch it assumes Cartesian coordinates), and `~wmq`

tells *lmchk* to write the symmetry lines file in direct coordinates (by default it writes the points in Cartesian coordinates). `n=21`

specifies that a line of unit length have 21 points; the number of points on each line is 21×length. `~lbl=LGXWGK`

tells the file maker to generate a file with lines **L-Γ**, **Γ-X**, **X-W**, **W-Γ**, and **Γ-K**.

For this exercise, just draw the most important lines, **L-Γ-X**. Do the following

```
lmchk ctrl.si --syml~n=21~mq~lblq:G=0,0,0,L=1/2,1/2,1/2,X=0,1/2,1/2,W=1/4,1/2,3/4,K=3/8,3/8,3/4~lbl=LGX
```

You should see *syml.si* was created with the following

```
# generated from text ~n=21~mq~lblq:G=0,0,0,L=1/2,1/2,1/2,X=0,1/2,1/2,W=1/4,1/2,3/4,K=3/8,3/8,3/4~lbl=LGX
18 0.5000000 0.5000000 0.5000000 0.0000000 0.0000000 0.0000000 L to G
21 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 G to X
```

Run *lmf* in band mode

```
lmf ctrl.si --band~fn=syml
```

*lmf* creates *bnds.si*, which contains the bands in a particular format the *plbnds* utility is designed to read. It create inputs for the *fplot* utility, which is the Questaal analog of *gnuplot*.

Do the following:

```
plbnds -range=-13,8 -fplot~scl=.7,.8~sh~ts=.5 -ef=0 -scl=13.6 -lbl=L,G,X bnds.si
fplot -f plot.plbnds
```

**-scl=13.6**Converts energies from atomic Rydberg units to eV**-range=-13,8**makes the energy bands over the window**-fplot~scl=.7,.8~sh~ts=.5**Tells*plbnds*not to make a postscript file directly, but to create a script for the*fplot*utility along with files containing the energy bands. It also sets the size of the canvas and the spacing between tic marks.**-ef=0**Shifts the energy origin so the Fermi level (which coincides with the valence band maximum here) falls at 0**-lbl=L,G,X**Labels for the points on the symmetry lines

*fplot* creates a postscript file, *fplot.ps*, which you can view with your favorite postscript viewer; e.g. `evince fplot.ps`

. The valence band maximum should at zero; and you should see the minimum point of the lowest conduction band on the **Γ-X** line close to X.