# 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

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

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

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.

VERS  LM:7 FP:7

HAM                                   # Parameters controlling hamiltonian
GMAX=   5.6                     # PW cutoff for charge density
AUTOBAS[PNU=1 LOC=1 MTO=4]      # Parameters defining basis

BZ                                    # Brillouin zone integration
NKABC= 4 4 4                    # 1 to 3 values

STRUC                                 # Crystal lattice information
FILE=site                       # Extract some or all of NBAS, ALAT, PLAT from file

SPEC                                  # Chemical and basis information of species
ATOM=Si       Z= 14  R= 2.221355  LMX=2 LMXA=3

SITE                                  # Read site data from auxiliary file
FILE=site


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:

TagExplanation
VERS_LMMajor version number for the entire Questaal package
VERS_FPMajor version number for the full potential code lmf
HAM_GMAXG-vector cutoff for the plane wave representation of smoothed Hankel envelope functions and the smooth density
HAM_AUTOBAS_PNUPNU is an $\ell$-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_LOCPZ is an $\ell$-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_MTOTells 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 $\ell$, and AUTOBAS_MTO=4 tells the input file reader to seek for as many as two envelope functions on a site for each $\ell$ (RSMH, EH, and RSMH2, EH2)
BZ_NKABCThe 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_FILETells 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_ATOMA 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_ZThe atomic number associated with SPEC_ATOM
SPEC_ATOM_RThe augmentation radius associated with SPEC_ATOM
SPEC_ATOM_LMXThe $\ell$ cutoff for the smooth-Hankel basis. Here $\ell=2$, which corresponds to a basis of spd orbitals (recall that there are two envelope functions for each $\ell$)
SPEC_ATOM_LMXAThe $\ell$ 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 $\ell$ in the basis. Si’s d orbitals ($\ell=2$) 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_FILETells 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 = 4 ... set nkabc = 4 4 4
Check k mesh for nkabc = 4 4 4 ...
BZMESH: 8 irreducible QP from 64 ( 4 4 4 )  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[=n]) and whether and how --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: 26 0 38 0
...


The table following this line shows each Si has 13 orbitals (9 spd orbitals for the first envelope functions RSMH,EH, 4 sp 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.161015  ec=0.204995
VBmax = 0.161015  CBmin = 0.204995  gap = 0.043980 Ry = 0.59813 eV
BZINTS: Fermi energy:      0.161015;   8.000000 electrons;  D(Ef):    0.000
Sum occ. bands:   -1.6827742  incl. Bloechl correction:    0.000000


The bandgap, 0.598 eV is about 0.01 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 ehfehk. 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.1349−2×−577.6408 = −0.8531 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 $\mathbf{R}({n})=n_{1}\mathbf{R}_{1}+n_{2}\mathbf{R}_{2}+n_{3}\mathbf{R}_{3}$, with $n_{1},\, n_{2},\, n_{3}$ integers.
$\mathbf{R}_{i}$ are primitive lattice vectors (three members of $\mathbf{R}(n)$ which make the smallest nonzero volume). See Wikipedia for more details.
In Questaal’s notation $\mathbf{P}_{i}$ are used for (dimensionless) primitive lattice vectors, with $\mathbf{R}_{i}=a\mathbf{P}_{i}$, $a$ is the lattice constant read as ALAT in the input file (atomic units).
The reciprocal lattice are the family of vectors $\mathbf{G}(m)$ defined from three primitive vectors, $\mathbf{G}({m})=m_{1}\mathbf{G}_{1}+m_{2}\mathbf{G}_{2}+m_{3}\mathbf{G}_{3}$. $\mathbf{G}_i$ is related $\mathbf{R}_j$ by $\mathbf{G}_i\cdot\mathbf{R}_j=2\pi\delta_{ij}$.
Questaal uses $\mathbf{Q}_{i}$ for reciprocal primitive lattice vectors, with $\mathbf{G}_{i}=(2\pi/a)\mathbf{Q}_{i}$.
The simplest way to relate $\mathbf{Q}_i$ and $\mathbf{P}_j$ is to construct 3×3 matrices $P = (\mathbf{P}_{1},\,\mathbf{P}_{2},\,\mathbf{P}_{3})$ and $Q = (\mathbf{Q}_{1},\,\mathbf{Q}_{2},\,\mathbf{Q}_{3})$.
Then $Q^T=P^{-1}$. 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:

POSi = X1 Pi,1  +  X2 Pi,2  +  X3 Pi,3   where Pi,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 $\mathbf{Q}_{1}/n_{k_1},\,\mathbf{Q}_{2}/n_{k_2},\,\mathbf{Q}_{3}/n_{k_3}$.   $n_{k_1},\,n_{k_2},\,n_{k_3}$ 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.

1) Short version of tutorial with one-liner command

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

blm init.si --gmax=6 --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=5.6) it should be essentially identical. The changes in output give a measure of how well converged the default value (5.6) 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  4×4×4  k mesh. The actual minimum occurs near k=(1,0,0), commonly referred to as the X point. The X point is on the  4×4×4  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.ext --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. gs 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.