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

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

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

lmf ctrl.si > out.lmfsc

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.

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.

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.

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

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 $\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:

POS_i = X₁ P_i,1  +  X₂ P_i,2  +  X₃ 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 $\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.

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

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 

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.