Annotated standard output, program lmf
Table of Contents
 Table of Contents
 Preliminaries
 Preprocessor’s transformation of the input file
 Display tags parsed in the input file
 Reading basis information from the basp file
 Header information
 Lattice information
 Symmetry and k mesh
 Augmentation parameters
 Interstitial mesh
 Counting the size of the basis
 Envelope function parameters and their G cutoffs
 Obtain an input density
 Potential and Matrix Elements
 Reading QSGW selfenergies+
 Parameters for local orbitals
 Brillouin zone integration
 Output density and update of augmentation sphere boundary conditions
 Total energy
 Forces
 New density
 End of selfconsistency loop
 Other Resources
Preliminaries
The output documented here is mostly taken from the lmf tutorial for PbTe. Some portions are adapted from other calculations, as will be indicated.
The standard output is organized by blocks. The sections below annotate each block approximately in the order they are made.
Preprocessor’s transformation of the input file
The input file is run through the preprocessor, which modifies the ctrl file before it is parsed for tags. Normally it does this silently. To see the effects of the preprocessor use lmf showp ...
Append showp
to the lmf command in the PbTe tutorial
$ lmf ctrl.pbte vnkabc=6 vgmax=7.8 showp
The box below compares side by side the original ctrl.pbte and its transformation by the preprocessor (the original file was edited slightly)
% const nit=10
% const met=5
% const so=0 nsp=so?2:1
% const lxcf=2 lxcf1=0 lxcf2=0
% const pwmode=0 pwemax=3
% const nkabc=0 gmax=0
VERS LM:7 FP:7 # ASA:7  VERS LM:7 FP:7
IO SHOW=f HELP=f IACTIV=f VERBOS=35,35  IO SHOW=f HELP=f IACTIV=f VERBOS=35,35
EXPRESS  EXPRESS
# Lattice vectors and site positions 
file= site  file= site

# Basis set 
gmax= {gmax}  gmax= 7.8
autobas[pnu=1 loc=1 lmto=5 mto=4 gw=0]  autobas[pnu=1 loc=1 lmto=5 mto=4 gw=0]

# Selfconsistency 
nit= {nit}  nit= 10
mix= B2,b=.3,k=7  mix= B2,b=.3,k=7
conv= 1e5  conv= 1e5
convc= 3e5  convc= 3e5

# Brillouin zone 
nkabc= {nkabc}  nkabc= 6
metal= {met}  metal= 5

# Potential 
nspin= {nsp}  nspin= 1
so= {so}  so= 0
xcfun= {lxcf},{lxcf1},{lxcf2}  xcfun= 2,0,0

#SYMGRP i*r3(1,1,1) r4x 
HAM  HAM
PWMODE={pwmode} PWEMIN=0 PWEMAX={pwemax}  PWMODE=0 PWEMIN=0 PWEMAX=3
FORCES={so==0} ELIND=0.7  FORCES=1 ELIND=0.7
SPEC  SPEC
ATOM=Pb Z= 82 R= 3.044814 LMX=3 LMXA=4  ATOM=Pb Z= 82 R= 3.044814 LMX=3 LMXA=4
ATOM=Te Z= 52 R= 3.028689 LMX=3 LMXA=3  ATOM=Te Z= 52 R= 3.028689 LMX=3 LMXA=3
Display tags parsed in the input file
After transformation by the preprocessor, lmf parses for tags and substitutes default values for tags it does not find. To see the value of tags lmf, whether parsed or defaults, use lmf show
or lmf show=2
. The latter causes lmf to stop after displaying tags, and is useful if you want to see whether lmf is doing what you expect. Using show
is useful if you want to record the input conditions in the output (be advised that the output is verbose).
Add show=2
to the lmf command from the PbTe tutorial:
$ lmf ctrl.pbte vnkabc=6 vgmax=7.8 show=2
The output is quite verbose so only a snippet from the SPEC category is shown
Tag Input cast (size,min,read,def) result
...
 Parameters for species data 
SPEC_SCLWSR opt r8 1, 1, *, 1 0
SPEC_OMAX1 opt r8v 3, 1, *, 3 0 0 0
SPEC_HFAC opt chr 1, 0, *
SPEC_HFAC_V opt lg 1, 1, *, 1 F
... Species 1
SPEC_ATOM reqd chr 1, 0, 1 Pb
SPEC_ATOM_Z reqd r8 1, 1, 1,  82
SPEC_ATOM_R reqd* r8 1, 1, 1,  3.04481
SPEC_ATOM_A opt r8 1, 1, *, 1 0.025
SPEC_ATOM_NR opt i4 1, 1, *, 1 497
SPEC_ATOM_RSMH opt r8v 4, 4, *, 4 0 0 0 0
SPEC_ATOM_EH opt r8v 4, 4, *, 
...
 Tags are optional except for SPEC_ATOM, SPEC_ATOM_Z, and SPEC_ATOM_R; the latter could have been supplied with equivalent information through a different tag (SPEC_ATOM_R/W or **SPEC_ATOM_R/A) in this case.
 Variables are read as scalars or vectors; integers, floatingpoint numbers, or strings. If strings have spaces you must enclose them in quotes or backets.
 Default values were substituted for the optional tags. In this case, the input file contained none of the optional tags: all values are taken from defaults, except for EH which was not assigned any value. lmf requires EH; but it will be read independently from the basp file.
There are two other special modes to help with managing the input.
$ lmf input
doesn’t attempt to read anything; instead, it prints out a (large) table of all the tags it would try to read, including a brief description of the tag, and then exits. See here for further description.
$ lmf help
performs a similar function for the command line arguments: it prints out a brief summary of arguments effective in the executable you are using. See the tutorial for further description.
Reading basis information from the basp file
After parsing the ctrl file, lmf may attempt to read basis set information from the basp file.
The basp file is automatically generated by lmfa. Tokens in EXPRESS_autobas or HAM_AUTOBAS control what is read from the basp file.
basp.pbte supplies basis information (parameters EH and RSMH defining the shape of the envelope functions, continuous principal quantum numbers P and information about local orbitals). The structure of the basp file is described on this page. When lmf reads from this file, it will print to standard output what it read :
rdctrl: reading basis parameters from file basp
ioorbp: read species Pb RSMH,EH RSMH2,EH2 P PZ
ioorbp: read species Te RSMH,EH RSMH2,EH2 P
reset nkaph from 1 to 3
You can also supply this information in the ctrl file. If data is present in both the ctrl and basp files, lmf decides on which to use depending on settings in EXPRESS_autobas or HAM_AUTOBAS.
To see a synopsis of tokens in AUTOBAS, invoke lmf input
and search for autobas. Tokens in AUTOBAS are documented here.
Header information
The header information presents a condensed synopsis of some key settings that are used in the calculation.
LMF: nbas = 2 nspec = 2 vn 7.11.i verb 35
special: forces
pot: XC:BH
float: float P LDAstyle
autoread: mto basis(4), pz(1), pnu(1)
bz: metal(5), tetra, invit
It tells you, for example that you are using the BarthHedin functional
Lattice information
This block prints informations about the lattice vectors and settings used in Ewald summations:
Plat Qlat
0.000000 0.500000 0.500000 1.000000 1.000000 1.000000
0.500000 0.000000 0.500000 1.000000 1.000000 1.000000
0.500000 0.500000 0.000000 1.000000 1.000000 1.000000
alat = 12.147006 Cell vol = 448.071898
LATTC: as= 2.000 tol=1.00E08 alat=12.14701 awald= 0.261
r1= 1.807 nkd= 87 q1= 5.403 nkg= 169
Note: When long, thin cells are used, or when APW’s are added to the basis set, some attention needs to be paid to the Ewald tolerance, input through EWALD_TOL.
Symmetry and k mesh
The block below shows symmetry operations it finds in the crystal, and the irreducible k mesh it obtains from the point group it is given:
SGROUP: 1 symmetry operations from 0 generators
SYMLAT: Bravais system is cubic with 48 symmetry operations.
SYMCRY: crystal invariant under 48 symmetry operations for tol=1e5
GROUPG: the following are sufficient to generate the space group:
i*r3(1,1,1) r4x
i*r3(1,1,1) r4x
MKSYM: found 48 space group operations ... includes inversion
BZMESH: 16 irreducible QP from 216 ( 6 6 6 ) shift= F F F
TETIRR: sorting 1296 tetrahedra ... 35 inequivalent ones found
Notes: (see also Additional Exercises).
The group operations were determined automatically from the given lattice. PbTe is cubic, with 48 group operations.
First the crystal system is determined; then the symmetry operations inconsistent with the basis are discarded.
Finally the 48 operations are distilled into a minimum set of generators (i*r3(1,1,1) r4x) that generate the entire space group.Generators comprise a minimum set of group operators. All possible products of them form a closed set of operations, which forms a group.
In this case there are no translations; all the group operators are pure point group operations. This is not true in general. For example, hcp Co has 24 space group operations. The generators defining its space group can be written as:
i*r3z::(1/3,1/3,1/2) r2z::(1/3,1/3,1/2) r2x
 You can also specify symmetry operations manually via the SYMGRP category. This is particularly useful when magnetic symmetry must be considered.
 The k mesh is specifed through the number of divisions along each of the three reciprocal lattice vectors, tag EXPRESS_nabc or BZ_NABC. In this case the 216 k points generated by a 6×6×6 mesh gets reduced by symmetry to 16 inequivalent points.
You can also specify whether the kmesh should pass through the origin or straddle it (see BZ_BZJOB).  Bloechl’s generalized tetrahedron method is used to perform Brillouin zone integrations. You can also use the generalized MethfesselPaxton sampling integration scheme or a sampling with a Fermi function, using BZ_N=1.
Augmentation parameters
The table below contains a synopsis of key parameters associated with augmentation spheres.
species data: augmentation density
spec rmt rsma lmxa kmxa lmxl rg rsmv kmxv foca rfoca
Pb 3.045 1.218 4 3 4 0.761 1.522 15 1 1.218
Te 3.029 1.211 3 3 3 0.757 1.514 15 1 1.211
 rmt: the augmentation radius
 rsma and kmxa: smoothing radius and polynomial order used to expand envelope functions around other sites.
 lmxa: the lcutoff for the augmentation of the smoothed Hankel envelope functions by partial waves. Because of the unique way augmentation is done in this method, properties converge much faster with lmxa than in standard augmented wave methods.
 lmxl: analogous to lmxa, but it controls the lcutoff of the charge density. lmxl defaults to lmxa; you can often make it smaller with minimal loss of accuracy. There is little efficiency gain, however, so in practice this feature is rarely used.
 rg: is concerned with adding local gaussian pseudocharges to manage the Hartree potential.
 rsmv, kmxv: concerned projections of mesh density onto local densities.
 foca, rfoca: allow for differing treatments of the core.
Interstitial mesh
The following block is concerned with the mesh used to represent the charge density, and to evaluate matrix elements of the (unaugmented) envelope functions. The spacing of the mesh is controlled by the G cutoff (7.8 for PbTe).
MSHSIZ: mesh has 18 x 18 x 18 divisions; length 0.477, 0.477, 0.477
generated from gmax = 7.8 a.u. : 3647 vectors of 5832 (62%)
GVLIST: gmax = 7.8 a.u. created 3647 vectors of 5832 (62%)
mesh has 18 x 18 x 18 divisions; length 0.477, 0.477, 0.477
SGVSYM: 126 symmetry stars found for 3647 reciprocal lattice vectors
Information about whether this mesh is sufficiently accurate is given in the table beginning with sugcut below.
Counting the size of the basis
How the basis is defined is described on this page.
In the table below, the size of the basis is presented. lmf doesn’t have downfolding capability, so the important numbers are those in the “low” column. Rows 1,2,3 indicate how many orbitals are connected respectively with the first Hankel envelope (EH), the second envelope (EH2), and local orbitals. The total basis (and hamiltonian rank) consists of 55 orbitals.
Makidx: hamiltonian dimensions Low, Int, High, Negl: 55 0 32 63
kappa Low Int High L+I L+I+H Neglected
1 32 0 9 32 41 9
2 18 0 23 18 41 9
3 5 0 0 5 5 45
all 55 0 32 55 87 63
suham : s 41 augmentation channels, 41 local potential channels Maximum lmxa=4
The last line refers to augmentation channels. An envelope of a particular l must be expanded around remote sites. The lcutoff for expanding tails of envelope functions centered elsewhere is lmxa, input though tag opeSPEC_ATOM_LMXA. For lmf, lmxa is usually reasonably converged if it is 3 for sp systems, 4 for d atoms and 6 for f shell elements.
Envelope function parameters and their G cutoffs
In the table envelope function parameters for each species is given, which defines their shape. Also shown is :
 gmax : the G cutoff that represents that envelope to the given tolerance
 last term : an estimate for the precision of the planewave repesentation
 cutoff : the number of plane waves with G<gmax.
sugcut: make orbitaldependent reciprocal vector cutoffs for tol=1.0e6
spec l rsm eh gmax last term cutoff
Pb 0 1.80 0.10 4.123 2.68E06 531
Pb 1 2.02 0.10 3.848 2.42E06 411
Pb 2* 2.03 0.10 4.013 1.37E06 531
Pb 3 2.03 0.10 4.193 1.54E06 537
Pb 0 1.80 0.90 4.123 2.68E06 531
Pb 1 2.02 0.90 3.848 2.42E06 411
Pb 2 2.03 0.90 4.013 1.37E06 531
Te 0 1.63 0.10 4.572 1.45E06 725
Te 1 1.71 0.10 4.575 1.52E06 725
Te 2* 2.02 0.10 4.037 1.63E06 531
Te 3 2.02 0.10 4.218 1.87E06 537
Te 0 1.63 0.90 4.572 1.45E06 725
Te 1 1.71 0.90 4.575 1.52E06 725
Te 2 2.02 0.90 4.037 1.63E06 531
Each envelope function must be expanded in plane waves in order to
 assemble matrix elements of the interstitial potential
 assemble the output charge density.
Both are assembled in reciprocal space. The number of plane waves needed for a particular orbital depends on how sharply peaked the function is, so the cutoff is orbitaldependent to allow for faster execution. gmax of any one orbital may safely be less than the global G cutoff (7.8 for PbTe); if it can, lmf will find a gmax for each orbital that meets a preset tolerance (10^{−6} in this case).
Otherwise it uses all the G vectors available to it, and appends a ‘*’ to the number indicating not enough orbitals are available to meet the specified tolerance. This is flagged by an asterisk appended to the cutoff, e.g.
Te 0 1.63 0.10 4.572 1.74E06 701*
If ‘last term’ is too high you can expect errors in the hamiltonian and you should increase gmax.
The tolerance defaults to 10^{−6}, but you can control it with tag HAM_TOL. 10^{−5} or smaller is usually safe.
Note: With APW’s also in the basis, it the overlap matrix becomes nearly singular. To stabilize the overlap matrix, you can set HAM_TOL to something close to machine precision, e.g. 10^{−16}.
At this stage the potential independent setup is complete.
Obtain an input density
In KohnSham theory, the potential is a functional of the density. lmf tries to read the density from the density restart file, rst.pbte.
Initially the density is not known and rst.pbte does not exist.
lmfp : no rst file ... try to overlap atomic densities
When lmf cannot read the density, it constructs a trial density by overlapping free atom densities. The true density is obtained itertively through the selfconsistency cycle.
lmf reads atomic densities from atm.pbte and overlaps them (Mattheis construction). This makes a reasonable guess for the density: the atomic potential is very deep, and the crystal potential is a relatively modest perturbation to it. Thus the true density is typically not so far from a superposition of atomic densities.
This table indicates what lmf is doing, and checks that atm.pbte is consistent the input file.
lmfp : no rst file ... try to overlap atomic densities
rdovfa: read and overlap freeatom densities (mesh density) ...
rdovfa: expected Pb, read Pb with rmt= 3.0448 mesh 497 0.025
rdovfa: expected Te, read Te with rmt= 3.0287 mesh 461 0.025
Next follows a table with some useful information about the overlapping process
ovlpfa: overlap smooth part of FA densities
site 1 spec 1 pos 0.0000 0.0000 0.0000 Qsmooth 7.659845
site 2 spec 2 pos 0.0000 0.0000 0.5000 Qsmooth 7.820498
total smooth Q = 15.480342
Free atom and overlapped crystal site charges:
ib true(FA) smooth(FA) true(OV) smooth(OV) local
1 12.535727 6.195571 12.881552 6.541397 6.340155
2 4.566605 6.387103 4.926059 6.746557 1.820498
Smooth charge on mesh: 15.480342
Sum of local charges: 4.519657
Total valence charge: 20.000000
Sum of core charges: 114.000000
Sum of nuclear charges: 134.000000
Homogeneous background: 0.000000
Deviation from neutrality: 0.000000
It is important that the system is charge neutral (last line).
If, on the other hand, the charge density has been saved in a restart file from a prior iteration lmf reads it. Usually restart files are saved in a binary form in a file rst.ext.
When lmf reads from a restart file it prints out a message similar to this:
iors : read restart file (binary, mesh density)
use from restart file: ef window, positions, pnu
ignore in restart file: *
The restart file contains:
 valence and core densities inside augmentation spheres for each nucleus.
 mesh density
 lattice vectors and site positions
 Information about the Fermi level and (in the sampling case) parameters controlling the Brillouin zone integration
 logarithmic derivative parameters pnu — boundary conditions for partial waves
Controlling what is read from the restart file
By default, parameters 35 will be read from rst.pbte and supersede preexisting values. In this tutorial lattice vectors and site positions are read from site.pbte, Brillouin zone integration parameters from category BZ in ctrl.pbte, pnu from basp.pbte.
Through commandline switch rs you can control what lmf reads from rst.pbte and what it keeps from input files. rs has five arguments: rs=#1,#2,#3,#4,#5; you can supply one to five numbers and default values are taken for the rest.
#1 and #2 mark whether and how to read and write the restart file. Use the remaining switches to suppress reading from the restart file:
#3=1 keep original site positions. This is necessary if you shift the atoms but retain this rst.ext
#4=1 keep original sampling parameters.
#5=1 keep original pnu.
Notes:
 If you use #5=1, a selfconsistent density will no longer be fully selfconsistent because the augmentation part of the basis set has changed. lmf automatically floats pnu to the band center of gravity, to make the basis set more accurate.
 The mesh density has a lump centered at each atom. If the site positions change, the lump will be centered in the wrong place, and the density very poor. lmf can partially compensate for this by adding an atomic density centered at the new position, and subtracting it at the file position. It will do this if you use #1=11. But if the shifts are large, it is better to discard the restart file and overlap atomic densities again.
 If the lattice has been rotated or sheared, the reader will detect this and will rotate or shear the mesh density accordingly. But the site densities are left unaltered, unless you tell the reader to rotate them. To do this use #1=101.
 If you are doing sampling integration and want to change sampling parameters, you must use #4=1.
Potential and Matrix Elements
Routine (fp/mkpot.f) makes the potential and relevant matrix elements.
LDA functionals
It was indicated in the header that you are using the LDA with the BarthHedin functional. This is controlled by EXPRESS_xcfun, which is an alias for HAM_XCFUN.. To perform a GGA calculation with Kieron Burke’s PBE functional, remove this line in ctrl.pbte
xcfun= {lxcf},{lxcf1},{lxcf2} # set lxcf=0 for libxc functionals
and add this to the HAM category:
XCFUN=3 GGA=3
When you run lmfa or lmf the header should now read
pot: XC:PW91+PBE(gga)
Note: you should start the calculation from the beginning with this change, since the atomic densities are made from a BarthHedin functional. The core densities made by lmfa are frozen and will carry over into the lmf calculation.
The smooth part of the potential is made first. This is necessary to determine the electrostatic potential on the augmentation sphere boundaries. It serves as a boundary condition for the solving Poisson’s equation inside each sphere. The output reads
Average es pot at rmt = 0.473380 avg sphere pot = 0.516004 vconst = 0.473380
smooth rhoves 41.429795 charge 15.480342
smooth rhoeps = 11.301007 rhomu = 14.720287 avg vxc = 0.680712
These few lines conceal a somewhat involved process. For the mesh density to yield the correct potential in the interstial, lumps of charge must be added to the mesh density. Provided the additional charge has multipole moments matching those obtained by the difference in the two local densities , multipoles of the (modified) mesh density will match those of the true density, and the electrostatic potential generated by it will be correct at each augmentation radius.
Run lmf with a high verbosity (pr50 quit=pot) and the following appears:
rhomom: ib ilm qmom Qval Qc Z
1 1 1.954651 6.929056 68.00 82.00
1 21 0.018533
1 25 0.015663
2 1 0.307248 1.089167 46.00 52.00
Smooth charges: Qmesh = 14.160111 Qgauss = 5.839889 corenuc = 20 tot = 0
lmf has found the multipole moments from the local densities ( is the charge, the l=0 multipole moment is ). It adds localized gaussian charges with smoothing radius rg. (Notice that rg is quite small, to ensure that the gaussian does not spill out of the augmentation sphere.)
Now the electrostatic potential at the augmentation radius can be constructed
site class ilm vval ves(rmax)
1 1 1 1.623863 0.458083
21 0.059347
25 0.050158
2 2 1 1.767044 0.498474
lmf shifts the total mesh potential by a constant so that the average potential on augmentation boundaries is close to zero:
Average es pot at rmt = 0.478171 avg sphere pot = 0.518433 vconst = 0.478171
Average es pot at MT boundaries after vconst shift
Site ves Site ves
1 0.020088  2 0.020302 
It is done order to put the Fermi level close to zero. This shift can be chosen differently; if you use oldvc the cellaveraged potential is set to zero instead.
From the local densities and electrostatic boundary conditions, the local (site) potentials can be constructed along with matrix elements of partial waves. Some limited output is shown in this table:
locpot:
site 1 z= 82.0 rmt= 3.04481 nr=497 a=0.025 nlml=25 rg=0.761 Vfloat=T
sm core charge = 0.029706 (sphere) + 0.000318 (spillout) = 0.030024
potential shift to crystal energy zero: 0.000065
site 2 z= 52.0 rmt= 3.02869 nr=461 a=0.025 nlml=16 rg=0.757 Vfloat=T
sm core charge = 0.32823 (sphere) + 0.005658 (spillout) = 0.333888
potential shift to crystal energy zero: 0.000023
Note the printout Vfloat=T. This indicates that the potential used integrate the partial waves is being updated to the current local potential. The potential used to make partial waves and is kept independently of the potential generated by the density. Initially this potential is taken from of the free atom; and normally it is overwritten (Vfloat=T). If you tell lmf to freeze it (HAM_FRZWF=t globally or SPEC_ATOM_FRZWF=f for one species) the basis set will not adapt to the potential.
Once again, many steps are concealed behind this limited output. You can get more information, e.g. matrix elements, using a higher print verbosity: the higher the verbosity, the more is output.
Run lmf with (pr70 quit=pot) and a great deal of information is printed out. A table of classical LMTO potential parameters is made. Particularly useful are the indicating the band center (c) and bandwith (srdel^{2}), and linearization energy enu:
l enu v c srdel qpar ppar
0 0.685722 1.220936 0.974210 0.311410 2.5442 2.037524
1 0.345247 0.961024 0.034414 0.322661 9.5614 3.768891
2 0.215079 0.351057 1.877979 0.404577 13.6181 3.646179
2(sc) 1.355835 1.700376 1.373127 0.018377 60.8790 0.594481
3 0.367178 0.903272 2.382595 0.347763 27.1688 7.102666
4 0.261796 0.813619 4.861503 0.414942 32.9588 13.704206
The following table presents (hamiltonian, potential, overlap) matrix elements of special linear combinations val and slo of and that have the property that value or slope vanishes at the MT boundary:
l valval valslo sloval sloslo
0 s 10.363095 7.303401 7.303401 7.075456
h 11.973814 0.299694 8.971198 0.531008
v 67.315135 84.761799 84.761799 111.419064
1 s 8.042350 4.234108 4.234108 3.159956
...
A number of other quantities are made, depending on print verbosity and computational conditions.
For verbosities above 50, electric field gradients are printed:
Electric Field Gradient:
Site Tensor axes esu/cm^2 V/m^2 eta line splitting
x10^13 x10^19 (mm/s)
1 0.2920.936 0.198 0.00 0.00 0.000 0.00000
0.771 0.1090.627 0.00 0.00
0.565 0.336 0.754 0.00 0.00
2 0.2920.936 0.198 0.00 0.00 0.000 0.00000
0.771 0.1090.627 0.00 0.00
0.565 0.336 0.754 0.00 0.00
These are directly measurable quantities. At least in , the electric field gradient on the site (measured to be ), is not well predicted by LDA, or LDA+U. QSGW, however, does quite well. See Fig.13 of this paper.
If spinorbit coupling is set (HAM_SO=t), matrix elements of are generated. Working from the Fe tutorial, invoke
$ lmf ctrl.fe vso=1 iactiv quit=pot pr55
These matrix elements of are printed out:
soprm: matrix elements for perturbation from spinorbit coupling
l <phi  phi> <dot  phi> <dot  dot>
up up 1 0.01903736 0.00053774 0.00001694
down down 1 0.01947919 0.00036545 0.00000870
up down 1 0.01925699 0.00036153 0.00001201
down up 1 0.01925699 0.00054370 0.00001201
up up 2 0.00559381 0.00229625 0.00097759
...
After completing the loop over sites a table is written for terms involving the total energy
Energy terms: smooth local total
rhoval*vef 36.158050 282.028989 318.187039
...
valence chg 14.160111 5.839889 20.000000
core charge 114.000000 0.000000 114.000000
and the following lines are written:
Charges: valence 20.00000 cores 114.00000 nucleii 134.00000
hom background 0.00000 deviation from neutrality: 0.00000
Is is important that the system be neutral.
Reading QSGW selfenergies+
If you have performed a QSGW calculation and read a selfenergy from file sigm.ext (or sigmrs.ext), will be added to the LDA potential.
Note: sigm.ext contains the difference, , so it is added directly to the KohnSham potential .
You will see an indication is read from the following standard output. It was taken from the Fe tutorial.
RDSIGM: read file sigm and create COMPLEX sigma(R) by FT ...
Sigm will be approximated by: diagonal Sigma for high and low states
Approximate sigma for energies E(lda)>2
For high states read constant sigii from sigm file
sigm file has 29 irreducible QP: nk = ( 8 8 8 ) shift= F F F
average SE read from sigm file: 0.122871 0.138308
hft2rs: make neighbor table for r.s. hamiltonian using range = 5 * alat
Notes:
 How lmf reads Σ^{0} and carries out an inverse Bloch sum is described in some detail in Section IIG of this Physical Review B paper.
 The cutoff E(lda)>2 is controlled by HAM_SIGP_EMAX. It affects the partitioning between calculated matrix elements of Σ^{0} and and the fixed ones above emax, as explained at the top of p165107 in the PRB paper (EMAX corresponds to E_{xccut2}).
 The average selfenergies (0.122871 spinup and 0.138308 spindown) above 2 Ry. When SIGP_MODE=4 (recommended), these numbers are used for Σ_{ii}^{0} for energies E > E_{F}+EMAX. In this mode the code prints out the line containing read constant sigii from sigm file.
 Range 5 for inverse Blochsummed Σ^{0}, in units of alat. If this number is not sufficiently large, the inverse Bloch transform is incomplete, and the program may stop with an error message.
This table:
Irr. qp for which sigma is calculated ...
BZMESH: 29 irreducible QP from 512 ( 8 8 8 ) shift= F F F
indicates the Brillouin zone mesh for Σ^{0}.
This table:
hft2rs created hrs: ndhrs=30 max Re(hrs) = 0.556 max Im(hrs) = 3.93e10
...
check FT s(R) against s(q) ... maximum error = 2.3e15 < tol (5e6)
...
check FT s(R) against s(q) ... maximum error = 7.8e15
...
indicates the precision in the construction of the inverse Bloch sum. It:
 checks the forward Bloch sum of the inverse Bloch summed Σ^{0} (2.3e15) to ensure that the error is less than HAM_RSSTOL. The combined inverse/forward sum should recover the original Σ^{0} to machine precision unless the inverse Bloch transform is inexact, as described above.
 Shows deviation (7.8e15) from original Σ^{0} after symmetrization by the group operations.
Parameters for local orbitals
Extended local orbitals must attach a smooth Hankel tail. A single smooth hankel function that fits both value and slope is found. The result displayed in a table:
Fit local orbitals to sm hankels, species Pb, rmt=3.044814
l type Pnu Eval K.E. Rsm Eh Q(r>rmt) Fit K.E.
2 low 5.934 1.568604 0.810887 1.05495 1.12869 0.02566 0.810887
Brillouin zone integration
The charge density and Fermi level must be found by integration over the Brillouin zone. The code makes two passes over the irreducible k points. The first is needed to obtain the Fermi level. E_{F} and kpoint weights are saved into a binary file wkp.pbte.
Start first of two band passes ...
bndfp: kpt 1 of 16, k= 0.00000 0.00000 0.00000
1.3108 1.3108 1.3076 1.3076 1.3076 1.0419 0.5426 0.2498 0.2498
0.2498 0.1728 0.1728 0.1728 0.2942 0.2942 0.2942
bndfp: kpt 11 of 16, k= 0.16667 0.16667 0.83333
1.3161 1.3103 1.3081 1.3058 1.3056 0.9310 0.6887 0.4314 0.3314
0.3190 0.1376 0.1914 0.2171 0.2802 0.5842 0.6122
BZWTS :  Tetrahedron Integration 
... only filled or empty bands encountered: ev=0.138629 ec=0.080593
VBmax = 0.138629 CBmin = 0.080593 gap = 0.058036 Ry = 0.78930 eV
BZINTS: Fermi energy: 0.138629; 20.000000 electrons; D(Ef): 0.000
Sum occ. bands: 18.1376921 incl. Bloechl correction: 0.000000
Saved qp weights ...
You can get more information by increasing the print verbosity.
PbTe is an insulator, both in real life and in the LDA. lmf automatically detects what appears to be an insulating state, in which case it prints out the highest occupied state and lowest unoccupied state it encountered and the “bandgap.”
Note 1: the “bandgap” is only the actual bandgap if the actual band edges coincide with one of the mesh of discrete k points used for integration. It is true in PbTe (both band edges fall at L) but not other cases, e.g. Si. See this tutorial.
Note 2: the band code may incorrectly determine that the system is an insulator. This can happen when band gaps are small and the k mesh is too coarse.
In a metal the situation is more complicated. The following table was extracted from the Fe tutorial (selfconsistent LDA level)
BZWTS :  Tetrahedron Integration 
...
BZINTS: Fermi energy: 0.000601; 8.000000 electrons; D(Ef): 14.389
Sum occ. bands: 1.3036417 incl. Bloechl correction: 0.001188
Mag. moment: 2.200354
The BZINTS table contains significant information, including the Fermi level (0.000601 Ry) and densityofstates there (14.389/Ry), and the singleparticle sum 1.3036417 Ry entering into the LDA total energy. The LDA magnetic moment/cell (2.200354 μ_{B}) is close to the experimental value. The Bloechl correction (0.001001 Ry) indicates how much the band sum changes when Bloechl weights are used in lieu of weights from the standard tetrahedron method. This is a measure of the convergence of the k mesh.
Integration weights and the METAL switch
Numerical quadrature is used to accumulate the output density or any property integrated over the Brillouin zone.
In insulators, each point in the full Brillouin zone gets equal weight; each point in the irreducible zone is weighted by the number of points in the full zone mapped to it. You can see these weights by running lmf at a high verbosity.
In metals, how the weights are made depends on the quadrature. See this page for a detailed discussion of the tetrahedron and sampling methods.
Add pr50 quit=ham to your invocation of lmf. You should see this table:
BZMESH: 16 irreducible QP from 216 ( 6 6 6 ) shift= F F F
Qx Qy Qz Multiplicity Weight
1 0.000000 0.000000 0.000000 1 0.009259
2 0.166667 0.166667 0.166667 8 0.074074
3 0.333333 0.333333 0.333333 8 0.074074
...
16 0.000000 0.333333 1.000000 12 0.111111
The weight of the first point (0.009259) is 2/216. The factor of 2 is for spin. In spin polarized calculations this factor gets reduced to 1.
If you know that the system is an insulator in advance, you can tell lmf to assume an insulator: in the input file replace % const met=5 with % const met=0. This will set tag EXPRESS_metal (an alias for BZ_METAL) to zero. In this mode it can accumulate the density on the first pass as well.
In metals, the weights depend on the Fermi level , which must be determined from the eigenvalues. Quadratures are of two types:
 the weight at k does not depend on values of neighboring k (sampling integration)
 the weight at k does depend on values of neighboring k (tetrahedron integration).
Questaal has several strategies to solve this chickenandegg problem, which you control through the BZ_METAL tag. The numbers below correspond to the values of BZ_METAL.
 System assumed to be an insulator; weights determined a priori. Only one pass is required.
 Eigenvectors are written to disk, in which case the integration for the charge density can be deferred until all the bands are obtained. (Note: This option is available only for the ASA. When accumulating just the spherical part of the charge, eigenvector data can be contracted over m, reducing memory demands.)
 Integration weights are read from file wkp.ext, which will have been generated in a prior band pass, e.g. from a prior cycle in the selfconsistent iteration. If wkp.ext is unavailable, the program will temporarily switch to METAL=3.
 Two band passes are made; the first generates only eigenvalues to determine E_{F}. It is slower than METAL=2, but it is more stable which can be important in difficult cases.
 Three distinct Fermi levels are assumed and weights generated for each. After E_{F} is determined, the actual weights are calculated by quadratic interpolation through the three points. The three Fermi levels are gradually narrowed to a small window around the true E_{F} in the course of the selfconsistency cycle. This scheme works for sampling integration where the kpoint weights depend on E_{F} only and not eigenvalues at neighboring k. You can also use this scheme in a mixed tetrahedron/sampling method: Weights for the band sum are determined by tetrahedron, but the charge density is integrated by sampling. It works better than straight sampling because the total energy is variational in the density.
 Like METAL=3 in that two passes are made but eigenvectors are kept in memory so there is no additional overhead in making the first pass. This is the recommended mode for lmf unless you are working with a large system and need to conserve memory.
The ASA implements METAL=0,1,2; the FP codes METAL=0,2,3,4,5. See Additional Exercises in the Fe tutorial.
Tetrahedron and sampling methods are explained and compared in some detail here.
Output density and update of augmentation sphere boundary conditions
A second pass over the irreducible kpoints is made with output similar to the first pass. After completion the output density is made and symmetrized. This table tells you how many electrons are inside each augmentation sphere:
mkrout: Qtrue sm,loc local
1 12.610810 6.902591 5.708219
2 5.275210 6.408075 1.132866
Then the logarithmic derivative parameters pnu are updated. Whether they are frozen or allowed to float depends on SPEC_ATOM_IDMOD. In this tutorial they are all allowed to float.
Make new boundary conditions for phi,phidot..
site 1 species 1:Pb
l idmod ql ebar pold ptry pfree pnew
0 0 1.700379 0.652781 6.895000 6.901582 6.500000 6.901582
1 0 0.846967 0.353483 6.809000 6.620518 6.250000 6.620518
2 0 0.121211 0.338472 6.250000 6.143261 6.147584 6.250000
2 0 sc 1.306275 5.934000 5.938843 5.147584 15.938843
3 0 0.032691 0.376447 5.183000 5.125635 5.102416 5.125635
...
The fractional parts of P is of order 0.9 for the 6s and 5d states, because they are deep. The 6d state is far above E_{F}, so its P is close to the free electron value. (For numerical stability, P is frozen for any l for which there is a corresponding deep local orbital. Thus for that state pnew is not updated to ptry.) The 6p state is near the Fermi level, and forms the predominant contribution to the bonding and band structure. See this page for an interpretation.
Total energy
The HarrisFoulkes energy is a functional of the input density and singleparticle sum. Information is printed out in the table below.
Harris energy:
...
rhoeps= 1094.806758 utot= 116742.136483 ehar= 55318.170567
The KohnSham total energy is also calculated:
KohnSham energy:
sumtv= 299.045855 sumtc= 62219.476115 ekin= 62518.521970
rhoep= 1094.805096 utot= 116741.865633 ehks= 55318.148758
This energy is a functional of n^{out} and V^{in}. It is a different functional than the HF functional, but the two should come together at the selfconsistent density.
Forces
In the same time total energies are calculated, forces are made. Since PbTe is cubic forces vanish by symmetry; this section refers to the output generated in the Bi_{2}Te_{3} tutorial.
The force evaluation has two parts. First a correction term to the forces is made. In contrast to the total energy, which is variational deviations of the input density relative to the selfconsistent density, n^{in}−n, thus δE ~ (n^{in}−n)^{2}, the forces are not. If it were known how the density shifts with the nucleus (which requires knowledge of the dielectric function), the forces could be made variational as well.
We can do almost as well by making an ansatz for how the density shifts. The simplest ansatz is to assume that there is a cloud of charge that shifts rigidly with nucleus. Assuming that this charge is the freeatomic density, a correction term can be devised which dramatically improves on the convergence of the forces with respect to deviations n^{in}−n. In many cases the error becomes almost variational, i.e. the error in the force to linear order in n^{in}−n becomes much smaller. At selfconsistency the correction term vanishes.
In the Bi_{2}Te_{3} tutorial the correction calculated from the starting density (Mattheis construction) reads:
Harris correction to forces: shift in freeatom density
ib deltan dVes deltan dVxc total
1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
2 0.00 0.00 605.58 0.00 0.00 23.90 0.00 0.00 629.48
3 0.00 0.00 605.58 0.00 0.00 23.90 0.00 0.00 629.48
4 0.00 0.00 292.27 0.00 0.00 9.54 0.00 0.00 301.81
5 0.00 0.00 292.27 0.00 0.00 9.54 0.00 0.00 301.81
shift forces to make zero average correction: 0.00 0.00 0.00
The forces are printed, including the correction mentioned above:
Forces, with eigenvalue correction
ib estatic eigval total
1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
2 0.00 0.00 79.59 0.00 0.00 66.61 0.00 0.00 12.98
3 0.00 0.00 79.59 0.00 0.00 66.61 0.00 0.00 12.98
4 0.00 0.00 6.78 0.00 0.00 19.25 0.00 0.00 12.46
5 0.00 0.00 6.78 0.00 0.00 19.25 0.00 0.00 12.46
Maximum Harris force = 13 mRy/au (site 3) Max eval correction = 629.5
In the last iteration when selfconsistency has been approximately achieved these tables become
Harris correction to forces: shift in freeatom density
ib deltan dVes deltan dVxc total
1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
2 0.00 0.00 0.39 0.00 0.00 0.00 0.00 0.00 0.39
3 0.00 0.00 0.39 0.00 0.00 0.00 0.00 0.00 0.39
4 0.00 0.00 0.77 0.00 0.00 0.19 0.00 0.00 0.58
5 0.00 0.00 0.77 0.00 0.00 0.19 0.00 0.00 0.58
shift forces to make zero average correction: 0.00 0.00 0.00
...
Forces, with eigenvalue correction
ib estatic eigval total
1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
2 0.00 0.00 638.27 0.00 0.00 627.55 0.00 0.00 10.71
3 0.00 0.00 638.27 0.00 0.00 627.55 0.00 0.00 10.71
4 0.00 0.00 278.13 0.00 0.00 262.20 0.00 0.00 15.93
5 0.00 0.00 278.13 0.00 0.00 262.20 0.00 0.00 15.93
Maximum Harris force = 15.9 mRy/au (site 5) Max eval correction = 0.6
At selfconsistency, the force on Te is 10.71 mRy/a.u.. With the correction term, the force from the initial guessed density is 12.98 mRy/a.u.. The error is thus of order 5 mRy/a.u., vastly smaller than it would be without the correction term (629.48 mRy/a.u). This term vanishes at selfconsistency, as it should.
The final forces give a measure of the error in the LDA predicting structure since site positions were taken from experiment. 10 mRy/a.u. is a fairly small force, implying that LDA lattice positions will be close to experimental ones.
To what extent the basis set affects the forces is taken up in Additional Exercises in the Bi_{2}Te_{3} tutorial.
Forces are used in molecular statics, and also to compute phonons.
New density
Questaal’s general strategy for guessing the selfconsistent density based on (n^{in},n^{out}) pairs, is described here.
It takes some linear combination of n^{in} and n^{out} from the current, and possibly prior iterations, based on settings in ITER_MIX.
lmf’s printout of this mixing looks like the following
mixing: mode=B nmix=2 beta=.3 elind=.531 kill=7
mixrho: sought 2 iter from file mixm; read 0. RMS DQ=2.17e2
charges: old new screened rms diff lin mix
smooth 15.480342 15.424647 15.424647 0.011907 15.463634
site 1 6.340155 5.708219 5.708219 0.022361 6.150574
site 2 1.820498 1.132866 1.132866 0.013213 1.614208
AMIX: nmix=0 mmix=8 nelts=6302 beta=0.3 tm=5 rmsdel=2.17e2
unscreened rms difference: smooth 0.011613 local 0.022361
screened rms difference: smooth 0.011907 local 0.022361 tot 0.021660
Notes:
 mode=B indicates that the Broyden method will be used when mixing prior iterations (in the first pass there are no prior iterations to mix)
 elind is the Fermi level entering into the model Lindhard function used to screen n^{out}−n^{in} (see Eq. (1)).
 beta=.3 is the mixing parameter, governing how much n^{out} is added to n^{in}, as explained in Eq. (4).
 RMS DQ=2.17e2 is a measure of the collective RMS change in n^{out}−n^{in}. This number is used to check convergence to selfconsistency. It is broken down into individual contributions in the table.
End of selfconsistency loop
At the end of the selfconsistency cycle the density is written to rst.pbte and a check is made for convergence. Checks in the first and second iterations read
iors : write restart file (binary, mesh density)
it 1 of 10 ehf= 55318.170567 ehk= 55318.148758
h nkabc=6 gmax=7.8 ehf=55318.1705672 ehk=55318.1487582
...
 BNDFP: begin iteration 2 of 10 
...
it 2 of 10 ehf= 55318.165774 ehk= 55318.156862
From last iter ehf= 55318.170567 ehk= 55318.148758
diffe(q)= 0.004793 (0.018916) tol= 0.000010 (0.000030) more=T
i nkabc=6 gmax=7.8 ehf=55318.1657745 ehk=55318.1568616
No check is made in the first iteration because there is no prior iteration to compare the change in total energy.
In subsequent iterations two checks are made: the change (0.004793) in total energy from the prior iteration and the RMS DQ (0.018916). Selfconsistency is reached when both checks fall below tolerances.. You can set tolerances with EXPRESS_conv and EXPRESS_convc, aliases for ITER_CONV and ITER_CONVC.
In the PbTe case convergence is reached in iteration 10, where the convergence checks read:
diffe(q)= 0.000000 (0.000005) tol= 0.000010 (0.000030) more=F
c nkabc=6 gmax=7.8 ehf=55318.1620974 ehk=55318.1620958
The last line prints out variables assigned on the command line (and variables in the ctrl file kept by the % save directive), total magnetic moment in magnetic calculations, and total energies from the HarrisFoulkes and KohnSham functionals. These functionals are different but they should approach the same value at selfconsistency.
This line is also written to file save.pbte; see here for further documentation.
Other Resources
This output was mostly taken from the basic lmf tutorial on a selfconsistent calculation for PbTe