Questaal Home

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


Here is a sample input file for the compound Bi2Te3 written for the lmf code.

    categories                  tokens

    VERS                    LM:7 FP:7
    HAM                     AUTOBAS[PNU=1 LOC=1 LMTO=3 MTO=1 GW=0]
    ITER                    MIX=B2,b=.3  NIT=10  CONVC=1e-5
    BZ                      NKABC=3  METAL=5  N=2 W=.01
                            NSPEC=2  NBAS=5  NL=4
                            PLAT= 1     0          4.0154392
                            -0.5   0.8660254  4.0154392
                            -0.5  -0.8660254  4.0154392
                            ATOM=Te         Z= 52  R= 2.870279
                            ATOM=Bi         Z= 83  R= 2.856141
                            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.


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.

Consider the Brillouin zone integration category. You plan to carry out the BZ integration using the Methfessel-Paxton sampling method. M-P integration has two parameters: polynomial order n and gaussian width w. Two tags are used to identify them: BZ_N and BZ_W; they are usually expressed in the input file as follows:

BZ  N=2     W=.01

This format style is the most commonly used because it is clean and easy to read; but it conceals the tree structure a little. The same data can equally be written:

BZ[ N=2     W=.01]

Now the tree structure is apparent: [..] delimits the scope of tag BZ.

Any tag that starts in the first column is a category, so any non-white character appearing in the first column automatically starts a new category, and also terminates any prior category. N= and W= mark tokens BZ_N and BZ_W.

Apart from the special use of the first column to identify categories, data is largely free-format, though there are a few mild exceptions. Thus:

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

all represent the same information.

Note: if two categories appear in an input file, only the first is used. Subsequent categories are ignored. Generally, only the first tag is used when more than one appears within a given scope.

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

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.


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:


is turned in to:


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.

Take for an example:

  lmchk --input

This switch tells the parser not to try and read anything, but print out information about what it would try to read. Several useful bits of information are given, including a brief description of each tag in the following format. A snippet of the output is reproduced below:

 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

The table tells you IO_VERBOS and IO_IACTIV are optional tags; default values are 35 and 0, respectively. A single integer will be read from the latter tag, and between one and five integers will be read from IO_VERBOS.

There is a brief synopsis explaining the functions of each. For these particular cases, the output gives alternative means to perform equivalent functions through command-line switches.

STRUC_FILE=fname is optional. Here fname is a character string: it should be the site file name fname.ext from which lattice information is read. If you do use this tag, other tags in the STRUC category (NBAS, PLAT, ALAT) may be omitted. Otherwise, STRUC_PLAT is required input; the parser requires 9 numbers.
The synopsis also tells you that you can specify the same information using EXPRESS_file=fname (see EXPRESS category below).

SPEC_ATOM_LMX is optional input whose default value depends on other input (in this case, atomic number).

SITE_ATOM_POS is required input in the sense that you must supply either it or SITE_ATOM_XPOS. The * in reqd* the information in SITE_ATOM_POS can be supplied by an alternate tag – SITE_ATOM_XPOS in this case.

Note: if site data is given through a site file, all the other tags in the SITE category will be ignored.

The cast (real, integer, character) of each tag is indicated, and also how many numbers are to be read. Sometimes tags will look for more than one number, but allow you to supply fewer. For example, BZ_NKABC in the snippet below looks for three numbers to determine the k-mesh, which are the number of divisions only each of the reciprocal lattice vectors. If you supply only one number it is copied to elements 2 and 3.

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.


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


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

GETQPl YFRead list of k-points from disk file qpts.ext. This is a special mode, and you normally would let the program choose its own mesh by specifying the number of divisions (see NKABC).
If token is not parsed, the program will attempt to parse NKABC.
NKABC1 to 3 i N The number of divisions in the three directions of the reciprocal lattice vectors. k-points are generated along a uniform mesh on each of these axes. (This is the optimal general purpose quadrature for periodic functions as it integrates the largest number of sine and cosine functions exactly for a specified number of points.)

The parser will attempt to read three integers. If only one number is read, the missing second and third entries assume the value of the first.

Information from NKABC, together with BZJOB below, contains specifications equivalent to the widely used “Monkhorst Pack” scheme. But it is more transparent and easier to understand. The number of k-points in the full BZ is the product of these numbers; the number of irreducible k-points may be reduced by symmetry operations.
PUTQPl YFIf T, write out the list of irreducible k-points to file qpts, and the weights for tetrahedron integration if available.
BZJOB1 to 3 i Y0Controls the centering of the k-points in the BZ:
0: the mesh is centered so that one point lies at the origin.
1: points symmetrically straddle the origin.

Three numbers are supplied, corresponding to each of the three primitive reciprocal lattice vectors. As with NKABC if only one number is read the missing second and third entries assume the value of the first.
I123l YTControls looping order when generating regular mesh of k-points.
T: inner loop is along QLAT(1); outer loop along QLAT(3)
F: inner loop is along QLAT(3); outer loop along QLAT(1)
METALilmf, lm, tbeY5Specifies how the weights are generated for Brillouin zone integration. For a detailed description, see this page.
The METAL token accepts the following:
0. System assumed to be an insulator; weights determined a priori.
1. Eigenvectors are written to disk, in which case the integration for the charge density can be deferred until all the bands are obtained.
2. Integration weights are read from file wkp.ext, which will have been generated in a prior band pass. If wkp.ext is unavailable, the program will temporarily switch to METAL=3.
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.
4. 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.
5. Like METAL=3 in which 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.
TETRAilmf,lm,tbeY1Selects BZ integration method.
0: Methfessel-Paxton sampling integration. Tokens NPTS, N, W, EF0, DELEF, DOS (see below) are relevant to this integration scheme.
1: tetrahedron integration, with Bloechl weights
Nilmf,lm,tbeY0Polynomial order for M-P sampling integration. (Not used with tetrahedron integration or for insulators).
 0: integration uses standard gaussian method.
>0: integration uses generalized gaussian functions, i.e. polynomial of order N × gaussian to generate integration weights.
−1: use the Fermi function rather than gaussians to broaden the δ-function. This generates the actual electron (fermi) distribution for a finite temperature.
Add 100: by default, if a gap is found separating occupied and unoccupied states, the program will treat the system as an insulator, even when METAL>0. To suppress this, add 100 to N (use −101 for Fermi distribution).
Wrlmf,lm,tbeY5e-3Case BZ_N≥0 :  broadening (Gaussian width) for Gaussian sampling integration (Ry).
Case BZ_N<0 :  kBT (Ry) where kB is the Boltzmann constant and T the temperature.
W is not used for insulators or when using tetrahedron integration.
EF0rlmf,lm,tbeY0Initial guess at Fermi energy. Used when  TETRA=0, or when  BZ_METAL=4 (which does not use the tetrahedron method for the density).
DELEFrlmf,lm,tbeY0.05Initial uncertainty in Fermi level for sampling integration.
Used when  TETRA=0, or when  BZ_METAL=4 (which does not use the tetrahedron method for the density). As the system approaches self-consistency this window is reduced.
ZBAKrlmf,lmY0Homogeneous background charge
SAVDOSilmf,lm,tbeY00: does not save dos on disk.
1: writes the total density of states on NPTS energy mesh points to disk file dos.ext.
2: Write weights to disk for partial DOS (does not work for lmf; in the ASA this occurs automatically).
4: Same as (2), but write weights m-resolved (ASA).
1. SAVDOS>0 uses DOS and NPTS tags also.
2. You may also cause lm or lmf to generate m-resolved dos the from command-line (see --pdos).
DOS2 r Y-1,0Energy window over which DOS accumulated (Ry). Needed either for sampling integration or if SAVDOS>0.
NPTSi Y1001Number of points in the density-of-states energy mesh used in conjunction with sampling integration. Needed either for sampling integration or if SAVDOS>0.
EFMAXrlmf,lm,tbeY2Only eigenvectors whose eigenvalues are less than EFMAX are computed; this improves execution efficiency.
NEVMXilmf,lm,tbeY0>0 : Find at most NEVMX eigenvectors.
=0 : program uses internal default.
<0 : no eigenvectors are generated (and correspondingly, nothing associated with eigenvectors such as density).

Caution: if you want to look at partial DOS well above the Fermi level (which usually comes out around 0), you must set EFMAX and NEVMX high enough to encompass the range of interest.
ZVALr Yall LDANumber of electrons to accumulate in BZ integration. Normally zval is computed by the program.
NOINVllmf,lm,tbeYFSuppress the automatic addition of the inversion to the list of point group operations. Usually the inversion symmetry can be included in the determination of the irreducible part of the BZ because of time reversal symmetry. There may be cases where this symmetry is broken:
e.g. when spin-orbit coupling is included or when the (beyond LDA) self-energy breaks time-reversal symmetry. In most cases, the program will automatically disable this addition in cases that it knows the symmetry is broken.
FSMOM2 rlmf,lmY0 0Set the global magnetic moment (collinear magnetic case). In the fixed-spin moment method, a spin-dependent potential shift Beff is added to constrain the total magnetic moment to value assigned by FSMOM. No constraint is imposed if this value is zero (the default).
Optional second argument #2 supplies an initial Beff. It is applied whether or not the first argument #1 is 0. If #1 ≠ 0, Beff is made consistent with it.
DMATKllmf,lmgfYFCalculate the density matrix. Implementation still not ready.
INVITllmf,lmYTGenerate eigenvectors by inverse iteration (this is the default). It is more efficient than the QL method, but occasionally fails to find all the vectors. When this happens, the program stops with the message:
DIAGNO: tinvit cannot find all evecs
If you encounter this message set INVIT=F.
EMESHrlmgf,lmpgY10,0,-1,…Parameters defining contour integration for Green’s function methods.
See also the GF documentation.
1. number of energy points n.
2. contour type:
 0: Uniform mesh of nz points: Real part of z between emin and emax
 1: Same as 0, but reverse sign of Im z
 10: elliptical contour
 11: same as 10, but reverse sign of Im z
 100s digit used for special modifications
 Add 100 for nonequil part using Im(z)=delne
 Add 200 for nonequil part using Im(z)=del00
 Add 300 for mixed elliptical contour + real axis to find fermi level
 Add 1000 to set nonequil part only.
3. lower bound for energy contour emin (on the real axis).
4. upper bound for energy contour emax, e.g. Fermi level (on the real axis).
5. (elliptical contour) eccentricity: ranges between 0 (circle) and 1 (line)
 (uniform contour) Im z.
6. (elliptical contour) bunching parameter eps : ranges between 0
 (distributed symmetrically) and 1 (bunched toward emax)
 (uniform contour) not used.
7. (nonequilibrium GF, lmpg) nzne = number of points on nonequilibrium contour.
8. (nonequilibrium GF, lmpg) vne = difference in fermi energies of right and left leads.
9. (nonequilibrium GF, lmpg) delne = Im part of E for nonequilibrium contour.
10 (nonequilibrium GF, lmpg) substitutes for delne when making the surface self-energy.
MULLitbeY0Mulliken population analysis. Mulliken population analysis is also implemented in lmf, but you specify the analysis with a command-line argument.

See Table of Contents


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

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


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

NITilmf, lmmc, tbeY maximum number of relaxation steps (molecular statics).
SSTAT[…] lm, lmgfY (noncollinear magnetism) parameters specifying how spin statics (rotation of quantization axes to minimze energy) is carried out.
SSTAT_MODEilm, lmgfN00: no spin statics or dynamics.
-1: Landau-Gilbert spin dynamics.
1: spin statics: quantization axis determined by making output density matrix diagonal.
2: spin statics: size and direction of relaxation determined from spin torque.
Add 10 to mix angles independently of P,Q (Euler angles are mixed with prior iterations to accelerate convergence).
Add 1000 to mix Euler angles independently of P,Q.
SSTAT_SCALEilm, lmgfN0(used with mode=2) scale factor amplifying magnetic forces.
SSTAT_MAXTilm, lmgfN0maximum allowed change in angle.
SSTAT_TAUilm, lmgfN0(used with mode=-1) time step.
SSTAT_ETOLilm, lmgfN0(used with mode=-1) Set tau=0 this iter if etot-ehf>ETOL.
MSTAT[…] lmf, lmmc, tbeY (molecular statics) parameters specifiying how site positions are relaxed given the internuclear forces.
MSTAT_MODEilmf, lmmc, tbeN00: no relaxation.
4: relax with conjugate gradients algorithm (not generally recommended).
5: relax with Fletcher-Powell alogirithm. Find minimum along a line; a new line is chosen. The Hessian matrix is updated only at the start of a new line minimization. Fletcher-Powell is more stable but usually less efficient then Broyden.
6: relax with Broyden algorithm. This is essentially a Newton-Raphson algorithm, where Hessian matrix and direction of descent are updated each iteration.
MSTAT_HESSllmf, lmmc, tbeNTT: Read hessian matrix from file, if it exists.
F: assume initial hessian is the unit matrix.
MSTAT_XTOLrlmf, lmmc, tbeY1e-3Convergence criterion for change in atomic displacements.
>0: criterion satisfied when xtol > net shift (shifts summed over all sites).
<0: criterion satisfied when xtol > max shift of any site.
0: Do not use this criterion to check convergence.

Note: When molecular statics are performed, either GTOL or XTOL must be specified. Both may be specified.
MSTAT_GTOLrlmf,lmmc,tbeY0Convergence criterion for tolerance in forces.
>0: criterion satisfied when gtol > “net” force (forces summed over all sites).
<0: criterion satisfied when xtol > max absolute force at any site.
0: Do not use this criterion to check convergence.

Note: When molecular statics are performed, either GTOL or XTOL must be specified. Both may be specified.
MSTAT_STEPrlmf, lmmc, tbeY0.015Initial (and maximum) step length.
MSTAT_NKILLilmf, lmmc, tbeY0 0: Never remove Hessian.
>0: Remove hessian after NKILL iterations.
<0: Remove hessian after -NKILL iterations, and also remove all memory of the hessian in the relaxation algorithm.
MSTAT_LASTRilmfY-1Controls how positions are set in restart file, on final exit from relaxation algorithm
−1: restore to position of minimum gradient
 0: Retain positions of last cycle.
 1: Use algorithm’s estimates for the next cycle (same pos written with --wpos).
MSTAT_PDEF=rlmf, lmmc, tbeY0 0 0 …Lattice deformation modes (not documented).
MD[…] lmmc, tbeY Parameters for molecular dynamics.
MD_MODEilmmcN00: no MD
1: NVE
2: NVT
3: NPT
MD_TSTEPrlmmcY20.671Time step (a.u.)
NB: 1 fs = 20.67098 a.u.
MD_TEMPrlmmcY0.00189999Temperature (a.u.)
NB: 1 deg K = 6.3333e-6 a.u.
MD_TAUPrlmmcY206.71Thermostat relaxation time (a.u.)
MD_TIMErlmmcN20671000Total MD time (a.u.)
MD_TAUBrlmmcY2067.1Barostat relaxation time (a.u.)

The following are specific to lmmag.

SDYN[…] N Subcategory setting parameters for spin dynamics (LL equations) with thermostat
SDYN_KTrN Temperature, in a.u.
SDYN_TSrN Time step, in a.u.
SDYN_TEQUrY0Equilibration time, a.u.
SDYN_TTOTrN Duration of total simulation, in a.u.
MMHAMcN Rules for micromagnetics hamiltonian
GD[…] Y Subcategory to set parameters for thermostat global daemons
GD_NTHERMiY3Number of thermostats
GD_MODETiY31,32,33Thermostat mode(s)
GD_CTrY1thermostat coefficient(s)
BSINT[…] Y Subcategory to set parameters for Bulirsch-Stoer integration
BSINT_NEWlYTStart new SD run, Bulirsch-Stoer integration
BSINT_TOLrN Tolerance in numerical integration
BSINT_TS0rY0Minimum time step in units of TS
BSINT_MXiY7Maximum order of rational function extrapolation
BSINT_MIiY11Maximum number of midpoint rules to invoke
BSINT_NSEQiY2Sequence of number of midpoint divisions

See Table of Contents


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

ASr Y2Controls the relative number of lattice vectors in the real and reciprocal space.
TOLr Y1e-8Tolerance in the Ewald sums. At times you may need to set this to a small value, like 10−12 — when the overlap matrix may have a small eigenvalue. Instances of this are when you use the PMT method or when calculating long superlattices.
NKDMXi Y800The maximum number of real-space lattice vectors entering into the Ewald sum, used for memory allocation.
Normally you should not need this token. Increase NKDMX if you encounter an error message like this one: xlgen: too many vectors, n=…
RPADr Y0Scale rcutoff by RPAD when lattice vectors padded in oblong geometries.

See Table of Contents


This category contains parameters defining the one-particle hamiltonian.

Portions of HAM are read by these codes:


GMAXrlmf, lmfgwdN G-vector cutoff used to create a regular mesh for smoothed Hankel envelope functions and the smooth density (Ry1/2). A uniform mesh with spacing between points in the three directions as homogeneous as possible, with G vectors |G| < GMAX.
This input is required; but you may omit it if you supply information with the FTMESH token.
FTMESHi1 [i2 i3]FPN The number of divisions specifying the uniform mesh density along the three lattice vectors. The second and third arguments default to the value of the first one, if they are not specified.
This input is used only if the parser failed to read the GMAX token.
NSPINiALLY11 for non-spin-polarized calculations.
2 for spin-polarized calculations.
NB: For the magnetic parameters below to be active, use NSPIN=2.
RELiALLY10: for nonrelativistic Schrödinger equation.
1: for scalar relativistic approximation to the Dirac equation.
2: for Dirac equation (ASA only).
11: compute cores with the Dirac equation (lmfa only).
SOiALLY00: no SO coupling.
1: Add L·S to hamiltonian. However, only the spin-diagonal part of the density is retained.
2: Add Lz·Sz only to the hamiltonian, so the spin channels remain distinct.
3: Like 2, but also L·S−LzSz is included perturbatively in the eigenvalues only and in a manner that preserves independence in the spin channels. This generates eigenvalues very close to L·S for a given potential, but the eigenfunctions are generated from H + Lz·Sz only. As a result the eigenfunctions (and then the density) remain spin-diagonal. There is some effect on the density, but the approximation seems to be rather good since the error on the eigenfunctions is of 2nd order in the perturbation.
4: Compute entire L·S perturbatively
11: Same as 1, but additionally decompose SO by site.

GW-based codes at present requires the spin channels to be kept separated and works with SO=2,3 only.
NONCOLlASAYFF: collinear magnetism.
T: non-collinear magnetism.
SS4 rASAY0Magnetic spin spiral, direction vector and angle.

Example: nc/test/ 1
BFIELDilm, lmfY0Tag to add external magnetic field, as explained in this tutorial
0: no external magnetic field applied.
1: add site-dependent constant external Zeeman field (requires NONCOL=T).
Fields are read from file bfield.ext.
2: add Bz·Sz only to hamiltonian.

fp/test/test.fp gdn
nc/test/ 5
BXCSCALilm, lmgfY0This tag provides an alternative means to add an effective external magnetic field in the LDA.
0: no special scaling of the exchange-correlation field.
1: scale the magnetic part of the LDA XC field by a site-dependent factor 1 + sbxci as described below.
2: scale the magnetic part of the LDA XC field by a site-dependent factor as described below.
This is a special mode used to impose constraining fields on rotations, used, e.g. by the CPA code.
Site-dependent scalings sbxci are read from file bxc.ext.
XCFUNcALLY2Specifies local part exchange-correlation functional by name, using names in the libxc library library. Selected two names from this list, prepending names with xc_. Names should be separated by a comma or space, e.g. “xc_GGA_X_PBE xc_GGA_C_PBE”. Be sure to enclose the string in quotes.
XCFUNiALLY2Specifies local part exchange-correlation functional by number.
1: Ceperly-Alder
2: Barth-Hedin (ASW fit)
3: PW91 (same as PBE)
4: PBE (same as PW91)

0,#2,#3: Specify functional in the libxc library by number: use exchange functional #2 and correlation functional #3.
GGAiALLY0Specifies gradient additions to exchange-correlation functional (not used when XCFUN=0,#2,#3).
0. No GGA (LDA only)
1. Langreth-Mehl
2. PW91
3. PBE
4. PBE with Becke exchange

This tutorial uses the PBE functional. To compare the internally coded PBE functional with libxc, try fp/test/test.fp te
PWMODEilmf, lmfgwdY0Controls how APWs are added to the LMTO basis.
1s digit:
0. LMTO basis only
1. Mixed LMTO+PW
2. PW basis only
Examples: fp/test/test.fp srtio3  and  fp/test/test.fp felz 4
10s digit:
0. PW basis fixed to (less accurate, but simpler)
1. PW basis symmetry-consistent, but basis depends on k.
Example:  fp/test/test.fp te
PWEMINrlmf, lmfgwdY0Include APWs with energy E > PWEMIN (Ry)
PWEMAXrlmf, lmfgwdY Include APWs with energy E < PWEMAX (Ry)
NPWPADilmf, lmfgwdY-1If >0, overrides default padding of variable basis dimension. Certain arrays have fixed dimension that must be at least as large as the rank of the hamiltonian. The APW basis is depends on k if PWMODE>10, so some padding must be added to this fixed dimesion to ensure that these arrays can accommodate any k. Normally the code will internally select a sensible default. In the event it is not large enough (the program will stop), you can enlarge the padding with this token.
RDSIGilmf, lmfgwd, lm, lmgfY0Controls how the QSGW self-energy Σ0 substitutes for the LDA exchange correlation functional.
Note: the GW codes store in file sigm.ext.
1s digit:
 0 do not read Σ0
 1 read file sigm.ext, if it exists, and add it to the LDA potential
 2 same as 1 but symmetrize sigm after reading
 Add 4 to retain only real part of real-space sigma
10s digit:
 0 simple interpolation (not recommended).
 1 approximate high energy parts of sigm by diagonal.
Optionally add the following (the same functionality using --rsig on the command line):
10000 to indicate the sigma file was stored in the full BZ (no symmetry operations are assumed).
20000 to use the minimum neighbor table (only one translation vector at the surfaces or edges; cannot be used with symmetrization).
40000 to allow mismatch between expected k-points and file values.
RSSTOLrALLY5e-6Max tolerance in Bloch sum error for real-space Σ0.

Σ0 is read in k-space and is immediately converted to real space by inverse Bloch transform.
The real space form is forward Bloch summed and checked against the original k-space Σ0.
If the difference exceeds RSSTOL the program will abort.
The conversion should be exact to machine precision unless the range of Σ0 is truncated.
You can control the range of real-space Σ0 with RSRNGE below.
RSRNGErALLY5Maximum range of connecting vectors for real-space Σ0 (units of ALAT).
NMTOiASAY0Order of polynomial approximation for NMTO hamiltonian.
KMTOrASAY Corresponding NMTO kinetic energies.
Read NMTO values, or skip if NMTO=0.
EWALDllmYFMake strux by Ewald summation (NMTO only).
VMTZrASAY0Muffin-tin zero defining wave functions.
QASAiASAY3A parameter specifying the definition of ASA moments Q0,Q1,Q2
0. band code accumulates Q1, Q2 from true energy moments of sphere charges (KKR style).
 Sphere code generates density from Q0× + Q2×.
 This (Methfessel convention) is approximate but decouples potential parameters from charges.
1. Sphere code generates density from Q0× + Q2×; thus Q0 is the sphere charge.
2. Q1,Q2 accumulated from and , rather than power moments (not applicable to lmgf, lmpg).
3. 1+2 (Standard conventions).
Add 4 to cause the sphere integrator to construct and by outward radial integration only.
PMINr,r,…ALLY0 0 0 …Global minimum in fractional part of the continuous principal quantum number .
Enter values for l=0,..lmx.
0: no minimum constraint.
# : with #<1, fractional part of .
1: use free-electron value as minimum.

Note: lmf always uses a minimum constraint, the free-electron value (or slightly higher if AUTOBAS_GW is set).
You can set the floor still higher with PMIN=#.
PMAXr,r,…ALLY0 0 0 …Global maximum in fractional part of the continuous principal quantum number . Enter values for l=0,..lmx.
0 : no maximum constraint.
#: with #<1, uppper bound of fractional P is #.
OVEPSrALLY0The overlap is diagonalized and the hilbert space is contracted, discarding the part with eigenvalues of overlap < OVEPS.
Especially useful with the PMT basis, where the combination of smooth Hankel functions and APWs has a tendency to make the basis overcomplete.
OVNCUTiALLY0This tag has a similar objective to OVEPS.
The overlap is diagonalized and the hilbert space is contracted, discarding the part belonging to lowest OVNCUT evals of overlap.
Supersedes OVEPS, if present.
TOLrFPY1e-6Specifies the precision to which the generalized LMTO envelope functions are expanded in a Fourier expansion of G vectors.
FRZWFlFPYFSet to T to freeze the shape of the augmented part of the wave functions. Normally their shape is updated as the potential changes, but with FRZWF=t the potential used to make augmentation wave functions is frozen at what is read from the restart file (or free-atom potential if starting from superposing free atoms).
This is not normally necessary, and freezing wave functions makes the basis slightly less accurate. However, there are slight inconsistencies when these orbitals are allowed to change shape. Notably the calculated forces do not take this shape change into account, and they will be slightly inconsistent with the total energy.
FORCESilmfY0Controls how forces are to be calculated, and how the second-order corrections are to be evaluated. Through the variational principle, the total energy is correct to second order in deviations from self-consistency, but forces are correct only to first order. To obtain forces to second order, it is necessary to know how the density would change with a (virtual) displacement of the core+nucleus, which requires a linear response treatment.
lmf estimates this change using one of ansatz:1.  the free-atom density is subtracted from the total density for nuclei centered at the original position and added back again at the (virtually) displaced position.
The core+nucleus is shifted and screened assuming a Lindhard dielectric response. You also must specify ELIND, below.
NFORCEilmfY-1Controls when to turn on calculation of forces, if they are to be calculated (see FORCES above). This switch can save time by deferring computation of forces until some iterations have been completed.
-1: Always compute forces
 0: Wait until self-consistency criterion has been met; then iterate once more to calculate forces
#>0: Calculate forces starting at iteration  #  or until self-consistency has been reached; then iterate at least once more with forces turned on.
ELINDrlmfY-1A parameter in the Lindhard response function, (the Fermi level for a free-electron gas relative to the bottom of the band). You can specify this energy directly, by using a positive number for the parameter. If you instead use a negative number, the program will choose a default value from the total number of valence electrons and assuming a free-electron gas, scale that default by the absolute value of the number you specify. If you have a simple sp bonded system, the default value is a good choice. If you have d or f electrons, it tends to overestimate the response.
Use something smaller, e.g. ELIND=-0.7.

ELIND is used in three contexts:
(1) in the force correction term; see FORCES= above.
(2) to estimate a self-consistent density from the input and output densities after a band pass.
(3) to estimate a reasonable smooth density from a starting density after atoms are moved in a relaxation step.
SIGP[…] lmf, lmfgwdY Parameters used to interpolate the self-energy . Used in conjunction with the GW package. See gw for description. Default: not used.
SIGP_MODEilmf, lmfgwdY4Specifies the linear function used for matrix elements of at highly-lying energies. High-lying states should be far enough away from the Fermi level that their effect should be small, and the result should depend very little on the choice of the constraint. By approximating for these states, one ensures that the LDA and quasiparticle eigenvectors for those states are the same.
0. constrain to be greater than .
1. constrain to be equal to .
2. constrain to be defined in the interval .
3. constrain as in SIGP_MODE=1. The difference between modes 1 and 3 are merely informational.
4. constrain to be a constant. Its value is calculated by the GW package and read from sigm.ext. This mode requires no information from the user. It is the recommended mode, available in version 7.7 or later.
SIGP_NMAXilmf, lmfgwdY0Integer specifying which of the highest self-energy matrix elements are to be approximated. States higher than SIGP_NMAX have the off-diagonal part of sigma stripped; unlike the low-lying states, the diagonal part of is constrained (see SIGP_MODE above). If SIGP_NMAX is lower or equal to 0, it is not used; see SIGP_EMAX below.
SIGP_EMAXrlmf, lmfgwdY2.0Alternative way to specify approximation of high-lying elements of the self-energy matrix. It is only used if SIGP_NMAX is lower or equal to 0, which case SIGP_EMAX is an energy cutoff: states above SIGP_EMAX are approximated.
SIGP_WEMAXrlmf, lmfgwdY0Smooth transition to energy cutoff in . energy window (EMAX,EMAX+abs(WEMAX)). Use WEMAX>0 for linear interpolation, WEMAX<0 for gaussian interpolation.
SIGP_NMINilmf, lmfgwdY0Integer specifying how many of the lowest-lying states are approximated by discarding the off-diagonal parts in the basis of LDA functions. If SIGP_NMIN is zero, no low-lying states are approximated.
SIGP_EMINrlmf, lmfgwdY0.0Alternative way to specify approximations of low-lying elements of the self-energy matrix. It is only used if SIGP_NMIN<0, which case SIGP_EMIN is an energy cutoff: states below SIGP_EMIN are approximated.
SIGP_Arlmf, lmfgwdY0.02Coefficient in the linear fit (see SIGP_MODE=0,…,3). If SIGP_MODE=4, SIGP_A is not used. In the linear constraints (SIGP_MODE=0,1) it is the constant coefficient; for SIGP_MODE=2, it is the lower bound. Note that its default value is a good estimate for Si.
SIGP_Brlmf, lmfgwdY0.06Coefficient in the linear fit (see SIGP_MODE=0,…,3). If SIGP_MODE=4, SIGP_B is not used. In the linear constraints (SIGP_MODE=0,1) it is the linear coefficient; for SIGP_MODE=2, it is the upper bound. Note that its default value is a good estimate for Si.
SIGP_EFITrlmf, lmfgwdY0Lower bound for the least squares fit required for a reasonable evaluation of the above coefficients SIGP_A and SIGP_B when SIGP_MODE=0,…,3. For SIGP_MODE<3, lmf will make a least-squares fit to for states higher than SIGP_EFIT. For SIGP_MODE=3, lmf will make a least-squares fit for states between SIGP_EFIT and SIGP_EMAX, which must be used if one is going to evaluate for states above some SIGP_EMAX. For the case SIGP_MODE<3 one must invoke lmf one the mesh of k-points for which the self-energy is known (there appear to be fewer problems with interpolation on that mesh). lmf accumulates the minimum, maximum, and least-squares fit for the for all the states above the cutoff. Look in the output for a line beginning with “hambls:”. Also, setting the verbosity above 45, lmf will print out the calculated for each of these states, together with the constrained value. lmf will write to file sigii.ext the data used to make the fit, and summarize the fit and the end of the file. If SIGP_MODE=4, SIGP_EFIT is not needed.
SXiASAY0Calculate screened exchange potential
0. Do nothing
1. Calculate SX self-energy Sigma
2. Calculate SX Sigma with onsite term W.
SXOPTScASAY-Options for SX, e.g. SXOPTS=rssig;nit=3
AUTOBAS[…] lmfa, lmf, lmfgwdY Parameters associated with the automatic determination of the basis set.
These switches greatly simplify the creation of an input file for lmf.
Note: Programs lmfa and lmf both use tokens in the AUTOBAS tag but they mean different things, as described below. This is because lmfa generates the parameters while lmf uses them.

Default: not used.
AUTOBAS_GWilmfaY0Set to 1 to tailor the autogenerated basis set file basp0.ext to a somewhat larger basis, better suited for GW.
AUTOBAS_GWilmfY0Set to 1 to float log derivatives a bit more conservatively — better suited to GW calculations.
AUTOBAS_LMTOilmfaY0lmfa autogenerates a trial basis set, saving the result into basp0.ext.
LMTO is used in an algorithm to determine how large a basis it should construct: the number of orbitals increases as you increase LMTO. This algorithm also depends on which states in the free atom carry charge.
Let lq be the highest which carries charge in the free atom.

There are the following choices for LMTO:
0. standard minimal basis; same as LMTO=3.
1. The hyperminimal basis, which consists of envelope functions corresponding those which carry charge in the free atom, e.g. Ga sp and Mo sd (this basis is only sensible when used in conjunction with APWs).
2. All up to lq+1 if lq<2; otherwise all up to lq.
3. All up to min(lq+1, 3).
For elements lighter than Kr, restrict ≤2.
For elements heavier than Kr, include to 3.
4. (Standard basis) Same as LMTO=3, but restrict ≤2 for elements lighter than Ar.
5. (Large basis) All up to max(lq+1,3) except for H, He, Li, B (use =spd).
Use the MTO token (see below) in combination with this one. MTO controls whether the LMTO basis is 1-κ or 2-κ, meaning whether 1 or 2 envelope functions are allowed per channel.
AUTOBAS_MTOilmfaY0Autogenerate parameters that control which LMTO basis functions are to be included, and their shape.

Tokens RSMH,EH (and possibly RSMH2,EH2) determine the shape of the MTO basis. lmfa will determine a reasonable set of RSMH,EH automatically (and RSMH2,EH2 for a 2-κ basis), fitting to radial wave functions of the free atom.

Note: lmfa can generate parameters and write them to file basp0.ext.
lmf can read parameters from basp.ext.
You must manually create basp.ext, e.g. by copying basp0.ext into basp.ext. You can tailor basp.ext with a text editor. Here are the following choices for MTO:
0: do not autogenerate basis parameters.
1 or 3 : 1-κ parameters with Z-dependent LMX.
2 or 4: 2-κ parameters with Z-dependent LMX. For lmfa 1 and 3 are equivalent, as are 2 and 4.
AUTOBAS_MTOilmf, lmfgwdY0Read parameters RSMH,EH,RSMH2,EH2 that control which LMTO basis functions enter the basis.

Once initial values have been generated you can tune these parameters automatically for the solid, using lmf with the –optbas switch; see here (or for a simple input file guide, here) and here.

The –optbas step is not essential, especially for large basis sets, but it is a way to improve on the basis without increasing the size.

Here are the following choices for MTO:

0 Parameters not read from basp.ext; they are specified in the input file ctrl.ext.
1 or 3: 1-κ parameters may be read from the basis file basp.ext, if they exist.
2 or 4: 2-κ parameters may be read from the basis file basp.ext, if they exist.
1 or 2: Parameters read from ctrl.ext take precedence over basp.ext.
3 or 4: Parameters read from basp.ext take precedence over those read from ctrl.ext.
AUTOBAS_PNUilmfaY0Autoset boundary condition for augmentation part of basis, through specification of the continuous principal quantum number .

0 do not make P
1 Find P for l < SPEC_LMXA from free atom wave function; save in basp0.ext.
AUTOBAS_PNUilmf, lmfgwdY0Autoset boundary condition for augmentation part of basis, through specification of the continuous principal quantum number .

0 do not attempt to read P from basp.ext.
1 Read P from basp.ext, for species which P is supplied.
AUTOBAS_LOCilmfa, lmf, lmfgwdY0Autoset local orbital parameters PZ, which determine which high-lying deep or states are to be included as local orbitals (P.Q.N. differs by ±1).
For deep states, LO can either be conventional or extended (they spill into the interstitial)
PZ can be read either from the ctrl file or the basp file.

Used by lmfa to control whether parameters PZ are to be sought:
0 (or 3): do not autogenerate PZ. Option 3 writes PZ to basp file contents of ctrl file.
1 (or 5) : autogenerate extended (or conventional) PZ. Nonzero values from ctrl file take precedence over basp file
2 (or 6) : same as 1 or 5, but contents of ctrl file are ignored.
Default: 0

Used by lmf and lmfgwd to control how PZ is read:
0 or 3 do not read read parameters PZ from basp file
1-3 read read parameters. 1: Nonzero values from ctrl file take precedence over basp file; otherwise basp takes precedence.
Adding 4 to this parameter has no effect.
Default: 1
AUTOBAS_RSMMXrlmfaY2/3sets an upper bound to LMTO smoothing radius RSMH, when autogenerating a basis set. Value is a multiple of the MT radius.
AUTOBAS_EHMXrlmfaYsets an upper bound to LMTO smoothed Hankel energy EH, when autogenerating a basis set. Default depends on whether AUTOBAS_GW is set.
AUTOBAS_ELOCrlmfaY-2 RyThe first of two criteria to decide which orbitals should be included in the valence as local orbitals. If the energy of the free atom wave function exceeds (is more shallow than) ELOC, the orbital is included as a local orbital.
AUTOBAS_QLOCrlmfaY0.005The second of two criteria to decide which orbitals should be included in the valence as local orbitals.
If the fraction of the free atom wave function’s charge outside the augmentation radius exceeds QLOC, the orbital is included as a local orbital.
AUTOBAS_PFLOATi1 i2lmf, lmfgwdy1 1Governs how the Pnu are set and floated in the course of a self-consistency cycle.
The 1st argument controls default starting values of P and lower bounds to P when it is floated.
0: Use pre-2002 (i.e. version 6) lower bound for P (lmf only).
1: Use defaults and float lower bound designed for LDA.
2: Use defaults and float lower bound designed for GW.
The 2nd argument controls how the band center of gravity (CG) is determined — used when floating P.
0: band CG is found by a traditional method.
1: band CG is found from the true energy moment of the density.

See Table of Contents


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.

MODEiASAY0Tells lmgf what function to perform. See also the Green’s function web page
0: do nothing.
1: self-consistent cycle.
10: Transverse magnetic exchange interactions J(q).
11: Read J(q) from disk and analyze results.
14: Longitudinal exchange interactions.
20: Transverse χ+− from ASA Green’s function.
21: Read χ from disk and analyze results.
20: Transverse χ++, χ−− from ASA Green’s function
Caution: Modes 14 and higher have not been maintained.
GFOPTScASAY ASCII string with switches governing execution of lmgf or lmpg. Use  ’;’  to separate the switches, e.g. GFOPTS=p3;padtol=1e-7 . Switches in GFOPTS are documented on the Green’s function web page.
DLMiALLY0Disordered local moments for CPA.
Governs self-consistency for both chemical CPA and magnetic CPA.
12 : normal CPA/DLM calculation: charge and coherent potential Ω both iterated to self-consistency.
32 : Ω alone is iterated to self-consistency.
BXY1ALLYF(DLM) Setting this switch to T generates a site-dependent constraining field to properly align magnetic moments.

In this context constraining field is applied by scaling the LDA exchange-correlation field.
The scaling factor is [1+bxc(ib)^2]1/2.

A table of bxc is kept for each site in the first column of file shfac.ext.
TEMPrALLY0(DLM) spin temperature.

See Table of Contents


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.

CODEilmfgwdY2This token tells what GW code you are creating input files for.
lmfgwd serves as a driver to several GW codes.
0. First GW version v033a5 (code is no longer maintained) .
2. Current version of GW codes .
1. Driver for the Julich spex code (not fully debugged or maintained).
NKABC1 to 3 ilmfgwdY k-mesh for GW. This token serves the same function for GW as BZ_NKABC does for lmf, and the input format is the same.
Analog to n1n2n3 in GWinput.
MKSIGigw suiteY3(self-consistent calculations only).
Controls the form of (the QSGW approximation to the dynamical self-energy , where refers to a matrix element of Σ between eigenstates n and n′, at energy E relative to EF.
Values of this tag have the following meanings.
0. do not make Σ0
1. Σ0 = Σnn (EF) if nn’, and Σnn(En) if n=n’: mode B, Eq.(11) in Phys. Rev. B76, 165106 (2007)
3. Σ0 = 1/2[Σnn (En) + Σnn (En)]: mode A, Eq.(10) in Phys. Rev. B76, 165106 (2007)
5. “eigenvalue only” self-consistency Σ0 = δnnΣnn‘ (En).
Analog to MKSIG in GWinput.
ECUTPBrlmfgwdY-6.5 RyFor core states with energy >ECUTPB, include this wave in `occ’ column of the core product basis setup
GCUTBrgw suiteY2.7G-vector cutoff for basis envelope functions as used in the GW package (Ry1/2).
Analog to QpGcut_psi in GWinput.
GCUTXrlmfgwdY2.2G-vector cutoff for interstitial part of two-particle objects such as the screened coulomb interaction (Ry1/2)..
Analog to QpGcut_cou in GWinput.
ECUTSrlmfgwdY2.5 Ry(for self-consistent calculations only). Maximum energy for which to calculate the described in MKSIG above.
This energy should be larger than HAM_SIGP_EMAX which is used to interpolate .
When generating a GWinput template, lmfgwd passes ECUTS+1/2 to the emax_sigm tag in the GWinput file.
NBANDigw suiteY6 
NIMEigw suiteY6Number of frequencies on the imaginary integration axis when making the correlation part of Σ.
Analog to niw in GWinput.
DELRErgw suiteY0.01, 0.1Frequency mesh parameters DW and OMG defining the real axis mesh in the calculation of Im .
The ith mesh point is given by: ωi = DW×(i−1) + [DW×(i−1)]2/OMG/2
Points are approximately uniformly spaced, separated by DW, up to frequency OMG, around which point the spacing begins to increase linearly with frequency.
Analog to dw and omg_c in GWinput.
Note: OPTICS category also has this tag, with the same meaning but applied to computing the macroscopic dielectric function.
DELTArgw suiteY-1e-4δ-function broadening for calculating χ0, in atomic units. Tetrahedron integration is used if DELTA<0.
Analog to delta in GWinput.
DELTAWrgw suiteY0.02(one-shot GW only) Width for finite difference in energy differentiation of Σ(ω) for Z factor.
Analog to deltaw in GWinput.
GSMEARrgw suiteY0.003Broadening width for smearing pole in the Green’s function when calculating Σ.
This parameter is sometimes important in metals, e.g. Fe.
Analog to esmr in GWinput. The tag is described in this manual
PBTOLrgw suiteY0.001Overlap criterion for product basis functions inside augmentation spheres. The overlap matrix of the basis of product functions generated and diagonalized for each . Functions with overlaps less than PBTOL are removed from the product basis.
Analog to the second line of the product basis section in GWinput.
QOFFPigw suiteY1Not documented
Q0PSCALErgw suiteY0.1Scaling factor for the offset-Gamma point used to handle the q→0 limit. Analog to Q0Pscale in the GWinput file
Q0PNORMrgw suiteYTSet length of all offset q to average length before symmetry reduction
MIXBETArgw suiteY1Mixing parameter beta for Anderson mixing of the self-energy. Analog to mixbeta in the GWinput file.
NMIXigw suiteY4Number of prior iterations to include in Anderson mixing of the self-energy. Analog to mixpriorit in the GWinput file.
TRSYMlgw suiteYTTurn on/off time-reversal symmetry. Analog to TimeReversal in the GWinput file.
VTFrgw suiteY2e-5Artificial screening added to render the bare coulomb interaction ). Analog to Vkappa in the GWinput file.
Hazard: Too small value can lead to negative eigenvalues in the screened coulomb interaction.
BSE_TDAlgw suiteYTMake the Tamm-Dancoff approximation in the BSE.
BSE_NVrgw suiteYallNumber of occupied states in BSE calculation of W. Analog to NumValStates in the GWinput file.
BSE_NCrgw suiteY4Number of unoccupied states in BSE calculation of W. Analog to NumCondStates in the GWinput file.
BSE_ETAr,…gw suiteY0.02,0.02Smearing for BSE W (Ry), linear interpolation between the two values. No longer used?

See Table of Contents


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 .

NKABC1 to 3 i Y Defines the k-mesh on each of 3 lattice vectors for DMFT driver. If not present, substitute BZ_NKABC.
NLOHI2 i Y-first and last eigenstates to include in projector, relative to EF
WLOHI2 r Y-(used only if NLOHI not found) lower, upper bound to frequency to include in projector, relative to EF
PROJi Y-DMFT projector type
KNORMi Y0How local projectors are normalized:
0: k-independent normalization
1:k-dependent normalization
BROADr N0.0025(for ω on real axis only) additional broadening of sigma, in eV
BETARr --Inverse temperature, in Ry−1
BETAKr --(read only if the preceding tag is missing) Inverse temperature, in K−1
BETAr --(read only if the preceding two tag are missing) Inverse temperature, in eV−1
NOMEGAi Y2000Number of points on frequency mesh

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.

Li N  quantum number defining this correlated subblock
QSPLITi Y2for compatibility with Haule. For now, QSPLIT should always be 2
SITESi, i, … N List of sites with this correlated block.
Note: you can use a negative number for a site index. The minus sign indicates that spins 1 and 2 are to be flipped. In the nonmagnetic case this should have no effect, but for a magnetic site, sites with negative indices are antiferromagnetic (same moment amplituded) to their counterparts with positive index.
SIDXDi, i, … N (diagonal Σ only)   list of (2+1) components of diagonal ΣDMFT(ω) to calculate (see tutorial).
Equal values imply equivalent elements, and 0 value implies matrix element is not calculated. List must contain contiguous numbers.
SIDXMi, i, … - (full Σ, read only if SIDXD is missing)   (2+1)2 components of ΣDMFT(ω) to calculate.
Read in (11, 12, 13, … 21, 22, 23, …) order.
SIDXA  - (full Σ, read only if SIDXM is missing)   populate all (2+1)2 elements of the Σ matrix.
UMODEi Y10specifies approximation for U. 1s, 10s, 100s digit are independent numbers. For now, only 10s digit is implemented.
1s digit:
 0: u(1) = Hubbard U, u(2) = Hubbard J
 1: u(1) = F0, u(2) = F2, u(3) = F4
 2: u(1) = screening length, Yukawa potential
10s digit:
 0: density-density
 1: full matrix U
100s digit:
 0: U is static
 1: U is dynamical
Alternatively, specify by strings separated by ~, one string for each of 1s digit (uj, slater, yukawa), 10s digit (density, full), 100s digit (static, dynamic). Thus UMODE=full~static is equivalent to UMODE=10.

See Table of Contents

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.


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

See Table of Contents


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

SHOW1allYFEcho lines as they are read from input file and parsed by the proprocessor.
Command-line argument --show provides the same functionality.
HELP1allYFShow what input would be sought, without attempting to read data.
Command-line argument --input provides the same functionality.
VERBOS1 to 3allY30Sets the verbosity. 20 is terse, 30 slightly terse, 40 slightly verbose, 50 verbose, and so on. If more than one number is given, later numbers control verbosity in subsections of the code, notably the parts dealing with augmentation spheres.
May also be set from the command-line: --pr#1[,#2]
IACTIV1allYFTurn on interactive mode. Programs will prompt you with queries, in various contexts.
May also be controlled from the command-line: --iactiv or --iactiv=no.
TIM1 or 2allY0, 0Prints out CPU usage of blocks of code in a tree format.
First value sets tree depth. Second value, if present, prints timings on the fly.
May also be controlled from the command-line: --time=#1[,#2]

See Table of Contents


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

It applies to all programs that iterate to self-consistency: lmlmflmmclmgflmpgtbelmfa.

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

NITiallY1Maximum number of iterations in the self-consistency cycle.
MIXcallY A string of mixing rules for mixing input, output density in the self-consistency cycle.
The syntax is given below.
See here for detailed description of the mixing.
CONVrallY1e-5Maximum energy change from the prior iteration for self-consistency to be reached.
See annotated lmf output.
CONVCrallY3e-5Maximum in the RMS difference in the density noutnin. See below.
UMIXrallY1Mixing parameter for density matrix; used with LDA+U
TOLUrallY0Tolerance for density matrix; used with LDA+U
NITUiallY0Maximum number of LDA+U iterations of density matrix
AMIXcASAY Mixing rules when extra degrees of freedom, e.g. Euler angles, are mixed independently. Uses the same syntax as MIX.
NRMIXi1 i2ASA, lmfaY80, 2Used when self-consistency is needed inside an augmentation sphere. This occurs when the density is determined from the momentsQ0,Q1,Q2 in the ASA; or in the free atom code, just Q0.
i1: max number of iterations
i2: number of prior iterations for Anderson mixing2 of the sphere density
Note: You will probably never need to use this token.

See Table of Contents


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
    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
     (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 . WINDOW | r1 r2 | OPTICS | N | 0 1 | Energy (frequency) window over which to calculate Im[ε(ω)].
    Im ε is calculated on a mesh of points .
    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:
    = 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:
    = WINDOW(1) + DW×(i−1) + [DW×(i−1)]2/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 χ2, 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
    2. approximate correction to 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


Portions of OPTIONS are read by these codes:


HF1lm, lmfYFIf T, use the Harris-Foulkes functional only; do not evaluate output density.
SHARM1ASA, lmf, lmfgwdYFIf T, use true spherical harmonics, rather than real harmonics.
FRZlallYF(ASA) If T, freezes core wave functions.
(FP) If T, freezes the potential used to make augmented partial waves, so that the basis set does not change with potential.
SAVVEC1lmYFSave eigenvectors on disk. (This may be enabled automatically in some circumstances)
Q=strncallY  Q=SHOW,  Q=ATOM,  Q=HAM,  Q=POT,  Q=BAND,  Q=DOS,  Q=RHO 
make the program stop at selected points without completing a full iteration.
SCRiASAY0Is connected with the generation or use of the q->0 ASA dielectric response function. It is useful in cases when there is difficulty in making the density self-consistent.
See here for documentation.
0. Do not screen qout−qin.
1. Make the ASA response function P0.
2. Use P0 to screen qout−qin and the change in ves.
3. 1+2 (lmgf only).
4. Screen qout−qin from a model P0.
5. Illegal input.
6. Use P0 to screen the change in ves only.
P0 and U should be updated every iteration, but this is expensive and not worth the cost. However, you can:
Add 10k to recompute intra-site contribution U every kth iteration, 0<k≤9.
Add 100k to recompute P0 every kth iteration (lmgf only).
 Examples: testing/test.scr and gf/test/ mnpt 6
ASA[…]rASAN Parameters associated with ASA-specific input.
ASA_ADNF1ASAYFEnables automatic downfolding of orbitals.
ASA_NSPH1ASAY0Set to 1 to generate l>0 contributions (from neighboring sites) to l=0 electrostatic potential
ASA_TWOCiASAY0Set to 1 to use the two-center approximation ASA hamiltonian
ASA_GAMMAiASAY0Set to 1 to rotate to the (orthogonal) gamma representation. This should have no effect on the eigenvalues for the usual three-center hamiltonian, but converts the two-center hamiltonian from first order to second order.
Set to 2 to rotate to the spin-averaged gamma representation.
The lm code does not allow downfolding with GAMMA≠0.
ASA_CCORllmYTIf F, suppresses the combined correction. By default it is enabled.

Note: NB: if any orbitals are downfolded, CCOR is automatically enabled.
ASA_NEWREPilmYFSet to 1 to rotate structure constants to a user-specified representation.
It requires special compilation to be effective
ASA_NOHYB1lmYFSet to 1 to turn off hybridization
ASA_MTCOR1lmYFSet to T to turn on Ewald MT correction
ASA_QMTrNCY0Override standard background charge for Ewald MT correction
Input only meaningful if MTCOR=T
ASA_FILEMiASAY0Flag to preserve Madelung matrix on disk, to conserve memory (useful only for very large systems).
0: keep in RAM
1: Write to disk and exit
2: never generate; always read from disk
3: Read from disk if available; if not, generate it and write to disk
RMINESrlmchkN1Minimum augmentation radius when finding new empty sites (--getwsr)
RMAXESrlmchkN2Maximum augmentation radius when finding new empty sites (--getwsr)
NESABCi,i,ilmchkN100Number of mesh divisions when searching for empty spheres (--getwsr)
NEPHDilmf,lmfgwdY0Controls the number of energy points for computing the energy derivative of the partial wave by finite difference.
0: lmf uses 2 points, lmfgwd uses 4 points (for historic compatibity)
2: use 2 points for both lmf and lmfgwd
4: use 4 points for both lmf and lmfgwd

See Table of Contents


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

It is read by lmpg and lmstr.

MODEiASAY 0: do nothing.
1: diagonal layer GF.
 Examples: pgf/test/test.pgf -all 5 and pgf/test/test.pgf -all 6
2: left- and right-bulk GF.
3: find k(E) for left bulk.
 Example: pgf/test/test.pgf 2
4: find k(E) for right bulk.
5: Calculate ballistic current.
 Example: pgf/test/test.pgf femgo
SPARSEiASAY00: Calculate G layer by layer using Dyson’s equation
 Example: pgf/test/test.pgf -all 5
1: Calculate G using LU decomposition
 Example: pgf/test/test.pgf -all 6
PLATLrASAN The third lattice vector of left bulk region
PLATRrASAN The third lattice vector of right bulk region
GFOPTScASAY ASCII string with switches governing execution of lmgf or lmpg. Use  ‘;’ to separate the switches.
Available switches:
p1 First order of potential function
p3 Third order of potential function
pz Exact potential function (some problems; not recommended)
Use only one of the above; if none are used, the code makes second order potential functions
idos integrated DOS (by principal layer in the lmpg case)
noidos suppress calculation of integrated DOS
pdos accumulate partial DOS
emom accumulate output moments; use noemom to suppress
noemom suppresss accumulation of output moments
sdmat make site density-matrix
dmat make density-matrix
frzvc do not update potential shift needed to obtain charge neutrality
‘padtol Tolerance in Pade correction to charge. If tolerance exceeded, lmgf will repeat the band pass with an updated Fermi level
omgtol (CPA) tolerance criterion for convergence in coherent potential
omgmix (CPA) linear mixing parameter for iterating convergence in coherent potential
nitmax (CPA) maximum number of iterations to iterate for coherent potential
lotf (CPA)
dz** (CPA)

See Table of Contents


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

FILEcallY Provides a mechanism to read site data from a separate file. File subs/iosite.f documents the syntax of the site file structure.
The reccommended (standard) format has the following syntax:

The first line should contain a ‘%’ in the first column, and a `version’ token vn=#.
Structural data (see category STRUC documentation) may also be included in this line. Each subsequent line supplies input for one site. In the simplest format, a line would have the following:
spid x y z
where spid is the species identifier (same information would otherwise be specified by token ATOM= below) and x y z are the site positions.

Examples: fp/test/test.fp er and fp/test/test.fp tio2

Bug: when you read site data from an alternate file, the reader doesn’t compute the reference energy.

Kotani format (documented here but no longer maintained). In this alternative format the first four lines always specify data read in the STRUC category; see FILE= in STRUC.
Then follow lines, one line for each site
ib iclass spid x y z
The first number is merely a basis index and should increment 1,2,3,4,… in successive lines. The second class index is ignored by these programs. The remaining columns are the species identifier for the site positions.

If SITE_FILE is missing, the following are read from the ctrl file:
ATOMcallN Identifies the species (by label) to which this atom belongs. It is a fatal error for the species not to have been defined.
ATOM_POSr1 r2 r3allN The basis vector (3 elements), in dimensionless Cartesian coordinates. As with the primitive lattice translation vectors, the true vectors (in atomic units) are scaled from these by ALAT in category STRUC.
NB: XPOS and POS are alternative forms of input. One or the other is required.
ATMOM_XPOSr1 r2 r3allN Atom coordinates, as (fractional) multiples of the lattice vectors.
NB: XPOS and POS are alternative forms of input. One or the other is required.
ATOM_DPOSr1 r2 r3allY0 0 0Shift in atom coordinates to POS
ATOM_RELAXi1 i2 i3allY1 1 1Relax site positions (lattice dynamics or molecular statics) or Euler angles (spin dynamics, ASA).
Three numbers correspond to , , Cartesian components. 0 constrains component not to move; 1 allows it to move.
ATOM_RMAXSrFPY Site-dependent radial cutoff for structure constants, in a.u.
ATOM_ROTcASAY Rotation of spin quantization axis at this site
ATOM_PLilmpgY0(lmpg) Assign principal layer number to this site

See Table of Contents


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:

SCLWSRrALLY0SCLWSR>0 turns on the automatic sphere resizer. It defaults to 0, which turns off the resizer.
The 10’s digit tells the resizer how to deal with resizing empty spheres; see this page.
OMAX1r1 r2 r3ALLY0.16, 0.18, 0.2Constrains maximum allowed values of sphere overlaps. This overlap is defined as , where and ae the two sphere radii, and is the bond length. See this page.
You may input up to three numbers, which correspond to atom-atom, atom-empty-sphere, and empty-sphere-empty-sphere overlaps respectively.
OMAX2r1 r2 r3ALLY0.4, 0.45, 0.5Constrains maximum allowed values of sphere overlaps defined as ; see this page.
Both constraints are applied.
WSRMAXrALLY0Imposes an upper limit to any one sphere radius
HFACrlmf, lmfgwdY0Supplies default scale factors that can parameterize the LDA hamiltonian. These values only supply input to file shfac when it does not yet exist. If shfac already exists these parameters are not used. For species-specific factors, edit shfac or use SPEC_ATOM_HFAC[B,L,V].
HFAC_Brlmf, lmfgwdY1Scale BXC by factor B. To invoke this scaling, also turn on HAM_BXCSCAL. Use SPEC_ATOM_HFACB for species-specific scalings.

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

ATOMcallN A character string (8 characters or fewer) that labels this species. This label is used, e.g. by the SITE category to associate a species with an atom at a given site.

The species ID also names a disk file with information about that atom (potential parameters, moments, potential and some sundry other information). More precisely, species are split into classes, the program differentiates class names by appending integers to the species label. The first class associated with the species has the species label; subsequent ones have integers appended.
 Example: testing/test.ovlp 3
ZrallN Nuclear charge. Normally an integer, but Z can be a fractional number. A fractional number implies a virtual crystal approximation to an alloy with some Z intermediate between the two integers sandwiching it.
RrallN The augmentation sphere radius, in atomic units. This is a required input for most programs:
choose one of R=, R/W= or R/A=. Read descriptions of the R/W AND R/A below for further remarks; also see this page for a more complete discussion on the choice of sphere radii.

lmchk can find sphere radii automatically. Invoke lmchk with -–getwsr.
You can also rescale as-given radii to meet constraints with the SCLWSR token.
R/WrallN R/W= ratio of the augmentation sphere radius to the average Wigner Seitz radius W.
W is the radius of a sphere such that (4πW3/3) = V/N, where V/N is the volume per atom.
Thus if all radii are equal with R/W=1, the sum of sphere volumes would fill space, as is usual in the ASA.

You must choose the radii so that the sum of sphere volumes () equals the unit cell volume V; otherwise results may become unreliable. The space-filling requirement means sphere may overlap quite a lot, particularly in open systems. If sphere overlaps get too large (>20% or so), accuracy becomes an issue. In such a case you should add “empty spheres” to fill space. Use lmchk to print out sphere overlaps. lmchk also has an automatic empty spheres finder, which you invoke with the -–findes switch; see here for a discussion.

Example: testing/test.ovlp 3

FP results are much less sensitive to the choice of sphere radii. Strictly, the spheres should not overlap, but because of lmf’s unique augmentation scheme, overlaps of up to 10% cause negligibly small errors as a rule.
(This does not apply to GW calculations!)
Even so, it is not advisable to let the overlaps get too large. As a general rule the -cutoff should increase as the sphere radius increases. Also it has been found in practice that self-consistency is harder to accomplish when spheres overlap significantly.
R/ArallN R/A = ratio of the aumentation sphere radius to the lattice constant
ArallY0.025Radial mesh point spacing parameter. All programs dealing with augmentation spheres represent the density on a shifted logarithmic radial mesh. The ith point on the mesh is . b is determined from the number of radial mesh points specified by NR.
NRiallYDepends on other inputNumber of radial mesh points
LMXiallYNL-1Basis -cutoff inside the sphere. If not specified, it defaults to NL−1
RSMHr,r,…lmf, lmfgwdY0Smoothing radii defining basis (a.u.), one radius for each .
RSMH and EH together define the shape of basis function in lmf.
To optimize, try running lmf with --optbas.
EHr,r,…lmf, lmfgwdY Hankel energies for basis (Ry), one energy for each . RSMH and EH together define the shape of basis function in lmf.
RSMH2r,r,…lmf, lmfgwdY0Basis smoothing radii, second group.
EH2r,r,…lmf, lmfgwdY Basis Hankel function energies, second group.
LMXAiFPYNL - 1Angular momentum -cutoff for projection of wave functions tails centered at other sites in this sphere.
Must be at least the basis -cutoff (specified by LMX=).
IDXDNiASAY1A set of integers, one for each -channel marking which orbitals should be downfolded.
0 use automatic downfolding in this channel.
1 leaves the orbitals in the basis.
2 folds down about the inverse potential function at
3 folds down about the screening constant alpha.
In the FP case, 1 includes the orbital in the basis; >1 removes it
KMXVilmfY15k cutoff for expansion of smoothed Hankel or Gaussian functions in polynomials PkL used in several contexts e.g. overlapping of free atomic densities in the Mattheis construction, and in the charge mixing algorithm.
KMXAilmf, lmfgwdY3k cutoff for expansion of envelope basis functions in polynomials PkL in augmentation spheres. The PkL are explained in Section 3.1 of Questaal’s methods paper, Ref 4. Smoothed Hankels are expanded in polynomials around other sites instead of Bessel functions as in the case of normal Hankels. Polynomials are parameterized by the polynomial order, KMXA, and the smoothing radius RSMA, explained next.
RSMArlmf, lmfgwdYR * 0.4Smoothing radius for projection of smoothed Hankel tails onto augmentation spheres. These functions are expanded in polynomials by integrating with Gaussians of radius RSMA at that site. RSMA very small reduces the polynomial expansion to a Taylor series expansion about the origin. For large KMXA the choice is irrelevant, but RSMA is best chosen that maximizes the convergence of smooth Hankel functions with KMXA.
Note: It has been found from experience that larger values of KMXA need relatively smaller values of RSMA (particularly important when the APW part of the PMT basis has a significant number of plane waves).

A value RSMA<0 is treated as (the negative of) scaling the default value.
LMXLilmf, lmfgwdYNL - 1Angular momentum -cutoff for explicit representation of local charge on a radial mesh.
RSMGrlmf, lmfgwdYR/4Smoothing radius for Gaussians added to sphere densities to correct multipole moments needed for electrostatics. Value should be as large as possible but small enough that the Gaussian doesn’t spill out significantly beyond the Radius of the Muffin-Tin (RMT).
LFOCAiFPY1Prescribes how the core density is treated.
0 confines core to within RMT. Usually the least accurate.
1 treats the core as frozen but lets it spill into the interstitial
2 same as 1, but interstitial contribution to vxc treated perturbatively.
RFOCArFPYR × 0.4Smoothing radius fitting tails of core density. A large radius produces smoother interstitial charge, but less accurate fit.
RSMFArFPYR/2Smoothing radius for tails of free-atom charge density.
Irrelevant except first iteration only (non-self-consistent calculations using Harris functional).
A large radius produces smoother interstitial charge, but somewhat less accurate fit.
RS3rFPY1Minimum allowed smoothing radius for local orbital
HCRrlmY Hard sphere radii for structure constants.
If token is not parsed, attempt to read HCR/R below
HCR/RrlmY0.7Hard sphere radii for structure constants, in units of R
ALPHArASAY Screening parameters for structure constants
DVrASAY0Artificial constant potential shift added to spheres belonging to this species
MIX1ASAYFSet to suppress self-consistency of classes in this species
IDMODiallY00 : floats Pl aka continuous principal quantum number to band center of gravity
1 : freezes
2 : freezes linearization energy .
lmf,lmfgwd only: Add 10 to orthogonalize partial waves to highest core level for a particular .
CSTRMX1allYFSet to T to exclude this species when automatically resizing sphere radii
GRP2iASA, lmfY0Species with a common nonzero value of GRP2 are symmetrized, independent of symmetry operations.
The sign of GRP2 is used as a switch, so species with negative GRP2 are symmetrized but with spins flipped (NSPIN=2)
FRZWF1FPYFSet to T to freeze augmentation wave functions for this species
IDUi[,i…]allY0LDA+U mode. IDU is a vector, with one number for each corresponding to s, p, d, f. A number signifies:
0. No LDA+U
1. LDA+U with Around Mean Field limit double counting
2. LDA+U with Fully Localized Limit double counting
3. LDA+U with mixed double counting.
4. U is treated as a potential, not coulomb interaction. A fixed potential U is added to both spin channels.
5. Same as 4, but U is a potential applied to the majority channel and J to the minority channel. Add 10 to apply U to both phi and phidot (generally recommended).
UHr[,r…]allY0Hubbard U for LDA+U (Ry). UH is a vector, with one number for each .
JHr[,r…]allY0Exchange parameter J for LDA+U (Ry). JH is a vector, with one number for each .
EREFrallY0Reference energy subtracted from total energy
AMASSrFPY Nuclear mass in a.u. (for nuclear dynamics)
C-HOLEclmf, lmY Channel for core hole. You can force partial core occupation.
Syntax consists of two characters, the principal quantum number and the second one of ‘s’, ‘p’, ‘d’, ‘f’ for the quantum number, e.g. ‘2s’
See Partially occupied core holes for description and examples.

Default: nothing
C-HQr[,r]allY-1 0First number specifies the number of electrons to remove from the channel specified by C-HOLE=.
Second (optional) number specifies the hole magnetic moment.
See Partially occupied core holes for description and examples.
Pr,r,…allY Starting values for Pl, aka “continuous principal quantum number”, one for each =0..LMXA

Default: taken from an internal table.
PZr,r,…FPY0Starting values for local orbital’s potential functions, one for each of =0..LMX. Setting PZ=0 for any means that no local orbital is specified for this .
The integer part of nonzero PZ must be either one less than P (semicore state) or one greater (high-lying state).
Add 10 to render the local orbital an extended local orbital (semicore states only)
Add 100 to replace the partial wave with a relativistic analog.
Qr,r,…allY Charges for each -channel making up the free-atom density

Default: taken from an internal table.
MMOMr,r,…ASA, lmfaY0Magnetic moments for each -channel making up the free-atom density
Relevant only for the spin-polarized case.
HFACBrFPY0Species-specific parameter overriding global SPEC_HFAC_B.
HSFITKr,rFPY0,0(JPO basis) radius beyond which kinetic energy should approach asymptotic value.

See Table of Contents


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

RMAXSrallY Radial cutoff for strux, in a.u.
If token is not parsed, attempt to read RMAXA, below
RMAXArallY Radial cutoff for strux, in units of the lattice constant.
If token is not parsed, attempt to read RMAX, below
RMAXrallY0The maximum sphere radius (in units of the average Wigner Seitz radius) over which neighbors will be included in the generation of structure constants.
This takes a default value and is not required input. It is an interesting exercise to see how much the structure constants and eigenvalues change when this radius is increased.
NEIGHBiFPY30Minimum number of neighbors in cluster
ENV_MODEiFP, lm, lmstrY0lmf: 1 turns on streened version of basis
lm: Type of envelope functions:
0 2nd generation
1 SSSW (3rd generation)
3 SSSW and val-lap basis
ENV_NELiFP, lm, lmstrY lmf: Number of screened envelope functions.
In the ASA context (lm): number of NMTOs (SSSW)
ENV_ELrFP, lm, lmstrN0lmf: Number of screened smooth Hankel functions.
In the ASA context (lm): NMTO (SSSW) energies, in a.u.
DELRXrASAY3Range of screened function beyond last site in cluster
TOLGrFPY1e-6Tolerance in =0 gaussians, which determines their range
RVL/RrallY0.7Radial cutoff for val-lap basis (this is experimental)
VLFUNiallY0Functions for val-lap basis (this is experimental)
0 G0 + G1
1 G0 + Hsm
2 G0 + Hsm-dot
MXNBRiASAY0Make lmstr allocate enough memory in dimensioning arrays for MXNBR neighbors in the neighbor table. This is rarely needed.
SHOW1lmstrYFShow strux after generating them
EQUIV1lmstrYFIf true, try to find equivalent neighbor tables, to reduce the computational effort in generating strux.
Not generally recommended
LMAXWilmstrY-1-cutoff for (optional) Watson sphere, used to help localize strux
DELRWrlmstrY0.1Range extending beyond cluster radius for Watson sphere
IINV_NIT=ilmstrY0Number of iterations
IINV_NCUTilmstrY0Number of sites for inner block
IINV_TOLrlmstrY0Tolerance in errors

*IINV parameters govern iterative solutions to screened strux

See Table of Contents


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.

BEGMOMiASAY1When true, causes program lm to begin with moments from which potential parameters are generated. If false, the potential parameters are used and the program proceeds directly to the band calculation.
FREE1ASAYFIs intended to facilitate a self-consistent free-atom calculation. When FREE is true, the program uses rmax=30 for the sphere radius rather than whatever rmax is passed to it; the boundary conditions at rmax are taken to be value=slope=0 (rmax=30 should be large enough that these boundary conditions are sufficiently close to that of a free atom.); subroutine atscpp does not calculate potential parameters or save anything to disk; and lm terminates after all the atoms have been calculated.
CNTROL1ASAYFWhen CONTRL=T, the parser attempts to read the “continuously variable principal quantum numbers” P and moments Q0,Q1,Q2 for each channel; see P,Q below.
ATOMcASAY Class label. P,Q (and possibly other data) is given by class.
Tokens following a class label and preceding the next class label belong to that class.
ATOM_P= and ATOM_QcASAY Read “continuously variable principal quantum numbers” for this class (P=…), or energy moments Q0,Q1,Q2 (Q=…). P consists of one number per channel, Q of three numbers (Q0,Q1,Q2) for each .

Note In spin
polarized calculations, a second set of parameters must follow the first, and the moments should all be half of what they are in non-spin polarized calculations.

In this sample input file for Si, P,Q is given as:
ATOM=SI P=3.5 3.5 3.5 Q=1 0 0 2 0 0 0 0 0
ATOM=ES P=1.5 2.5 3.5 Q=.5 0 0 .5 0 0 0 0 0

One electron is put in the Si s orbital, 2 in the p and none in the d, while 0.5 electrons are put in the s and p channels for the empty sphere. All first and second moments are zero. This rough guess produces a correspondingly rough potential.

You do not have to supply information here for every class; but for classes you do, you must supply all of (P,Q0,Q1,Q2). Data read in START supersedes whatever may have been read from disk.

Remarks below provide further information about how P,Q is read and printed.
RDVES1ASAYFRead Ves(RMT) from the START category along with P,Q
ATOM_ENUrASAY Linearization energies
Sample START category

The following is taken for the distribution’s test of La2 Cu O4.

        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 channel) the corresponding first moments Q1 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


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.

FILEcallY Read structural data (ALAT, NBAS, PLAT) from an independent site file. The file structure is documented here; see also this tutorial. Note: EXPRESS_file performs the same function as STRUC_FILE, and supersedes STRUC_FILE if it is present.
NBASiallN† Number of the sites in the primitive unit cell.
NSPECiallY Number of atom species
ALATrallN† A scaling, in atomic units, of the lattice and basis vectors
DALATrallY0is added to ALAT. It can be useful in contexts certain quantities that depend on ALAT are to be kept fixed (e.g. SPEC_ATOM_R/A) while ALAT varies.
PLATr,r,…allN† (dimensionless) primitive translation vectors
SLATr,r,…lmscellN Superlattice vectors
NLiallY3Sets a global default value for -cutoffs cut = NL−1. NL is used for both basis set and augmentation cutoffs.
SHEARr,r,r,rallY Enables shearing of the lattice in a volume-conserving manner.
If SHEAR=#1,#2,#3,#4,  #1,#2,#3=direction vector;  #4=distortion amplitude.
Example: SHEAR=0,0,1,0.01
distorts a lattice in initially cubic symmetry to tetragonal symmetry, with 0.01 shear.
ROTcallY Rotates the lattice and basis vectors, and the symmetry group operations by a unitary matrix.
Example: ROT=z:pi/4,y:pi/3,z:pi/2 generates a rotation matrix corresponding to the Euler angles α=π/4, β=π/3, γ=π/2. See this document for the general syntax.
Lattice and basis vectors, and point group operations (SYMGRP) are all rotated.
The direction of rotation is such that a +π/2 rotation around z transforms (x,y) into (−y,x).
DEFGRDr,r,…allY A 3×3 matrix defining a general linear transformation of the lattice and basis vectors.
STRAINr,r,…allY A sequence of six numbers defining a general distortion of the lattice and basis vectors.
ALPHArallN Amount of Voigt strain.

†Information may be obtained from a site file

See Table of Contents


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.



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


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


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:


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

where and 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. is the exchange or correlation GGA enhancement factor which introduce semilocal effects and can be choosen using the token GGA. 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 parameterbeta.
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 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  nin and output density  nout to make a trial guess  n*  for the self-consistent density (see for example Chapter 9 in Richard Martin’s book1. Questaal uses two independent techniques to accelerate convergence to the self-consistency condition  noutnin. First, the quantities are mixed making use of a model for the dielectric function. Second, multiple (nin, nout) 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 (Xin,Xout) 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 Broyden3 and Anderson2 methods. The string belonging to ITER_MIX should begin with one of


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                           &larr; iteration 2
   tj:-0.39609  -0.44764                 &larr; iteration 3
   tj:-0.05454   0.01980                 &larr; iteration 4
   tj: 0.24975
   tj: 0.48650

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 Xout and less of Xin. If you see successive iterations where tj is significantly negative (less than −1), increase beta.

Broyden mixing3 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 (XoutXin) 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=#][;...]  or

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 Vin is a functional of the input density nin. For a given Vin, eigenfunctions and nout can be generated. Self-consistency is reached when a particular nin=n* is found for for which n*(in)=n*(out). In Green’s function theory, nout is a functional of Gin, which contains more information than nin, but in any case the self-consistency requires the condition nout=nin. 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 n0 represented on a uniform mesh, true site densities n1, and finally site densities n2 that are a one-center expansions of n0. 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 n1+n2, keeping n1n2 fixed, has little effect on the total potential.

    In brief, each of the three density components n0, n1, n2 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 n0, n1, n2 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 [ dependent radial functions] × [spherical harmonics (the real form of true spherical harmonics Ylm)]. Radial functions are tabulated numerically on a shifted logarithmic mesh. For the site densities n1, n2 the most sensitive part is the =0 radial function as the total charge is a function of it. The =0 component of whatever combination of n1, n2 is mixed is always included in the main density in some form. Depending on locm0 (see below), the >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 n0 is also included in the main density.

    The output densities mixed here are not n0, n1, n2, 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 nout 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

    1. the =0 part of r1’ and r2’, and rsm. The residual is the >0 part of r1’ and r2’
    2. the =0 part of r1’ only, and rsm. The residual is the >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
    3. the =0 part of r2’ only, and rsm. The residual is the >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.
    4. the =0 part of r1’ and r2’, and rsm, and >0 part of r1’ and r2’, approximated as Generalized Gaussian functions or their polynomial parts . The residual is the difference between the numerically tablulate representation and the compacted form of the >0 parts of r1’ and r2’.
    5. 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.
    6. Like mode 3, except all of r1’ is contained in the residual, as in mode 2

    locm1: scales r1’

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

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

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

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

    1. The =0 radial density is used as is.
    2. A Fourier transform is taken, and the Fourier transform is A/B mixed. (For linear mixing locm3=0 and locm3=1 are equivalent)
    3. 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

    1. The mesh density is included in the A/B mixing
    2. 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 paper3 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: .
    elind<0: Use the free-electron gas value, scaled by elind. The default value is −1.

  • 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 scheme2, 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


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

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

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

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