Annotated standard output, program lmf

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
- Controlling what is read from the restart file
Potential and Matrix Elements
- LDA functionals
Reading QSGW self-energies+
Parameters for local orbitals
Brillouin zone integration
- Integration weights and the METAL switch
Output density and update of augmentation sphere boundary conditions
Total energy
Forces
New density
End of self-consistency 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)

# Autogenerated from init.pbte using:				     |	 ------------------ LMF - ref: 9c0462a3 ------------------
# blm --brief --ctrl=ctrl init.pbte --nk~ins			     <

# Variables entering into expressions parsed by input		     <
% const nit=10							     <
% const met=5							     <
%cconst so nsp=2						     <
% const so=0 nsp=so?2:1						     <
#%const lxcf=0 lxcf1=101 lxcf2=130    # PBE functional from libxc    <
% const lxcf=2 lxcf1=0 lxcf2=0        # built-in B-H functional.  Us <
% const pwmode=0 pwemax=3             # Use pwmode=1 or 11 to add AP <
% const nk1=4 nk2=nk1 nk3=nk2 gmax=8.2 beta=0.3			     <
% const conv=1e-5 convc=3e-5          # Convergence parameters in ch <

VERS  LM:7 FP:7 # ASA:7						     |	VERS  LM:7 FP:7
IO    SHOW=f HELP=f IACTIV=f VERBOS=35,35  OUTPUT=*			IO    SHOW=f HELP=f IACTIV=f VERBOS=35,35  OUTPUT=*

#SYMGRP i*r3(1,1,-1) r4x					     <

HAM                                   # Parameters controlling hamil |	HAM
      GMAX=   {gmax}                  # PW cutoff for charge density |	      GMAX=   7.8
      AUTOBAS[PNU=1 LOC=1 MTO=4 LMTO=4 GW=0]  # Parameters defining  |	      AUTOBAS[PNU=1 LOC=1 MTO=4 LMTO=4 GW=0]
      FORCES= {so==0} NFORCE=0 ELIND=-0.7			     |	      FORCES= 1 NFORCE=0 ELIND=-0.7
      NSPIN=  {nsp}                   # 2 for spin polarized calcula |	      NSPIN=  1
      SO=     {so}                    # 1 turns on spin-orbit coupli |	      SO=     0
      XCFUN=  {lxcf},{lxcf1},{lxcf2}  # set lxcf=0 for libxc functio |	      XCFUN=  2,0,0
								     |
ITER                                  # Management of iterations to  |	ITER
      NIT=    {nit}                   # Maximum number of iterations |	      NIT=    10
      MIX=    B2,b={beta},k=7         # Charge density mixing parame |	      MIX=    B2,b=.3,k=7
      CONV=   {conv}                  # Convergence tolerance (energ |	      CONV=   1e-5
      CONVC=  {convc}                 # tolerance in RMS (output-inp |	      CONVC=  3e-5
								     |
BZ                                    # Brillouin zone integration   |	BZ
      METAL=  {met}                   # Management of k-point integr |	      METAL=  5
      NKABC=  {nk1},{nk2},{nk3}       # 1 to 3 values		     |	      NKABC=  4,4,4

STRUC                                 # Crystal lattice information  |	STRUC
      FILE=site                       # Extract some or all of NBAS, |	      FILE=site

SPEC                                  # Chemical and basis informati |	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

SITE                                  # Read site data from auxiliar |	SITE
      FILE=site								      FILE=site

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, floating-point 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.

See Table of Contents

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  verb 35
 special:  forces
 pot:      XC:BH
 float:    float P LDA-style
 autoread: mto basis(4), pz(1), pnu(1)
 bz:       metal(5), tetra, invit

It tells you, for example that you are using the Barth-Hedin 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.00e-8   alat=12.14701   awald= 0.261   lmax=6
         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=1e-5
 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](/docs/numerics/rotations/#(/docs/numerics/rotations/#specifying-space-groups-in-questaal) 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 k-mesh 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 Methfessel-Paxton sampling integration scheme or a sampling with a Fermi function, using BZ_N=-1.

See Table of Contents

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 l-cutoff 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 l-cutoff 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 20 x 20 x 20 divisions; length 0.429, 0.429, 0.429
         generated from gmax = 8.2 a.u. : 4279 vectors of 8000 (53%)

 GVLIST: gmax = 8.2 a.u. created 4279 vectors of 8000 (53%)
         mesh has 20 x 20 x 20 divisions; length 0.429, 0.429, 0.429
 SGVSYM: 144 symmetry stars found for 4279 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 :  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 l-cutoff 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 plane-wave repesentation
cutoff : the number of plane waves with G<gmax.

 sugcut:  make orbital-dependent reciprocal vector cutoffs for tol=1.0e-6
 spec      l    rsm    eh     gmax    last term   cutoff
  Pb       0    1.80  -0.10   4.123    2.68E-06     531
  Pb       1    2.02  -0.10   3.848    2.42E-06     411
  Pb       2*   2.03  -0.10   4.013    1.37E-06     531
  Pb       3    2.03  -0.10   4.193    1.54E-06     537
  Pb       0    1.80  -0.90   4.123    2.68E-06     531
  Pb       1    2.02  -0.90   3.848    2.42E-06     411
  Pb       2    2.03  -0.90   4.013    1.37E-06     531
  Te       0    1.63  -0.10   4.572    1.45E-06     725
  Te       1    1.71  -0.10   4.575    1.52E-06     725
  Te       2*   2.02  -0.10   4.037    1.63E-06     531
  Te       3    2.02  -0.10   4.218    1.87E-06     537
  Te       0    1.63  -0.90   4.572    1.45E-06     725
  Te       1    1.71  -0.90   4.575    1.52E-06     725
  Te       2    2.02  -0.90   4.037    1.63E-06     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 orbital-dependent 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⁻⁶ 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.74E-06     701*

If ‘last term’ is too high you can expect errors in the hamiltonian and you should increase gmax.

The tolerance defaults to 10⁻⁶, but you can control it with tag HAM_TOL. 10⁻⁵ 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⁻¹⁶.

At this stage the potential independent setup is complete.

See Table of Contents

Obtain an input density

In Kohn-Sham 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 self-consistency cycle.

lmf reads atomic densities from atm.pbte and overlaps them (Mattheis construction). (Note that lmfa generated atm.pbte in a prior step.) A Mattheis construction 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 free-atom 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.659844
 site   2  spec  2  pos  0.0000  0.5000  0.0000  Qsmooth 7.820498
 total smooth Q = 15.480342

 Free atom and overlapped crystal site charges  spec 1  site 1:
   ib    true(FA)    smooth(FA)  true(OV)    smooth(OV)    local
    1   12.535727    6.195571   12.881552    6.541397    6.340155

 Free atom and overlapped crystal site charges  spec 2  site 2:
   ib    true(FA)    smooth(FA)  true(OV)    smooth(OV)    local
    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 3-5 will be read from rst.pbte and supersede pre-existing 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 command-line 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.

Switches for --rs=#1,#2,#3,#4,#5 mean the following:

#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 self-consistent density will no longer be fully self-consistent 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.

See Table of Contents

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 Barth-Hedin 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 Barth-Hedin 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.473394  avg sphere pot = 0.516004   vconst = -0.473394
 smooth rhoves     41.429795   charge    15.480342
 smooth rhoeps =  -11.301006   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 $n_1 - n_2$ , 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:

            1     1    1.954652    6.929061    68.00    82.00
            1    21    0.018533
            1    25    0.015664
            2     1   -0.307248   -1.089166    46.00    52.00

lmf has found the multipole moments from the local densities ( $\mathbf{Qval}$ is the charge, the l=0 multipole moment is $\mathbf{qmom} = Y_{00}\times \mathbf{Qval}$ ). 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.623903    0.458095
              21   -0.059350
              25   -0.050160
   2     2     1    1.767092    0.498487

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.478184  avg sphere pot = 0.518433   vconst = -0.478184
 Average es pot at MT boundaries after vconst shift
 Site    ves        Site    ves
   1  -0.020089  |    2   0.020304  |

It is done order to put the Fermi level close to zero. This shift can be chosen differently; if you use --oldvc the cell-averaged 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

 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

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 $\phi$ and $\dot{\phi}$ 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²), and linearization energy enu:

 l        enu         v           c          srdel        qpar        ppar
   -0.685734   -1.220935   -0.974210    0.311409      2.5442    2.037535
   -0.345260   -0.961024    0.034413    0.322662      9.5614    3.768910
    0.215078   -0.351058    1.877978    0.404577     13.6181    3.646178
 2(sc) -1.355836   -1.700377   -1.373128    0.018376     60.8791    0.594481
   -0.367191   -0.903273    2.382599    0.347764     27.1687    7.102701
   -0.261797   -0.813620    4.861502    0.414942     32.9588   13.704207

The following table presents (hamiltonian, potential, overlap) matrix elements of special linear combinations val and slo of $\phi$ and $\dot{\phi}$ that have the property that value or slope vanishes at the MT boundary:

  l         val-val     val-slo     slo-val     slo-slo
  0   s   10.370596   -7.396703   -7.396703    7.271920
      h  -12.303585    0.102866    9.373758   -0.869622
      v  -67.729851   86.547317   86.547317 -115.586449
  1   s    8.280043   -4.050478   -4.050478    2.712195
 ...

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.292-0.936 0.198      0.00      0.00      0.000      0.00000
        0.771 0.109-0.627      0.00      0.00
        0.565 0.336 0.754      0.00      0.00

   2    0.292-0.936 0.198      0.00      0.00      0.000      0.00000
        0.771 0.109-0.627      0.00      0.00
        0.565 0.336 0.754      0.00      0.00

These are directly measurable quantities. At least in $CuAlO_2$ , the electric field gradient on the $Cu$ site (measured to be $\pm 10.6 \cdot 10^{21} V/cm^2$ ), is not well predicted by LDA, or LDA+U. QSGW, however, does quite well. See Fig.13 of this paper.

If spin-orbit coupling is set (HAM_SO=t), matrix elements of $\lambda L \cdot S$ are generated. Working from the Fe tutorial, invoke

$ lmf ctrl.fe -vso=1 --iactiv  --quit=pot --pr55

These matrix elements of $\lambda L \cdot S$ are printed out:

 soprm:  matrix elements for perturbation from spin-orbit 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.

See Table of Contents

Reading QSGW self-energies+

If you have performed a QSGW calculation and read a self-energy $\Sigma^0$ from file sigm.ext (or sigmrs.ext), $\Sigma^0$ will be added to the LDA potential.
Note: sigm.ext contains the difference, $\Sigma^0{-}V_{xc}^\mathrm{LDA}$ , so it is added directly to the Kohn-Sham potential $V_{xc}^\mathrm{LDA}$ .

You will see an indication $\Sigma^0$ 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 Σ⁰ 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 Σ⁰ 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 self-energies (0.122871 spin-up and 0.138308 spin-down) above 2 Ry. When SIGP_MODE=4 (recommended), these numbers are used for Σ_ii⁰ 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 Bloch-summed Σ⁰, 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 Σ⁰.

This table:

 hft2rs created hrs:  ndhrs=30  max Re(hrs) = 0.556  max Im(hrs) = 3.93e-10
 ...
 check FT s(R) against s(q) ... maximum error = 2.3e-15 < tol (5e-6)
 ...
 check FT s(R) against s(q) ... maximum error = 7.8e-15
 ...

indicates the precision in the construction of the inverse Bloch sum. It:

checks the forward Bloch sum of the inverse Bloch summed Σ⁰ (2.3e-15) to ensure that the error is less than HAM_RSSTOL. The combined inverse/forward sum should recover the original Σ⁰ to machine precision unless the inverse Bloch transform is inexact, as described above.
Shows deviation (7.8e-15) from original Σ⁰ 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 k-point 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 (self-consistent 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 density-of-states there (14.389/Ry), and the single-particle sum -1.3036417 Ry entering into the LDA total energy. The LDA magnetic moment/cell (2.200354 $\mu_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 $E_F$ , 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 chicken-and-egg 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 self-consistent 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 self-consistency cycle. This scheme works for sampling integration where the k-point 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.

See Table of Contents

Output density and update of augmentation sphere boundary conditions

A second pass over the irreducible k-points 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 Harris-Foulkes energy is a functional of the input density and single-particle sum. Information is printed out in the table below.

 Harris energy:
 ...
 rhoeps=   -1094.806758     utot= -116742.136483    ehar=  -55318.170567

The Kohn-Sham total energy is also calculated:

 Kohn-Sham 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ⁱⁿ. It is a different functional than the HF functional, but the two should come together at the self-consistent 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₂Te₃ 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 self-consistent density, nⁱⁿ−n, thus δE ~ (nⁱⁿ−n)², 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 free-atomic density, a correction term can be devised which dramatically improves on the convergence of the forces with respect to deviations nⁱⁿ−n. In many cases the error becomes almost variational, i.e. the error in the force to linear order in nⁱⁿ−n becomes much smaller. At self-consistency the correction term vanishes.

In the Bi₂Te₃ tutorial the correction calculated from the starting density (Mattheis construction) reads:

 Harris correction to forces: shift in free-atom density
  ib         delta-n dVes             delta-n 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 self-consistency has been approximately achieved these tables become

 Harris correction to forces: shift in free-atom density
  ib         delta-n dVes             delta-n 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 self-consistency, 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 5mRy/a.u., vastly smaller than it would be without the correction term (629.48 mRy/a.u). This term vanishes at self-consistency, 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₂Te₃ tutorial.

Forces are used in molecular statics, and also to compute phonons.

See Table of Contents

New density

Questaal’s general strategy for guessing the self-consistent density based on (nⁱⁿ,n^out) pairs, is described here.

It takes some linear combination of nⁱⁿ 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.17e-2
 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.17e-2
 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ⁱⁿ (see Eq. (1)).
beta=.3 is the mixing parameter, governing how much n^out is added to nⁱⁿ, as explained in Eq. (4).
RMS DQ=2.17e-2 is a measure of the collective RMS change in n^out−nⁱⁿ. This number is used to check convergence to self-consistency. It is broken down into individual contributions in the table.

End of self-consistency loop

At the end of the self-consistency 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). Self-consistency 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 Harris-Foulkes and Kohn-Sham functionals. These functionals are different but they should approach the same value at self-consistency.

This line is also written to file save.pbte; see here for further documentation.

See Table of Contents

Other Resources

This output was mostly taken from the lmf tutorial on a self-consistent calculation for PbTe