Detailed lmf tutorial: selfconsistent LDA calculation for PbTe
This tutorial carries out a selfconsistent densityfunctional calculation for PbTe using the lmf code. It has a purpose similar to the basic tutorial on Si but provides much more detail. See also the Fe tutorial, an LDA+QSGW calculation for a ferromagnetic metal.
It synchronizes with an ASA tutorial on the same system, enabling a comparison of the ASA and full potential methods, and forms a starting point for other tutorials, e.g. on optics.
nano init.pbte
blm init.pbte #makes template actrl.pbte and site.pbte
cp actrl.pbte ctrl.pbte
Free atomic density and basis parameters
lmfa ctrl.pbte #use lmfa to make basp file, atm file and to get gmax
cp basp0.pbte basp.pbte #copy basp0 to recognised basp prefix
lmfa ctrl.pbte #remake atomic density with updated valencecore partitioning
lmf ctrl.pbte vnkabc=6 vgmax=7.8
Table of Contents
 1. Building the input file
 2. How the input file is organized
 3. The EXPRESS category
 4. Determining what input an executable seeks
 5. Initial setup: free atomic density and parameters for basis
 6. Remaining Inputs
 7. Self consistency
 Other Resources
 FAQ
 Additional exercises
Preliminaries
Some of the basics are covered in the basic lmf tutorial for Si. It is easier to read but less detailed. See also the the tutorial on building input files.
The standard outputs from running this tutorial are explained in the annotation of lmfa output and annotation of lmf output. Many details omitted here are given in the annotated outputs.
Executables blm, lmchk, lmfa, and lmf are required and are assumed to be in your path.
1. Building the input file
(See also the tutorial on building input files).
PbTe crystallizes in the rocksalt structure with lattice constant a = 6.428 Å. You need the structural information in the box below to construct the main input file, ctrl.pbte. Start in a fresh working directory and cut and paste the box’s contents to init.pbte.
LATTICE
ALAT=6.427916 UNITS=A
PLAT= 0.0000000 0.5000000 0.5000000
0.5000000 0.0000000 0.5000000
0.5000000 0.5000000 0.0000000
SITE
ATOM=Pb X= 0.0000000 0.0000000 0.0000000
ATOM=Te X= 0.5000000 0.5000000 0.5000000
The primitive lattice vectors are in row format (the first row contains the x, y and z components of the first lattice vector and so forth). In the SITE section, the atom type and coordinates are shown. X= specifies the site coordinates. They are specified in “direct” representation, i.e., as fractional multiples of lattice vectors PLAT. You can also use Cartesian coordinates; instead of X= you would use POS= (see additional exercises below). Positions in Cartesian coordinates are in units of ALAT, like the lattice vectors.
Note : You can also import structural data from other sources. See this page for options.
Use the blm tool as in the box below to create the input file (ctrl.pbte) and the site file (site.pbte):
blm init.pbte
cp actrl.pbte ctrl.pbte
Note : If you are preparing for a later QSGW calculation, use blm gw init.pbte
. See this page for documentation of blm’s commandline switches.
2. How the input file is organized
In this tutorial, blm is used in “standard” mode. (The basic tutorial creates a simpler file with blm express init.si
). Standard mode makes limited use of the preprocessing capabilities of the Questaal input system : it uses algebraic variables which can be modified on the command line. Thus lmf vnit=10 ...
sets variable nit to 10 before doing anything else. Generally:
 Lines which begin with ‘#’ are comment lines and are ignored. (More generally, text following a `#’ in any line is ignored).
 Lines beginning with ‘%’ are directives to the preprocessor.
The beginning of the ctrl file generated by blm should look like the following:
# Variables entering into expressions parsed by input
% const nit=10
% const met=5
% const so=0 nsp=so?2:1
% const lxcf=2 lxcf1=0 lxcf2=0 # for PBE use: lxcf=0 lxcf1=101 lxcf2=130
% const pwmode=0 pwemax=3 # Use pwmode=1 or 11 to add APWs
% const nkabc=0 gmax=0
% const tells the preprocessor that it is declaring one or more variables. nit, met, etc, used in expressions later on. The parser interprets the contents of brackets {…} as algebraic expressions: The contents of {…} is evaluated and the numerical result is substituted for it. Expression substitution works for input lines proper, and also in the directives.
For example this line
metal= {met} # Management of kpoint integration weights in metals
becomes
metal= 5
because met is a numerical expression (admittedly a trivial one). It evaluates to 5 because met is declared as an algebraic variable and assigned value 5 near the top of the ctrl file. The advantage is that you can do algebra in the input file, and you can also reassign values to variables from the command line, as we will see shortly.
Lines corresponding to actual input are divided into categories and tokens within the categories.
A category begins when a character (other than % or #) occurs in the first column. Each token belongs to a category; for example in box below IO contains three tokens, SHOWMEM, IACTIV and VERBOS :
IO SHOWMEM=f
IACTIV=f VERBOS=35,35
(Internally, a complete identifier (aka tag) would be IO_IACTIV=, though it does not appears in that form in the ctrl file.)
This link explains the structure of the input file in more detail.
3. The EXPRESS category
blm normally includes an EXPRESS category in ctrl.pbte.
EXPRESS
# Lattice vectors and site positions
file= site
# Basis set
gmax= {gmax} # PW cutoff for charge density
autobas[pnu=1 loc=1 lmto=5 mto=4 gw=0]
Tags in the EXPRESS category are effectively aliases for tags in other categories, e.g. EXPRESS_gmax corresponds to the same input as HAM_GMAX. If you put a tag into EXPRESS, it will be read there and ignored in its usual location; thus in this instance adding GMAX to the HAM category would have no effect.
The purpose of EXPRESS is to simplify the input file, collecting the most commonly used tags in one place.
4. Determining what input an executable seeks
Executables accept input from two primary streams : tags in the ctrl file and additional information through commandline switches. Each executable reads its own particular set, though most executables share many tags in common.
Usually an input file contains only a small subset of the tags an executable will try to read; defaults are used for the vast majority of tags.
There are four special modes designed to facilitate managing input files. For definiteness consider the executable lmfa.
$ lmfa input
$ lmfa help
$ lmfa showp
$ lmfa show  lmfa show=2
input
puts lmfa in a special mode. It doesn’t attempt to read anything; instead, it prints out a (large) table of all the tags it would try to read, including a brief description of the tag, and then exits.
See here for further description.
help
performs a similar function for the command line arguments: it prints out a brief summary of arguments effective in the executable you are using.
See annotated lmfa output for further description.
showp
reads the input through the preprocessor, prints out the preprocessed file, and exits.
See the annotated lmf output for a comparison of the pre and postprocessed forms of the input file in this tutorial.
show
tells lmfa to print out tags as it reads them (or the defaults it uses).
It is explained in the annotated lmf output.
5. Initial setup: free atomic density and parameters for basis
 The basp file
 In order for lmf carry out a selfconsistent calculation for the crystal, it requires information about the basis set, and an initial trial density.
lmfa computes densities for free atoms, which are overlapped to form an initial trial density (called the Mattheis construction).
lmfa will also automatically generate a basis set, and save it in a file called basp.ext. What is stored in this file depends on the contents of HAM_AUTOBAS, but typically three pieces of information are kept:
 Parameters PZ specifying local orbitals
 Parameters (RSMH,EH) defining the shape of the smooth Hankel envelope functions
 Parameters P defining the boundary conditions of partial waves on augmentation boundaries
As will be described in more detail in the pages below, lmfa prepares the following.

Make a selfconsistent atomic density for each species.

Fit the density outside the augmentation radius. lmf needs this information to overlap atomic densities for an initial trial density.
Information about the augmented and interstitial parts of the density are written to file atm.pbte. 
Provide a reasonable estimate for the Gaussian smoothing radius r_{s} and hankel energy ε) that fix the shape of the smooth Hankel envelope functions for l=0, 1,…. The l cutoff is determined internally, depending on the setting of HAM_AUTOBAS_LMTO.
These parameters are written to file basp0.pbte as RSMH and EH. 
Provide a reasonable estimate for boundary conditions that fix linearization energies, parameterized by the logarithmic derivative parameter P_{l}, aka the “continuous principal quantum number.”
These parameters are written to basp0.pbte as P. 
Decide on which shallow cores should be included as local orbitals.
Local orbitals are written basp0.pbte as nonzero values of PZ. 
Supply an estimate for the interstitial density plane wave cutoff GMAX.
lmfa will provide all of this information automatically. It will write atomic density information to atm.pbte and basis set information to template basp0.pbte. The Questaal suite reads from basp.pbte, but lmfa writes to basp0 to avoid overwriting a file you may want to preserve. You can edit basp.pbte and customize the basis set.
As a first step, do:
lmfa ctrl.pbte #use lmfa to make basp file, atm file and to get gmax
cp basp0.pbte basp.pbte #copy basp0 to recognised basp prefix
The output is annotated in some detail here. It begins with a header:
LMFA: nbas = 2 nspec = 2 vn 7.11.i verb 35
special
pot: XC:BH
autogen: mto basis(4), pz(1), pnu(1) Autoread: pz(1)
The pot line says that lmfa the potential will be made from the BarthHedin functional. To use a GGA, see here.
The autogen line says that lmfa will make the basis set information (points 35 outlined above).
The next few sections amplify on these three points. Point 6 is discussed here.
5.1 Local orbitals
Part of lmfa’s function is to identify local orbitals that extend the linear method. Linear methods are reliable only over a limited energy window; certain elements may require an extension to the linear approximation for accurate calculations. This is accomplished with local orbitals. lmfa will automatically look for atomic levels which, if certain criteria are satisfied it designates as a local orbital, and includes this information in the basp0 file. The annotated lmfa output explains how lmfa analyzes core states for local orbitals.
Inspect basp.pbte. Note in particular this text connected with the Pb atom:
PZ= 0 0 15.934
(The same information can be supplied in the input file, through SPEC_ATOM_PZ.)
lmfa is suggesting that the Pb 5d state is shallow enough that it be included in the valence. Since this state is far removed from the Fermi level, we would badly cover the Hilbert space spanned by Pb 6d state were we to use Pb 5d as the valence partial wave. (In a linear method you are allowed to choose a single energy to construct the partial wave; it is usually the “valence” state, which is near the Fermi level.)
This problem is resolved with local orbitals : these are partials wave at an energy far removed from the Fermi level. The three numbers following PZ correspond to specifications for local orbitals in the s, p, and d channels. Zero indicates “no local orbital;” there is only a d orbital here.
15.934 is actually a compound of 10 and the “continuous principal quantum number” 5.934. The 10’s digit tells lmf to use an “enhanced” local orbital as opposed to the usual variety found in most densityfunctional codes. Enhanced orbitals append a tail so that the density from the orbital spills into the interstitial. You can specify a “traditional” local orbital by omitting the 10, but this kind is more accurate, and there is no advantage to doing so.
The continuous principal quantum number (5.934) specifies the number of nodes and boundary condition. The large fractional part of P is large for core states, typically around 0.93 for shallow cores. lmfa determines the proper value for the atomic potential. In the selfconsistency cycle the potential will change and lmf will update this value.
lmfa automatically selects the valencecore partitioning; the information is given in basp.pbte. You can set the partitioning manually by editing this file.
Note: highlying states can also be included as local orbitals; they improve on the Hilbert space far above the Fermi level. In the LDA they are rarely needed and lmfa will not add them to the basp.pbte. But they can sometimes be important in GW calculations, since in contrast to the LDA, unoccupied states also contribute to the potential.
5.2 Valencecore partitioning of the free atomic density
After basp.pbte has been modified, you must run lmfa a second time:
lmfa ctrl.pbte
This is necessary whenever the valencecore partitioning changes through the addition or removal of a local orbital. Even though lmfa writes the atomic to atm.pbte, this file will change when partitioning between core and valence will change with the introduction of local orbitals, as described next. This is because core and valence densities are kept separately.
Relativistic core levels
Normally lmfa determines the core levels and core density from the scalar Dirac equation. However there is an option to compute the core levels from the full Dirac equation.
Tag HAM_REL controls how the Questaal package manages different levels of relativistic treatment. Run lmfa input
and look for HAM_REL. You should see:
HAM_REL opt i4 1, 1 default = 1
0 for nonrelativistic Schrödinger equation
1 for scalar relativistic Schrödinger equation
2 for Dirac equation (ASA only for now)
10s digit 1: compute core density with full Dirac equation
10s digit 2: Like 1, but neglect coupling (1,2) pairs in 4vector
Set HAM_REL=11 to make lmfa calculate the core levels and core density with the full Dirac equation.
You might want to see the core level eigenvalues; they can shift significantly relative to the scalar Dirac solution. Also, l is no longer a good quantum number so there can be multiple eigenvalues connected with the scalar Dirac l. To see these levels, invoke lmfa with a sufficiently high verbosity. In the present instance insert HAM REL=11 into ctrl.pbte and do
lmfa pr41 ctrl.pbte
You should see the following table:
Dirac core levels:
nl chg <ecore(S)> <ecore(D)> <Tcore(S)> <Tcore(D)> nre
1s 2 6461.412521 6461.420614 9160.575645 9160.568216 439
ec(mu) 6461.420614 6461.420614
2s 2 1154.772794 1154.777392 2201.484620 2201.485036 473
ec(mu) 1154.777392 1154.777392
3s 2 277.137428 277.136313 700.148783 700.160432 501
ec(mu) 277.136313 277.136313
4s 2 62.683976 62.678557 231.671152 231.686270 531
ec(mu) 62.678557 62.678557
5s 2 10.589828 10.580503 60.826909 60.833608 567
ec(mu) 10.580503 10.580503
2p 6 990.094400 1001.984462 1702.510726 1772.365432 475
ec(mu) 948.389636 1109.174115 948.389636 1109.174115 948.389636 948.389636
3p 6 229.993746 232.623198 568.649082 585.156080 505
ec(mu) 220.667558 256.534478 220.667558 256.534478 220.667558 220.667558
4p 6 47.246014 47.902771 184.751871 189.523363 537
ec(mu) 44.969950 53.768412 44.969950 53.768412 44.969950 44.969950
5p 6 6.300710 6.422904 43.507054 44.670581 577
ec(mu) 5.869706 7.529300 5.869706 7.529300 5.869706 5.869706
3d 10 182.032939 182.146340 501.452676 502.171493 509
ec(mu) 179.091564 186.728504 179.091564 186.728504 179.091564 186.728504 179.091564 186.728504 179.091564 179.091564
4d 10 29.432703 29.453418 150.979227 151.198976 545
ec(mu) 28.796634 30.438595 28.796634 30.438595 28.796634 30.438595 28.796634 30.438595 28.796634 28.796634
5d 10 1.566638 1.562069 23.907636 23.945913 605
ec(mu) 1.485638 1.676716 1.485638 1.676716 1.485638 1.676716 1.485638 1.676716 1.485638 1.485638
4f 14 9.755569 9.751307 117.412788 117.457023 569
ec(mu) 9.592725 9.962749 9.592725 9.962749 9.592725 9.962749 9.592725 9.962749 9.592725 9.962749 9.592725 9.962749 9.592725 9.592725
qcore(SR) 78.000000 qcore(FR) 78.000000 rho(rmax) 0.00000
sum ec : 25841.9031 (SR) 25934.9233 (FR) diff 93.0203
sum tc : 48113.1010 (SR) 48677.3220 (FR) diff 564.2210
The scalar Dirac Pb 5d eigenvalue (1.566638 Ry) gets split into 6 levels with energy 1.485638 Ry and four with 1.676716 Ry. The mean (1.56207 Ry) is close to the scalar Dirac value. In the absence of a magnetic field a particular l will split into two distinct levels with degeneracies 2l and 2l+2, respectively.
The bottom part of the table shows how much the free atom’s total energy changes as a consequence of the fully relativistic Dirac treatment.
5.3 Automatic determination of basis set
Some details of the basis set (envelope functions, augmentation, local orbitals) are explained in this tutorial.
lmfa loops over each species, generating a selfconsistent density.
Given a density and corresponding potential, lmfa will construct some estimates for the basis set, namely the generation of envelope function parameters RSMH and EH (and possibly RSMH2 and EH2, depending on the setting of HAM_AUTOBAS_MTO), analyzing which cores should be promoted to local orbitals, and reasonable estimates for the boundary condition of the partial wave.
 Envelope functions
 The envelope functions (smoothed Hankel functions) are characterized by RSMH and EH. RSMH is the Gaussian “smoothing radius” and approximately demarcates the transition between shortrange behavior, where the envelope varies as $r^l$, and asymptotic behavior where it decays exponentially with decay length $1/\kappa_l = 1/\sqrt{\varepsilon_l}$, where $\varepsilon_l$ is one of the EH. lmfa finds an estimate for RSMH and EH by fitting them to the “interstitial” part of the atomic wave functions (the region outside the augmentation radius).
 Fitting the smooth Hankel function to the numerically tabulated exact function is usually quite accurate. For Pb, the error in the energy (estimated from the single particle sum) is 0.00116 Ry — very small on the scale of other errors. The fitting process is described in more detail in the annotated lmfa output.
 lmf requires RSMH and EH. Those generated by lmfa are reasonable, but unfortunately not optimal choices for the crystal, as explained in the annotated lmfa output. You can change them by hand, or optimize them with lmf’s optimizing function,
opt
. To make an accurate basis, a second envelope function is added through RSMH2 and EH2. (lmfa automatically does this, depending on the setting of HAM_AUTOBAS_MTO). Alternatively you can add APW’s to the basis. For a detailed discussion on how to select the basis set, see this tutorial.
 Note 1: The envelope functions $H_{\alpha\mathbf{R}L}(\mathbf{r})$ [$\alpha{=}1$ for the (RSMH,EH) group, $\alpha{=}2$ for the (RSMH2,EH2) group], are augmented by partial waves in augmentation spheres. Thus the lmf basis set $\chi_{\alpha\mathbf{R}L}(\mathbf{r})$ consists of augmented smooth Hankel functions.
 Note: 2: The new Jigsaw Puzzle Orbital basis is expected significantly improve on the accuracy of the existing Questaal basis. High quality envelope functions are automatically constructed that continuously extrapolate the accurate augmented partial waves smoothly into the interstitial; the kinetic energy of the envelope functions are continuous across the augmentation boundary.
 Local orbitals
 lmfa searches for core states which are shallow enough to be treated as local orbitals, using the core energy and charge spillout of the augmentation radius (rmt) as criteria; see annotated lmfa output.
 When it was run for the first time, lmfa singled out the Pb 5d state, using information from the table below taken from lmfa’s standard output. Once local orbitals are specified lmfa is able to appropriately partition the valence and core densities. This is essential because the two densities are treated differently in the crystal code. Refer to the annotated lmfa output for more details.
Find local orbitals which satisfy E > 2 Ry or q(r>rmt) > 5e3
l=2 eval=1.569 Q(r>rmt)=0.0078 PZ=5.934 Use: PZ=15.934
l=3 eval=9.796 Q(r>rmt)=3e8 PZ=4.971 Use: PZ=0.000
 Boundary conditions
 The free atomic wave function satisfies the boundary condition that the wave function decay as r→∞. Thus, the value and slope of this function at rmt are determined by the asymptotic boundary condition. This boundary condition is needed for fixing the linearization energy of the partial waves in the crystal code. lmfa generates an estimate for this energy and encapsulates it into the “continuous principal quantum number”, saved as P in basp0.pbte (normally P will updated in the selfconsistency cycle).
Refer to the annotated lmfa output for more details.
5.4 Fitting the interstitial density
lmfa fits valence and core densities to a linear combination of smooth Hankel functions. This information will be used to overlap freeatomic densities to obtain a trial starting density. This is explained in the annotated lmfa output.
5.5 Estimate for GMAX
After looping over all species lmfa writes basis information to basp0.pbte, atomic charge density data to file atm.pbte, and exits with the following printout:
FREEAT: estimate HAM_GMAX from RSMH: GMAX=4.3 (valence) 7.8 (local orbitals)
This is the G cutoff EXPRESS_gmax or HAM_GMAX needed to determine the interstitial mesh spacing. Two values are printed, one determined from the shape of valence envelope functions (4.3) and, if local orbitals are present the largest value found from their shape, as explained in the annotated lmfa output.
6. Remaining Inputs
k mesh
We are almost ready to carry out a selfconsistent calculation. Try the following:
lmf ctrl.pbte
lmf stops with this message:
Exit 1 bzmesh: illegal or missing kmesh
We haven’t yet specified a k mesh. You must supply it yourself since there are too many contexts to supply a sensible default value. In this case a kmesh of $6\times 6\times 6$ divisions is adequate. With your text editor change nkabc=0 in the ctrl file to nkabc=6, or alternatively assign variable nkabc on the command line using vnkabc=6 (which is what this tutorial will do).
Charge density mesh
We also haven’t specified the G cutoff for the density mesh. blm does not determine this parameter automatically because it is sensitive to the selection of basis parameters, which local orbitals are included. lmfa conveniently supplies that information for us, based in the shape of envelope functions it found. In this case the valence G cutoff is quite small (4.3), but the Pb 5d local orbital is a much sharper function, and requires a larger cutoff (7.8). You must use the larger of the two.
Note: if you change the shape of the envelope functions you must take care that gmax is large enough. This is described in the lmf output below.
Change variable gmax=0 in the ctrl file, or alternatively add a variable to the command line (vgmax=7.8), as we do in the next section. Or, you can run blm again, with commandline argument gmax=7.8.
7. Self consistency
Carry out a selfconsistent calculation as follows:
lmf ctrl.pbte vnkabc=6 vgmax=7.8
lmf will iterate up to 10 iterations. The cycle is capped to 10 iterations because of the following lines in ctrl.pbte, which before and after preprocessing read:
before preprocessing after preprocessing
% const nit=10 
... 
EXPRESS  EXPRESS
... 
nit= {nit}  nit= 10
Initialization steps
lmf begins with some initialization steps. Each step is explained in more detail in the annotated lmf output.
 Read basis set parameters from basp.pbte. This information can also be given via the ctrl file, depending on the settings of EXPRESS_autobas.
 Informational printout about computing conditions, lattice structure, and atomic parameters such as augmentation radii and lcutoffs
 Automatic determination of crystal symmetry
 Setup for the Brillouin zone integration
 Construction of the mesh for interstitial density and potential
 Assemble and display information about the size and constituents of the basis set
 Read or assemble an input density
 A QSGW potential Σ^{0} may be read in.
Selfconsistent cycle
Each iteration of the selfconsistency cycle begins with
 BNDFP: begin iteration 1 of 10 
...
 BNDFP: begin iteration 2 of 10 
...
One iteration consists of the following steps. The standard output is annotated in some detail here.
 Construct the potential and matrix elements.
 Interstitial and local parts of the potential are made.
 Partial waves $\phi$ and $\dot{\phi}$ are integrated from the potential subject to the boundary conditions.
 Matrix elements of the partial waves (kinetic energy, potential energy, overlap) are assembled for the KohnSham hamiltonian.
 Matrix elements of the interstitial potential $\langle\chi_iV\chi_j\rangle$. for envelope functions $\chi_i$.
This is sufficient to make the KohnSham hamiltonian. Other matrix elements may be made depending on circumstances, matrix elements for optics or for spinorbit coupling (HAM_SO=t).
 Makes an initial pass through the irreducible k points in the Brillouin zone to obtain the Fermi level and obtain integration weights for each band and k point into a binary file wkp.pbte. In general until the Fermi level is known, the weights assigned to each eigenfunction are not known, so the charge density cannot be assembled. How labor is divided between the first and second pass depends on BZ_METAL. See here for further discussion.
 Makes a second pass to accumulate the output mesh and local densities. For the latter essential information is retained as coefficients of the local density matrix (a compact form).
 Assembles the local densities from the density matrix.
 Symmetrizes the density.
 Finds new logarithmic derivative parameters pnu by floating them to the band centerofgravity
 Computes the HarrisFoulkes and KohnSham total energies and forces.

Mixes the input and output densities to form a new trial density. A segment of the output reads:
mixrho: sought 2 iter from file mixm; read 0. RMS DQ=2.17e2
DQ=2.17e2 is a measure of the root mean square deviation n^{out}−n^{in}. At selfconsistency this number should be small.
 Checks for convergence.
lmf should converge to selfconsistency in 10 iterations.
At the end of the selfconsistency cycle the density is written to rst.pbte
iors : write restart file (binary, mesh density)
and a check is made for convergence. No check is made in the first iteration because there is no prior iteration to compare the change in total energy. The second iteration reads:
↓ ↓ diffe(q)= 0.004793 (0.018916) tol= 0.000010 (0.000030) more=T i nkabc=6 gmax=7.8 ehf=55318.1657745 ehk=55318.1568616
Two checks are made: against the change (0.004793) in total energy and the RMS DQ (0.018916). When both checks fall below tolerances selfconsistency is reached. In this case it occurs in iteration 10, where the convergence check reads:
↓ ↓ diffe(q)= 0.000000 (0.000005) tol= 0.000010 (0.000030) more=F c nkabc=6 gmax=7.8 ehf=55318.1620974 ehk=55318.1620958
The first line prints out the change in HarrisFoulkes energy relative to the prior iteration and some norm of RMS change in the charge density n^{out}−n^{in} (see arrows), followed by the tolerances required for selfconsistency.
The last line prints out variables specified on the command line, and total energies from the HarrisFoulkes and KohnSham functionals. Theses are different functionals but they should approach the same value at selfconsistency. The c at the beginning of the line indicates that this iteration achieved selfconsistency with the tolerances specified. See the annotated output for more details.
Other Resources

Click here to see annotated standard output from lmfa, and here to see annotated standard output from lmf.

An input file’s structure, and features of the programming language capability, is explained in some detail here. The full syntax of categories and tokens can be found in the input file manual.

This tutorial more fully describes some important tags the lmf reads, and this one presents alternative ways to build input files from various sources such as the VASP POSCAR file.

This tutorial more fully explains the lmf basis set. There is a corresponding tutorial on the basics of a selfconsistent ASA calculation for PbTe. A tutorial on optics can be gone through after you have finished the present one.

This document gives an overview of some of lmf’s unique features and capabilities.

The theoretical formalism behind the lmf is described in detail in this book chapter: M. Methfessel, M. van Schilfgaarde, and R. A. Casali, ``A fullpotential 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, 114147. H. Dreysse, ed. (SpringerVerlag, Berlin) 2000.
FAQ

How does lmf iterate to selfconsistency?
It mixes the input density n^{in} with output density n^{out} generated by lmf, to construct a new input density n^{in}. This process is repeated until n^{out}=n^{in} (within a specified tolerance). The actual mixing algorithm can be quite involved; see this page.

The gap is small and Pb is a heavy element. Doesn’t spinorbit coupling affect the band structure?
Yes, it does. The bandgap will change significantly when spinorbit coupling is added.

The LDA is supposed to underestimate bandgaps. But the PbTe bandgap looks pretty good. Why is that?
This turns out to be largely an accident. If spin orbit coupling is included, the bandgap appears to be pretty good, but in fact levels L_{6}^{+} and L_{6}^{−} that form the valence and conduction band edges are inverted in the LDA. See Table I of this paper. As the paper notes, they are well described in QSGW.

How do you know where the band edges are?
PbTe is has a quite simple band structure with high symmetry. It’s a good bet that the band edges are on highsymmetry lines. But in general the position of band edges can be quite complex. A slightly more complicated case is Si. See this tutorial.

Is there an easy way to calculate effective masses?
Yes, once you know where the band edge is. See this tutorial.

The augmentation l cutoff in the ctrl file is 4 for Pb and 3 for Te. LAPW methods typically require much higher l cutoffs (6 or 8) to be well converged. Why is this? Is it a worry?
Questaal uses a different kind of augmentation. It a difference in the basis set (smoothed Hankels vs LAPWs; Questaal in fact has an LAPW basis that may be combined with the augmented smooth Hankels), but because of the different way augmentation is carried out. Both kinds of augmentation converge to the same answer in the limit of large l, but Questaal’s augmentation is much more rapidly convergent. It can be incorporated into APW basis sets as well. See this tutorial.
Additional exercises

Try selfconsistent calculations with the Pb 5d in the valence as a local orbital. Repeat the calculation but remove the PZ part from basp.pbte.

Specify symops manually.

Turn on spin orbit coupling and observe how the band structure changes.

Try rotations

kconvergence. Try BZ_BZJOB //: {::comment}
Edit This Page