This tutorial demonstrates how to set up and run band calculation in the Atomic Spheres Approximation (ASA), using the lm program. It:
1. Building the Input file
This tutorial explains an a sample calculation for Si, using this input file as a template. You may choose to create some or all of your own input file from scratch; or you can just copy doc/ASAsamples/ctrl.si to your current working directory. The top-level directory is assumed to be "~/lm." It is assumed that executables lmchk, lm, lmstr, etc are in your path.
Before going through these steps, read how the input file is structured, including the categories and tokens that organize it. You should also read the ASA documentation.
CONST a0=.5292Category CONST is special. It does not supply any tokens for input; its only purpose is to declare symbolic variables for use elsewhere.
You supply information defining the lattice structure in the STRUC category through ALAT and PLAT; supply site position data in the SITE category. Each site must be associated with a chemical element, or species; you supply information about the "personality" of species in the SPEC category.
Si has one species and two atoms/cell, but because the ASA makes shape approximations to the potential, its use becomes problematic when dealing with open structures. To compensate two artificial "empty" sites (atomic number 0) are added to the basis. Thus for purposes of this calculation, the basis has four atoms: two Si atoms and two others just to fill the interstitial. Specify the number of sites in STRUC_NBAS.
Note that each species label (following ATOM) in the SITE category must match a label in the SPEC category. This is how the code determines the chemical character of element at each site.
These programs will splits species into symmetry-equivalent classes. In the ASA codes (lm, lmgf, lmpg) labels are assigned to each class. Class labels are derived from species labels by appending numbers to the species labels. They cannot exceed 8 characters, so usually you should limit the species labels to 6 or so characters (fewer if you have many atoms).
The atomic number Z= is required for each species.
You must choose a sphere radius, choosing one of R=, R/W= or R/A=. The choice of sphere radii is the Achilles heel in augmented wave methods. This is particularly true in the ASA because one must choose a balance between the following competing needs:
It is difficult for the user to choose radii that satisfy (1), but a switch is available in program lmchk to determine these radii automatically, as described HEREXX. (The automatic input file generator, program blm uses this same algorithm to supply sphere radii.)
(2) and (4) are diametrically opposed. In the ASA it is essential the the sum-of-sphere volumes equals the cell volume (i.e. the interstitial volume be zero) which means empty spheres are needed for all but close-packed structures. In the FP case one still prefers to make the interstitial as small as possible because the augmented wave basis is better than the envelope (interstitial) basis.
The sample input file makes a simple choice: all the spheres have the same radii, namely the average Wigner-Seitz radius W; radii are chosen as R/W=1.
OPTIONS ASA[ CCOR=F ADNF=F]Brackets [..] are ncessary because OPTIONS_ASA_CCOR and OPTIONS_ASA_ADNF are tokens nested to the 3rd level. CCOR tells the ASA not to include the ``combined correction'' term (not generally recommended). ADNF=T would turn on a feature to automatically search for orbitals to downfold (remove from the basis). All these tags (the OPTIONS category itself, and OPTIONS_ASA_CCOR and OPTIONS_ASA_ADNF optional; see see below.
ITER also holds the content of the MIX token, MIX=A,b=.8, a string. This string is passed charge mixer, which mixes the input density nin (some trial density which generates the potential), with the output density nout generated by the potential. The aim is to find the self-consistent density, where nout= nin. Some kind of mixing is usually necessary to stabilize convergence. 'A' tells the mixer to use the Anderson mixing scheme; it mixes 80% of output charge with 20% of input charge because b=.8. The ASA density is completely specified by the logarithmic derivatives at the sphere boundary, via parameters P, and the multiple moments Q0,Q1,Q2. Theu the parameters mixed in the ASA are (P,Q): four numbers for each l channel and class. In the FP implementation, the full charge density is mixed. Other programs, e.g. self-consistent empirical tight-binding, mix still other quantities. But the mixing commands are mostly common to all programs.
There are many other choices in this string. Unfortunately difficult to make a general specification about the best mixing procedure to converge P and Q to self-consistency (in the full-potential case, it mixing the sphere and interstitial densities) If you have problems, The first thing to do is down the mixing parameter `b'. You can also use Broyden mixing (MIX=B) rather than Anderson mixing; indeed this is recommended. For difficult cases, see linear-response-asa.html.
2. Running programs in the ASA package
Before carrying out a self-consistent calculation, run program lmchk for some preliminary checks. This program's main function is to checks sphere packing and overlaps, but it has a number of other features.
Use --showp to see the operation of the preprocessor:
lmchk si --showp (or equivalently) lmchk ctrl.si --showpIn this case what is printed on the screen is nearly identical to the input file, except commented lines are not printed. This file makes no use of the preprocessor so no lines are modified.
Next, use --show=2 to see what data is read. This switch tells the program to print out what the parser reads, or infers from defaults, and exits without further processing.
lmchk ctrl.si --show=2Many tokens are printed out; note in particular NSPIN token was read. There should be a line HAM_NSPIN in the output:
Token Input cast (size,min,read,def) result ... HAM_NSPIN opt i4 1, 1, 1, 0 1
Program lmchk checks sphere packing, and has the following options:
lmchk --shell siwhich produces this output for the sample ctrl.si.
lmchk --mino~z --wpos=newpos file-ext ...
lmchk --rpos=newpos file-ext ...
lmchk siVerify output that the sum-of-sphere volumes equals the total volume.
lmchk --inputand you can see what values lmchk is actually using (including for tokens you didn't specify) by invoking
lmstr siNotice in the output that the each site has 27 neighbors, which corresponds to three shells. The number of neighbors depends on the choice of RMAX=, which is 3.2 for this example. If we hadn't specified anything, the program defaults to RMAX= is 2.7 which for bcc packing here will generate 15, or two shells of neighbors. This isn't too bad, for total energies but is a little crude. Few neighbors correspond to a low Fourier series cutoff in reciprocal space, and thus don't pick up rapid oscillations in the bands. More neighbors are also needed for reasonable description of higher-lying bands.
lmstr si --pr41
lm siThe output shown here was doctored slightly to incorporate some hyperlinks. Note the following points:
BZWTS : --- Tetrahedron Integration --- Est E_f Window Tolerance n(E_f) -0.148571 -0.150937 -0.148426 0.002511 14.958994 -0.148573 -0.148577 -0.148552 0.000025 14.819478 -0.148573 -0.148573 -0.148573 0.000000 14.820559 BZINTS: Fermi energy: -0.148573; 41.000000 electrons Sum occ. bands: -24.813178, incl. Bloechl correction: -0.006528The tetrahedron method is generally the most accurate way to perform Brillouin zone integrations. This code incorporates Blochl corrections so that the integration error scales as something like h3 (h = k-point spacing) rather than h2 of the standard tetrahedron integration method. But band crossings may occur which messes up the integration. There is a telltale sign of a band crossing, namely the integrated number of electrons doesn't come out to an exact integer. There is alternative generalized generalized sampling integration method developed by Methfessel and Paxton. For M/P integration you must choose N and W by trial and error---but you get a good idea by looking at the bands at the Fermi surface. For steep bands use W about 0.05 - 0.1 and for flat bands, use W about 0.01 to 0.05 and then adjust N between 1 and 4 or 5. See this page for a detailed comparison of the two methods.
3. Generating energy bands and densities-of-states
cp startup/syml.fcc ./syml.sisyml.si has rows such as:
41 .5 .5 .5 0 0 0 L to Gamma (Lambda) ... 0 0 0 0 0 0 0Each line of text specifies a line segment in k-space. In the first line of this file the segment runs from k=(1/2,1/2,1/2) to 0 in units of 2π/a (L to Γ). '41' specifies the number of k points are generated in this segment. At present the remainder of the line is not used. The last line tells the line reader not to read any more lines.
To generate the ASA energy bands, do:
lm si --band:fn=symlThe bands are generated and saved in file bands.si.
If you have the plbnds program installed from the FPLOT package, you can plot the bands by invoking one of the following
echo -15 15 5 10 | plbnds -scl=13.6 -ef=0 siwhich creates directly a postscript file ps.si, or better:
echo -15 15 5 10 | plbnds -fplot -scl=13.6 -ef=0 siwhich creates an fplot command, in file plot.plnds. Use fplot to create a postscript file:
fplot -f plot.plbndswhich creates a postscript file ps.dat.
The energy bands generator has an optional feature that enables you to highlight a particular orbital character in the energy bands. See color weights below.
When making DOS, you should run lm once more with many k-points, and also use tetrahedron integration. (Note: check that BZ_SAVDOS=T). You can reset the switches in ctrl.si tetrahedron or sampling (it doesn't matter how you made it self-consistent). For DOS plotting, W must be about 1/10th the distance between peaks (N = 1) or larger for N > 1. The defaults are N=0, W=.005; N=0 always gives conventional Gaussian smearing with Gaussian width, W.
Instead of altering the ctrl file, we take advantage of the ability to alter the value of variables from command-line arguments. Since CONST doesn't change the value of variables already declared, command-line declarations take precedence over the CONST entries. Also, because we don't want the program to change the potential, we use interactive mode, and make it stop after the bands have been completed. Finally, we turn down the verbosity, since there will be so many k-points
lm si -vnk=16 -vmet=1 --iactiv --pr30You will encounter lines beginning with QUERY: Press return to each query until you reach this line:
QUERY: beta (def=0.8)?which is asking you if you want to change the mixing parameter. Just type q and the program will stop.
*Aside it is interesting to note that RMS DQ=3e-4, which is the error we made because we only used 4 divisions of k-points to make the potential self-consistent. This is a relatively small error; but in metals the errors would be significantly larger: 4 divisions is not usually sufficient in such a case.
The ASA program has the ability to use files generated as a byproduct of the band pass to immediately generate DOS resolved by orbital and class. This is done using lmdos. Invoke
lmdos si -vnk=16 -vmet=1 --iactiv --pr30You will be queried with a prompt:
Enter npts (def=501), emin and emax (def=-1,0):enter 1001 -1 .3 RETURN, or do the whole thing in one go:
echo 1001 -1 .3 / | lmdos si -vnk=16 -vmet=1 --iactiv --pr30This generates partial partial DOS in 6 channels: three for the s,p,d channels on the Si and three for the s,p,d channels on the "empty" site (effectively the interstitial).
If you have installed the FPLOT package, you can use pldos to create pictures of DOS from this file. For example,
echo 8 7 / | pldos -fplot '-lst=1,2,3;4,5,6' dos.sicreates a file dosp.dat with two columns containing partial dos (one column combines dos 1,2,3 ---the Si dos and the other combines dos 4,5,6---the Es dos), and a file plot.pldos which holds the fplot command that creates a postscript file for this DOS.
fplot -disp -pr10 -f plot.doscreates and displays postscript file ps.dat. You can see from the DOS that there is a direct gap near the Fermi level of about 0.04 Ry, which is the LDA gap for Si. Note: there is a facility to draw two DOS in a single panel, one above the 'zero' and one below. This is particulary convenient in spin polarized cases when you want to compare the majority and minority DOS. Example: in file dos.dat the majority DOS are in channels 1,3,5 (atom 1) and 7,9,11 (atom 2), and the minority DOS are in channels 2,4,6 (atom 1) and 8,10,12 (atom 2), invoke, e.g.
echo 8 7 / | pldos -fplot '-lst=1,3,5;7,9,11' '-lst2=2,4,6;8,10,12' dos.dat
Starting with v6.16, the energy bands maker will allow you to generate
weights together with the bands themselves. With this option, each
energy levels is assigned a corresponding weight which can be used by
a graphics package to highlight orbital character associated with
them. The fplot graphics package, for example, can
read these weights and to draw bands with continuously varying color
fixed by the weights.
The orbital character is defined through the Mulliken weights, defined as follows. Consider a basis of orthonormal orbitals, with orbital component i of the nth eigenvector, zin. The inner product
Σi (zin)† zinis unity. This is just a statement that the wave function is normalized, or that the eigenstate carries unit charge. The unit norm can be resolved into individual orbital contributions (Mulliken decomposition), which gives insight into the orbital character of a particular state. By supplying the orbital-list above, the weight assigned to each state is the fraction of the total norm comes from the specified list of orbitals. This weight is written into the energy bands file.
Makidx: basis arranged in downfolding order: ib low intermed high .. offH .. 1 spdf (16) (0) g (9) 0 0 0 k2 sd (6) (0) pfg (19) 16 0 9 k3 d (5) (0) (0) 22 0 28 2 spdf (16) (0) g (9) 27 0 28 k2 p (3) (0) sdfg (22) 43 0 37 k3 s (1) (0) (0) 46 0 59 3 spd (9) (0) fg (16) 47 0 59 k2 (0) (0) spdfg (25) 56 0 75 4 spd (9) (0) fg (16) 56 0 100 k2 (0) (0) spdfg (25) 65 0 116 Makidx: hamiltonian dimensions Low, Int, High, Negl: 65 0 141 94 kappa Low Int High L+I L+I+H Neglected 1 50 0 50 50 100 0 2 9 0 91 9 100 0 3 6 0 0 6 6 94 all 65 0 141 65 206 94 Orbital positions in hamiltonian, resolved by l: Site Spec Total By l ... 1 Ga 1:27 1:1(s) 2:4(p) 5:9(d) 10:16(f) 17:17(s) 18:22(d) 23:27(d) 2 As 28:47 28:28(s) 29:31(p) 32:36(d) 37:43(f) 44:46(p) 47:47(s) 3 EA1 48:56 48:48(s) 49:51(p) 52:56(d) 4 EC1 57:65 57:57(s) 58:60(p) 61:65(d)This particular hamiltonian has three kinds of orbitals per l channel ("κ's"), so there are potentially three groups of s orbitals, three groups of p orbitals, three groups of d orbitals, and so on. They are labelled '1', 'k2', 'k3' in the top table, one set for each site. ('k3' orbitals are local orbitals in lmf.) The left column under offH tabulates the number of orbitals preceding a particular block. The 'low' block contains all basis functions in the actual hamiltonian.
The Orbital Positions table offers the most convenient way to find orbitals for color weights for color weights. The As atom, for example, occupies columns 28..47 in the hamiltonian. To assign a color to orbitals associated with As, invoke lmf with the --band switch, modified by the list of orbitals you want to highlight, i.e. --band~col=orbital-list~.... For example:
lmf --band~col=28:47~syml ...will read information about what lines to plot from file syml.ext, and assign a color weight to orbitals 28..47. For the general syntax of `orbital-list', see Syntax of Integer Lists. Supposing you wanted to assign a color to just the As p orbitals. (The valence band maximum is pure p character; the conduction band minimum is pure s character.) As either of the above tables show, there are two sets of As p orbitals (no local orbitals in this example). Invoke lmf with:
lmf --band~col=29:31,44:46~syml ...lmf (or lm) will place into the bands band file bnds.ext both the eigenvalues and a a corresponding set of weights which can be used by a graphics package to highlight band features connected with them. For example, the fplot plotting package will use these weights to color the energy bands according to the weight. It is a very useful way to pick out a particular orbital character in the energy bands.
lmstr gas lm gasLet's select pick out the orbitals of As p character. To see where there are in the hamiltonian, do
lm gas --pr51 --quit=hamThese tables should appear in the output:
Makidx: basis arranged in downfolding order: ib low intermed high .. offH .. 1 sp (4) d (5) + 0 0 0 2 sp (4) d (5) + 4 5 0 3 sp (4) d (5) + 8 10 0 4 sp (4) d (5) + 12 15 0 Makidx: hamiltonian dimensions Low, Int, High, Negl: 16 20 0 0 kappa Low Int High L+I L+I+H Neglected - 16 20 0 36 36 0 Orbital positions in hamiltonian, resolved by l: Site Spec Total By l ... 1 GA 1:4 1:1(s) 2:4(p) 2 AS 5:8 5:5(s) 6:8(p) 3 E1 9:12 9:9(s) 10:12(p) 4 E2 13:16 13:13(s) 14:16(p)The second atom is the As atom. Orbitals 6,7,8 are the As p orbitals. Create a symmetry-lines file or copy one:
cp startup/syml.fcc ./syml.gasThen generate the bands with
lm gas --band~col=6,7,8~fn=symlIf you have the FPLOT graphics installed, you can draw a picture with
echo -14,10,5,10 | plbnds -fplot -ef=0 -scl=13.6 -lt=1,col=1,.2,.1,colw=.1,.2,1 -lbl=L,G,X,W,G bnds.gas fplot -f plot.plbndsThe blue lines are valence bands, The two heavy hole bands px and py character between 0 and -3 eV. The valence band are apparently carried mostly by As p states. The light hole band of pz character is p-like only at Γ ; most of this valence band (running between 0 and -7 eV) is purple, reflecting the hybridization of the As p state with other states.
Σi (zni)-1 zinThis changes nothing except that it need not be true that each individual orbital contributes a positive weight to the norm. The sum of all contributions is still unity.
For a ready-made example using lmf with spin orbit coupling, run the Co test case in the standard distribution. If lm is the top-level directory where the codes are installed, run
~/lm/fp/test/test.fp co 1This test has two sets of spd orbitals in the basis. The --band switch reads:
--band~col=5:9,dup=9~col2=18+5:18+9,dup=9~fn=symlOrbitals 5-9 are the Co d orbitals of the first set (corresponding to EH...). String 'dup=#' is an extension of the standard integer-list syntax. It means: take the list generate so far, and replicate each element adding # to it. In this example orbitals 14:18 are the d bands of the EH2 channel. The basis consists of 18 orbitals, which is doubled to 36 with spin-orbit coupling. Thus
col=5:9,dup=9 and col2=18+5:18+9,dup=9correspond respectively to the two pairs of d orbitals of the first spin and of the second spin. Thus the first color weight singles out all spin-up states of d character, the second spin-down d states.
A figure with showing bands with colors corresponding to Mulliken weights,
can be found at the top of the energy bands documentation
. Instructions for drawing the figure using fplot are also there.
The following outlines a more complex example. Suppose we are studying Ba3ZnTa2O9, and we want to know which bands are predominantly Ta character and which are predominantly Zn character. Suppose
lmf --pr51 --quit=bandyields this table:
Orbital positions in hamiltonian, resolved by l: Site Spec Total By l ... 1 Ba2 1:19 1:1(s) 2:4(p) 5:9(d) 10:16(f) 17:19(p) 2 Ba 20:38 20:20(s) 21:23(p) 24:28(d) 29:35(f) 36:38(p) 3 Ba 39:57 39:39(s) 40:42(p) 43:47(d) 48:54(f) 55:57(p) 4 Ta 58:78 58:58(s) 59:61(p) 62:66(d) 67:67(s) 68:70(p) 71:75(d) 76:78(p) 5 Ta 79:99 79:79(s) 80:82(p) 83:87(d) 88:88(s) 89:91(p) 92:96(d) 97:99(p) 6 Zn 100:114 100:100(s) 101:103(p) 104:108(d) 109:109(s) 110:114(d) 7 O1 115:127 115:115(s) 116:118(p) 119:123(d) 124:124(s) 125:127(p) 8 O1 128:140 128:128(s) 129:131(p) 132:136(d) 137:137(s) 138:140(p) 9 O1 141:153 141:141(s) 142:144(p) 145:149(d) 150:150(s) 151:153(p) 10 O2 154:166 154:154(s) 155:157(p) 158:162(d) 163:163(s) 164:166(p) 11 O2 167:179 167:167(s) 168:170(p) 171:175(d) 176:176(s) 177:179(p) 12 O2 180:192 180:180(s) 181:183(p) 184:188(d) 189:189(s) 190:192(p) 13 O2 193:205 193:193(s) 194:196(p) 197:201(d) 202:202(s) 203:205(p) 14 O2 206:218 206:206(s) 207:209(p) 210:214(d) 215:215(s) 216:218(p) 15 O2 219:231 219:219(s) 220:222(p) 223:227(d) 228:228(s) 229:231(p) suham : 375 augmentation channels, 375 local potential channelsInvoke
lmf --band~col=58:99~col2=100:114~syml ...to assign the first color weights to the two Ta atoms, and the second color weight to Zn.