# 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 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.

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 plane-wave 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 1-kappa parameters with Z-dependent LMX
2 or 4 2-kappa parameters with Z-dependent 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 **lmfa**-specific 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.

#### 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

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 5*s*, 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 5*d* to be a local orbital. When the 5*d* local orbital specified (by **PZ= 0 0 15.934**) in *basp.pbte*, **lmfa** will include the 5*d* 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 5*d* 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.

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 5
*d*,**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 5*d* state is too shallow to be treated as a core.

#### 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**.

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 single-particle sum.

The following table displays this information for each *l* that carries electrons. **rmt** is the augmentation radius.

```
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
```

**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 (self-consistent) 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 QS*GW*calculations the envelope functions need to be reasonably short range so that the self-energy 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 plane-wave 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 already-existing *basp* file, depending on how **HAM_AUTOBAS_LOC** is set. **lmfa** identifies local orbitals by from the highest-lying core level. If a core level of a particular *l* satisfies one of two criteria it will convert this core state into an extended local orbital:

- the charge spilling out of the MT radius exceeds
**QLOC** - the eigenvalue exceeds
**ELOC**

Near the beginning of execution **lmfa** prints a header with **QLOC** and **ELOC**:

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

You can modify default values with **HAM_AUTOBAS_QLOC** or **HAM_AUTOBAS_ELOC**.

For each species and each *l* within the species **lmfa** will inform you about local orbitals it finds. In this case it prints out

```
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
```

*Note:*: If **HAM_AUTOBAS_LOC=0**, **lmfa** will not look for LO; if it is set to unity, it searches for LO in channels where you have already set; if it is set to 2, your specifications from the *basp* or *ctrl*. file are ignored and all LO are determined by the internal algorithm.

For any extended LO that are given (either explicitly by you or by its internal algorithm), **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 8.2
```

*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**.

#### 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
```

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 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*.

### 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 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*^{in} 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.