Building input files for lmf

This tutorial shows how to build an input file (ctrl.ext) for the main DFT code, lmf. The easiest way is to use the input file maker blm. Its primary input is the crystal structure (lattice and basis vectors), and it builds a (nearly) complete input file from the structural information.

This tutorial serves in part as an accompaniment to the blm documentation, explaining its application through examples. It demonstrates different ways to feed structural information to blm:

Preliminaries

This tutorial assumes you have blm in your path. Additionally, poscar2init, lmchk, lmscell, poscar2site and cif2site may be required for some sections. Also, it uses some files from the Questaal repository. This tutorial assumes ~/lm is the top-level directory where you keep it.

blm is an input file generator. Usually it reads from init.ext to build the input file (ctrl.ext). You can create the ctrl file from scratch without blm, but most people find going through blm to be easier and faster. You can customize the file once it’s been made.

blm writes the input file to actrl.ext to avoid overwriting ctrl.ext you might already have. (You can specify otherwise: blm has many command line switches; --ctrl lets you name the file).

Usually, but depending on the switches you use, blm splits off the structural information and writes it to a separate file, site.ext. The Questaal package can read information from either the ctrl file or a site file, or some combination, but site files are very convenient because you can change the structural information without modifying the ctrl file.

blm writes the input file to actrl.ext to avoid overwriting ctrl.ext you might already have. (You can specify otherwise: blm has many command line switches; --ctrl lets you name the file).

The init file structure is documented in the blm documentation, and explained here through examples.

1. Importing Crystal Structure

1.1 Supplying Crystal Structure by hand

Classic books like Wyckoff’s Structure of Crystals provide information about crystal structure in a compact form.

Bi2Te3 (and also Bi2Se3) has the space group R3M. R3M is a rhombohedral lattice with hexagonal axes, which has two independent dimensions: the length a=4.3835Å of the basal plane and c=30.487Å of the third (c) axis normal to the basal plane.

Bi2Te3 contains 3 Te and 2 Bi atoms in the primitive unit cell; the two Bi and two of the Te and are equivalent by symmetry. Thus coordinates need be specified for only three atoms. Wykcoff writes them as (0,0,0) and (0,0,0.788) for two inequivalent Te, and (0,0,0.4) for one of the Bi. These coordinates are expressed in terms fractional multiples of the lattice vectors in the conventional cell.

This completely specifies the lattice structure. Enter this information into file init.ext as follows.

SITE_C denotes site positions expressed as fractional multiples of the conventional cell. See below for alternatives

We can supply equivalent information through the lattice vectors. The Bravais lattice R3M tells us about the point group operations as well as the lattice vectors. There are 12 rotations in all, which can be generated from the generators I R3Z MY. These operations tell blm the positions of symmetry-equivalent Te and Bi atoms. To get all five atoms you must either give blm the generators of the group, or specify the positions of all five atoms in the unit cell. The init file below supplies the generators.

Typically when the lattice is designated by a space group as Wyckoff does, the accompanying positions are expressed in terms of the conventional cell (C=). In this mode you must supply the lengths (e.g. A) and angles (e.g. ALPHA) of the lattice vectors that are dictated by symmetry.

When instead the lattice is specified directly through lattice vectors PLAT, site positions are usually expressed in as multiples of the primitive cell (X=), as was done in the second init file, or in Cartesian coordinates (POS=). In this mode you must also supply the lattice constant ALAT (not the same as A, as you can see by comparing the two init files), and possibly the generators of the point group, as discussed above.

Lines which begin with ‘#’ are comment lines and are ignored. (More generally, text following a `#’ in any line is ignored).

Lines beginning with ‘%’ are directives to the preprocessor. Directives can perform various functions similar to a normal programming language, such as assigning variables, evaluating expressions, conditionally readings some lines, and repeated loops over sections of input.

Quantities in brackets {…} are algebraic expressions parsed by the preprocessor, evaluated and the result the contents of brackets replaced by the result.

SPCGRP=R-3M tells blm that you are specifying the lattice by its space group. R3M is entered as R-3M. You can also refer to this group by its number: SPCGRP=166.

A={a} C={c} are the two lengths of this hexagonal lattice. After parsing by the preprocessor, this line gets transformed into A=4.3835 C=30.487.

UNITS=A tells blm that lengths A and C are in Å units. By default blm uses atomic units.

The SITE category operates in a tree structured syntax in much the same form as the ctrl file.

Symbols Te and Bi tell blm that the atoms are Tellurium and Bismuth, with atomic numbers 52 and 83. You can use any symbol for the species name, but if you don’t use a standard one you must specify the atomic number, e.g. write ATOM=A Z=52.

1.2 Importing a CIF file

Crystallographic Information Files (CIF files for short) is a standard text file format for representing crystallographic information, whose standards are set by the International Union of Crystallography. If you have a CIF file, you can automatically make either an init.ext file, or a site.ext file. The tools in this package do not read CIF files directly, but parse the output of the cif2cell program (version 1.1.0). cif2cell is a very versatile tool, freely available on the web.

To import data from a CIF file you need cif2cell installed and cif2init in your path (cif2init should be automatically compiled with this package). The steps are:

1. Run cif2cell (without any special switches) and capture the output in a file, e.g. cif2cell.out . 2. Run cif2init to generate an init file (called simply ` init ‘). 3. Rename init to init.ext and use the blm tool.

Example : create an init file for the orthorhombic form of BaTiO3:

cp testing/cif2cell.batio3 .
cif2init cif2cell.batio3

Note:cif2cell.batio3 was obtained by running cif2cell on the CIF file supplied with the cif2cell-1.1.0 distribution: cif2cell-1.1.0/cifs/BaTiO3_orthorhombic.cif.

If you prefer to make your own input file but import structure information via a site file, use cif2site in place of cif2init. To summarize:

cif2init creates an init file from the output of cif2cellcif2site creates a site file from the output of cif2cell

1.3 Importing a POSCAR file

VASP is a very popular electronic structure program which can store structural information in the POSCAR files. If you want to import data from such a file in a form suitable for this program suite, use one of the following commands:

poscar2init creates an init file from a VASP POSCAR fileposcar2site creates a site file from a VASP POSCAR file

If you already have an input file, use poscar2site to translate structural data from a POSCAR file into site. If you want an input file instead, use poscar2init; then run blm.

Example: create ctrl file for Zn3As2 from a POSCAR file:

~/lm/testing/test.blm 2

The script has the following steps

poscar2init > init.zn3as2
blm zn3as2 --fixpos:tol=1e-6 > out.zn3as2

which creates an input file template (actrl.zn3as2) from a POSCAR file in two steps.

Example: create a site file for the Kesterite Cu2ZnSnS4:

testing/test.blm 3

It creates file site from file POSCAR.

2. Running blm

blm will create a template file for the main input file of the Questaal package, ctrl.ext. ext is a name you select; we will use bi2te3 corresponding to the material. blm actually generates actrl.ext, prepending the a so as not to overwrite any file named ctrl.ext.

Almost all programs in this package require the input file ctrl.ext. You can do an entire calculation starting only with this file; but often you supply other files. For example, symmetry line points for plotting energy bands are read from a separate file (e.g. syml.bi2te3). Structural data is typically split off into a separate file (site.bi2te3 in this tutorial); blm will create a site file by default.

2.1 init file

Copy either the first or second form of init files from one of the boxes above into init.bi2te3.

As explained in the blm documentation, init files and ctrl files are styled the same way: data is divided into categories (LATTICE and SITE) — tags which begin in the first column, and tokens (e.g. SPCGRP=, A=, X=) — tags belonging to a particular category.

2.2 ctrl file

To create a template input file invoke blm:

blm bi2te3

Now take a look at the standard output. blm finds the lattice vectors from the crystal symmetry (if you had supplied the lattice vectors, it would have found the symmetry). In this case, the symmetry is given; the corresponding point group has 12 symmetry operations, and two extra atoms had to be added to make the basis compatible with the group. The output snippet below indicates what blm found.

blm automatically determined augmentation sphere radii, which it accomplishes by attempting to find spheres with equal potentials on each sphere surfaces (as well as it can). If you already have an input file, you can run lmchk with –getwsr to determine radii for you (it uses the same algorithm as blm). Particularly in polar compounds, this algorithm probably does a better job than you can do by hand, and it is recommended that you use the radii it finds, or some scaled version of them.

You must specify a valid k mesh for Brillouin zone integration; see here.

You must define a basis set which can be done manually or automatically, as described next.

blm creates a family of tags belonging to AUTOBAS to enbable other programs to automatically find a basis set for you. We will use this tag, which sets up a standard minimal basis:

AUTOBAS[PNU=1 LOC=1 LMTO=3 MTO=3 GW=0]

(Note that you must modify actrl.bi2te3 a little; the default gives [.. LMTO=4 MTO=4], which makes a double kappa basis.

lmfa calculates wave functions and atomic densities for free atoms. It also has a mode that automatically generates information for basis sets, using tokens in AUTOBAS to guide it. This information is written to a file basp0.ext. AUTOBAS specifies set of conditions that enable lmfa to automatically build the basis set for you, including specification of envelope function parameters RSMH and EH. Alternatively you can define parameters such as EH and RSMH basis by hand, as described in this tutorial.

Tokens in AUTOBAS tell lmfa to do the following:

PNU=1 Calculate logarithmic derivative parameter for the free atom; save P parameters in file basp0._ext_ .
Nothing about P is written if PNU=0 .
LOC=1 Look for shallow cores to be explicitly treated as valence electrons, using [local orbitals](/docs/code/fpoverview/#local-orbitals).
Shallow cores that meet specific criteria are identified and written to basp0._ext_ as PZ=.
No search is made if LOC=0 .
LMTO=3 Pick a default choice about the size of basis. LMTO=3 is a standard minimal basis.
Run lmfa --input and look for HAM_AUTOBAS_LMTO to see what other choices there are.
Note that **lmfa** will pick some defaults for the _l_-cutoff
e.g. _spd_ or _spdf_ depending on the value of LMTO.
MTO=1 Choose 1-kappa basis set (single orbital per _l_ channel).
For better quality calculations, it is recommended you use MTO=2.
GW=0 If GW=1, tailor basis for a GW calculation, e.g. changing the
criteria for including shallow cores as valence, and the size of basis.

These tokens thus define some reasonable default basis for you. lmfawritesbasp0.ext. This file is never read, but lmf will readbasp.ext and use this information when assembling the basis set. The two files have the same structure and you can copy basp0.ext to basp.ext. What lmfa generates is not cast in concrete. You are free to adjust the parameters to your liking, e.g. add a local orbital or remove one from the basis.

The AUTOBAS tokens tell lmf what to read from basp.ext. It uses tokens in a manner similar, but not identical, to the way lmfa uses them:

PNU=1 Read parameters P for all species present in basp._ext_
LOC=1 LOC=1 or 2 tell **lmf** to read local orbital parameters PZ.
Since these parameters may also be specified by the input file,
LOC=1 tells lmf to give precedence to parameters specified by ctrl file
LOC=2 tells lmf to give precedence to parameters specified by basp.
LMTO= is not used by **lmf**.
MTO=1 RSMH and EH may also be specified by the input file
LMTO=1 or 3 tells lmf to read 1-kappa parameters specified by basp
LMTO=2 or 4 tells lmf to read 2-kappa parameters specified by basp
LMTO=1 or 2 tells lmf that parameters in the ctrl file take precedence
LMTO=2 or 4 tells lmf that parameters in the basp file take precedence
GW=0 If GW=1, tune basis for a GW calculation: log derivative parameters P
are floated a little differently in the self-consistency cycle.
They are weighted to better represent unoccupied states, at a slight cost
to their representation of occupied states.

3. Checking sphere overlaps

Sphere overlaps can be checked using lmchk. To do this copy the template actrl.bi2te3 to the input file and run lmchk:

cp actrl.bi2te3 ctrl.bi2te3
lmchk bi2te3

By default, blm makes the spheres as large as possible without overlapping, as the output shows. In this case the Bi and Te radii are nearly the same.

4. Making the atomic density

Make the free atom density. If you did not do so already copy actrl.bi2te3 to the input file (changing [.. LMTO=4 MTO=4] to [.. LMTO=3 MTO=3]) and invoke lmfa:

cp actrl.bi2te3 ctrl.bi2te3
lmfa bi2te3

The primary purpose of lmfa is to generate a free atom density. A secondary purpose is to supply additional information about the basis set in an automatic way. All of this information can be supplied manually in the input file, but the autogenerated input file supplies a minimum amount of information. lmfa generates basp0.bi2te3 which contains

Every species gets one line. This file specifies a basis set consisting of spdf orbitals on Te sites, and spdf orbitals on Bi sites, and a local 5d orbital on Bi. The contents of this file are explained above; see also RSMH and EH and PZ.

Note: Remember that lmf reads from basp.ext, not basp0.ext.

3.1 Automatically finding deep states to include as valence electrons

The partitioning between valence and core states is something that requires a judgement call. lmfa has made a default choice for you: the output shows that for Bi, lmfa selected the 6s, 6p, 6d, 5f states, populating them with charges 2, 3, 0, 0. Note that the total sphere charge is Q=0. You can override the default, e.g. choose the 5d over the 6d with SPEC_ATOM_P; override the l channel charges with SPEC_ATOM_Q.

As was explained earlier, when HAM_AUTOBAS_QLOC is set lmfa will look for shallow core levels below 6s, 6p, 6d, 5f states, and as this table shows lmfa selected the 5d orbital which is normally a core state, to be included as a local orbital so that the usual 6d state and the 5d state are simultaneously included in the basis. Even though the 5d state is fairly deep (the output shows it lies at −2 Ry), the criterion of having a charge density outside the smoothing radius greater than 3×10−3 was met. (Use HAM_AUTOBAS_ELOC and HAM_AUTOBAS_QLOC to change these criteria.) lmfa supplies information about this to basp0.bi2te3, in the form PZ=0 0 15.936 (no local orbitals for s or p states). The 0.936 is significant: it tells lmf what boundary condition to use for the 5d radial function.

3.2 Automatically finding linearization energies

Because HAM_AUTOBAS_P is set, lmfa save estimates for logarithmic derivative parameters P into basp0.ext. As is well known from elementary quantum mechanics there is a relation between the energy of a wave function and its logaritmic derivative at some radius. This information is supplied through the parameters P.

lmfa calculates P for the free-atom potential. Since this potential is not so far removed from the crystal potential, these parameters can find the “band center” reasonably well for each partial wave l. In any case, these are only estimates; they normally get “floated” in the self-consistency cycle.

3.3 Automatically finding envelope function parameters

Finally, the input file contains AUTOBAS[MTO=1]. This causes lmfa to envelope function parameters RSMH and EH (RSMH is the most important of the two) that fit the free atom wave functions well, and save the result in basp0.bi2te3. Unfortunately, what is optimum for the free atom is not optimum for the crystal, but the parameters are a reasonable starting point. These parameters are important, as they determine the quality of the basis. Later we discuss ways to optimize them, or improve the basis quality by adding APWs. blm cannot automatically determine every required input from the structural data. But for the following reasons: lmf will not run properly as the situation now stands

The basis set has to be specified. lmfa autogenerated parameters for a basis set, and saved the result to basp0.bi2te3. A decision must be made whether to use lmfa’s choices for RSMH, EH, P, and PZ, to supply your own, or to modify lmfa’s choice. In any case, you can do it in one of two ways:

The atm file was created by lmfa without prior knowledge that the 5d local orbital is to be included as a valence state (via a local orbital). Thus it incorrectly partitioned the core and valence charge. You must do one of the following:

Remove PZ=0 0 15.936 from basp.bi2te3. It will no longer be treated as a valence state. Removing it means the remaining envelope functions are much smoother, which allows you to get away with a coarser mesh density, as described below. Whether you need it or not depends on the context.

Copy basp0.bi2te3 to basp.bi2te3 and run lmfa over again:

cp basp0.bi2te3 basp.bi2te3
lmfa bi2te3

With the latter choice lmfa operates a little differently from before as can be seen by comparing new output with the old. Initially the Bi 5d was part of the core; now is included as part of the valence.

blm does not by default assign any value to the plane wave cutoff for the interstitial density. lmf reads this information through HAM_GMAX. It is a required input; but blm does not pick a value because its proper choice depends on the smoothness of the basis. lmfa will determine a suggested value for HAM_GMAX for you. In the present instance, when the usual 6s, 6p, 6d, 5f states are included lmfa recommends GMAX=4.4 as can be seen by inspecting the first lmfa run. In the second run it recommends GMAX=4.1 from the valence states alone (as before), but because of the 5d state lmfa recommends that GMAX=8.1. The 5d state is strongly peaked at around the atom, and requires more plane waves to represent reasonably, even a smoothed version of it, than the other states. The difference between 8.1 and 4.4 is substantial, and it reflects the additional computational cost of including deep core-like states in the valence. This is the all-electron analog of the “hardness” of the pseudopotential in pseudopotential schemes. If you want high-accuracy calculations (especially in the GW context), you will need to include these states as valence. This particular choice of local orbital is rather overkill for LDA calculations however. If you eliminate the Bi 5d local orbital you can set GMAX=4.4 and significantly speed up the execution time.

blm assigns the initial k-point mesh to zero. Note the following lines in actrl.bi2te3:

BZ_NKABC governs the mesh of k-points. An appropriate choice will depend strongly on the context of the calculation and the sytem of interest; the density-of-states at the Fermi level; whether Fermi surface properties are important; whether you want optical properties as well as total energies well described; the precision you need; the integration method, and so on. Any automatic formula can be dangerous, so blm will not choose a default for you. In this case, a 4×4×4 mesh works well. Use your text editor to change nk=0 to nk=4. Alternatively, supply –nk=.. to blm on the command line, as was done in this tutorial.

Note that as generated, ctrl.bi2te3 will reflect METAL=5. Using METAL=5 with the tetrahedron integration is the recommended way to handle Fermi surface integration in metals. See this tutorial for some discussion.

Input files for GW

GW calculations demand more of the basis set because unuoccupied states are important. To set up a job in preparation for a GW calculation, invoke blm as :

blm --gw bi2te3

Compare actrl.bi2te3 generated with the –gw switch to one without. One important difference will be that the default basis parameters are modified because AUTOBAS becomes:

AUTOBAS[PNU=1 LOC=1 LMTO=5 MTO=4 GW=1]

The basis is similar to LMTO=4 but EH has been set a little deeper. This helps the QS_GW_ implementation interpolate between k-points. The larger basis makes a minor difference to the valence bands; but the conduction bands change, especially the higher in energy you go.