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.
mkdir si && cd si #create working directory and move into it nano init.si #create init file using lines from box below blm init.si --express #use blm tool to create actrl and site files cp actrl.si ctrl.si #copy actrl to recognised ctrl prefix lmfa ctrl.si #get basp file, atm file and gmax estimate cp basp0.si basp.si #copy basp0 to recognised basp prefix sed -i.bak 's/gmax= *NULL/gmax=5.0/' ctrl.si #assign value for gmax sed -i.bak 's/nkabc= *NULL/nkabc=4/' ctrl.si #assign value for k-point messh lmf ctrl.si > out.lmfsc #make self-consistent
Alternatively, you can run the following one-liner to get the same result. This assumes you have already created the init file. See Additional Exercises for more details.
blm si --express --nk=4 --gmax=5 && cp actrl.si ctrl.si && lmfa si && cp basp0.si basp.si && lmf si > out.lmfsc
This tutorial assumes you have downloaded the repository (see on 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. Brief installation instructions are given on this page (which assumes your build directory is ~/lmbuild and you want to install them in ~/bin).
The source code for all Questaal executables can be found here.
To begin with, create a new working directory and move into it. Here we will call it “si”.
mkdir si && cd si
We will start by creating a file called init.si. init files in Questaal contains structural information in a format that is recognised by the input file maker, the blm code. blm will translate the init file an input file ctrl.si (analogous to the POSCAR file in VASP).
Create a file named init.si and copy the following lines to it:
LATTICE ALAT=10.26 PLAT= 0.00000000 0.50000000 0.50000000 0.50000000 0.00000000 0.50000000 0.50000000 0.50000000 0.00000000 # pos means cartesian coordinates, units of alat SITE ATOM=Si POS= 0.00000000 0.00000000 0.00000000 ATOM=Si POS= 0.25000000 0.25000000 0.25000000
ALAT specifices the lattice constant, which is in atomic units by default (Questaal code 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 information can be found in the Additional Exercises.
In order to run the DFT program lmf, you need an input file along with structural information. The blm utility reads the init file as input and creates a template input file actrl.si and structure file site.si. The Questaal package recognises certain names as file types such as
- ctrl the main input file
- site a file containing structural information, which can (but need not) be split off from the ctrl file
and many more, which depend on the context. blm writes to actrl.si rather than ctrl.si to prevent overwriting of an existing ctrl file. Extension “si” labels the material; it is something you choose. Most files read or generated by Questaal programs will append this extension.
Run blm as shown below and then copy the template file actrl.si to ctrl.si, which is the name of the main input file recognised by most codes in the Questaal package. The
--express switch tells blm to make a particularly simple input file; we will see more complicated examples in later tutorials.
blm init.si --express cp actrl.si ctrl.si
The start of the blm output shows some structural and symmetry information. Further down, the “makrm0:” part gives information about creating the augmentation spheres, both silicon atoms were assigned spheres of radii 2.22 Bohr. Now open up the site file and you can see it contains the lattice constant and lattice vectors in the first line. The other terms in the first line are just standard settings and a full explanation can be found in the online page for the site file. The second line is a comment line and the subsequent lines contain the atomic species labels and coordinates. Note that running blm produces a new actrl.si and site.si file each time.
Next take a look at the input file ctrl.si.
The first few lines are just header information, then you have a number of basic parameters for a calculation. We won’t talk about these values now but a full description is provided on the ctrl file page. Defaults are provided by blm for most of the variables except gmax and nkabc, which are left as “NULL”. nkabc specifies the k mesh and there is no sensible way for blm to select a default value, as it depends on the bandgap or details of the Fermi surface, and also the precision needed for the physical property of interest. gmax specifies how dense a space mesh is used to represent the interstitial charge density and this depends on the basis set as described below. Take a look at the last line, it contains information about the different atoms in the system (here we only have silicon) and their associated augmentation spheres.
We now need a basis set and based on this we can get an estimate for gmax. This can done automatically by using the lmfa tool. lmfa has a dual purpose: first, it calculates free atom wavefunctions and densities. The latter will be used later to make a trial density to start the crystal calculation (a superposition of atomic densities is called the Mattheis construction).
Second, it provides estimates for gmax and basis set parameters. Questaal uses an all-electron approach, which partitions space into two kinds of regions: an augmentation sphere part (around each atom) and an interstitial part (region between spheres). In the augmentation spheres, the basis set is composed of local atomic functions, tabulated numerically. In the interstitial regions, the basis consists of analytic, smooth-Hankel functions, characterised by their smoothing radius (RSMH) and Hankel energies (EH). The number of interstitial functions and their parameters (RSMH and EH) are defined in the basp file (the size of which is determined by parameters in the ctrl file).
Again the output shows some structural information and then details about finding the free atom density, basis set fitting and, at the end an estimate of gmax is printed out. Note that the Barth-Hedin exchange-correlation functional is used, as indicated by “XC:BH”, this was specified by “xcfun=2” in the ctrl file (the default). We won’t go into more detail now, but a full description can be found on the lmfa page. One thing to note is a recommended value for gmax given towards the end: “GMAX=5.0”. Now that we have a gmax value, open up the ctrl file and change the default NULL value to 5.0.
Or simply use the in-line sed editor:
sed -i.bak 's/gmax= *NULL/gmax=5.0/' ctrl.si
Check the contents of your working directory and you will find two new files atm.si and basp0.si. File atm.si contains the free atom densities calculated by lmfa, that will be used in the Mattheis construction. File basp0.si is the template basis set file; the standard basis set name is basp and the extra 0 is appended to avoid overwriting. Take a look at basp0.si and you will see that it contains basis set parameters that define silicon’s smooth Hankel functions. Changing these values would change their functional form, but lmfa does a reasonable job (also later on parameters can be automatically optimized, if desired) so we will leave them as they are. Copy basp0.si to basp.si, which is the name lmf recognises for the basis set file.
cp basp0.si basp.si
The second unknown parameter is the k-mesh. A 4×4×4 k mesh is sufficient for Si. Set this value with your text editor by simply changing “nkabc=NULL” to “nkabc=4” (4 is automatically extended to the second and third mesh dimensions, you could equivalently use “nkabc=4 4 4”).
Or simply use the in-line sed editor:
sed -i.bak 's/nkabc= *NULL/nkabc=4 4 4/' ctrl.si
We now have everything we need to run an all-electron, full-potential DFT calculation. This is done using the lmf program. Double check that you have specified the k mesh (nkabc) and a gmax value and then run the following command. 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
The first part of the output is similar to what we’ve seen from the other programs. 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. Towards the end of the output, 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 20”. At the end of each iteration the ehf and ehk total energies are printed and a check is made for self-consistency. The two parameters conv and convc in the ctrl file specify, respectively, are the self-consistency tolerances for the total energy and root mean square (RMS) change in the density. 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.
Further down the Fermi energy and band gap values, and other key bits of information are reported in the Brillouin zone integration section. You should find something similar to the output snippet below.
BZWTS : --- Tetrahedron Integration --- ... only filled or empty bands encountered: ev=0.185509 ec=0.229539 VBmax = 0.185509 CBmin = 0.229539 gap = 0.044029 Ry = 0.59880 eV BZINTS: Fermi energy: 0.185509; 8.000000 electrons; D(Ef): 0.000 Sum occ. bands: -1.4864280 incl. Bloechl correction: 0.000000
To see how the density 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).
grep 'DQ' out.lmfsc grep 'ehk=-' out.lmfsc
You can also check how the bandgap changes as iterations proceed to self-consistency by grepping out.lmfsc for “gap”.
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 to look at are those to generate energy band structures, and density-of-states, or calculate a mechanical property such as the optical mode frequency.
This tutorial on PbTe also covers the basic self-consistency cycle, in a bit 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.
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; the formalism behind the method is described in Comp. Phys. Comm. 249, 107065 (2020).
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:
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 there. 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 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.
1) Short version of tutorial with one-liner command
The entire tutorial can be run with the following one-liner:
blm si --express --nk=4 --gmax=5 && cp actrl.si ctrl.si && lmfa si && cp basp0.si basp.si && lmf 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. First of all, notice that here we use
blm si instead of
blm ctrl.si. blm accepts either the full init file name or just the extension. Next, try running the above command without the --express switch. Compare this new ctrl file to the original one and you will see that the file is different, and provides more options. Note also that values for nkabc and gmax were supplied as switches (--) to blm. This removes the need to edit the ctrl file by hand.
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
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 of cif conversion, can be found in 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 --express --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.