# Properties of the lmf basis set

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

Make an input file:

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 Slater-Koster form,

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; 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 Pl 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 l-cutoff, e.g. spd or spdf depending on the value of LMTO.

MTO=1  Choose 1-kappa 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 1-kappa parameters specified by basp
LMTO=2 or 4 tells lmf to read 2-kappa 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 self-consistency 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, 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:

 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: 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
$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 ninn 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 Ylm. 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

##### 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 Pkl 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 brute-force minimization of the Harris-Foulkes energy. Starting from the self-consistent 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=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


Each cycle carries out non self-consistent calculations to get the Harris-Foulkes 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.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  The orbitals that mad the most difference were the Te 5p (1.6811.617) and Bi 5p (1.8671.658). Optimized parameters where the total energy decreased by more than etol (5e-4) 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 self-consistent 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 self-consistency.
• 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]


    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  Run lmf to self-consistency $ 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 self-consistency 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=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 ei(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 Emax=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  #### 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.