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
 5. Estimating the planewave 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.
Add input
to the lmfa and lmf commands
$ lmfa ctrl.pbte input > out.lmfa
$ lmf ctrl.pbte input > out.lmf
The output is quite verbose so only some differences are highlighted.
This appears in the lmf output
 Parameters for hamiltonian 
gmax reqd* r8 1, 1
(alias for HAM_GMAX)
Energy cutoff for planewave mesh (Ry)
 If preceding token is not parsed, attempt to read the following:
lmfa does not read EXPRESS_gmax (alias for HAM_GMAX). For this reason you don’t need to specify it when running lmfa, but do need to when running lmf. Similarly, lmfa does not read BZ_NKABC (another tag lmf requires), and many other tags that are optional in lmf, e.g. HAM_SO.
lmfa looks for this tag, but lmf does not:
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)
Two others lmfa looks for but lmf does not are SPEC_ATOM_Q and SPEC_ATOM_MMOM. These affect the construction of the atomic charge density.
The following tag is read by both codes, but the meaning differs. From the output of lmfa:
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 1kappa parameters with Zdependent LMX
2 or 4 2kappa parameters with Zdependent LMX
...
From the output of lmf:
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
...
In the first case the token affects how the basis parameters are generated; in the second how the basis parameters are read. The same applies to HAM_AUTOBAS_PNU and HAM_AUTOBAS_LOC; also HAM_AUTOBAS_PFLOAT have different meanings.
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.
Consider the output of lmfa help
.
First appears a list of command line options available in most Questaal codes, a portion of which is shown below. They are described in more detail here.
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
pr#1[,#2...] Set the verbosity (stack) to values #1,#2, ...
...
vnam=expr Define numerical variable "nam"; set to result of 'expr'
cnam=strn Define character variable "nam"; set to 'strn'
...
Then follow a synopsis of lmfaspecific options:
noopt Suppress optimization of s.m. Hankel basis
norscnst In optimization of s.m. Hankel basis, do not constrain rsm < rmt
plotwf Writes atomic radial wave functions to disk files
dumprho Writes out the density for each atom to out.ext
basp Turns on autofind EH,RSMH (better to use HAM_AUTOBAS)
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.
Selfconsistent 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(i1)}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 selfconsistent charge density should be the same in either case, but the valencecore partitioning is different.
lmfa starts with a crude guessed density and after 55 iterations converges to the selfconsistent 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.614E05 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.
The following is taken from standard output:
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
 eval is the eigenvalue
 node at is the largest radius for which the wave function has a node
 max at is the radius where the wave function has a maximum value
 c.t.p is the classical turning point
 rho(r>rmt) is the charge spilling outside the augmentation radius. For Pb 5d, 0.007786 electrons spill out: this is on the ragged edge of whether it needs to be included as a local orbital (see Additional Exercises in the tutorial).
Note: for GW calculations the Pb 5d state is too shallow to be treated as a core.
Generating basis information
From the selfconsistent 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.
Fitting RSMH and EH to the numerically derived wave functions can be readily accomplished. (It is a nonlinear procedure.) lmfa actually does it twice. First it fits the occupied levels only, for which boundary conditions are known, and at the same time it computes a variational estimate for the energy. In perturbation theory this differs from the exact value computed from numerical wave functions as the difference in the singleparticle sum.
The following table displays this information for each l that carries electrons. rmt is the augmentation radius.
Optimise freeatom 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
 Rsm and Eh are the optimum values for RSMH and EH in the free atom
 StiffR and stiffE are the sensitivity of the total energy to changes in Rsm and Eh
 Eval is the expectation value of the approximate wave function (numerical for r<rmt envelope for r>rmt)
 Exact is the exact eigenvalue for the given (selfconsistent) potential.
 Pnu is the logarithmic derivative of the exact atomic wave function at rmt, converted into the continuous principal quantum number
Notes:
 A basis of augmented partial waves for r<rmt and the envelopes for r>rmt is the atomic analog of the crystal code.
 Eval, determined from the variational principle, should be less binding than the Exact value
 The error in the total energy from augmented partial waves, is given in the last line.
 The fit is evidently very good. Unfortunately, it is less good in the crystal because crystal eigenfunctions consist of tails from envelopes centered on multiple atoms. So, even if the potential were the same as atomic potential there (in reality it changes a little) the crystal eigenfunction will not be perfect at rmt. The new Jigsaw Puzzle Orbital basis is designed to deal exactly with these issues.
After the initial fit, a second fit is carried out, this time with all the partial waves. Unoccupied states may not be bound, so Eh is set to a default value. This information is displayed in the following table.
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

Empty states are not necessarily bound, in which case they need not decay exponentially as r→∞. Nevertheless the smooth Hankel basis set requires decaying functions. For these states EH is set to a default value.
Note: For QSGW calculations the envelope functions need to be reasonably short range so that the selfenergy can be smoothly interpolated between k points. Usually EH should be −0.3, or deeper for atoms where the potential is deep (second row or column VI, VII, and VII). 
The key new piece of information is Gmax. The maximum value of this number is used in estimating the planewave cutoff for the density mesh.
Local orbitals may already be specified when lmfa begins execution. Information about them is supplied by PZ in the input file or an alreadyexisting basp file (depending on how HAM_AUTOBAS_LOC is set)
In this event lmfa will try to fit the tail of this orbital to a smooth Hankel for r>rmt, in a manner similar to its fit of the envelope functions. For the Pb atom, the fit reads
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 7.8
Notes:
 The eigenvalues can be quite deep with the wave function steeply decaying. Left unconstrained Rsm can become quite small. This makes for a very sharply peaked function. In such a case the weight of the wave function (measured by Q(r>rmt)) is small and the fit does not need to be especially accurate. To protect against functions too sharply peaked, a constraint is placed on lower bound Rsm. Most of the time the default is ok, but if not the eigenvalue will deviate significantly from the exact one. You can set the lower bound with tag SPEC_ATOM_RS3. Be careful not to allow it to become too small; it makes representation of the interstitial density difficult and numerically unstable.
 Using the fit Rsm and Eh lmfa can estimate how many plane waves are needed to accurately represent the intersttial part of this function. That number is Gmax. See the section on GMAX.
Find local orbitals which satisfy E > 2 Ry or q(r>rmt) > 5e3
l=2 eval=1.569 Q(r>rmt)=0.0078 PZ=5.934 Use: PZ=15.934
l=3 eval=9.796 Q(r>rmt)=3e8 PZ=4.971 Use: PZ=0.000
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 onecenter 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
The valence density was fit to a linear combination of 6 smooth Hankel functions, with varying energies but a fixed smoothing radius rsm=1.51434. The fit is carried out subject to the constraint that the integrated charge is exact (and the integrated magnetic moment, in the spin polarized case).
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
The core density was fit to a single smooth Hankel function. The energy (14.401) was determined by the fitting procedure.
5. Estimating the planewave cutoff GMAX
The charge density is represented on a uniform mesh of points. The spacing between points or equivalently the planewave 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.3 (valence) 7.8 (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.
FAQ
 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 HelmanFeynman theorem.
Around a selfconsistent point (the independent variable is the density, in KohnSham theory) the change in total energy from a perturbation is the change in the singleparticle 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 singleparticle sum is a pretty good estimate of the change in total energy.
There are reasons why this is approximately true even when not selfconsistent. It is exact for onebody model Hamiltonians (i.e. that that have no electronelectron interactions). Matthew Foulkes showed that densityfunctional theory can be cast in the form similar to model hamiltonians, with the total energy a sum of the singleparticle sum and an extra term that is functional of the input density n^{in} only. The extra term vanishes at selfconsistency in response to loworder perturbation, and often doesn’t change by a large amount even if the density isn’t so different from a selfconsistent one.
Edit This Page