ASA tutorial (v7.10)

Purpose

This tutorial demonstrates how to set up and run band calculation in the Atomic Spheres Approximation (ASA), using the lm program. It:

1. describes how to build the input file.

2. explains the output of a self-consistent calculation.

3. shows how to generate energy bands and partial DOS.
This link provides a corresponding tutorial for the full-potential LDA program lmf. Much of the structure is the same, so you are advised to go through this tutorial first.

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.

This completes the creation of the input file.

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 --showp
In 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=2
Many 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:

Once created, lmchk and other programs can read these positions as well:
  • Invoke lmchk to check sphere overlaps.
       lmchk si 
    Verify output that the sum-of-sphere volumes equals the total volume.
    Note that you can see what tokens lmchk looks for by invoking
       lmchk --input 
    and you can see what values lmchk is actually using (including for tokens you didn't specify) by invoking
       lmchk --show 

  • Invoke lmstr to generate real-space structure constants.
       lmstr si 
    Notice 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.

    Setting STR SHOW=t has the effect of displaying some portion of the structure constants; the amount of detail displayed depends on the verbosity.

    By increasing verbosity, lmstr prints out more information. For example, setting the verbosity to 41 or higher causes lmstr to print out a neighbor table. The verbosity in the ctrl file is set to 40. You can change the VERBOS= token in the file, or override it using a command line argument:
       lmstr si --pr41 

  • Invoke lm to generate a self-consistent potential
       lm si 
    The output shown here was doctored slightly to incorporate some hyperlinks. Note the following points:
    1. At the beginning of the output there is a brief epitome of some important parameters and program switches, e.g. ``nonrel'' and ``no-ccor''.

    2. No SYMGRP category was supplied, so the program found its own symmetry group operations.

    3. The BZ mesh generated 10 k-points from the specified number of divisions (4x4x4) and the BZJOB=1 switch. The latter causes the k-mesh to be offset from the Γ point (mesh does not contain k=0)

    4. Note the potential parameters PPAR generated in the sphere branch. Parameter C is the band center of gravity; you see that it is (in this first iteration) rather deep for the Si s, about -0.7 Ry, near zero for the Si p and about 1.5 Ry for the Si d. This is as expected: the bottom of the valence bands are Si s-like; the valence bands are Si p-like, and the conduction bands are Si d like. It is interesting to note how much these parameters have changed when self-consistency is reached; see the corresponding data for the last iteration.
      Parameter E&nu is the linearization energy, which is connected to the choice of P. P, and consequently E&nu, change quite a lot in the iterations towards self-consistency. Because SPEC IDMOD=0, P was allowed to float so that enu falls at the center-of-gravity of the occupied parts of the bands. This corresponds to making the first energy moment Q1=0.

    5. Next comes a table of average potentials at each MT boundary. This data is averaged to construct the ASA muffin-tin zero, VMTZ.

    6. In the band pass some energy bands are printed out. The amount of data printed depends on the verbosity. Because BZJOB=1, the no k-point lies at Γ where the valence-band maximum lies.

    7. Next comes the determination of the Fermi level. Because Si is an insulator, this branch is particularly simple. For a metal, you will get something like the following (depending on what integration method and verbosity you select)
       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.006528
      
      The 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.
      Also printed is the sum of occupied bands, which is needed together with the potential to compute the total energy. The Bloechl correction is of interest because it shows how much the energy changes from straight tetrahedron (converges as h2) and the Bloechl-corrected tetrahedron (which converges faster than h2)

    8. Next comes the printout of the Harris-Foulkes total energy. delsev should be small when near self-consistency. (Note: in this special case the Harris-Foulkes total energy is not generated in the first iteration. This is an artifact of the initial condition where the moments and therefore empty sphere potential is exactly zero, which deceives lm into thinking that double-counting terms were never generated. It is a buglet not worth pursuing, since it is not connected to any physical condition.)

    9. Immediately following the output moments Q0,Q1,Q2 are shifted so that the  P's  can float to the band centers of gravity.

    10. The output moments are mixed with the starting moments to form a new set of moments, which with luck will closer to the self-consistent numbers; also the RMS deviation between the two is printed.

    11. A new pass is made through the sphere branch, creating a new potential and set of potential parameters.

    12. Some information summarizing the total energy and how close this pass was to achieving the tolerances for self-consistency is printed. (See note above about 1st iteration Harris-Foulkes energy.) This ends the first iteration. The cycle is repeated until self-consistency is reach, or the maximum number of iterations is encountered.

    13. Invoke lmctl to extract the self-consistent P,Q. (This step is optional). Edit the file  log.si  and paste the results to the bottom of the  ctrl  file.
      *Caution  take care with what P,Q the program uses. Since the  START  category is set up as  BEGMOM=T CTRL=T, it will always read any  P,Q  it finds in  START, and make a new potential from those values. Since our specification of  P,Q  have been commented out, we don't have any worries; the program will use whatever it finds on the disk.

    3. Generating energy bands and densities-of-states

    1. Energy bands are typically plotted along specified symmetry lines, which information you supply through a symmetry-line file. Rather than create one, we will just copy one from the startup directory to syml.si.
        cp startup/syml.fcc ./syml.si
                
      syml.si has rows such as:
         41  .5 .5 .5     0  0 0                L to Gamma   (Lambda)
         ...
          0    0 0 0  0 0 0
                
      Each 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=syml 
      The 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 si 
      which creates directly a postscript file ps.si, or better:
        echo -15 15 5 10 | plbnds -fplot -scl=13.6 -ef=0 si 
      which creates an fplot command, in file plot.plnds. Use fplot to create a postscript file:
        fplot -f plot.plbnds 
      which 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.

    2. To generate the total density-of-states (DOS), you don't have to do anything special. If BZ_SAVDOS=T, the total DOS is automatically written to disk every band pass. We will generate the partial DOS here. It is particularly simple in the ASA; in the FP code you need to set some extra switches.

      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 --pr30 
      You 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 --pr30 
      You 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 --pr30 
      This 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).
      *Aside  if you add the command-line argument --dos:totdos, then lmdos would combine the partial DOS to make only the total DOS.
      Either way, the DOS are output in file  dos.si.

      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.si 
      creates 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.dos 
      creates 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 

    Use of color to highlight orbital character in energy bands or DOS  

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

    To see how the orbitals are ordered, run the energy band program (lm or lmf) with rather high verbosity (>51) and look for the tables following Makidx.
    Here is a sample output for GaAs run using lmf. This case contains two sites, and two additional sites with floating orbitals.
     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.

    To see this feature illustrated, copy doc/ASAsamples/ctrl.gas to your current working directory. Verify that the potential is self-consistent:
      lmstr gas
      lm gas
    
    Let's select pick out the orbitals of As p character. To see where there are in the hamiltonian, do
      lm gas --pr51 --quit=ham
    
    These 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.gas 
    Then generate the bands with
      lm gas --band~col=6,7,8~fn=syml
    
    If 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.plbnds
    
    The 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.

    Another useful application is to distinguish between majority and minority bands in the spin-orbit coupled case. Suppose the basis consists of 80 orbitals. With SO coupling, the basis is doubled to 160 orbitals. Using option ~col=1:80 in the --band switch assigns whatever weight the first spins contribute to the eigenvalue. The color of the energy bands as generated by fplot will be related to its spin character.

    One final point: the LMTO basis is not orthogonal, so the inner product defining the norm is a little different from the orthogonal case. It is:
      Σi (zni)-1 zin
    
    This 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.

    Pairs of color weights: examples using lmf

    Both lm and lmf will generate two sets of color weights.

    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 1
    
    This 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=syml
    
    Orbitals 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=9
    
    correspond 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=band
    
    yields 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 channels
    
    Invoke
        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.