Annotated standard output, program lmfa

Table of Contents

• 1. Preprocessor’s transformation of the input file
• 2. How lmf and lmfa read input
• 3. Header, lattice, and symmetry blocks
• 4. Loop over species
  • Self-consistent density
    • Radial mesh
    • Parameters that specify the charge density
  • Valence and core wave functions
  • Generating basis information
  • Fitting the charge density outside the augmentation radius
• 5. Estimating the plane-wave cutoff GMAX
• FAQ

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 organised by blocks. The sections below annotate each block approximately in the order they are made. Some parts are similar to the lmf output; in places where they are similar a cursory treatment is given here; the reader is referred to that page.

Note: This document is annotated for an output vebosity of 35 (medium verbosity) Raising the verbosity causes lmfa to print additional information. Verbosity 31 is a little terse; 41 is a little verbose while 55 is quite verbose. 80 can give you a headache.

1. 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 lmfa --showp ... The result is very similar to lmf --showp ..., which is documented here.

2. How lmf and lmfa read input

To see what tags lmfa will look for, use lmfa --input. This web page explains the function of --input.

The result is similar to lmf produces. However lmfa parses fewer tags; moreoever, some tags they both parse have different meanings.

$ lmfa ctrl.pbte --input > out.lmfa
$ lmf  ctrl.pbte --input > out.lmf

 --- Parameters for hamiltonian ---
 gmax                   reqd*  r8       1,  1
   (alias for HAM_GMAX)
   Energy cutoff for plane-wave mesh (Ry)
 - If preceding token is not parsed, attempt to read the following:

 autobas_lmto           opt    i4       1,  1     default = 0
   (alias for HAM_AUTOBAS_LMTO)
   Controls lmfa's autogeneration of LMTO basis parameters (RSMH,EH,RSMH2,EH2)

 autobas_mto            opt    i4       1,  1     default = 0
   (alias for HAM_AUTOBAS_MTO)
   Controls lmfa's autogeneration of LMTO basis parameters (RSMH,EH,RSMH2,EH2)
   0 do not autogenerate basis parameters
   1 or 3 1-kappa parameters with Z-dependent LMX
   2 or 4 2-kappa parameters with Z-dependent LMX
   ...

 autobas_mto            opt    i4       1,  1     default = 0
   (alias for HAM_AUTOBAS_MTO)
   Autoset basis:  controls what part of MTO basis is read from basis file basp
   0       No parameters are read from basp
   ...

--help performs a function similar to --input but for the command line arguments: it prints out a brief summary of arguments effective in the executable you are using.

 --help         Print this message, and quit
 --input        List categories, tokens, and data program expects, and quit
 --show         Print control file after parsing by preprocessor,
                and echo input data as read from the control file
 --showp        Same as --show, but quit after input parsed
 --iactiv       (--no-iactiv) Turn on (off) interactive mode
                This switch overrides input file setting
 --pr#1[,#2...] Set the verbosity (stack) to values #1,#2, ...
 --time=#1[,#2] Print timing info to # levels (#1=summary; #2=on-the-fly)
 --shomem       Print memory allocation of some arrays (equivalent to IO_SHOMEM=1)

 --rpos=filnam  After reading input file, read site positions from "filnam"
 --rdpos=filnam Add contents of "filnam" to site positions.  For options, see web page
 --shorten[opts]Shorten site positions.  For options, see web page
 --fixlat       Adjust lattice vectors and point group ops, attempt to
                render them internally consistent
 --fixpos[:tol=#] Adjust positions slightly, rendering them
                as consistent as possible with the symmetry group
 --nosym        Suppress symmetry operations

 -vnam=expr     Define numerical variable "nam"; set to result of 'expr'
 -cnam=strn     Define character variable "nam"; set to 'strn'

 --tbeh         Take sm-Hankel energies EH(1) and EH(2) from STR_ENV_EL
 --noopt        Suppress optimization of s.m. Hankel basis
 --nopz         Suppress search for local orbitals
 --norscnst     In optimization of s.m. Hankel basis, do not constrain rsm < rmt
 --basp[~opts]  Turns on autofind EH,RSMH (better to use HAM_AUTOBAS)
                Options are :  ctrl  eh=#  eh2=#  rsmmx  incrlmx
 --basfile=fn   write the autogenerated basis set parameters to file fn instead of default basp0
 --usebasp      Remake the core density after the search for local orbitals; write basis parameters to file basp
 --plotwf       Writes atomic radial wave functions to disk files
 --dumprho      Writes out the density for each atom to out.ext
 --getallloc    Look for local orbitals (better to use HAM_AUTOBAS)

After transformation by the preprocessor, lmfa parses for tags and substitutes default values for tags it does not find. To see the value of tags lmfa, whether parsed or defaults, use lmfa --show or lmfa --show=2. The latter causes lmfa to stop after displaying tags, and is useful if you want to see whether lmfa 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).

The action of --show is similar to what occurs with lmf; see see annotation of lmf output for further discussion.

3. Header, lattice, and symmetry blocks

The header information presents a condensed synopsis of some key settings that are used in the calculation. It is similar to the header table produced by lmf.

The next blocks print information about the lattice vectors and settings used in Ewald summations. This is not relevant for lmfa; but it is printed out anyway and is identical to lmf output; similarly for the following block on symmetry and k mesh.

4. Loop over species

lmfa begins a loop over each species, and performs the following steps.

Self-consistent density

If no local orbitals have been specified, lmfa’s printout for the first atom (Pb) begins with:

 Species Pb:  Z=82  Qc=78  R=3.044814  Q=0
 mesh:   rmt=3.044814  rmax=47.629088  a=0.025  nr=497  nr(rmax)=607
  Pl=  6.5     6.5     6.5     5.5     5.5
  Ql=  2.0     2.0     0.0     0.0     0.0

The first line shows the atomic number, number of core charges (core levels are assumed to be filled), augmentation radius and net sphere charge. The next line shows the augmentation radius and cutoff radius of the free atom and radial mesh parameters.

Radial mesh

Wave functions are integrated on a uniform radial mesh of points. All the Questaal codes use a shifted logarithmic radial mesh: point i has a radius $$r_i = b[e^{a(i-1)}-1]$$

Mesh points are linearly spaced by ab near the nucleus. For r_i large compared to a, the mesh points are spaced exponentially (equally spaced on a log scale, spacing a).

Three numbers are used to specify the mesh: augmentation radius, the number of points inside it, and the spacing parameter a (b can be determined from them). These can be specified in the input file as SPEC_ATOM_R, SPEC_ATOM_NR,  SPEC_ATOM_A though usually you can rely on default values.

The free atom calculation doesn’t need to know about the augmentation radius, but it is needed for atm.pbte where the augmentation and interstitial parts are kept separate. More details about the how the Questaal codes perform radial integrations can be found here.

Parameters that specify the charge density

The Pl in the table are the logarithmic derivative parameters, (also called continuous principal quantum numbers) for s, p, d, … valence orbitals. For free atoms the fractional part is not relevant since there is a boundary condition that the wave function decay exponentially as r→∞). The integer parts of Pl are important because they define what states are valence. All states with smaller principal quantum number, e.g. the 5s, are treated as core. The Ql are corresponding charges.

Pl and Ql are specified by  SPEC_ATOM_P  and  SPEC_ATOM_Q. Neither Pl nor Ql are required inputs: lmfa will use default values from a lookup table for whatever is missing.

As described in the tutorial, lmfa finds the Pb 5d to be a local orbital. When the 5d local orbital specified (by PZ= 0 0 15.934) in basp.pbte, lmfa will include the 5d in the valence. The printout then reads:

 Species Pb:  Z=82  Qc=68  R=3.044814  Q=0
 mesh:   rmt=3.044814  rmax=47.629088  a=0.025  nr=497  nr(rmax)=607
  Pl=  6.5     6.5     5.5     5.5     5.5
  Ql=  2.0     2.0     10.0    0.0     0.0

Qc is smaller by 10 and the 5d are included in the valence with 10 electrons (Pl is reduced to 5.5 and Ql becomes 10).

The Ql and the boundary condition at r→∞ are sufficient to completely determine the charge density. The total self-consistent charge density should be the same in either case, but the valence-core partitioning is different.

lmfa starts with a crude guessed density and after 55 iterations converges to the self-consistent one.

  iter     qint         drho          vh0          rho0          vsum     beta
    1   82.000000   2.667E+04      410.0000    0.4078E+03     -164.7879   0.30
   55   82.000000   4.614E-05     1283.9616    0.3612E+08     -309.4131   0.30

Valence and core wave functions

In this block information about the eigenvalues of the valence and core states it finds along with some additional information, such as what fraction of the state falls outside the augmentation radius.

 valence:      eval       node at      max at       c.t.p.   rho(r>rmt)
   6s      -0.91143         1.014       1.961       2.814     0.168779
   6p      -0.27876         1.185       2.643       4.790     0.524423
   5d      -1.56879         0.523       1.073       2.252     0.007786
...

 core:        ecore       node at      max at       c.t.p.   rho(r>rmt)
   1s   -6465.77343         0.000       0.010       0.022     0.000000
   2s   -1155.91509         0.020       0.057       0.090     0.000000
...
   5p      -6.31315         0.486       0.882       1.314     0.000052

Generating basis information

From the self-consistent density and potential, lmfa will build some basis set information which is written to template basp0.pbte. It supplies

• Envelope function (parameters EH and RSMH) defining the shape of the envelope functions.
• Boundary condition at the augmentation radius, encapsulated as logarithmic derivative parameter P.
  P is also called the “continuous principal quantum number”.
• Information about local orbitals, indicated as PZ in the input.

These quantities are supplied in the input file as SPEC_ATOM_EH  SPEC_ATOM_RSMH,  SPEC_ATOM_P,  and SPEC_ATOM_PZ.
RSMH, EH, P, and PZ are also saved in basp0.pbte. lmf may read these parameters from basp.pbte, depending on settings in HAM_AUTOBAS.

 Optimise free-atom basis for species Pb, rmt=3.044814
 l  it    Rsm      Eh     stiffR   stiffE      Eval      Exact     Pnu    Ql
 0  50   1.794  -0.698      88.6    108.3   -0.91141  -0.91143    6.89   2.00
 1  26   2.026  -0.161     182.6    936.8   -0.27824  -0.27876    6.81   2.00
 eigenvalue sum:  exact  -2.38037    opt basis  -2.37931    error 0.00106

 Make LMTO basis parms for species Pb to lmxb=3, rmt=3.0448  vbar=0
 l  it    Rsm       Eh        Eval      Exact     Pnu    Ql   Gmax
 0  30   1.803   -0.706    -0.91141  -0.91143    6.89   2.00   3.9
 1  18   2.024   -0.160    -0.27825  -0.27876    6.81   2.00   3.6
 2   1   2.030+  -0.100+   -0.11512   0.01352    6.24   0.00
 3   1   2.030+  -0.100+    0.20219   0.02051    5.18   0.00

 Autogenerated Pnu:  6.895 6.809 6.250 5.183 5.089

 FREEAT: G tolerance for envelope functions : 1e-6
 Search criteria for local orbitals:  QLOC: qistl > 5e-3  ELOC: eig < -2 Ry

 Search for possible local orbitals
 LO l=2: qistl=0.00779 exceeds QLOC, autogen PZ=15.93359 use PZ=15.93359
 LO l=2: e=-1.57 exceeds ELOC, autogen PZ=15.93359 use PZ=15.93359

 Fit local orbitals to sm hankels, species Pb, rmt=3.044814
 l   Rsm    Eh     Q(r>rmt)   Eval      Exact     Pnu     K.E.   fit K.E.  Gmax
 2  1.041 -1.083   0.00792  -1.56878  -1.56879   5.934  -0.8111  -0.8111    8.2

Fitting the charge density outside the augmentation radius

lmfa retains two densities on the numerical mesh for points inside the augmentation sphere: the core and valence densities.

In addition it must keep some information about the density outside. This density will become part of the interstitial density in the crystal. Thus it must be represented on the interstitial charge density mesh and one-center expansions of it made to include the contribution to neighboring sites from this atom’s density.

Both can be readily accomplished if the density is represented as a linear combination of smooth Hankel functions. lmfa fits the numerical density to a linear combination of such functions; the smoothing radii, energies, and fit coefficients are stored in atm.pbte.

 tailsm: fit tails to 6 smoothed hankels, rmt= 3.02869, rsm= 1.51434
 HNSMFT:  85 points in interval  3.02869  24.73301;  q=  1.433395
 E:    -1.00000    -2.00000    -4.00000    -6.00000    -9.00000    -15.0000
 C:    -0.10708    13.99477    -208.816    3199.570    -32652.0    736827.1
        r          rho         fit         diff
    3.028689    0.013779    0.013777    0.000002
    3.888922    0.003349    0.003348    0.000000
    4.993484    0.000541    0.000541    0.000000
    6.411769    0.000054    0.000054    0.000000
    8.232883    0.000003    0.000003    0.000000
    q(fit):     1.433395    rms diff:   0.000002
    fit: r>rmt  1.433395   r<rmt  6.387103   qtot  7.820498
    rho: r>rmt  1.433395   r<rmt  4.566605   qtot  6.000000

Also the core density spills out into the interstitial. Rather than renormalizing the core density inside the augmentation sphere, it is allowed spill out and included in the interstitial density. lmfa fits the core tail to a single smoothed Hankel function.

 coretail: q=0.00215, rho(rmt)=0.00566.  Fit with Hankel e=-14.401  coeff=639.|
      r            rhoc          fit
    3.028689    0.01975240    0.01975240
    3.347222    0.00650048    0.00651738
    3.792904    0.00136878    0.00136091
    4.297927    0.00023338    0.00022687
    4.870194    0.00003132    0.00002930
    5.518656    0.00000320    0.00000283
    6.253461    0.00000024    0.00000020
    7.086104    0.00000001    0.00000001

5. Estimating the plane-wave cutoff GMAX

The charge density is represented on a uniform mesh of points. The spacing between points or equivalently the plane-wave cutoff GMAX used to represent the density in reciprocal space is governed by how sharply peaked the envelope functions are.

The cutoff that will reasonably represent the expansion of envelope functions is determined separately for valence orbitals and local orbitals. Both maximum values for all species are printed in the table below.

 FREEAT:  estimate HAM_GMAX from RSMH:  GMAX=4.6 (valence)  8.2 (local orbitals)

lmfa is now complete. Parmeters needed for both the augmentation and interstitial parts of the atomic density are saved in atm.pbte for each species, and basis information parameters are saved in basp0.pbte.

See Table of Contents

FAQ

1. Why can the total energy change be calculated as a sum of eigenvalues, when fitting an epproximate wave function to the exact one?

This follows from the Helman-Feynman theorem.

Around a self-consistent point (the independent variable is the density, in Kohn-Sham theory) the change in total energy from a perturbation is the change in the single-particle sum to lowest order, because to first order the wave function is stationary. This is an exact statement.

For many properties, e.g. magnetocrystalline anisotropy, for the effect of spin orbit coupling on the total energy, the change in single-particle sum is a pretty good estimate of the change in total energy.

There are reasons why this is approximately true even when not self-consistent. It is exact for one-body model Hamiltonians (i.e. that that have no electron-electron interactions). Matthew Foulkes showed that density-functional theory can be cast in the form similar to model hamiltonians, with the total energy a sum of the single-particle sum and an extra term that is functional of the input density nⁱⁿ only. The extra term vanishes at self-consistency in response to low-order perturbation, and often doesn’t change by a large amount even if the density isn’t so different from a self-consistent one.