Detailed lmf tutorial: self-consistent LDA calculation for PbTe

This tutorial carries out a self-consistent density-functional 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.

1. Building the input file

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.

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

% 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 k-point 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 re-assign values to variables from the command line, as we will see shortly.

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 command-line 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.

--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 post-processed 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

In order for lmf carry out a self-consistent 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:

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.

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 density-functional 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 self-consistency cycle the potential will change and lmf will update this value.

lmfa automatically selects the valence-core partitioning; the information is given in basp.pbte. You can set the partitioning manually by editing this file.

Note: high-lying 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 Valence-core 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 valence-core 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:

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

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.

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.

The envelope functions (smoothed Hankel functions) are characterized by RSMH and EH. RSMH is the Gaussian “smoothing radius” and approximately demarcates the transition between short-range behavior, where the envelope varies as rlr^lr​l​​, and asymptotic behavior where it decays exponentially with decay length 1/κl=1/−εl1/\kappa_l = 1/\sqrt{-\varepsilon_l}1/κ​l​​=1/√​−ε​l​​​​​, where εl\varepsilon_lε​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αRL(r)H_{\alpha\mathbf{R}L}(\mathbf{r})H​αRL​​(r) [α=1\alpha{=}1α=1 for the (RSMH,EH) group, α=2\alpha{=}2α=2 for the (RSMH2,EH2) group], are augmented by partial waves in augmentation spheres. Thus the lmf basis set χαRL(r)\chi_{\alpha\mathbf{R}L}(\mathbf{r})χ​αRL​​(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.

lmfa searches for core states which are shallow enough to be treated as local orbitals, using the core energy and charge spill-out of the augmentation radius (rmt) as criteria; see annotated lmfa output.

When it was run for the first time, lmfasingled 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.

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 self-consistency 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 free-atomic 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:

6. Remaining Inputs

k mesh

We are almost ready to carry out a self-consistent calculation. Try the following:

lmf ctrl.pbte

lmf stops with this message:

Exit -1 bzmesh: illegal or missing k-mesh

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 k-mesh of 6×6×66\times 6\times 66×6×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 command-line argument-gmax=7.8.

This is sufficient to make the Kohn-Sham hamiltonian. Other matrix elements may be made depending on circumstances, matrix elements for optics or for spin-orbit 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).

Two checks are made: against the change (0.004793) in total energy and the RMS DQ (0.018916). When both checks fall below tolerances self-consistency is reached. In this case it occurs in iteration 10, where the convergence check reads:

The first line prints out the change in Harris-Foulkes energy relative to the prior iteration and some norm of RMS change in the charge density nout−nin (see arrows), followed by the tolerances required for self-consistency.

The last line prints out variables specified on the command line, and total energies from the Harris-Foulkes and Kohn-Sham functionals. Theses are different functionals but they should approach the same value at self-consistency. The c at the beginning of the line indicates that this iteration achieved self-consistency with the tolerances specified. See the annotated output for more details.

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 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, 114-147. H. Dreysse, ed. (Springer-Verlag, Berlin) 2000.

FAQ

How does lmf iterate to self-consistency?

It mixes the input density nin with output density nout generated by lmf, to construct a new input density nin. This process is repeated until nout=nin (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 spin-orbit coupling affect the band structure?

Yes, it does. The bandgap will change significantly when spin-orbit 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 L6+ and L6− 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 high-symmetry lines. But in general the position of band edges can be quite complex. A slightly more complicated case is Si. 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 self-consistent 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.