This tutorial covers the basics for running a self-consistent DFT calculation for silicon. The goal is to introduce you to the different file types and the basics of running the code.
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 nano ctrl.si #set k-point mesh and real space mesh 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
Executables blm, lmfa, and lmf are required and are assumed to be in your path. 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 that contains structural information in a format that is recognised by the code (analogous to the POSCAR file in VASP). Create a file called 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 specficies the atom types and their positions in cartesian coordinates (indicated by POS) in units of the lattice constant. 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 for input file and “site” for structure file). blm writes 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 you need. gmax specifies how fine a real space mesh is used for 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 also calculates free atom wavefunctions and densities. The latter will be used later to make a trial density for the crystal calculation; the free atom wavefunctions are used to fit the basis set parameters. Note that in the all-electron approach, space is partitioned 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). There are a few more subtleties to the code’s basis set but we will leave further details to the theory pages and the input file pages. Run the following command:
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.
Check the contents of your working directory and you will find two new files atm.si and basp0.si . The file atm.si contains the free atom densities calculated by lmfa. 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\times4\times4$ k mesh is sufficient for Si. Set this value with your text editor by simply changing “nkabc=NULL” to “nkabc=4” (4 is automatically used for each mesh dimension, you could equivlaently use “nkabc=4 4 4”).
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 Brillouin zone integration 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.185518; 8.000000 electrons; D(Ef): 0.000 Sum occ. bands: -1.4864297 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 this book chapter.
Below is a list of frequently asked questions. Please get in contact if you have other questions.
- 1) How does blm determine the augmentation sphere radii?
- Overlaps free atom densities and adjusts the radii to make the potential as similar as possible on teh surface of each sphere.
- 2) What is the log file?
- The log file log.si keeps a compact record of key outputs in the current directory. In successive runs, data is appended to the log file.
- 3) What is the Harris-Foulkes energy?
- It is a functional of the input density, rather than the output density. At self-consistency it should be the same as the standard Kohn-Sham functional. The Harris-Foulkes functional tends to be more stable, and like the Kohn-Sham functional, it is stationary at the self-consistent density. But it is not necessarily a minimum there. See this paper by M. Foulkes and R. Haydock.
- 4) What is the relation between Cartesian coordinates and coordinates expressed as fractions of lattice vectors?
- The connection between positions expressed in Cartesian coordinates POS and fractional multiples of the lattice vectors X is:POSi = X1Pi,1 + X2 Pi,2 + X3 Pi,3 where Pi,j is Cartesian component i of lattice vector j.
The same expression applies in reciprocal if the primitive reciprocal lattice vectors are subsituted for P.
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\times4\times4$ 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\times4\times4$ 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-point 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.