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
 Tutorial
 Additional exercises
nano init.bi2te3
blm init.bi2te3 #makes template actrl.bi2te3 and site.bi2te3
nano actrl.bi2te3 #change to autobas[pnu=1 loc=1 lmto=3 mto=1 gw=0]
cp actrl.bi2te3 ctrl.bi2te3
Free atomic density and basis parameters
lmfa ctrl.bi2te3 #use lmfa to make basp file, atm file and to get gmax
nano basp0.bi2te3 #remove PZ=.. from file
cp basp0.bi2te3 basp.bi2te3 #copy basp0 to recognised basp prefix
... to be finished
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 SlaterKoster 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. (68) 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^{l1}e^{\kappa\,r}$, $\kappa^2{=}\varepsilon$) and the smoothing radius $r_s$ demarcates the approximate point of transition from Gaussianlike 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 closepacked. 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 selfconsistency 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 highlying local orbitals (one principal quantum number above the valence partial waves). These are usually less relevant in DFT; however, highlying 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:

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. 
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 highlying (far above the linearization energy) or deep, to include highlying (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.

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 allelectron 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://cstwww.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=R3M
UNITS=A
A={a} C={c}
SITE
ATOM=Te X=0 0 0
ATOM=Te X=0 0 {uTe}
ATOM=Bi X=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 planewave cutoff GMAX; see the PbTe tutorial.
 You must specify a valid k mesh for Brillouin zone integration; see the PbTe tutorial.
 You must define a basis set which can be done manually or automatically, as described next.
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 actrl.bi2te3 a little; the default gives [… LMTO=5 MTO=4], which makes a double kappa basis.
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.
Tokens in AUTOBAS tell lmfa to do the following:
 PNU=1 Calculate the logarithmic derivative parameter P_{l} for the free atom.
 Calculated parameters are saved in file basp0.ext as P=…. Nothing about P is written if PNU=0.
 LOC=1 Look for shallow cores to be explicitly treated as valence electrons, as local orbitals.
 Shallow cores that meet specific criteria are identified and written to basp0.ext as PZ=…. No search is made if LOC=0
 LMTO=3 Pick a default choice about the size of basis. LMTO=3 is a standard minimal basis.
 Run
lmfa input
and look for HAM_AUTOBAS_LMTO to see what other choices there are.lmfa will pick some defaults for the lcutoff, e.g. spd or spdf depending on the value of LMTO.
 MTO=1 Choose 1kappa basis set (single orbital per l channel).
 Small basis for fast calculations. For better quality calculations, it is recommended you use MTO=2 or MTO=4.
 GW=0 Create a setup for an LDA calculation.
 If GW=1, tailor basis for a GW calculation, selecting stricter criteria for including shallow cores as valence, and the size of basis.
These tokens thus define some reasonable default basis for you. lmfa writes basp0.ext. This file is never read, but lmf will read basp.ext and use this information when assembling the basis set. The two files have the same structure and you can copy basp0.ext to basp.ext. What lmfa generates is not cast in concrete. You are free to adjust the parameters to your liking, e.g. add a local orbital or remove one from the basis.
The AUTOBAS tokens tell lmf what to read from basp.ext. It uses tokens in a manner similar, but not identical, to the way lmfa uses them:
 PNU=1 Read parameters P for all species present in basp.ext.
 If PNU=0, these parameters will not be read.
 LOC=1 tells lmf to read local orbital parameters PZ.
 Since these parameters may also be specified by the input file,
LOC=1 tells lmf to give precedence to parameters specified by ctrl file
LOC=2 tells lmf to give precedence to parameters specified by basp.
LMTO= is not used by lmf.
 MTO=1 controls what part of RSMH and EH is read from basp.bi2te3.
 LMTO=1 or 3 tells lmf to read 1kappa parameters specified by basp
LMTO=2 or 4 tells lmf to read 2kappa parameters specified by basp
LMTO=1 or 2 tells lmf that parameters in the ctrl file take precedence
LMTO=2 or 4 tells lmf that parameters in the basp file take precedence  GW=0 tunes the basis for an LDA calculation
 If GW=1, tune basis for a GW calculation. For example log derivative parameters P are floated a little differently in the selfconsistency cycle. They are weighted to better represent unoccupied states, at a slight cost to their representation of occupied states.
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
By default, blm makes the spheres as large as possible without overlapping. In this case the Bi and Te radii are nearly the same.
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.
3. The atomic density and basis set parameters
If you did not do so already copy actrl.bi2te3 to ctrl.bi2te3 and change [… LMTO=4 MTO=4] → [… LMTO=3 MTO=1]).
Invoke lmfa:
$ lmfa bi2te3
The primary purpose of lmfa is to generate a free atom density. A secondary purpose is to supply additional information about the basis set in an automatic way. All of this information can be supplied manually in the input file, but the blm provides a minimum amount of information. lmfa generates basp0.bi2te3 which contains
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.936
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.
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, singlekappa 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:
AUTOBAS  lmfa writes, lmf reads 
MTO  RSMH,EH (and RSMH2,EH2 if double kappa basis) 
P  P 
LOC  PZ 
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:

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

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.
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 corelike states in the valence. This is the allelectron analog of the “hardness” of the pseudopotential in pseudopotential schemes. If you want highaccuracy 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 kpoint 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 kpoints. An appropriate choice will depend strongly on the context of the calculation and the sytem of interest; the densityofstates 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 selfconsistent calculation for that: the kconvergence 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 HarrisFoulkes 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.
5. Selfconsistent 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 selfconsistent 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 selfconsistency in 9 iterations.
The HarrisFoulkes and KohnSham total energies are ehf=126808.3137885 and ehk=126808.2492178 from the Mattheis construction. At selfconsistency they come together: ehf=126808.2950696 and ehk=126808.2950608 (the large integer part comes almost entirely from the core states.)
The selfconsistent value is 18 mRy less binding than the HarrisFoulkes energy of the Mattheis construction, and 46 mRy more binding than the corresponding KohnSham energy. That the two initial functionals bracket the selfconsistent 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 planewave expansions of envelope functions as shown this table:
sugcut: make orbitaldependent reciprocal vector cutoffs for tol=1.0e6
spec l rsm eh gmax last term cutoff
Te 0 1.61 0.89 4.603 3.31E06 1635*
Te 1 1.68 0.29 4.662 5.08E06 1635*
Te 2* 1.91 0.10 4.273 1.15E06 1539
Te 3 1.91 0.10 4.471 1.71E06 1635*
Bi 0 1.67 0.84 4.441 1.29E06 1635*
Bi 1 1.87 0.21 4.183 1.08E06 1411
Bi 2* 1.90 0.10 4.297 1.37E06 1539
Bi 3 1.90 0.10 4.497 2.06E06 1635*
gmax=4.4 isn’t quite large enough to push the error below tolerance (1.0e6), 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.
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^{in}−n in the guessed density from the selfconsistent 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 onecenter expansions in spherical harmonics Y_{lm}. LMXA is a cutoff that truncates the onecenter 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
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 onecenter 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
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 bruteforce minimization of the HarrisFoulkes energy.
Starting from the selfconsistent density generated previously do the following:
$ lmf ctrl.bi2te3 vgmax=4.4 vnkabc=3 optbas > out
This will loop through RSMH in each species, minimizing the total energy for each one by one. It does not vary EH because the total energy is less sensitive to it. lmf will print out a table of the parameters it will optimize:
LMFOPB: optimizing energy wrt 8 parameters, etol=5e4:
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
Each cycle carries out non selfconsistent calculations to get the HarrisFoulkes total energy. After completing lmf prints out
LMFOPB: before optimization ehf=126808.2951. After optimization ehf=126808.3092
Exit 0 Optimization complete. See file basp2
By tuning the parameters the energy decreased by 0.0141 Ry. To see which orbitals contributed the most, do:
$ grep LMFOPB out  grep optimized
You should see
LMFOPB: var #1 optimized to 1.602 (7 iter): ehf=126808.295097, delta ehf=2.38e5
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.77e4
LMFOPB: var #5 optimized to 1.605 (5 iter): ehf=126808.301, delta ehf=9.12e4
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.11e4
LMFOPB: var #8 optimized to 1.679 (6 iter): ehf=126808.309198, delta ehf=3.42e4
The orbitals that mad the most difference were the Te 5p (1.681 → 1.617) and Bi 5p (1.867 → 1.658).
Optimized parameters where the total energy decreased by more than etol (5e4) were modified; the optimized parameters are written to basp2.bi2te3
You can copy basp2.bi2te3 to basp.bi2te3, but there is little point in optimizing any but the Te 5p and Bi 6p. Instead we use a text editor and modify basp.bi2te3 in those channels:
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
Run lmf again
$ lmf ctrl.bi2te3 vgmax=4.4 vnkabc=3 quit=band
and confirm that ehf comes out close (126808.3086945) to the optimized value (126808.309198).
You can at this point make the density selfconsistent again.
Note: Before continuing, be sure to restore the original basp file.
In sum,
 without optimizing the basis, you should achieve total energy ehf=126808.2950696 at selfconsistency.
 With the optbas switch it should become ehf=126808.306962.
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 (doublekappa)
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
Run lmf to selfconsistency
$ rm f mixm.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 selfconsistency you should find that the total energy converges to ehf=126808.313403 Ry.
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=1d16
EWALD TOL=1d16
One problem with the PMT method is that the smooothHankel 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 doublekappa 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 selfconsistent doublekappa 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 2kappa 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 '(x2ehf)*1000' etot.bi2te3
The figure shows that the doublekappa basis (additional 45 orbitals) stabilizes the singlekappa spdf basis by about 18 mRy, and that the same can be accomplished with APWs with a planewave cutoff of 2 Ry (additional 4862 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: 343370 APWs.
Note that the APW basis is generally less efficient than the LMTO basis. To reach a precision comparable to the 2kappa basis (125 orbitals) starting from the 1kappa 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
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 kpoints. 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.
Additional exercises

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.

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.

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.

At selfconsistency, 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_{2}Te_{3} 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) 1kappa (optimized) 10.6 15.7 2kappa 7.5 13.1 1kappa(opt)+ APW(pwemax=2) 8.0 13.3 1kappa(opt)+ APW(pwemax=6) 8.0 13.4
Edit This Page