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.

Table of Contents

• 1. Building the input file
• 2. How the ctrl file is organized
• 3. The EXPRESS category
• 4. Determining what input an executable seeks
• 5. Initial setup: free atomic density and parameters for basis
  • 5.1 Local orbitals
    • Relativistic core levels
  • 5.3 Automatic determination of envelope function parameters
  • 5.4 Fitting the interstitial density
  • 5.5 Estimate for GMAX
• 6. Remaining Inputs
  • k mesh
  • Charge density mesh
• 7. Sanity Checks
• 8. Autogenerate basis parameters using blm
• 9. Self consistency
  • Initialization steps
  • Self-consistent cycle
  • Self-consistency in the LDA
• Other Resources
• FAQ
• Additional exercises

cat << END > 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
END
blm --nk=6 --gmax=8.2 --ctrl=ctrl init.pbte

lmfa ctrl.pbte
cp basp0.pbte basp.pbte
lmfa ctrl.pbte

lmfa ctrl.pbte --usebasp

lmf ctrl.pbte -vnkabc=6 -vgmax=8.2

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 for lmf and the general tutorial on generating 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.

A corresponding tutorial for PbTe in the ASA approximation can be found here.

Executables blm, lmchk, lmfa, and lmf are required and are assumed to be in your path.

Note: This tutorial provides many specifics which are useful only if you want to understand the logic of lmf’s program flow in detail. After creating the init file (Section 1), you can skip directly to Section 8.

1. Building the input file

The input file is constructed with the input file maker, blm, documented here. See also the general tutorial on building input files and the one specific to lmf.

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

(Take care that the LATTICE and SITE tags are in the leftmost column of the file.)

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 a tom type and coordinates are shown. X= specifies the site coordinates. They are specified as fractional multiples of lattice vectors PLAT (sometimes called “direct” representation). 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

Or just simply

blm init.pbte --ctrl=ctrl

The output of blm is explained in the introductory tutorial.

Note : If you are preparing for a later QSGW calculation, use blm --gw init.pbte. For quick reference to the command-line switches blm recognizes, try blm --h. This page offers more complete documentation.

2. How the ctrl file is organized

In this tutorial, blm uses a “standard” mode to create the input file ctrl.pbte. (The basic tutorial creates a simpler form using blm --simple 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 in the input file:

• which begin with  #  are comment lines and are ignored. (More generally, text following a  #  in any line is ignored).
• beginning with  %  are directives to the preprocessor. They are not part of the input file after it has been preprocessed.

The beginning of the ctrl file generated by blm should look like the following:

% const nit=10 conv=1e-5 convc=3e-5   # convergence parameters
% const beta=0.3                      # charge mixing parameters
%cconst so nsp=2
% const so=0 nsp=so?2:1
#%const lxcf=0 lxcf1=101 lxcf2=130    # PBE functional from libxc
% const lxcf=2 lxcf1=0 lxcf2=0        # built-in B-H functional.  Use lxcf=1 for Ceperly-Alder
% const pwmode=0 pwemax=3             # to add APWs, use pwmode=1 or 11
% const met=5 nkabc=0 gmax=0          # BZ-related variables

% 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.

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, and documents the categories and tokens it knows about.

3. The EXPRESS category

blm may include 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.

To create an input file with an EXPRESS category, use blm --express. On the other hand EXPRESS category is logically somewhat confusing. Here we use blm --brief, which does not add EXPRESS.

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.

$ lmfa --input
$ lmfa --help
$ lmfa --showp
$ lmfa --show | lmfa --show=2

Command-line switches for lmfa and other programs are documented on this page. In particular:

--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 specific to 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.

See Table of Contents

5. Initial setup: free atomic density and parameters for basis

lmfa is an initialization code with a dual purpose: to create free-atomic densites as a framework to make a starting trial density for a density functional calculation (stored in file atm.ext), and to generate default parameters defining the basis set (stored in file basp.ext). The atm file must be generated by lmfa; while basp.ext is a small, editable file which you can create on your own. This file’s purpose is to specify information for the basis set independently of the ctrl file: whether data is read from basp or the ctrl file depends on  HAM_AUTOBAS.

The atm file

lmfa generates self-consistent densities for free atoms. These are overlapped to form an initial crystal trial density (called the Mattheis construction). The following is stored in atm.ext:

• valence electron density on a mesh inside the augmentation spheres
• core electron density on a mesh inside the augmentation spheres
• The potential that generates this density
• fit to the interstitial density as a linear combination of smooth Hankel functions

The basp file

lmfa will also automatically generate a basis set, and save it in a file called basp0.ext (see callbox.) What is stored in basp depends on the contents of HAM_AUTOBAS, but typically three pieces of information are kept:

• Parameters P defining the boundary conditions of partial waves on augmentation boundaries
• Parameters (RSMH,EH) defining the shape of the smooth Hankel functions that form the envelope functions of the basis set
• Parameters PZ specifying local orbitals

As will be described in more detail in the pages below, lmfa prepares the following.

1. Make a self-consistent atomic density for each species.

2. 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.

3. 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.

4. 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.

5. Decide on which shallow cores should be included as local orbitals.
   Local orbitals are written basp0.pbte as nonzero values of  PZ.

6. 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  verb 35
 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 Barth-Hedin functional. To use a GGA, see here.

The autogen line says that lmfa will make the basis set information (points 3-5 outlined above).
The next few sections amplify on these three points. Point 6 is discussed here.

See Table of Contents

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.9336

(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, i.e. the state nearest 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 in this case.

15.9336 is actually a compound of 10 and the “continuous principal quantum number” 5.9336. 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 the extended kind is more accurate, and there is no advantage to doing so.

The continuous principal quantum number (5.9336) 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 less important and lmfa will not add them to basp.pbte. (Caveat: even at the LDA level local orbitals for transition metal d states have a modest effect, and local orbitals for 4f states can be important.) But they can sometimes be important in GW calculations, since in contrast to the LDA, unoccupied states also contribute to the potential.

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.

Alternatively, you can cause lmfa to remake the valence-core partitioning after the LO have been found, and write directly to file basp.pbte.

lmfa ctrl.pbte --usebasp

It saves effort since you don’t need to copy basp0.pbte to basp.pbte or run lmfa twice. On the other hand it gives you less flexibility since you no longer control which orbitals are to be treated as LO.

To see how to control and monitor lmfa’s search for local orbitals, see the annotated lmfa output.

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 4-vector

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

rm basp.pbte
lmfa --pr41 ctrl.pbte

You should see the following table:

 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 envelope function parameters

Some details of the basis set (envelope functions, augmentation, local orbitals) are explained in this tutorial.

lmfa loops over each species, generating a self-consistent 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 short-range 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.

Even while the free-atomic wave functions are very accurate, they are not optimal for the crystal. Tails of envelope functions from remote sites are not designed to solve the Schrödinger equation locally, and they contaminate the solution locally.

lmf requires RSMH and EH. Those generated by lmfa are reasonable, but unfortunately not optimal choices for the crystal, as just noted and further 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 APWs 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 spill-out 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) > 5e-3
 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)=3e-8  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 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:

 FREEAT:  estimate HAM_GMAX from RSMH:  GMAX=4.6 (valence)  8.2 (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.

See Table of Contents

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. Information is supplied through BZ_NKABC, which value is governed by variable nkabc in this input file. See this page for more details on how the Brillouin mesh is constructed, and how nkabc affects BZ_NKABC.

In this case a k-mesh of $6\times 6\times 6$ divisions is more than adequate. With your text editor or the sed command-line editor change nkabc=0 in the ctrl file to nkabc=6, e.g.

    sed -i s/nkabc=0/nkabc=6/ ctrl.pbte

or alternatively run blm again with the switch --nk=6; or alternatively assign variable nkabc on the command line using  -vnkabc=6 (which is what this tutorial will do).

Note: You can let blm selected a default value for this parameter by using --nk~ins for insulators or --nk~met for metals. However, be advised it is a somewhat hazardous thing to do. In this case blm selects a 4×4×4 mesh.

Charge density mesh

We also haven’t specified the G cutoff for the density mesh. blm will this for you if you use the --simple or --brief options, but by default blm will not set this value itself 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 (8.2). You must use the larger of the two.

Warning: 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=8.2), as we do in the next section. Or, you can run blm again, with command-line argument  -gmax=8.2.

See Table of Contents

7. Sanity Checks

This step shows how perform some sanity checks on the input. It is optional, but always advisable when you are getting started with a new system. Run lmf with a high verbosity, quitting after hamiltonian parameters are generated.

lmf ctrl.pbte -vnkabc=6 -vgmax=8.2 --pr55 --quit=ham

The output should confirm :

• the k mesh is 6×6×6 divisions, which generates 16 points on the irreducible wedge

   BZMESH: 16 irreducible QP from 216 ( 6 6 6 )  shift= F F F
  

• gmax=8.2 yields a charge density mesh of 20×20×20 divisions

   MSHSIZ: mesh has 20 x 20 x 20 divisions; length 0.429, 0.429, 0.429
  

  and this mesh is accurate enough to converge the plane-wave expansion of the smooth Hankel basis to about 10⁻⁶ (estimated error is printed in table sugcut)

• The basis consists of 30 Pb orbitals (spdfspd+LO) and 25 Te orbitals (spdfspd)

Also notice the table  Orbital positions ... . It locates the position of each orbital in the Hamiltonian, which is very useful when bands are to be plotted with color weights.

8. Autogenerate basis parameters using blm

Sections 5 and 6 carefully explained how to construct parameters for basis set, and charge density and Brillouin zone integration.

blm can perform these steps automatically. This saves you from knowing the details, but the downside is that the details can sometimes be important. For most purposes, it is sufficient and simpler to let blm perform these steps for you. The following pair of instructions

blm init.pbte --brief --ctrl=ctrl --nk~ins
lmfa pbte

will carry out all the steps of sections 5 and 6 automatically, creating a complete setup, building the ctrl file, the basp file, the site file, and the atm file, without any need for manual intervention.

The only difference is that blm selects a default k mesh with 4×4×4 divisions, rather than 6×6×6. Replacing --nk~ins with --nk~6 should generate an identical setup.

9. Self consistency

Carry out a self-consistent calculation as follows:

lmf ctrl.pbte -vnkabc=6 -vgmax=8.2

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                     |
...                                |
ITER                                  # Management of iterations to self-consistency
      NIT=    {nit}                   # Maximum number of iterations

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 l-cutoffs
• 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 Σ⁰ may be read in.

Self-consistent cycle

Each iteration of the self-consistency 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 Kohn-Sham hamiltonian.
  • Matrix elements of the interstitial potential $\langle\chi_i|V|\chi_j\rangle$. for envelope functions $\chi_i$.

  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. In the metallic case, the structure of the Fermi surface is very important for many properties. One simple property is the electronic contribution to the specific heat; see this tutorial for a demonstration within Questaal.
• 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 center-of-gravity
• Computes the Harris-Foulkes and Kohn-Sham 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.16e-2
  

  DQ=2.17e-2 is a measure of the root mean square deviation n^out−nⁱⁿ. At self-consistency this number should be small. This tutorial provides some additional explanation for how lmf mixes charge densites as it iterates to self-consistency.

• Checks for convergence.

lmf should converge to self-consistency in 10 iterations.

At the end of the self-consistency 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.004805 (0.018562)    tol= 0.000010 (0.000030)   more=T
i nkabc=6 gmax=8.2 ehf=-55318.1657749 ehk=-55318.1568595

Two checks are made: against the change (0.004805) in total energy and the RMS DQ (0.018562). When both checks fall below tolerances self-consistency 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=8.2 ehf=-55318.1620975 ehk=-55318.1620958

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 n^out−nⁱⁿ, 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. There is also a brief discussion in the Fe tutorial about monitoring convergence in the metallic case, where the presence of the Fermi surface makes Brillouin zone integrations more challenging.

Self-consistency in the LDA

See Table of Contents

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 self-consistent 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 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). See also 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

1. How does lmf iterate to self-consistency?

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

2. 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.

3. 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₆⁺ and L₆⁻ 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.

4. 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.

5. Is there an easy way to calculate effective masses?

   Yes, once you know where the band edge is. See this tutorial.

6. The augmentation $\ell$ cutoff in the ctrl file is 4 for Pb and 3 for Te. LAPW methods typically require much higher $\ell$ 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

1. 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. Be sure take note of the callout boxes for the atm file and the rst file.

2. Specify symops manually.

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

4. Try rotating the lattice and confirm that the total energy and energy bands are unchanged.

5. k-convergence. Try varying the number of k points or using staggering the mesh with BZ_BZJOB.

See Table of Contents