Properties of the lmf basis set

This tutorial describes the lmf basis set, and various kinds of cutoffs that affect convergence in the basis.

Table of Contents

• Preliminaries
• Definition of the basis set
  • Envelope functions
  • Augmentation
  • Local orbitals
  • Hilbert space and rank of the lmf Hamiltonian
• Tutorial
  • 1. Building the input file
  • 2. Checking sphere overlaps
  • 2a. General Remarks on choosing sphere radii
  • 3. The atomic density and basis set parameters
  • 4. GMAX and NKABC
    • 4.1 Setting GMAX
    • 4.2 Setting NKABC
  • 5. Self-consistent LDA calculation, small basis
    • 5.1 Convergence in the charge density
    • 5.2 Convergence in G cutoff
    • 5.3 Convergence in the forces
    • 5.4 Convergence in LMXA
    • 5.5 Convergence in KMXA
  • 6. Optimizing the basis set
    • 6.1 Tuning the envelope function parameters
    • 6.2 Increasing the number of envelope functions
    • 6.3 Adding APWs : the PMT method
  • Modifying the input file for GW
• Additional exercises

blm init.bi2te3
cp actrl.bi2te3 ctrl.bi2te3
sed -i s/nkabc=0/nkabc=3/ ctrl.bi2te3
sed -i s/gmax=0/gmax=4.4/ ctrl.bi2te3
sed -i 's/autobas\[pnu=1 loc=1 mto=4 lmto=4\]/autobas\[pnu=1 loc=1 lmto=3 mto=1 gw=0\]/' ctrl.bi2te3

lmchk bi2te3

lmfa bi2te3
mv basp0.bi2te3 basp.bi2te3
lmfa bi2te3

lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=2
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=3
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=4
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=5
lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=6
cat save.bi2te3 | vextract h nkabc ehf | tee dat

rm -f mixm.bi2te3 
rm -f save.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3

lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
lmf ctrl.bi2te3 -vgmax=6 -vnkabc=3 --quit=band

lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp ctrl.bi2te3 ctrl.bak
cat ctrl.bak | sed s/LMXA=./LMXA=6/ > ctrl.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp ctrl.bak ctrl.bi2te3

lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp ctrl.bi2te3 ctrl.bak
cat ctrl.bak | sed 's/LMXA=/KMXA=5 LMXA=/' > ctrl.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp ctrl.bak ctrl.bi2te3

rm -f mixm.bi2te3 
rm -f save.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --optbas
cp basp.bi2te3 basp.bak
cp basp2.bi2te3 basp.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
cp basp.bak basp.bi2te3

sed -i 's/autobas\[pnu=1 loc=1 lmto=3 mto=1 gw=0\]/autobas\[pnu=1 loc=1 lmto=4 mto=4 gw=0\]/' ctrl.bi2te3
lmfa bi2te3
cp basp.bi2te3 basp.bak
cp basp0.bi2te3 basp.bi2te3
sed -r -i "s/PZ=(([ ]{1,})?[0-9]+(\.[0-9]{1,})?){3}//" basp.bi2te3
lmfa bi2te3
rm -f mixm.bi2te3
rm -f rst.bi2te3
lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3

rm -rf out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=0 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=1 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=2 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=3 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=4 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=5 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=6 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=7 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf

cat save.bi2te3 | vextract i pwemax ehf ehk | tee etot.bi2te3
fplot -frme 0,.8,0,.5 -frmt th=3,1,1 -tmy 2.5 -vehf=-126808.313403 -s circ -ord '(x2-ehf)*1000' etot.bi2te3

blm --gw bi2te3

Preliminaries

The input file structure is briefly described in this lmf tutorial for PbTe, which you may wish to go through first.

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

Definition of the basis set

The lmf basis set $\chi_{\alpha\mathbf{R}L}(\mathbf{r})$ consist of smooth Hankel functions envelope functions $H_{\alpha\mathbf{R}L}(\varepsilon,r_s;{\mathbf{r}})$, which get augmented by solutions of the radial wave equation in augmentation spheres (partial waves). An additional set of local orbitals may be added to make the basis more complete in the augmentation region.

Envelope functions

The lmf basis set begins with smooth envelope functions. These functions are smooth Hankel functions, centered on an atom R. They are defined by smoothing radius $r_s$ and an energy $\varepsilon$. Smooth Hankels have the Slater-Koster form, $$H_{L}(\varepsilon,r_s;{\mathbf{r}}) = h_l(\varepsilon,r_s;r) Y_L(\hat\mathbf{r})$$

Here L is a compound index denoting the l and m angular momentum quantum numbers; the $h_l(\varepsilon,r_s;r)$ are radial functions defined by Eqns. (6-8) on this page, and the $Y_L(\hat\mathbf{r})$ are spherical harmonics.

The $H_{L}$ are convolutions of ordinary Hankel functions of energy $\varepsilon$ and Gaussian functions of smoothing radius $r_s$. These two parameters define the shape of the envelope: $\varepsilon$ controls the asymptotic decay for large r ($H_L \sim r^{-l-1}e^{-\kappa\,r}$, $\kappa^2{=}-\varepsilon$) and the smoothing radius $r_s$ demarcates the approximate point of transition from Gaussian-like shape at small r ($H_L \sim r^l$) to asymptotic behavior.

Each site has its own family of $H_{L}$. While it would be possible to have any number of $H_L$ on a particular site, each with unique values of $\varepsilon$ and $r_s$, in practice the Questaal codes allow two types ($\alpha{=}1,2$) of $H_{L}$ for a particular site R and angular momentum l. For the first envelope ($\alpha{=}1$), you define $(r_s,\varepsilon)$ pair through parameters RSMH and EH: Each l gets its own RSMH and EH. You can elect to limit the basis to one envelope function for a particular l, or choose a second envelope. The second set ($\alpha{=}2$) is optional: parameters are defined through RSMH2 and EH2. A single $H_l$ (“single kappa”) has the advantage of making for a small basis, but its accuracy is limited unless the system is fairly close-packed. At least for low l, i.e. for orbitals whose atomic levels are below or not far above the Fermi level, it is recommended that the “double kappa” basis (two functions per l) be used.

This flexibility in choosing is both a blessing and a curse. A blessing because the flexibility allows for more variational freedom, keeping the basis at low rank while maintaining high accuracy. But it is a curse because of the added burdens imposed on the user to determine the parameters. Usually you can allow the atom program lmfa to automatically generate parameters for basis set.

Thus we can identify the entire family of envelopes by $H_{\alpha\mathbf{R}L}$, labeling whether the first or second envelope, the site, and the L. Note that for each L, there are 2l+1 basis functions.

Augmentation

Each smooth envelope is “augmented” in a sphere around each nucleus by partial waves, which are numerical solution of the Schrodinger equation in the sphere up to an augmentation radius (specified by tag R in the species category). For each l (and m for a given l) up to a cutoff lmxa, the $H_{\alpha\mathbf{R}L}$ is replaced (“augmented”) by a linear combination of partial waves $\phi_l$ and $\dot\phi_l$. These partial waves form the Hilbert space of essentially exact solutions to the Schrodinger equation, to linear order in energy around some linearization energy $\varepsilon_\nu$. This energy is normally allowed to float in the course of the self-consistency cycle, to minimize the linearization error. See also the description of partial waves on this page. Suitable linear combinations of $\phi_l Y_L$ and $\dot\phi_l Y_L$ are taken in each sphere so as to make the augmented $H_{\alpha\mathbf{R}L}$ continuous and differentiable.

Local orbitals

It may be the linear method is not sufficient and that a third partial wave is required to make the basis complete over a sufficiently wide energy window. A conventional local orbital is realized by integrating the Schrodinger equation at another energy either far above or far below the energy used to make $\phi$ and $\dot\phi$, to obtain another partial wave $\phi_{z}$. A suitable combination of $\phi$ and $\dot\phi$ is subtracted to make the value and slope of the local orbital vanish at the augmentation radius. With such a construction it need not have an interstitial part. An “extended” local orbital (suitable for semicore states far below the Fermi level) is constructed by attaching a single smooth Hankel function, whose energy and smoothing radius are varied to match value and slope of the partial wave. You specify the local orbital through tag PZ in the species category. lmfa can automatically find deep local orbitals (one principal quantum number below the valence partial waves) that may be needed to make the basis reasonably complete. You may also select high-lying local orbitals (one principal quantum number above the valence partial waves). These are usually less relevant in DFT; however, high-lying d local orbitals on transition metals were necessary to obtain good agreement with other codes in the DeltaCodes validation exercise. Note that addition of local orbitals increases the rank of the Hamiltonian.

Hilbert space and rank of the lmf Hamiltonian

The Hilbert space of the lmf basis $\chi_{\alpha\mathbf{R}L}(\mathbf{r})$ then consists of the following:

1. In the interstitial, smooth Hankel functions $H_{\alpha\mathbf{R}L}(\mathbf{r}){=}H_{\alpha{L}}(\mathbf{r{-}R})$. They can be expanded in plane waves to make matrix elements of the potential.
   Note: For a given set $\alpha{=}1$ or $\alpha{=}2$ you must include all l’s up to some basis cutoff, i.e. l=0,1,…lmxb.  lmxb need not be the same for the $\alpha{=}2$ set of envelopes as for the first.

2. In augmentation sphere R up to the augmentation l cutoff lmxa, the pair of partial waves ($\phi_{\mathbf{R}l}(r)Y_L(\hat\mathbf{r})$, $\dot\phi_{\mathbf{R}l}(r)Y_L(\hat\mathbf{r})$). The dominant partial wave is $\phi_{\mathbf{R}l}(r)$; mostly it attaches on to smoothed Hankels centered at R. $\dot\phi_{\mathbf{R}l}(r)Y_L(\hat\mathbf{r}$ mostly attaches on to the (tails) of smooth Hankels $H_{\alpha\mathbf{R'}L}(\mathbf{r})$ centered at other sites and makes a smaller contribution. Finally there may possibly be local orbitals $\phi_{\mathbf{R}zl}(r)Y_L(\hat\mathbf{r})$ in some l channels. Local orbitals may be high-lying (far above the linearization energy) or deep, to include high-lying (semi)core levels in the valence. This extra set is specified through parameter PZ in the main input (ctrl file) file or the basp file.

3. In augmentation sphere R above lmxa, the tails $H_{\alpha\mathbf{R'}L}(\mathbf{r})$ of smoothed Hankels at sites other than R. The special manner in which augmentation is done enables the lmf basis set to converge very rapidly with lmxa — much faster than traditional all-electron basis sets. See below.

The total rank of the hamiltonian is then the number of $H_{\alpha\mathbf{R}L}$ you specify on all sites, plus any local orbitals specified. The size of the basis (number of first kappa, second kappa, and local orbitals) is printed in a table in lmf’s standard output.

Tutorial

1. Building the input file

This step is essentially identical to the first step in the PbTe tutorial. An abbreviated version is presented here.

Cut and paste the following into init.bi2te3.

# from http://cst-www.nrl.navy.mil/lattice/struk/c33.html
# Bi2Te3 from Wyckoff
% const a=4.3835 c=30.487 uTe=0.788 uBi=0.40
LATTICE
   SPCGRP=R-3M
   UNITS=A
   A={a} C={c}
SITE
    ATOM=Te C=0     0    0
    ATOM=Te C=0     0   {uTe}
    ATOM=Bi C=0     0   {uBi}

This tutorial explains how the input files init.ext and ctrl.ext are structured.

To create the skeleton input file invoke blm:

$ blm init.bi2te3
$ cp actrl.bi2te3 ctrl.bi2te3

There 2 Bi atoms and 3 Te atoms in the unit cell. Two of the Te atoms are symmetry equivalent.

This template will not work as is; three essential pieces of information which blm does not supply are missing, to rectify this:

• You must specify plane-wave cutoff GMAX. This can be obtained by running the lmfa tool as seen in the PbTe tutorial. For the purpose of this tutorial we will use the value GMAX=4.4. A more detailed description of GMAX is given below.
• You must specify a valid k mesh for Brillouin zone integration. The correct k mesh size is dependent on your problem, you can find an example in the PbTe tutorial. For this tutorial we will use a k mesh of $6\times 6\times 6$ divisions, nkabc=3. A more detailed description of nkabc is given below.
• You must define a basis set which can be done manually or automatically, as described next.

Note: if you are using the GMAX and k-mesh values listed above, you can use the following commands to automatically edit your ctrl.bi2te3 file:

$ sed -i s/nkabc=0/nkabc=3/ ctrl.bi2te3
$ sed -i s/gmax=0/gmax=4.4/ ctrl.bi2te3

In either case, to edit these values you simply need to find the gmax and nkabc variables in the % const section of your ctrl.bi2te3 file and modify them to the desired value. The sed command above does this for you.

blm creates a family of tags belonging to AUTOBAS to enable other programs to automatically find a basis set for you. We will use this tag, which sets up a standard minimal basis:

autobas[pnu=1 loc=1 lmto=3 mto=1 gw=0]

This is an alias for the tag in the HAM category

AUTOBAS[PNU=1 LOC=1 LMTO=3 MTO=1 GW=0]

Note that you must modify ctrl.bi2te3 a little; the default gives [… LMTO=5 MTO=4], which makes a double kappa basis. Once again, you can use the following command to perform this replacement automatically:

$ sed -i 's/autobas\[pnu=1 loc=1 mto=4 lmto=4\]/autobas\[pnu=1 loc=1 lmto=3 mto=1 gw=0\]/' ctrl.bi2te3

lmfa calculates wave functions and atomic densities for free atoms. It also has a mode that automatically generates information for basis sets, using tokens in AUTOBAS to guide it. This information is written to a file basp0.ext. AUTOBAS specifies set of conditions that enable lmfa to automatically build the basis set for you, as described below. Parameters are generated, but you can modify them as you like.

See Table of Contents

2. Checking sphere overlaps

Sphere overlaps can be checked using lmchk. To do this copy the template actrl.bi2te3 to the input file and run lmchk:

$ lmchk bi2te3

You should see the site positions for all the atoms:

   Site    Class            Rmax        Hcr                 Position
    1      1  Te          2.870279    2.009195    0.00000    0.00000    0.00000
    2      3  Te2         2.870279    2.009195   -0.50000   -0.86603    1.46162
    3      3  Te2         2.870279    2.009195    0.50000    0.86603   -1.46162
    4      2  Bi          2.856141    1.999299    0.50000    0.86603    0.80309
    5      2  Bi          2.856141    1.999299   -0.50000   -0.86603   -0.80309

and a table of overlaps

 ib jb  cl1     cl2        Pos(jb)-Pos(ib)      Dist  sumrs   Ovlp    %    summt   Ovlp   %
  1  4  Te      Bi        2.391 -4.142  3.841  6.134  5.726  -0.41  -6.6   4.008  -2.13 -34.7
  1  5  Te      Bi       -2.391 -4.142 -3.841  6.134  5.726  -0.41  -6.6   4.008  -2.13 -34.7
  2  4  Te2     Bi       -2.391 -4.142 -3.149  5.726  5.726   0.00   0.0*  4.008  -1.72 -30.0
  3  5  Te2     Bi        2.391 -4.142  3.149  5.726  5.726   0.00   0.0*  4.008  -1.72 -30.0

 OVMIN, 107 pairs:  fovl = 0   <ovlp> = 0%   max ovlp = 0%

By default, blm makes the spheres as large as possible without overlapping. In this case the Bi and Te radii are nearly the same. and at least some spheres touch without overlapping.

The packing fraction is printed

 Cell volume= 1141.20380   Sum of sphere volumes= 492.34441 (0.43143)

Generally larger packing fractions are better because the augmentation partial waves are quite accurate. 0.43 is a fairly good packing fraction.

2a. General Remarks on choosing sphere radii

Choosing sphere radii is an important business in an augmented-wave method. Ideally results should not depend on the choice, and in well converged cases this is nearly true. But it can never be exactly because as the radii change, the basis set of some portion of the volume is represented differently (partial waves inside spheres up to a finite l cutoff, and a limited famly of envelope functions in the interstitial). Particularly when basis sets are small, there is some dependence on radii and a poor choice can lead to poor results. As a minimum you should be aware these radii are parameters, and that there may be some dependence on their choice.

Sphere radii are chosen by balancing between making them on the one hand as large as possible: the augmentation basis is of high quality and partial waves are accurate and more complete than envelope functions. This argues for large radii, but l convergence of the augmented part is slower with larger radii, and also there is a geometry violation when spheres overlap. (Overlap is defined as $(r_1{+}r_2{-}d)/d$, where $r_1$ and $r_2$ ae the two sphere radii, and $d$ is the bond length.) Both considerations argue for smaller radii.

Finally there is the question of how to select the relative size of radii of different species.

Fortunately blm can determine relative sphere radii in a nearly optimal way, which it does by constructing an approximate potential and varying radii to make the potentials on the augmentation boundaries as nearly equal as possible. This makes the potential as close to a constant as possible, which improves the ability of the envelope functions to represent the interstitial potential and also reduces the burden on the interstial matrix element maker. Thus is is a good idea to let blm find the sphere radii for you, or if you are starting from a ctrl file, use lmchk --getwsr to supply suggested radii. You can tell blm or lmchk what maximum overlap you will tolerate (blm accepts a command-line switch –omax, and lmchk reads tag OMAX1 in the SPEC category).

What to consider when selecting upper bounds to overlaps depend on the context. It is useful to keep in mind the following obervations. Further observations, particularly in an ASA context, can be found on this page.

• ASA codes lm, lmgf, lmpg: Sphere overlaps are necessary because the ASA has no interstitial — thus the approximations in the ASA cannot be systematically reduced. Care must be taken that they are not too large; sometimes this necessitates using empty spheres. bcc lattices (overlap 14%) and fcc lattices (overlap 11%) are good. Overlaps of 18% or 19% can be necessary for some bonds in complex cases, e.g. incommensurate interfaces.

• Full potential code lmf, traditional basis set: Strictly speaking there should be no geometry violation. All the same lmf has a special three-component augmentation (see Sec. 3.6. of Ref. 1) that allows it to tolerate sphere overlaps up to about 10% without incurring significant error. It has been found empirically however, that convergence to self-consistency becomes more delicate when overlaps become too large.

• Full potential code lmf, Jigsaw Puzzle Orbital basis set: The JPO basis is still undergoing development, but when it implemented it should provide a marked improvement over the traditional one. We anticipate that the sphere radii can be chosen rather small without significant loss of precision, e.g. −10% overlap. At present this is only speculation.

• GW code driven by lmfgwd: The GW suite does not have the same elegant augmentation as the LDA. Therefore you should observe the constraint of no overlaps fairly carefully. (1% seems to be acceptable).

See Table of Contents

3. The atomic density and basis set parameters

If you did not do so already copy actrl.bi2te3 to ctrl.bi2te3 and change [… MTO=4  LMTO=4] → [… MTO=1  LMTO=3]).

$ cat actrl.bi2te3 | sed 's/mto=4 lmto=4/mto=1 lmto=3/' > ctrl.bi2te3

Invoke lmfa:

$ lmfa bi2te3

The primary purpose of lmfa is to generate a free atom density. A secondary purpose but nevertheless a very useful one, is to supply additional information about the basis set in an automatic way. All of this information that blm omits can be supplied manually in the input file, but once an input file is available, lmfa can automatically supply most of the key missing information for a minimal working input file. For example it generates a basp file basp0.bi2te3m which contains information about the shape of the basis:

BASIS:
Te RSMH= 1.615 1.681 1.914 1.914 EH= -0.888 -0.288 -0.1 -0.1 P= 5.901 5.853 5.419 4.187
Bi RSMH= 1.674 1.867 1.904 1.904 EH= -0.842 -0.21 -0.1 -0.1 P= 6.896 6.817 6.267 5.199 5.089 PZ= 0 0 15.9362

Every species gets one line. This file specifies a basis set consisting of spdf orbitals on Te sites, and spdf orbitals on Bi sites, and a local 5d orbital on Bi. See the PbTe tutorial for further description of these parameters. This choice is generally a good choice, but it optimized for the atom, not the crystal. There is an automatic procedure to optimise envelope function parameters RSMH and EH if you wish to do so.

Note: Remember that lmf reads from basp.ext, not basp0.ext.

The partitioning between valence and core states is something that requires a judgement call. lmfa has made a default choice: from basp0.bi2te3 you can see that lmfa selected the 6s, 6p, 6d, 5f states, populating them with charges   2, 3, 0, 0, making the total sphere charge Q=0. You can override the default, e.g. choose the 5d over the 6d with SPEC_ATOM_P; override the l channel charges with SPEC_ATOM_Q.

lmfa does the following to find basis set parameters:

• Automatically finds deep states to include as valence electrons.
• Selects sphere boundary condition for partial waves
• Finds envelope function parameters

The process is essentially the same as described in the PbTe tutorial; it is described in some detail there and in the annotated lmfa output. The main difference is that in this case a smaller, single-kappa basis was specified (LMTO=3); lmfa makes (RMSH,EH) instead of the double kappa (RMSH,EH; RMSH2,EH2). Later we will improve on the basis by adding APW’s.

With your text editor insert lines from basp0.bi2te3 in the SPEC category of the ctrl file, viz

  ATOM=Te         Z= 52  R= 2.870279  LMX=3  LMXA=3
    RSMH= 1.615 1.681 1.914 1.914 EH= -0.888 -0.288 -0.1 -0.1 P= 5.901 5.853 5.419 4.187
  ATOM=Bi         Z= 83  R= 2.856141  LMX=3  LMXA=4
    RSMH= 1.674 1.867 1.904 1.904 EH= -0.842 -0.21 -0.1 -0.1 P= 6.896 6.817 6.267 5.199 5.089 PZ= 0 0 15.936

Alternatively make lmf read these parameters from basp.bi2te3. Copy basp0.bi2te3 to basp.bi2te3, and modify it as you like. basp.ext is read after the main input file is read, if it exists. According to which of following tokens is present, their corresponding parameters will be be read from the basp file, superseding prior values for these contents:

If this information is supplied in both the ctrl file and the basp file, values of MTO and LOC tell lmf which to use. In this tutorial we will work just with the basp file.

$ cp basp0.bi2te3 basp.bi2te3

The atm file was created by lmfa without prior knowledge that the 5d local orbital is to be included as a valence state (via a local orbital). Thus it incorrectly partitioned the core and valence charge. You must do one of the following:

1. Remove PZ=0 0 15.936 from basp.bi2te3. It will no longer be treated as a valence state. Removing it means the remaining envelope functions are much smoother, which allows you to get away with a coarser mesh density. Whether you need it or not depends on the context: with GW, for example, this state is a bit too shallow to be treated with Fock exchange only (which is how cores are handled in GW).

2. Copy basp0.bi2te3 to basp.bi2te3 and run lmfa over again:

$ cp basp0.bi2te3 basp.bi2te3
$ lmfa bi2te3

With the latter choice lmfa operates a little differently from the first pass. Initially the Bi 5d was part of the core (similar to the Pb 5d in the Pbte tutorial; now it is included in the valence.

See Table of Contents

4. GMAX and NKABC

4.1 Setting GMAX

blm makes no attempt to automatically assign a value to the plane wave cutoff for the interstitial density. This value determines the mesh spacing for the charge density. lmf reads this information through HAM_GMAX or EXPRESS_gmax, or HAM_FTMESH. It is a required input; but blm does not automatically pick a value because its proper choice depends on the smoothness of the basis. In general the mesh density must be fine as the most strongly peaked orbital in the basis.

lmfa makes RSMH and EH and can determine an appropriate value for HAM_GMAX based on them. In the present instance, when the usual 6s, 6p, 6d, 5f states are included lmfa recommends GMAX=4.4 as can be seen by inspecting the first lmfa run. If you keep the 5d in the valence and run lmfa second time you will see this output

 FREEAT:  estimate HAM_GMAX from RSMH:  GMAX=4.4 (valence)  8.1 (local orbitals)

The valence states still require only GMAX=4.4 but the 5d state is strongly peaked but GMAX=8.1 for the local orbitals. The 5d state is strongly peaked at around the atom, and requires more plane waves to represent reasonably, even a smoothed version of it, than the other states. The difference between 8.1 and 4.4 is substantial, and it reflects the additional computational cost of including deep core-like states in the valence. This is the all-electron analog of the “hardness” of the pseudopotential in pseudopotential schemes. If you want high-accuracy calculations (especially in the GW context), you will need to include these states as valence. Including the Bi 5d is a bit of overkill for LDA calculations however. If you eliminate the Bi 5d local orbital you can set GMAX=4.4 and significantly speed up the execution time. It is what this tutorial does.

4.2 Setting NKABC

blm assigns the initial k-point mesh to zero. Note the following lines in ctrl.bi2te3:

% const nkabc=0 gmax=0
...
# Brillouin zone
  nkabc=  {nkabc}                  # 1 to 3 values

Note: nkabc is simultaneously a variable and a tag here. This can be somewhat confusing. The expression  {nkabc}  gets parsed by the preprocessor and is turned into the value of variable nkabc (see how nit gets turned into 10 in the PbTe tutorial), so after preprocessing, the argument following tag is a number.

EXPRESS_nkabc (alias BZ_NKABC) governs the mesh of k-points. An appropriate choice will depend strongly on the context of the calculation and the sytem of interest; the density-of-states at the Fermi level; whether Fermi surface properties are important; whether you want optical properties as well as total energies well described; the precision you need; the integration method, and so on. Any automatic formula can be dangerous, so blm will not choose an operational default.

The most reliable way determine the appropriate mesh is to vary nkabc and monitor the convergence. We don’t need a self-consistent calculation for that: the k-convergence will not depend strongly whether the potential is converged or assembled from free atoms.

Do the following (assuming no 5d local orbital)

$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=2
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=3
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=4
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=5
$ lmf ctrl.bi2te3 -vgmax=4.4 --quit=band -vnkabc=6

The meaning of --quit=band is described here.

Total energies are written to the save file save.bi2te3. It should read:

h gmax=4.4 nkabc=2 ehf=-126808.2361717 ehk=-126808.1583957
h gmax=4.4 nkabc=3 ehf=-126808.3137885 ehk=-126808.2492178
h gmax=4.4 nkabc=4 ehf=-126808.3168406 ehk=-126808.2505156
h gmax=4.4 nkabc=5 ehf=-126808.3165536 ehk=-126808.2497121
h gmax=4.4 nkabc=6 ehf=-126808.3164058 ehk=-126808.2494041

You can use the vextract tool to conveniently extract a table of the Harris-Foulkes total energy as a function of nkabc

$ cat save.bi2te3 | vextract h nkabc ehf | tee dat

You can plot the data, or just see by inspection that the energy is converged to less than a mRy with 4×4×4 k mesh and about 0.1 mRy with a 5×5×5 k mesh. A detailed analysis of k point convergence is given here.

See Table of Contents

5. Self-consistent LDA calculation, small basis

In the following we will use nkabc=3, though after convergence is reached we might consider increasing it a little. Before doing the calculation you can use your text editor to set nkabc=3 and gmax=4.4, or continue to set assign values from the command line.

The density can be made self-consistent with

$ rm -f mixm.bi2te3 save.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3

5.1 Convergence in the charge density

You should lmf reach self-consistency in 9 iterations.

The Harris-Foulkes and Kohn-Sham total energies are ehf=-126808.3137885 and ehk=-126808.2492178 from the Mattheis construction. At self-consistency they come together: ehf=-126808.2950696 and ehk=-126808.2950608 (the large integer part comes almost entirely from the core states.)

The self-consistent value is 18 mRy less binding than the Harris-Foulkes energy of the Mattheis construction, and 46 mRy more binding than the corresponding Kohn-Sham energy. That the two initial functionals bracket the self-consistent result, and that the HF is generally closer to the final result than the HK functional, is typical behavior. (The HF functional is generally more stable.)

5.2 Convergence in G cutoff

The G cutoff, 4.4u Ry^1/2 yields precisions in plane-wave expansions of envelope functions as shown this table:

 sugcut:  make orbital-dependent reciprocal vector cutoffs for tol=1.0e-6
 spec      l    rsm    eh     gmax    last term   cutoff
  Te       0    1.61  -0.89   4.603    3.31E-06    1635*
  Te       1    1.68  -0.29   4.662    5.08E-06    1635*
  Te       2*   1.91  -0.10   4.273    1.15E-06    1539
  Te       3    1.91  -0.10   4.471    1.71E-06    1635*
  Bi       0    1.67  -0.84   4.441    1.29E-06    1635*
  Bi       1    1.87  -0.21   4.183    1.08E-06    1411
  Bi       2*   1.90  -0.10   4.297    1.37E-06    1539
  Bi       3    1.90  -0.10   4.497    2.06E-06    1635*

gmax=4.4 isn’t quite large enough to push the error below tolerance (1.0e-6), but it is pretty close. You can check how well the total energy is converged by increasing gmax:

$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band
$ lmf ctrl.bi2te3 -vgmax=6 -vnkabc=3 --quit=band

and comparing ehf in the last two lines of save.bi2te3. You should find that the energy is converged to about 0.1 mRy.

See Table of Contents

5.3 Convergence in the forces

The exact forces are the change in total energy with respect to displacement of the nucleus. As they are calculated in the Questaal package, some correction terms are added to accelerate convergence with respect to deviations nⁱⁿ−n in the guessed density from the self-consistent one. See the annotation of lmf output.

Note: The forces are not exact derivatives of the total energy. This is because the change in shape of the augmented partial waves $\phi$ and $\dot{\phi}$ is not taken into account as a nucleus displaces. The effect is usually very small. Forces should be exactly consistent with the energy if the shape of the partial waves is frozen, however, which you can do with HAM_FRZWF.

To what extent the basis set affects the forces is taken up in the Additional Exercises.

5.4 Convergence in LMXA

In an augmented wave method, envelope functions centered a different site is must be expanded around a local site, as one-center expansions in spherical harmonics Y_lm. LMXA is a cutoff that truncates the one-center expansion to a finite l = LMXA. The input file reads:

SPEC
  ATOM=Te         Z= 52  R= 2.870279  LMX=3  LMXA=3
  ATOM=Bi         Z= 83  R= 2.856141  LMX=3  LMXA=4

LMXA=3 and LMXA=4 are very low l cutoffs for an augmented wave method. They are about half of what is needed for a reasonably converged calculation with traditional augmentation. However, the Questaal suite has a unique augmentation procedure that converges very rapidly with l. Indeed, it can be seen as a justification for pseudopotential methods that limit the pseudopotential to very low l, e.g. l=2.

It is a useful exercise to see how well converged the total energy is using the default values LMXA=3 and LMXA=4. First, run lmf with the unaltered ctrl file:

$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band

Set both LMXA to 6: and run lmf again

$ cp ctrl.bi2te3 ctrl.bak
$ cat ctrl.bak | sed s/LMXA=./LMXA=6/ > ctrl.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band

When the restart file is read you should see

  site   1, species Te      : augmentation lmax changed from 3 to 6
  site   1, species Te      : inflate local density from nlm= 16 to 49

Compare the last two lines of save.bi2te3. You should be able to confirm that the energy change is 0.25 mRy, By playing around with the two LMXA you can confirm that increasing the Te LMXA to 4, nearly all the error is eliminated.

Before continuing, be sure to restore the original ctrl file

$ cp ctrl.bak ctrl.bi2te3

See Table of Contents

5.5 Convergence in KMXA

Ordinary Hankel functions can be expanded in Bessel functions around a remote site. This follows from the fact that both are solutions of the same second order differential equation, one being regular at the origin and the other being regular at ∞. Smooth Hankel functions are more complicated: the one-center expansion has no such elementary function, but it can be written as a linear combination of Laguerre polynomials P_kl of half integer order; see Eq. (12.17) in this J. Math. Phys. paper.

The polynomials are cut off at a fixed order given by KMXA. blm doesn’t write KMXA to the template file; instead a default value is used, which is written to this table

 species data:  augmentation                           density
 spec       rmt   rsma lmxa kmxa      lmxl     rg   rsmv  kmxv foca   rfoca
 Te       2.870  1.148    3    3         3  0.718  1.435    15    1   1.148
 Bi       2.856  1.142    4    3         4  0.714  1.428    15    1   1.142

You can check the convergence by in KMXA by supplying an input for it. First run lmf with the unmodified file:

$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band

Set KMXA to 5 and run lmf again:

$ cp ctrl.bi2te3 ctrl.bak
$ cat ctrl.bak | sed 's/LMXA=/KMXA=5 LMXA=/' > ctrl.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band

When the restart file is read you should see an indication that KMXA has been increased:

  site   1, species Te      : augmentation kmax changed from 3 to 5
  site   2, species Te      : augmentation kmax changed from 3 to 5

Compare the last two lines of save.bi2te3. You should be able to confirm that the energy change is 0.15 mRy.

Note: Before continuing, be sure to restore the original ctrl file.

$ cp ctrl.bak ctrl.bi2te3

See Table of Contents

6. Optimizing the basis set

6.1 Tuning the envelope function parameters

The envelope function shape generated by lmfa is tailored to the atom. They are very good basis sets for that case, but less so for the crystal.

Ideally a basis of the size could be constructed that yields energies converged almost as well as we accomlish for the free atom. This is difficult to do in practice, though it is the aim of the Jigsaw Puzzle Orbitals basis.

The best we can with the existing lmf basis is to optimize the energy by tuning RSMH and EH. This cannot be done analytically, but switch --optbas causes lmf to perform this optimization by brute-force minimization of the Harris-Foulkes energy.

$ rm -f mixm.bi2te3 
$ rm -f save.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3

$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --optbas > out

 LMFOPB:  optimizing energy wrt 8 parameters, etol=5e-4:
   spec       l  type     start
   1:Te       0   rsm     1.615
   1:Te       1   rsm     1.681
   1:Te       2   rsm     1.914
   1:Te       3   rsm     1.914
   2:Bi       0   rsm     1.674
   2:Bi       1   rsm     1.867
   2:Bi       2   rsm     1.904
   2:Bi       3   rsm     1.904

 LMFOPB:  before optimization ehf=-126808.2951.  After optimization ehf=-126808.3092
 Exit 0 Optimization complete.  See file basp2

$ grep LMFOPB out | grep optimized

 LMFOPB: var #1 optimized to 1.602 (7 iter):  ehf=-126808.295097, delta ehf=-2.38e-5
 LMFOPB: var #2 optimized to 1.617 (5 iter):  ehf=-126808.298989, delta ehf=-0.00392
 LMFOPB: var #3 optimized to 1.778 (6 iter):  ehf=-126808.299166, delta ehf=-1.77e-4
 LMFOPB: var #5 optimized to 1.605 (5 iter):  ehf=-126808.301, delta ehf=-9.12e-4
 LMFOPB: var #6 optimized to 1.658 (7 iter):  ehf=-126808.307945, delta ehf=-0.00695
 LMFOPB: var #7 optimized to 2.245 (5 iter):  ehf=-126808.308856, delta ehf=-9.11e-4
 LMFOPB: var #8 optimized to 1.679 (6 iter):  ehf=-126808.309198, delta ehf=-3.42e-4

$ cp basp.bi2te3 basp.bak

BASIS:
 Te RSMH= 1.615 1.617 1.914 1.914 EH= -0.888 -0.288 -0.1 -0.1 P= 5.901 5.853 5.419 4.187
 Bi RSMH= 1.674 1.658 1.904 1.904 EH= -0.842 -0.21 -0.1 -0.1 P= 6.896 6.817 6.267 5.199 5.089

$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3 --quit=band

$ cp basp.bak basp.bi2te3

In sum,

• without optimizing the basis, you should achieve total energy ehf=-126808.2950696 at self-consistency.
• With the --optbas switch it should become ehf=-126808.306962.

See Table of Contents

6.2 Increasing the number of envelope functions

lmfa uses the lmto tags to make default choices for the size of basis. The tutorial used a relative minimal basis, lmto=3.

    autobas[pnu=1 loc=1 lmto=3 mto=1 gw=0]

With your editor, modify this line to read

    autobas[pnu=1 loc=1 lmto=4 mto=4 gw=0]

This will increase the basis, most notably include two envelope functions per l channel (double-kappa)

Run lmfa again and note how the basp file changed:

$ lmfa bi2te3
$ diff basp0.bi2te3 basp.bi2te3

Remove the PZ=… as before and copy basp0.bi2te3 to basp.bi2te3. Save the original file in a backup.

$ nano basp0.bi2te3
$ cp basp.bi2te3 basp.bak
$ cp basp0.bi2te3 basp.bi2te3

Since we modified the basp file we need to re-run lmfa:

$ lmfa bi2te3

Run lmf to self-consistency

$ rm -f mixm.bi2te3
$ rm -f rst.bi2te3
$ lmf ctrl.bi2te3 -vgmax=4.4 -vnkabc=3

Note that size of the basis has increased: With lmto=3 there are 80 orbitals

 Makidx:  hamiltonian dimensions Low, Int, High, Negl: 80 0 18 27
 suham :  98 augmentation channels, 98 local potential channels  Maximum lmxa=4

which increases to 125 orbitals with lmto=5:

 Makidx:  hamiltonian dimensions Low, Int, High, Negl: 125 0 71 54
 kappa   Low   Int   High  L+I  L+I+H  Neglected
   1      80     0    18    80    98      27
   2      45     0    53    45    98      27
  all    125     0    71   125   196      54
 suham :  98 augmentation channels, 98 local potential channels  Maximum lmxa=4

At self-consistency you should find that the total energy converges to ehf=-126808.313403 Ry.

See Table of Contents

6.3 Adding APWs : the PMT method

You can also use augmented plane waves (APWs) to improve on the basis set. The combination of smooth Hankel functions with APW’s is called the PMT basis set. Adding APW’s provides a systematic way of making the basis of envelope functions complete.

These tags in the HAM category control how many APWs are added:

  PWMODE={pwmode} PWEMIN=0 PWEMAX={pwemax}  # For APW addition to basis

To use APW’s some tolerances in the ctrl file should be tightened. If you invoke blm with --pmt these extra tolerances are included automatically. That is the easiest way to make the changes but here we just add include tolerances HAM_TOL and EWALD_TOL by hand. Insert two lines at the end of the HAM category:

HAM
...
      FORCES={so==0} ELIND=-0.7
      TOL=1d-16
EWALD TOL=1d-16

One problem with the PMT method is that the smoooth-Hankel and APW basis set span nearly the same hilbert space. If you add too many plane waves the overlap matrix does not remain positive definite. Tightening the tolerances above helps, and so does increasing the fineness of the density mesh, as these points are also used for the PW expansion of the basis.

Note: as a last resort, you can stabilize the overlap with the HAM_OVEPS switch, but it is better to reduce the rank of the LMTO or APW basis sets.

PWMODE=11 adds APWs in a standard way: envelope functions of the form e^i(k+G)·r are added to the smooth Hankel basis. You must specify the PW cutoffs. Here you can specify both the lower and upper bounds since the smooth Hankels will efficiently cover most of the lower part of the space.

With these extra considerations, run lmf as follows

rm -f out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=0 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=1 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=2 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=3 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=4 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=5 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=6 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf
lmf ctrl.bi2te3 -vpwmode=11 -vpwemax=7 -vgmax=8 -vnkabc=3 --quit=band >> out.lmf

and observe how the total energy decreases with pwemax. Use the vextract tool:

cat save.bi2te3 | vextract i pwemax ehf ehk | tee etot.bi2te3

Draw a picture of the total energy relative to the double-kappa value. The fplot command shown in the box below will generate a postscript file; the figure actually shown has been elaborated a little. The red circle shows the self-consistent double-kappa result of Sec. 6.2. The light grey line follows the PMT procedure as above, but taking for the LMTO part a smaller, optimized spd single kappa basis (see Additional exercises). Numbers in parenthesis are the number of LMTO’s. The red circle is the 2-kappa basis without APWs; the purple is the LAPW basis without LMTOs.

$ fplot -frme 0,.8,0,.5 -frmt th=3,1,1 -tmy 2.5 -vehf=-126808.313403 -s circ -ord '(x2-ehf)*1000' etot.bi2te3

The figure shows that the double-kappa basis (additional 45 orbitals) stabilizes the single-kappa spdf basis by about 18 mRy, and that the same can be accomplished with APWs with a plane-wave cutoff of 2 Ry (additional 48-62 APWs). By further increasing the APW cutoff, another 5 mRy can be gained. For most purposes this extra gain is not important. Note that the APW basis with E_max=7 Ry is quite large: 343-370 APWs.

Note that the APW basis is generally less efficient than the LMTO basis. To reach a precision comparable to the 2-kappa basis (125 orbitals) starting from the 1-kappa spd basis, about 160 APWs are needed, or about 200 orbitals all told. The power of the PMT method can be compared against a straight LAPW basis. About 300 APWs are needed to achieve the convergence of the double kappa basis. See purple line in the Figure.

To see how many orbitals the APW basis adds, do:

$ grep ndham out.lmf

See Table of Contents

Modifying the input file for GW

GW calculations demand more of the basis set because unoccupied states are important. To set up a job in preparation for a GW calculation, invoke blm as :

$ blm --gw bi2te3

Compare actrl.bi2te3 generated with the --gw switch to one without. One important difference will be that the default basis parameters are modified because AUTOBAS becomes:

AUTOBAS[PNU=1 LOC=1 LMTO=5 MTO=4 GW=1]

The basis is similar to LMTO=4 but EH has been set a little deeper. This helps the QSGW implementation interpolate between k-points. The larger basis makes a minor difference to the valence bands; but the conduction bands change, especially the higher in energy you go.

Note also that the LMTO f orbitals are more efficient in converging the total energy per extra orbital added than plane waves are.

See Table of Contents

Additional exercises

1. Reduce the smooth Hankel basis set of Section 5 by eliminating the f orbitals. This reduces the basis set to 9 orbitals/atom. You should find that the total energy is ehf= -126808.271025 without --optbas, and ehf=-126808.280117 with it.

2. Add APWs to this basis, and observe that the total energy converges to the same value. Note also that the LMTO f orbitals are more efficient in converging the total energy than plane waves are. You should get something similar to the grey line in the Figure of Section 6.3.

3. Try an all LAPW basis : use pwmode=12. You should get something similar to the purple line in the Figure of Section 6.3. Results start to get reasonable around pwemax=4.5.

4. At self-consistency, the forces should reach convergence for given basis set, as described in the annotated lmf output. Forces should be derivatives of the total energy wrt nuclear displacements. As the basis set becomes complete, the forces should stop changing. Confirm that this is the case by varying the basis set. Bi₂Te₃ has 5 atoms; two of the Te atoms have one force equal and opposite, and the two Bi another force, also equal and opposite. You should find something similar to the following:

   Basis	F(Te)	F(Bi)
   1-kappa (optimized)	10.6	15.7
   2-kappa	7.5	13.1
   1-kappa(opt)+ APW(pwemax=2)	8.0	13.3
   1-kappa(opt)+ APW(pwemax=6)	8.0	13.4

See Table of Contents