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

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.

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 numbersP 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 :

The group operations were determined automatically from the given lattice. PbTe is cubic, with 48 group operations. First the crystal system is determined; then the symmetry operations inconsistent with the basis are discarded. Finally the 48 operations are distilled into a minimum set of generators (i*r3(1,1,-1) r4x) that generate the entire space group.

Generators comprise a minimum set of group operators. All possible products of them form a closed set of operations, which forms a group.

In this case there are no translations; all the group operators are pure point group operations. This is not true in general. For example, hcp Co has 24 space group operations. The generators defining its space group can be written as:

i*r3z::(1/3,-1/3,-1/2) r2z::(1/3,-1/3,1/2) r2x

You can also specify symmetry operations manually via the SYMGRP category. This is particularly useful when magnetic symmetry must be considered.

The k mesh is specifed through the number of divisions along each of the three reciprocal lattice vectors, tag EXPRESS_nabc or BZ_NABC. In this case the 216 k points generated by a 6×6×6 mesh gets reduced by symmetry to 16 inequivalent points. You can also specify whether the k-mesh should pass through the origin or straddle it (see BZ_BZJOB).

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.

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

Counting the size of the basis

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.

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

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−6 in this case).

Otherwise it uses all the G vectors available to it, and appends a ‘*’ to the number indicating not enough orbitals are available to meet the specified tolerance. This is flagged by an asterisk appended to the cutoff, e.g.

Te 0 1.63 -0.10 4.572 1.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−6, but you can control it with tag HAM_TOL. 10−5 or smaller is usually safe.

Note: With APW’s also in the basis, it the overlap matrix becomes nearly singular. To stabilize the overlap matrix, you can set HAM_TOL to something close to machine precision, e.g. 10−16.

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). This makes a reasonable guess for the density: the atomic potential is very deep, and the crystal potential is a relatively modest perturbation to it. Thus the true density is typically not so far from a superposition of atomic densities.

This table indicates what lmf is doing, and checks that atm.pbte is consistent the input file.

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.

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

Potential and 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

These few lines conceal a somewhat involved process. For the mesh density to yield the correct potential in the interstial, lumps of charge must be added to the mesh density. Provided the additional charge has multipole moments matching those obtained by the difference in the two local densities , multipoles of the (modified) mesh density will match those of the true density, and the electrostatic potential generated by it will be correct at each augmentation radius.

Run lmf with a high verbosity (--pr50 --quit=pot) and the following appears:

lmf has found the multipole moments from the local densities ( is the charge, the l=0 multipole moment is ). It adds localized gaussian charges with smoothing radius rg. (Notice that rg is quite small, to ensure that the gaussian does not spill out of the augmentation sphere.)

Now the electrostatic potential at the augmentation radius can be constructed

Note the printout Vfloat=T. This indicates that the potential used integrate the partial waves is being updated to the current local potential. The potential used to make partial waves and is kept independently of the potential generated by the density. Initially this potential is taken from of the free atom; and normally it is overwritten (Vfloat=T). If you tell lmf to freeze it (HAM_FRZWF=t globally or SPEC_ATOM_FRZWF=f for one species) the basis set will not adapt to the potential.

Once again, many steps are concealed behind this limited output. You can get more information, e.g. matrix elements, using a higher print verbosity: the higher the verbosity, the more is output.

Run lmf with (--pr70 --quit=pot) and a great deal of information is printed out. A table of classical LMTO potential parameters is made. Particularly useful are the indicating the band center (c) and bandwith (srdel2), and linearization energy enu:

These are directly measurable quantities. At least in , the electric field gradient on the site (measured to be ), is not well predicted by LDA, or LDA+U. QSGW, however, does quite well. See Fig.13 of this paper.

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

Reading QSGW self-energies+

If you have performed a QSGW calculation and read a self-energy from file sigm.ext (or sigmrs.ext), will be added to the LDA potential.Note:sigm.ext contains the difference, , so it is added directly to the Kohn-Sham potential .

You will see an indication is read from the following standard output. It was taken from the Fe tutorial.

The cutoff E(lda)>2 is controlled by HAM_SIGP_EMAX. It affects the partitioning between calculated matrix elements of Σ0 and and the fixed ones above emax, as explained at the top of p165107 in the PRB paper (EMAX corresponds to Exccut2).

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 Σii0 for energies E > EF+EMAX. In this mode the code prints out the line containing read constant sigii from sigm file.

Range 5 for inverse Bloch-summed Σ0, in units of alat. If this number is not sufficiently large, the inverse Bloch transform is incomplete, and the program may stop with an error message.

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

checks the forward Bloch sum of the inverse Bloch summed Σ0 (2.3e-15) to ensure that the error is less than HAM_RSSTOL. The combined inverse/forward sum should recover the original Σ0 to machine precision unless the inverse Bloch transform is inexact, as described above.

Shows deviation (7.8e-15) from original Σ0 after symmetrization by the group operations.

Parameters for local orbitals

Extended local orbitals must attach a smooth Hankel tail. A single smooth hankel function that fits both value and slope is found. The result displayed in a table:

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. EF and k-point weights are saved into a binary file wkp.pbte.

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)

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μ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:

The weight of the first point (0.009259) is 2/216. The factor of 2 is for spin. In spin polarized calculations this factor gets reduced to 1.

If you know that the system is an insulator in advance, you can tell lmf to assume an insulator: in the input file replace % const met=5 with % const met=0. This will set tag EXPRESS_metal (an alias for BZ_METAL) to zero. In this mode it can accumulate the density on the first pass as well.

In metals, the weights depend on the Fermi level , which must be determined from the eigenvalues. Quadratures are of two types:

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 EF. 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 EF 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 EF in the course of the self-consistency cycle. This scheme works for sampling integration where the k-point weights depend on EF only and not eigenvalues at neighboring k. You can also use this scheme in a mixed tetrahedron/sampling method: Weights for the band sum are determined by tetrahedron, but the charge density is integrated by sampling. It works better than straight sampling because the total energy is variational in the density.

Like METAL=3 in that two passes are made but eigenvectors are kept in memory so there is no additional overhead in making the first pass. This is the recommended mode for lmf unless you are working with a large system and need to conserve memory.

The ASA implements METAL=0,1,2; the FP codes METAL=0,2,3,4,5. See Additional Exercises in the Fe tutorial.

Tetrahedron and sampling methods are explained and compared in some detail here.

Output density and update of augmentation sphere boundary conditions

A second pass over the irreducible 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:

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 EF, 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.

This energy is a functional of nout and Vin. 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 Bi2Te3 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, nin−n, thus δE ~ (nin−n)2, the forces are not. If it were known how the density shifts with the nucleus (which requires knowledge of the dielectric function), the forces could be made variational as well.

We can do almost as well by making an ansatz for how the density shifts. The simplest ansatz is to assume that there is a cloud of charge that shifts rigidly with nucleus. Assuming that this charge is the free-atomic density, a correction term can be devised which dramatically improves on the convergence of the forces with respect to deviations nin−n. In many cases the error becomes almost variational, i.e. the error in the force to linear order in nin−n becomes much smaller. At self-consistency the correction term vanishes.

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 5 mRy/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 Bi2Te3 tutorial.

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:

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.