The Input File (CTRL)

This web page explains the structure of the main input file (called the ctrl file^†). You can automatically generate a template for this file from bare structural and chemical data using the blm utility.

^†Note: that the full name of the input file is ctrl.ext; you supply the extension on the command line. (If the extension is omitted, dat is used.) The same extension will be tacked onto names of most other files read or generated by the codes.

Input can also be supplied through a parallel input stream, namely the command line switches. Switches are flagged by command-line arguments beginning with - or --. They serve many purposes: some switches apply to all executables, others are specific to one or a few of them. Command-line arguments can also modify the contents of the input file described on this page: variables can be assigned from the command-line before the input file is parsed. Switches are documented in the command-line documentation for most executables; also any executable provides a brief summary of most switches it recognizes if you run it with --help, e.g. lmf --help.

As explained below, data is identified by a label called a token. A token is part of a tag, which is the full label with multiple parts, in the logical structure of a tree. The top-level or first part (branch), we denote as category; the last is the token, and data to be read in immediately follows the token. The tag does not itself appear in the input file, only the branches as such as the category and token as explained below. This web page documents the contents of each token, organized by category, in the same way the ctrl file is structured. A more detailed description of the syntax can be found in the input file manual. Note also that the reader does not parse lines directly as read from the ctrl file. It is first passed through a preprocessor, which can modify the contents of the input. See here for complete documentation of the preprocessor’s syntax.

Table of Contents

• 1. Input File Structure
  • Introduction
  • Tags, Categories and Tokens
  • Preprocessor
• 2. Help with finding tokens
• 3. The EXPRESS category
• 4. Input File Categories
  • Preliminaries
  • BZ
  • CONST
  • DYN
  • EWALD
  • HAM
  • GF
  • GW
  • DMFT
  • HEADER
  • IO
  • ITER
  • OPTICS
  • OPTIONS
  • PGF
  • SITE
  • SPEC
  • STR
  • START
    • Sample START category
  • STRUC
  • SYMGRP
    • Keywords in the SYMGRP category
  • VERS
  • Notes on gradient corrected functionals
  • Interactive mode
    • Simulacrum of interactive mode
    • Interactive mode with MPI
  • ITER_MIX
  • The ITER_MIX tag and how to use it
  • References

1. Input File Structure

Introduction

Here is a sample input file for the compound Bi₂Te₃ written for the lmf code.

    categories                  tokens


    VERS                    LM:7 FP:7
    HAM                     AUTOBAS[PNU=1 LOC=1 LMTO=3 MTO=1 GW=0]
                            GMAX=8.1
    ITER                    MIX=B2,b=.3  NIT=10  CONVC=1e-5
    BZ                      NKABC=3  METAL=5  N=2 W=.01
    STRUC
                            NSPEC=2  NBAS=5  NL=4
                            ALAT=4.7825489
                            PLAT= 1     0          4.0154392
                            -0.5   0.8660254  4.0154392
                            -0.5  -0.8660254  4.0154392
    SPEC
                            ATOM=Te         Z= 52  R= 2.870279
                            ATOM=Bi         Z= 83  R= 2.856141
    SITE
                            ATOM=Te         POS=  0.0000000   0.0000000   0.0000000
                            ATOM=Te         POS= -0.5000000  -0.8660254   1.4616199
                            ATOM=Te         POS=  0.5000000   0.8660254  -1.4616199
                            ATOM=Bi         POS=  0.5000000   0.8660254   0.8030878
                            ATOM=Bi         POS= -0.5000000  -0.8660254  -0.8030878

Each element of data follows a token. The token tells the reader what the data signifies.

Each token belongs to a category. In the sample input file above VERS, HAM, ITER, BZ, STRUC, SPEC, SITE are categories that organize the input by topic. Any text that begins in the first column is a category.

The full identifier is called the tag and it has the logical structure of a tree. The tag’s trunk (or top-level) is the category and the last is the token, e.g. GMAX associated with HAM and PLAT associated with STRUC. After the token comes the data to be parsed. In most cases category and token comprise comprise the entire tag, e.g. BZ_METAL. Thus the category groups tags into themes, the token identifies a particular type of data within the theme. Sometimes a tag has three branches, e.g. HAM_AUTOBAS_LOC.

Note: input files described here (ctrl.ext) can be automatically constructed from init files using the blm utility. init files and ctrl files are structured with categories and tokens in essentially the same way. For a another description of categories and tokens, see the init file documentation.

Tags, Categories and Tokens

The input file offers a very flexible free format: tags identify data to be read by a program, e.g.

    W=.01

reads a number (.01) from token W=. In this case W= belongs to the BZ category, so the full tag name is BZ_W.

A category holds information for a family of data, for example BZ contains parameters associated with Brillouin zone integration. The entire input system has at present a grand total of 18 categories, though any one program uses only a subset of them.

BZ  N=2     W=.01

BZ[ N=2     W=.01]

BZ  N=2
    W=.01
BZ  W=.01   N=2
BZ[ W=.01   N=2]

Usually the tag tree has only two levels (category and token) but not always. For example, data associated with atomic sites must be supplied for each site. In this case the tree has three levels, e.g. SITE_ATOM_POS. Site data is typically represented in a format along the following lines:

SITE    ATOM=Ga  POS=  0   0   0    RELAX=T
        ATOM=As  POS= .25 .25 .25
        ATOM=...
        ...
END

The scope of SITE starts at “SITE” and terminates just before “END”. There will be multiple instances of the SITE_ATOM tag, one for each site. The scope of the first instance begins with the first occurrence of ATOM and terminates just before the second:

ATOM=Ga  POS=  0   0   0    RELAX=T

And the scope of the second SITE_ATOM is

ATOM=As  POS= .25 .25 .25

Note that ATOM simultaneously acts like a token pointing to data (e.g. Ga) and as a tag holding tokens within it, in this case SITE_ATOM_POS and (for the first site) SITE_ATOM_RELAX.

Some tags are required; others are optional; still others (in fact most) may not be used at all by a particular program. If a code needs site data, SITE_ATOM_POS is required, but SITE_ATOM_RELAX is probably optional, or not read at all.

Note: this manual contains a more careful description of the input file’s syntax.

Preprocessor

Input lines are passed through a preprocessor, which provides a wide flexibility in how input files are structured. The preprocessor has many features in common with a programming language, including the ability to declare and assign variables, evaluate algebraic expressions; and it has constructs for branching and looping, to make possible multiple or conditional reading of input lines.

For example, supposing through a prior preprocessor instruction you have declared a variable range, and it has been assigned the value 3. This line in the input file:

    RMAX={range+1/4}

is turned in to:

    RMAX=3.25

The preprocessor treats text inside brackets {…} as an expression (usually an algebraic expression), which is evaluated and rendered back as an ASCII string. See this annotated lmf output for an example.

The preprocessor’s programming language makes it possible for a single file to serve as input for many materials systems in the manner of a database; or as documentation. Also you can easily vary input conditions in a parameteric fashion.

Other files besides ctrl.ext are first parsed by the preprocessor — files for site positions, Euler angles for noncollinear magnetism are read through the preprocessor, among others.

2. Help with finding tokens

Seeing the effect of the preprocessor

The preprocessor can act in nontrivial ways. To see the effect of the preprocessor, use the --showp command-line option.
See this annotated output for an example.

Finding what tags the parser seeks

It is often the case that you want to input some information but don’t know the name of the tag you need. Try searching this page for a keyword.

You can list each tag a particular tool reads, together with a synopsis of its function, by adding --input to the command-line.
Search for keywords in the text to find what you need.

  lmchk --input

 Tag                    Input   cast  (size,min)
 ------------------------------------------

 IO_VERBOS              opt    i4v      5,  1     default = 35
   Verbosity stack for printout.
   May also be set from the command-line: --pr#1[,#2]
 IO_IACTIV              opt    i4       1,  1     default = 0
   Turn on interactive mode.
   May also be controlled from the command-line:  --iactiv  or  --iactiv=no
...
 STRUC_FILE             opt    chr      1,  0
   (Not used if data read from EXPRESS_file)
   Name of site file containing basis and lattice information.
   Read NBAS, PLAT, and optionally ALAT from site file, if specified.
   Otherwise, they are read from the _ctrl_{:.path} file.
...
 STRUC_PLAT             reqd   r8v      9,  9
   Primitive lattice vectors, in units of alat
...
 SPEC_ATOM_LMX          opt    i4       1,  1     (default depends on prior input)
   l-cutoff for basis
...
 SITE_ATOM_POS          reqd*  r8v      3,  1
   Atom coordinates, in units of alat
 - If preceding token is not parsed, attempt to read the following:
 SITE_ATOM_XPOS         reqd   r8v      3,  1
   Atom coordinates, as (fractional) multiples of the lattice vectors

BZ_NKABC               reqd   i4v      3,  1
  (Not used if data read from EXPRESS_nkabc)
  No. qp along each of 3 lattice vectors.
  Supply one number for all vectors or a separate number for each vector.

Command-line options

--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. A more complete description of general-purpose command line options can be found on this page.
See this annotated lmfa output for an example.

Displaying tags read by the parser

To see what is actually read by a particular tool, run your tool with --show=2 or --show.
See the annotated lmf output for an example.

These special modes are summarized here.

3. The EXPRESS category

Section 3 provides some description of the input and purpose of tags in each category.

There is one special category, EXPRESS, whose purpose is to simplify and streamline input files. Tags in EXPRESS are effectively aliases for tags in other categories, e.g. reading EXPRESS_gmax reads the same input as HAM_GMAX.

If you put a tag into EXPRESS, it will be read there and any tag appearing in its usual location will be ignored. Thus including GMAX in HAM would have no effect if gmax is present in EXPRESS.

EXPRESS collects the most commonly used tags in one place. There is usually a one-to-one correspondence between the tag in EXPRESS and its usual location. The sole exception to this is EXPRESS_file, which performs the same function as the pair of tags, STRUC_FILE and SITE_FILE. Thus in using EXPRESS_file all structural data is supplied through the site file.

4. Input File Categories

This section documents the tokens read from the input file, arranged by category. Remember that each executable reads its only tokens specific to it.

Preliminaries

Note: The tables below list the input systems’ tokens and their function. Tables are organized by category.

• The Arguments column refers to the cast belonging to the token (“l”, “i”, “r”, and “c” refer to logical, integer, floating-point and character data, respectively)
• The Program column indicates which programs the token is specific to, if any
• The Optional column indicates whether the token is optional (Y) or required (N)
• The Default column indicates the default value, if any
• The Explanation column describes the token’s function.

See Table of Contents

BZ

Category BZ holds information concerning the numerical integration of quantities such as energy bands over the Brillouin Zone (BZ). The LMTO programs permit both sampling and tetrahedron integration methods. Both are described in bzintegration, and the relative merits of the two different methods are discussed. As implemented both methods use a uniform, regularly spaced mesh of k-points, which divides the BZ into microcells as described here. Normally you specify this mesh by the number of divisions of each of the three primitive reciprocal lattice vectors (which are the inverse, transpose of the lattice vectors PLAT); NKABC below.

These tokens are read by programs that make hamiltonians in periodic crystals (lmf,lm,lmgf,lmpg,tbe). Some tokens apply only to codes that make energy bands, (lmf,lm,tbe).

See Table of Contents

CONST

This category enables users to declare variables in algebraic expressions. The syntax is a string of declarations inside the category, e.g:

    CONST   a=10.69 nspec=4+2

Variables declared this way are similar to, but distinct from variables declared for the preprocessor, such as

    % const nbas=5

In the latter case the preprocessor makes a pass, and may use expressions involving variables declared by e.g. “% const nbas=5” to alter the structure of the input file.

Variables declared for use by the preprocessor lose their definition after the preprocessor completes.

The following code segment illustrates both types:

% const nbas=5
CONST   a=10.69 nspec=4
STRUC    ALAT=a  NSPEC=nspec NBAS={nbas}

After the preprocessor compiles, the input file appears as:

CONST   a=10.69 nspec=4
STRUC    ALAT=a  NSPEC=nspec NBAS=5

When the  CONST  category is read (it is read before other categories), variables  a  and  nspec  are defined and used in the  SPEC  category.

See Table of Contents

DYN

Contains parameters for molecular statics and dynamics. For a tutorial with molecular statics, see this page.

The following are specific to lmmag.

See Table of Contents

EWALD

Category EWALD holds information controlling the Ewald sums for structure constants entering into, e.g. the Madelung summations and Bloch summed structure constants (lmf). Most programs use quantities in this category to carry out Ewald sums (exceptions are lmstr and the molecules code lmmc).

See Table of Contents

HAM

This category contains parameters defining the one-particle hamiltonian.

Portions of HAM are read by these codes:

lm,  lmfa,  lmfgwd,  lmfgws,  lmf,  lmgf,  lmpg,  lmdos,  lmchk,  lmscell,  lmstr,  lmctl,  lmmc,  tbe,  lmmag.

See Table of Contents

GF

Category GF is intended for parameters specific to the Green’s function code lmgf. and is read by that code. See the Green’s function web page and also the Introductory Tutorial for lmgf.

See Table of Contents

GW

Category GW holds parameters specific to GW calculations, particularly for the GW driver lmfgwd. Most of these tokens will supply values for tags in the general category of the GWinput template when lmfgwd generates it (--job -1 --classicGWinput). However, in contrast to Questaal’s classical GW implementation, most of these tags can be read from the ctrl file; thus --classicGWinput works but it not recommended and is rarely needed. Any tags present in GWinput will take precedence over input read from the ctrl file.

See this page for a more detailed description of tags in the GWinput’s general category, and how tags there correspond to the ones documented below.

See Table of Contents

DMFT

Category DMFT holds parameters for the inteface to the DMFT, particularly for the DMFT driver lmfdmft. Unless otherwise specified, the only code reading tags from the DMFT category is lmfdmft .

The following tokens are read for each inequivalent correlated subblock. Data sandwiched between successive occurences of token BLOCK within DMFT apply to one DMFT correlated block.

See Table of Contents

HEADER

This category is optional, and merely prints to the standard output whatever text is in the category. For example:

HEADER  This line and the following one are printed to
        standard output whenever a program is run.
NEXT

Alternately:

HEADER [ In this form only two lines reside within the
        category delimiters,]
        and only two lines are printed.

See Table of Contents

IO

This optional category controls what kind of information, and how much, is written to the standard output file.

See Table of Contents

ITER

The ITER category contains parameters that control the requirements to reach self-consistency.

It applies to all programs that iterate to self-consistency: lm,  lmf,  lmmc,  lmgf,  lmpg,  tbe,  lmfa.

A detailed discussion can be found at the end of this document.

See Table of Contents

OPTICS

Optics functions available with the ASA extension packages OPTICS.

It is read by lm and lmf.

Token | Arguments | Program | Optional | Default | Explanation

• | - | - | - | - | - MODE | i | OPTICS | Y | 0 | 0: make no optics calculations
  1: generate linear $\varepsilon_2$
  20: generate second harmonic ε
   Example: optics/test/test.optics sic
  The following cases (MODE<0) generate joint or single density-of-states.
  Note: MODE<0 works only with LTET=3 described below.
  −1: generate joint density-of-states
   Examples:
   (ASA) optics/test/test.optics --all 4
   (FP) fp/test/test.fp zbgan
  −2: generate joint density-of-states, spin 2
   Example:optics/test/test.optics fe 6
  −3: generate up-down joint density-of-states
  −4: generate down-up joint density-of-states
  −5: generate spin-up single density-of-states
   Example: optics/test/test.optics --all 7
  −6: generate spin-dn single density-of-states LTET | i | OPTICS | Y | 0 | 0: Integration by Methfessel-Paxton sampling
  1: standard tetrahedron integration
  3: enhanced tetrahedron integration
  Note: In the metallic case, states near the Fermi level must be treated with partial occupancy. LTET=3 is the only scheme that handles this properly.
  It was adapted from the GW package and has extensions, e.g. the ability to handle non-vertical transitions $k^{occ} \ne k^{unocc}$. WINDOW | r1 r2 | OPTICS | N | 0 1 | Energy (frequency) window over which to calculate Im[ε(ω)].
  Im ε is calculated on a mesh of points $\omega_i$.
  The mesh spacing is specified by NPTS or DW, below. NPTS | i | OPTICS | N | 501 | Number of mesh points in the energy (frequency) window. Together with WINDOW, NPTS specifies the frequency mesh as:
  $\omega_i$ = WINDOW(1) + DW×(i−1)
  where DW = (WINDOW(2)−WINDOW(1))/(NPTS−1)
  Note: you may alternatively specify DW below. DW | r1 [r2] | OPTICS | Y | | Frequency mesh spacing DW[,OMG]. You can supply either one argument, or two.
  If one argument (DW) is supplied, the mesh will consist of evenly spaced points separated by DW.
  If a second argument (OMG) is supplied, points are spaced quadratically as:
  $\omega_i$ = WINDOW(1) + DW×(i−1) + [DW×(i−1)]²/OMG/2
  Spacing is approximately uniform up to frequency OMG; beyond which it increases linearly.
  Note: The quadratic spacing can be used only with LTET=3. FILBND | i1 [i2] | OPTICS | Y | 0 no. electrons | i1[,i2] occupied energy bands from which to calculate ε using first order perturbation theory, without local fields.
  i1 = lowest occupied band
  i2 = highest occupied band (defaults to no. electro ns) EMPBND | i1 [i2] | OPTICS | Y | 0 no. bands | i1[,i2] empty energy bands from which to calculate ε using first order perturbation theory, without local fields.
  i1 = lowest unoccupied band
  i2 = highest unoccupied band (defaults to no. bands) PART | i | OPTICS | Y | 0 | Resolve ε or joint DOS into band-to-band contributions, or by k.
  Result is output into file popt.ext.
  0. No decomposition
  1. Resolve ε or DOS into individual (occ,unocc) contributions
   Example: optics/test/test.optics ogan 5
  2. Resolve ε or DOS by k
   Example: optics/test/test.optics --all 6
  3. Both 1 and 2
  Add 10 to write popt as a binary file. CHI2[..] | | lm | Y | | Tag containing parameters for second harmonic generation.
  Not calculated unless tag is parsed.
   Example: optics/test/test.optics sic CHI2_NCHI2 | i | lm | N | 0 | Number of direction vectors for which to calculate χ₂, i.e. the nonlinear susceptibility tensor. CHI2_AXES | i1, i2, i3 | lm | N | | Direction vectors for each of the NCHI2 sets ESCISS | r | OPTICS | Y | 0 | Scissors operator (constant energy added to unoccupied levels, in Ry) ECUT | r | OPTICS | Y | 0.2 | Energy safety margin for determining (occ,unocc) window.
  lmf will attempt to reduce the number of (occ,unocc) pairs by restricting, for each k, transitions that contribute to the response, i.e. to those inside the optics WINDOW.
  The window is padded by ECUT to include states outside, but near the edge of the window.
  States outside window may nevertheless make contribution, e.g. because they can be part of a tetrahedron that does contribute.
  If you do not want lmf to restrict the range, use ECUT<0. NMP | i | OPTICS | Y | BZ_N | If present, supersedes BZ_N for the optics energy integration W | r | OPTICS | Y | BZ_W | If present, supersedes BZ_W for energy integration entering into the dielectric function MEFAC | i | OPTICS | Y | 0 | Contribution from nonlocal self-energy to velocity operator.
  1. include $d\Sigma/dk$
  2. approximate correction to $\mathrm{Im}\,\epsilon$ using ratio of QP to LDA eigenvalues. (Approximation is exact if LDA and QP eigenvalues are the same). FFMT | i | OPTICS | Y | 0 | Governs formatting of optics file
  0. fortran F format
  1. fortran E format IQ | i1, i2, i3 | OPTICS | Y | 0 | q vector for JDOS(q), in multiples of qlat/BZ_NKABC ESMR | r | OPTICS | Y | 0.05 | Energy smearing width for determining (occ,unocc) window. States are excluded for which occ<EF-ESMR or unocc>EF+ESMR. ALLTRANS | l | OPTICS | Y | F | Do not limit allowed transitions to occ<EF-ESMR and unocc>EF+ESMR FERMI | r | OPTICS | Y | NULL | If not NULL, supersede calculated Fermi level with given value when calculating dielectric function. IMREF | r1 r2 | OPTICS | Y | NULL | If not NULL, quasi-Fermi levels for occ and unocc states (nonequilibrium optics) KT | r | OPTICS | Y | - | Temperature for Fermi functions (Ry). Used when NMP<0. DELRE | r | gw suite | Y | 0.002, 0.08 | Equivalent to DW above, but used in the optics branch of the GW code. EPSQ | r,… | gw suite | Y | 0,…| List of q points when calculating the RPA dielectric function using the GW code. BSE | r | GW | Y | | Parameters associated with BSE-specific input for optics BSE_ETA | r | GW | Y | | Broadening for BSE dielectric function. The first parameter is use when computing the dielectric function for finite q; both are used when computing ε at q=0. BSE_QVECS| r,… | gw suite | Y | 1,0,0 | Direction vectors for q=0 BSE optics BSE_EMESH| r,… | gw suite | Y | 0,1,1e-3 | Energy mesh for q=0 BSE optics. Starting and ending energy, and mesh spacing.

See Table of Contents

OPTIONS

Portions of OPTIONS are read by these codes:

lm,  lmfa,  lmfgwd,  lmfgws,  lmf,  lmmc,  lmgf,  lmdos,  lmstr,  lmctl,  lmpg,  tbe.

See Table of Contents

PGF

Category PGF concerns calculations with the layer Green’s function program lmpg.

It is read by lmpg and lmstr.

See Table of Contents

SITE

Category SITE holds site information. As in the SPEC category, tokens must read for each site entry; a similar restriction applies to the order of tokens. Token ATOM= must be the first token for each site, and all tokens defining parameters for that site must occur before a subsequent ATOM=.

See Table of Contents

SPEC

Category SPEC contains species-specific information. Because data must be read for each species, tokens are repeated (once for each species). For this reason, there is some restriction as to the order of tokens. Data for a specific species (Z=, R=, R/W=, LMX=, IDXDN= and the like described below) begins with a token ATOM=;  input of tokens specific to that species must precede the next occurence of ATOM=.

The following tokens apply to the automatic sphere resizer:

The following tokens are input for each species. Data sandwiched between successive occurences of ATOM apply to one species.

See Table of Contents

STR

Category STR contains information connected with real-space structure constants, used by the ASA programs. It is read by lmstr, lmxbs, lmchk, and tbe.

*IINV parameters govern iterative solutions to screened strux

See Table of Contents

START

This category is specific to the ASA. It controls whether the code starts with moments P,Q or potential parameters; also P,Q may be input in this category. It is read by lm, lmgf, lmpg, and tbe.

Sample START category

The following is taken for the distribution’s test of La₂ Cu O₄.

START   BEGMOM=T  CNTROL=T RDVES=t
        ATOM=LA       P=  6.3055046  6.3000000  5.2308707
                      Q=  0.4770507  0.0000000  0.0610692
                          0.9882047 -0.3905638  0.2327244
                          2.0252993  0.0000000  0.1272500
        ATOM=CU       P=  4.6331214  4.3438861  3.8947075
                      Q=  0.4910799  0.0000000  0.0974578
                          0.6087341  0.0000000  0.1140513
                          9.4164169  0.0000000  0.2018023
        ATOM=OX       P=  2.8833091  2.8438183  3.1896353
                      Q=  1.6741779  0.0000000  0.0653497
                          4.2304006  0.0000000  0.1036699
                          0.0404676  0.0000000  0.0023966
        ATOM=OX2      P=  2.8840328  2.8447249  3.1806967
                      Q=  1.6660490  0.0000000  0.0257208
                          4.1318836  0.0000000  0.0365535
                          0.0083512  0.0000000  0.0003608

Notes on parsing P and Q

In the ASA, knowledge of P and Q is sufficient to completely determine the ASA density. Several ways are available to read these important quantities. The parser returns (P,Q) as a set according the following priorities:

• The (P,Q) set is read from the disk, if supplied, (along possibly with other quantities such as potential parameters El, C, Δ, γ.) One file is created for each class that contains this data and other class-specific information. Some or all of the data may be missing from the disk files. Alternatively, you may read these data from a restart file rsta.ext, which if it exists contains data for all classes in one file. The program will not read this data by default; use --rs=1 to have it read from the rsta file. To write class data to rsta, use --rs=,1 ( must be 0 or 1)

• If START_CONTRL=T, (P,Q) (and possibly other quantities) are read from START for classes you supply (usually all classes). Data read from this category supersedes any that might have been read from disk. If class data read from either of these sources, the input system returns it. For classes where none is available the parser will pick a default:

  • If data from a different class but in the same species is available, use it.

  • Otherwise use some preset default values for (P,Q).

After a calculation finishes you can run lmctl to read (P,Q) from disk and format it in a form ready to insert into the START category,e.g.

ATOM=SI       P=  3.8303101  3.7074067  3.2545634
              Q=  1.1694276  0.0000000  0.0297168
                  1.8803181  0.0000000  0.0489234
                  0.1742629  0.0000000  0.0063520
ATOM=ES       P=  1.4162942  2.2521617  3.1546386
              Q=  0.2873686  0.0000000  0.0129888
                  0.3485430  0.0000000  0.0165416
                  0.1400664  0.0000000  0.0055459

Thus all the information needed to generate a self-consistent ASA density can be embedded in the ctrl file.

Because the P’s float to the band center-of gravity (i.e. center of gravity of the occupied states for a particular site and $\ell$ channel) the corresponding first moments Q₁ vanish. P’s are floated by default since it minimizes the linearization error.

Caution: Sometimes it is necessary to override this default: If the band CG (of the occupied states) is far removed from the natural CG of a particular channel, you must restrict how far P can be shifted to the band CG. In some cases, allowing P to float completely will result in “ghost bands”.

The high-lying Ga 4d state is a classic example. To restrict P to a fixed value, see SPEC_ATOM_IDMOD.

In such cases, you want to pick the fractional part of P to be small, but not so low as to cause problems (about 0.5 for s orbitals and 0.15 for d orbitals; see here).

See Table of Contents

STRUC

By default structural information is read through the ctrl file. But some of the essential data can be read in multiple ways, in particular from site file. Questaal has utilities that will import this information from other formats such as cif files.

†Information may be obtained from a site file

See Table of Contents

SYMGRP

Category SYMGRP provides symmetry information; it helps in two ways. First, it provides the relevant information to find which sites are equivalent, this makes for a simpler and more accurate band calculations. Secondly, it reduces the number of k-points needed in Brillouin zone integrations.

Normally you don’t need SYMGRP; the program is capable of finding its own symmetry operations. However, there are cases where it is useful or even necessary to manually specify them. For example when including spin-orbit coupling or noncollinear magnetism where the symmetry group isn’t only specified by the atomic positions. In this case you need to supply extra information.

You can use SYMGRP to explicitly declare a set of generators from which the entire group can be created. For example, the three operations R4X, MX and R3D are sufficient to generate all 48 elements of cubic symmetry.

Unless conditions are set for noncollinear magnetism and/or SO coupling, the inversion is assumed by default as a consequence of time-reversal symmetry.

A tag describing a generator for a point group operation has the form O(nx,ny,nz) where O is one of M, I or Rj, or E, for mirror, inversion j-fold rotation and identity operation, respectively. nx,ny,nz are a triplet of indices specifying the axis of rotation. You may use X, Y, Z or D as shorthand for (1,0,0), (0,1,0), (0,0,1), and (1,1,1) respectively. You may also enter products of rotations, such as I*R4X.

Thus

SYMGRP  R4X MX R3D

specifies three generators (4-fold rotation around x, mirror in x, 3-fold rotation around (1,1,1)). Generating all possible combinations of these rotations will result in the 48 symmetry operations of the cube.

To suppress all symmetry operations, use

SYMGRP E

In the ASA, owing to the spherical approximation to the potential only the point group is required for self-consistency.

But in general you must specify the full space group. The translation part gets appended to rotation part in one of the following forms:  :(x1,x2,x3)  or alternatively  ::(p1,p2,p3)  with the double ‘::’. The first defines the translation in Cartesian coordinates in units of ALAT, second in crystal coordinates.

These two lines (taken from testing/ctrl.cr3si6) provide equivalent specifications:

SYMGRP   r6z:(0,0,0.4778973) r2(1/2,sqrt(3)/2,0)
SYMGRP   r6z::(0,0,1/3)      r2(1/2,sqrt(3)/2,0)

Keywords in the SYMGRP category

SYMGRP accepts, in addition to symmetry operations the following keywords:

• find tells the program to determine its own symmetry operations. Thus:

  SYMGRP find
  

  amounts to the same as not incuding a SYMGRP category in the input at all

  You can also specify a mix of generators you supply, and tell the program to find any others that might exist. For example:

  SYMGRP r4x find
  

  specifies that 4-fold rotation be included, and  find  tells the program to look for any additional symops that might exist.

• AFM: For certain antiferromagnets, certain translation operations exist provided the rotation/shift is accompanied by a spin flip. Say a translation of (-1/2,1/2,1/2)a restores the crystal structure, but all atoms after translation have opposite spin. Specify this symmetry with:

  SYMGRP   ...       AFM::-1/2,1/2,1/2
  

  This operation is used only by lmf.

• SOC or SOC=2: Tells the symmetry group generator to exclude operations that do not preserve the z axis. This is used particularly for spin-orbit coupling where the crystal symmetry is reduced (z is the quantization axis). SOC=2 is like SOC but allows operations that preserve z or flip z to −z. This works in some cases.

  Note: This keyword is only active when the two spin channels are linked, e.g. SO coupling or noncollinear magnetism.

• GRP2 turns on a switch that can force the density among inequivalent classes that share a common species to be averaged. In the ASA codes the density is spherical and the averaging is complete; in the FP case only the spherical part of the densities can be averaged. This helps sometimes with stabilizing difficult cases in the path to self-consistency. You specify which species are to be averaged with the SPEC_ATOM_GRP2 token.

  GRP2 averages the input density; GRP2=2 averages the output density; GRP2=3 averages both the input and the output density.

• RHOPOS turns on a switch that forces the density positive at all points.
  You can also accomplish this with the command-line switch --rhopos.

See Table of Contents

VERS

This category is used for version control. As of version 7, the input file must have the following tokens for any program in the suite:

VERS LM:7

It tells the input system that you have a v7 style input file.

For a particular program you need an additional token to tell the parser that this file is set up for that program. Thus your VERS category should read:

VERS LM:7 ASA:7              for lm, lmgf or lmpg
VERS LM:7 FP:7               for lmf or lmfgwd
VERS LM:7 MOL:3              for a molecules codes such as lmmc
VERS LM:7 TB:9               for the empirical tight-binding tbe
                             and so on.

Add version control tokens for whatever programs your input file supports.

See Table of Contents

Notes on gradient corrected functionals

The semilocal exchange-correlation potential is defined as $$v_{j}^{GGA}(n,\mathbf{r}) = v_{j}^{HEG}(n,\mathbf{r}) * F_{j}^{GGA}(n,\mathbf{r}) + \int d\mathbf{r}' \varepsilon^{HEG}_j(n,\mathbf{r}') \frac{\partial F_{j}^{GGA}(n,\mathbf{r}')}{\partial s(n,\mathbf{r}')} \frac{\partial s(n,\mathbf{r}')}{\partial n(\mathbf{r})} \quad\quad \text{with } j=x,c$$

where $v_{j}^{HEG}[n](r)$ and $\varepsilon^{HEG}_j[n](r)$ are the exchange or correlation potentials and energy densities, respectively, for the Homogeneous Electron Gas (HEG), which are determined by the options XCFUN=1 or XCFUN=2. $F_{j}^{GGA}[n](r)$ is the exchange or correlation GGA enhancement factor which introduce semilocal effects and can be choosen using the token GGA. $s(n,r)$ is the reduced gradient. Then, the final exchange-correlation potential/functional is determined using two different tokens, i.e. one for the local (HEG) part and the other one for the semilocal (GGA) correction.

Note that hybrid and meta-GGA exchange-correlation functional schemes are not implemented in the QUESTAAL package. As consequence, only LDA and GGA libXC functionals can be used.

Interactive mode

For the most part Questaal codes are designed to run without receiving information through standard input. The various editors are an exception (though even editor instructions can be run in batch mode; see e.g. dynamical self-energy editor tutorial).

It is often convenient to have some interactive facility, e.g. whether to limit the number of iterations in a self-consistency cycle. The Questaal codes have a interactive mode, which you can turn on with token IO_IACTIV in the ctrl file, or on the command-line with the  --iactive  switch.

For example, if you run lm or lmf interactively you will be prompted with

QUERY: beta (def=0.3)?

and the program will wait for input. To see what your options are, enter  ? <RET> You should see

 (A)bort  (I)active  (V)erb  (C)pu  (T)iming  (W)ork  (S)et value
 QUERY: beta (def=0.3)?

Here it is asking if you want to modify the existing value for the charge mixing parameter  beta.
Enter one of:

• a   Program aborts execution
• i   Toggles interactive mode
• v  #  Sets verbosity to  #
• c   Prints out CPU usage so far
• t   Turns on timing printout
• w   Not used now
• s  #  Sets parameter to  #

If you enter

s .5 <RET>

in this instance, the program will modify its value for beta to 0.5 and prompt you again. If you don’t want to make any (further) changes, just enter  <RET>.

The most commonly changed parameter is the number of iterations, called  maxit. You can increase or decrease it; if you decrease maxit below the current iteration number, the program will stop.

Simulacrum of interactive mode

You can enter interactive mode instructions that are normally read by the standard input by entering them into a file iact.ext.

The executable first looks for that file, reads its contents, and executes its instructions before prompting you. It will perform the instruction, e.g. set the value of the parameter, without turning on the true interactive mode. However if you do turn it on, e.g. put  i  into iact.ext, the program will revert to true interactive mode and prompt you for instructions.

There is an important difference in the normal and simulacrum operations when setting a parameter. In the latter case, you must tell the program which parameter. Do this by naming the parameter after the  s. Thus iact.ext would contain a line like

s maxit 3

Interactive mode with MPI

Normal interactive mode is not available when running with multiple processors. The simulacrum mode does work, but only for a subset of parameters. Most importantly,  maxit is one parameter that is read; thus you can adjust the number of iterations a job will do after execution starts.

slatsm/query.f contains the source code controlling this mode.

See Table of Contents

ITER_MIX

ITER_MIX is a token in the ITER category that controls how Questaal codes iterate to self-consistency in the charge density. Its contents is a string consisting of mixing options separated by a delimiter, as described here.

Questaal codes follow the usual procedure of mixing a linear combination of input density  nⁱⁿ and output density  n^out to make a trial guess  n^*  for the self-consistent density (see for example Chapter 9 in Richard Martin’s book¹. Questaal uses two independent techniques to accelerate convergence to the self-consistency condition  n^out→nⁱⁿ. First, the quantities are mixed making use of a model for the dielectric function. Second, multiple (nⁱⁿ, n^out) pairs (taken from prior iterations) can be used to accelerate convergence. They are typically taken in combination and the contents of ITER_MIX control options for both kinds of approaches.

Both ideas are explained on this page. See also the ASA NbFe superlattice tutorial.

The ITER_MIX tag and how to use it

In practice, mixing proceeds as described on this page, but additionally by combining multiple instances of (Xⁱⁿ,X^out) pairs. Each iteration produces a new pair, and methods exist to use this information and take linear combinations of them from both current and prior iterations, that is a better choice than any one pair. You can choose between Broyden³ and Anderson² methods. The string belonging to ITER_MIX should begin with one of

  MIX=An
  MIX=Bn

which tells the mixer which scheme to use. slatsm/amix.f describes the mathematics behind the Anderson scheme. n is the maximum number of prior iterations to include in the mix. As programs proceed to self-consistency, they dump prior iterations to disk, to read them the next time through. Data is I/O to mixm.ext.

The Anderson scheme is particularly simple to monitor. How much of δX from prior iterations is included in the final mixed vector is printed to stdout as parameter tj, e.g.

   tj: 0.47741                           ← iteration 2
   tj:-0.39609  -0.44764                 ← iteration 3
   tj:-0.05454   0.01980                 ← iteration 4
   tj: 0.24975
   tj: 0.48650
   tj:-1.34689

In the second iteration, one prior iteration was mixed; in the third and fourth, two; and after that, only one. (When the normal matrix picks up a small eigenvalue the Anderson mixing algorithm reduces the number of prior iterations).

Consider the case when a single prior iteration was mixed.

• If tj=0, the new X is entirely composed of the current iteration. This means self-consistency is proceeding in an optimal manner.
• If tj=1, it means that the new X is composed 100% of the prior iteration. This means that the algorithm doesn’t like how the mixing is proceeding, and is discarding the current iteration. If you see successive iterations where tj is close to (or worse, larger than) unity, you should change something, e.g. reduce beta.
• If tj<0, the algorithm thinks you can mix more of X^out and less of Xⁱⁿ. If you see successive iterations where tj is significantly negative (less than −1), increase beta.

Broyden mixing³ uses a more sophisticated procedure, in which it tries to build up the Hessian matrix. It usually works better but has more pitfalls than Anderson. Broyden has an additional parameter,  wc, that controls how much weight is given to prior iterations in the mix (see below).

Remember also that (X^out − Xⁱⁿ) is screened. In a simple metal, the Lindhard function pretty well describes the actual dielectric function, and tj should be small, as see in this tutorial. As for the dielectric function in the ASA, codes (lm, lmgf, lmpg) offer two options:

1. A rough ε is obtained from eigenvalues of the Madelung matrix (OPTIONS_SCR=4).
2. The q=0 discretized polarization is explicitly calculated (OPTIONS_SCR=11 generates it and OPTIONS_SCR=2 uses it; see OPTIONS_SCR).

The general syntax is for ITER_MIX is

An[,b=beta][,locm=locm][,b2=b2][,bv=betv][,n=nit][,w=w1,w2][,idpmx=#][,fn=name][,k=nkill][,elind=#][,elin2=#][;...]  or
Bn[,b=beta][,locm=locm][,b2=b2][,bv=betv][,wc=wc][,idpmx=#][,n=nit][,w=w1,w2][,fn=name][,elind=#][,elin2=#][,k=nkill]

The options are described below. They are parsed in routine subs/parmxp.f. Parameters (b, wc, etc.) may occur in any order.

• An or Bn:  maximum number of prior iterations to include in the mix (the mixing file may contain more than n prior iterations).
  n=0 implies linear mixing. Default: B2.

• b=beta:  the mixing parameter beta in Eq. 4 above. Default: 0.3.

• b2=b2:  (lmf only). mixing parameter for the “residual” density (see locm, next).

• locm=locm:  (lmf only). Options to control how the charge is mixed. locm is an integer where each digit has an independent meaning. In the discussion below, let locm0 be the 1s digit of locm, locm1 the 10s digit, and so on. The mixing algorithm is sophisticated and it has many branches to help stabilize convergence, which you control with locm. A full explanation is somewhat involved; see fp/mixrho.f for more details, and also this page.

  In density-functional theory the potential Vⁱⁿ is a functional of the input density nⁱⁿ. For a given Vⁱⁿ, eigenfunctions and n^out can be generated. Self-consistency is reached when a particular nⁱⁿ=n^* is found for for which n^*(in)=n^*(out). In Green’s function theory, n^out is a functional of Gⁱⁿ, which contains more information than nⁱⁿ, but in any case the self-consistency requires the condition n^out=nⁱⁿ. The mixing scheme generates an estimate for n^* that meets this condition.

  Questaal carries a three-fold representation of the charge density with three components: the smooth density n₀ represented on a uniform mesh, true site densities n₁, and finally site densities n₂ that are a one-center expansions of n₀. One unfortunate complication of Questaal’s implementation is that it complicates charge mixing to reach self-consistency: the three components must all be mixed. By default they are all mixed together as one large array using Anderson or Broyden (or linear) mixing. Most of the time the default algorithm suffices, but there are cases where it does not. This is in part because the system is overdetermined: for example a change in the sum of site densities n₁+n₂, keeping n₁−n₂ fixed, has little effect on the total potential.

  In brief, each of the three density components n₀, n₁, n₂ is partitioned into a “main” part and a “residual.” The “main” density is the part for which the potential is most sensitive: these portions of n₀, n₁, n₂ are grouped together and mixed with an Anderson/Broyden mixing scheme (shortened here to “A/B mix”). The “residual” is what is left over. It is linearly mixed with mixing parameter b2, as distinct from the mixing parameter beta used for the main density. For the meaning of beta or b2 in a linear mixing scheme, see Eq. (4) on this page.

  Site densites are expanded spherical harmonics, specifically a sum of [$\ell$ dependent radial functions] × [spherical harmonics $\Upsilon_{lm}(\hat{\mathbf{r}})$ (the real form of true spherical harmonics Y_lm)]. Radial functions are tabulated numerically on a shifted logarithmic mesh. For the site densities n₁, n₂ the most sensitive part is the $\ell$=0 radial function as the total charge is a function of it. The $\ell$=0 component of whatever combination of n₁, n₂ is mixed is always included in the main density in some form. Depending on locm0 (see below), the $\ell$>0 components may be included in the residual, or as part of the main density in a compacted form. Normally the smooth or mesh density n₀ is also included in the main density.

  The output densities mixed here are not n₀, n₁, n₂, but ones screened by the model Lindhard function. (See here for how the change (output-input) density can be screened by the Lindhard function.) These densities are a function of n^out and the Lindhard parameter, elind. In this section we will designate the three input densities r1in, r2in, rsmin, and the three screened output densities r1out, r2out, rsmout, while r1, r2, rsm apply to either kind. For mixing, r1,r2 are locally transformed into r1’ and r2’, depending on locm2 (see below).

  locm0: controls what combinations of r1, r2, rsm are A/B mixed

  0. the $\ell$=0 part of r1’ and r2’, and rsm. The residual is the $\ell$>0 part of r1’ and r2’
  1. the $\ell$=0 part of r1’ only, and rsm. The residual is the $\ell$>0 part of r1’ and all of r2’.
     Note: because r2’ is linearly mixed while rsm is not, this mode may not preserve charge neutrality
  2. the $\ell$=0 part of r2’ only, and rsm. The residual is the $\ell$>0 part of r2’ and all of r1’.
     Note: this is sensible when locm2=0 because r1’=r1+r2, which has little effect on the potential.
  3. the $\ell$=0 part of r1’ and r2’, and rsm, and $\ell$>0 part of r1’ and r2’, approximated as Generalized Gaussian functions $G_{kL}$ or their polynomial parts $P_{kL}$. The residual is the difference between the numerically tablulate representation and the compacted form of the $\ell$>0 parts of r1’ and r2’.
  4. Like mode 3, except all of r2’ is contained in the residual, as in mode 1
     Note: this is sensible when locm2=0 because r1’=r1+r2, which has little effect on the potential.
  5. Like mode 3, except all of r1’ is contained in the residual, as in mode 2

  locm1: scales r1’

  0. no scaling
  1. not used
  2. scale r1’ by 1/(1+(r/rg)²) [use with locm0=3]
  3. scale r1’ by gaussian exp(-(r/rg)²) [use with locm0=3]
  4. copy r2in’ to r2out’ ⇒ no mixing this channel
  5. copy r1in’ to r1out’ ⇒ no mixing this channel

  locm2: transformations r1,r2 into r1’,r2’

  0. r1’ = r1+r2 = r1’, r2’ = r1−r2 (generally preferred because r1+r2 has little on the total potential)
  1. No transformation

  locm3: transformations of the $\ell$=0 parts of the radial densities, entering into A/B mix

  0. The $\ell$=0 radial density is used as is.
  1. A Fourier transform is taken, and the Fourier transform is A/B mixed. (For linear mixing locm3=0 and locm3=1 are equivalent)
  2. The output density assumes the same shape as the input density but is scaled by a constant factor, determined by the ratio of input to output charges. The difference between the true rout and the scaled rin is included in the residual density. In this mode, the mesh density is similarly scaled by a constant factor.

  locm4: How the smooth mesh density is mixed

  0. The mesh density is included in the A/B mixing
  1. The mesh density is included in the residual.

• n=nit:  the number of iterations to use mix with this set of parameters before passing on to the next set. After the last set is exhausted, it starts over with the first set.

• name=fn:  mixing file name (mixm is the default). Must be eight characters or fewer.

• k=nkill:  kill mixing file after nkill iterations. This is helpful when the mixing runs out of steam, or when the mixing parameters change. Default: 7.

• wc=wc:  (Broyden only) that controls how much weight is given to prior iterations in estimating the Jacobian. wc=1 is fairly conservative. Choosing wc<0 assigns a floating value to the actual wc, proportional to −wc/rms-error. This increases wc as the error becomes small. wc defaults to −1 if it is not specified. See Johnson’s paper³ for the definition of  wc.

• w=w1,w2:  (spin-polarized calculations only) The up- and down- spin channels treated separately, the sum (up+down), i.e. charge and difference (up-down), i.e. spin densities are mixed. The two combinations are weighted by w1 and w2 in the mixing, more heavily emphasizing the more heavily weighted. As special cases, w1=0 freezes the charge and mixes the magnetic moments only while w2=0 freezes the moments and mixes the charge only.

• idpmx=#:  (spin-polarized calculations only, used with w=…) if  #>0, the charge and spin densities are mixed independently of each other. Additionally you can use the 10s digit to specify that one or another channel be mixed with linear mixing only. In these cases the mixing beta for the charge and spin channels become beta·w1 and beta·w2.
  idpmx=11:  Linearly mix the charge channel, Anderson/Broyden mix the spin channel
  idpmx=21:  Linearly mix the spin channel, Anderson/Broyden mix the charge channel

• elind=elind:  The Fermi energy entering into the Lindhard dielectric function: $\mathbf{elind} = k_F^2$.
  elind<0: Use the free-electron gas value, scaled by −elind. The default value is −1.

• elin2=elin2:  Performs same function as elind=elind: elind=elind is used for smooth density and elin2=elin2 for local densities.
  If not specified, elin2 defaults to elind2.

• wa:  (ASA only) weight for extra quantities included with P,Q in the mixing procedure. For noncollinear magnetism, includes the Euler angles.

• r=expr:  continue this block of mixing sequence until rms error < expr.

  Example:  MIX=A4,b=.2,k=4
  uses the Anderson mixing scheme², killing the mixing file each fourth iteration. The mixing  beta  is 0.2.

  You can string together several rules. One set of rules applies for a certain number of iterations; followed by another set.
  Rules are separated by a “ ; ”.

  Example:  MIX=B10,n=8,w=2,1,fn=mxm,wc=11,k=4;A2,b=1
  does 8 iterations of Broyden mixing, followed by Anderson mixing. The Broyden iterations weight the (up+down) double that of (up-down) for the magnetic case, and iterations are saved in a file which is deleted at the end of every fourth iteration. wc is 11. beta assumes the default value. The Anderson rules mix two prior iterations with beta=1.

See Table of Contents

References

¹ R. M. Martin, Electronic Structure, Cambridge University Press (2004).

² D. G. Anderson. Iterative procedures for nonlinear integral equations. J. Assoc. Comput. Mach., 12:547–560, 1965

³ D. D. Johnson, Phys. Rev. B 38, 12807 (1988).

⁴ (Questaal’s methods paper) Dimitar Pashov, Swagata Acharya, Walter R. L. Lambrecht, Jerome Jackson, Kirill D. Belashchenko, Athanasios Chantis, Francois Jamet, Mark van Schilfgaarde, Questaal: a package of electronic structure methods based on the linear muffin-tin orbital technique, Comp. Phys. Comm. 249, 107065 (2020).