Plotting bands using lmf and plbnds

Plotting bands using lmf is very easy. Once a calculation has been converged, lmf is restarted using the “–band” option, which calculates the eigenvalues at $k$ along a predefined path. Questaal comes with several tools for simplifying the plotting bands and producing simple figures, but the bnds.extn file has a simple format and is easily plotted otherwise.

We will look at the cubic perovskite strontium titanate and repeat the steps to setup the calculation starting from a cif file:

cp /home/vol05/tmp15/sto_icsd_80872.cif .
cif2cell sto_icsd_80872.cif > out
cif2init out
mv init init.sto
blm sto
cp actrl.sto ctrl.sto

Edit the ctrl file and increase the number of iterations to ~20, and maybe choose the libxc’s PBE implementation (XCFUN=0 101 130).

lmfa --usebasp sto | tee log_lmfa

The option –usebasp causes lmfa to write to the basp file directly. lmfa suggests a GMAX of 10.0 a.u., add this to the ctrl file also. A k-mesh of at least 8^3 is reasonable (feel free to check): either edit the ctrl file or invoke lmf with the option -vnkabc=8.

lmf reports in the first output lines that it is indeed using the PBE-GGA:

 pot:
 libxc     GGA X (Perdew, Burke & Ernzerhof) + GGA C (Perdew, Burke & Ernzerhof)

Converging the calculation should take 10 seconds or so.

Simple band plot

Create a file named qp.sto containing the following:

0.5 0.0 0.0  0.0 0.0 0.0 # X--G
0.0 0.0 0.0  0.5 0.5 0.0 # G--M
0.5 0.5 0.0  0.5 0.5 0.5 # M--R
0.5 0.5 0.5  0.0 0.0 0.0 # R--G
 0.0 0.0 0.0  0.0 0.0 0.0 # terminates

Each line of the file describes a section of a path in $k$ . The first column is the number of points in each segment followed by Cartesian coordinates of the start and end points. lmchk is capable of setting this up automatically, including proportionate scaling of the divisions.

mpirun -np 24 lmf --band sto | tee log_lmf-bands

Produces bnds.sto, an ASCII tabulation of the eigenvalues for $k$ along the path. –band takes various options, eg, –band~fn=syml causes syml.sto to be read instead of qp.sto.

Plotting using the plbnds utility

The plbnds tool reads the bnds.sto file and prepares instructions for the fplot plotting program (part of Questaal) that can generate postscript figures more or less automatically.

echo -8  10   8   6 |plbnds -fplot -ef=0 -scl=13.6057 -lbl=X,G,M,R,G bnds.sto
      a   b   c   d              e     f    g            h

The meaning of the different options are:

a: ymin, here in eV, of the figure to be produced
b: ymax
c: width in cm of the figure
d: height (these could alternatively be provided interactively)
e: prepare output for fplot
f: shift $E_f$ to zero in the plot
g: scale energies to eV
h: labels for the plot – these should correspond to the coordinates in qp.sto

Finally, fplot can be executed to give a postscript file:

fplot -f plot.plbnds
mv fplot.ps basic_bands_sto.ps

Which should produce a completely acceptable band structure plot: copy it to the local machine to have a look.

fplot is quite an advanced package for general plotting tasks – fplot.
Documentation for plbnds is here : plbnds.

Colouring band plots according to Mulliken decomposition

The bnd file can also be written with one or two weights for each eigenvalue, corresponding to the weight of different basis functions in the eigenfunction. In order to use this, we need to find the position in the Hamiltonian of the basis functions we are interested in – this is printed by lmf when the output verbosity is higher than normal:

lmf --quit=ham --pr60 sto | tee log_hampos

yields:

 Orbital positions in hamiltonian, resolved by l:
 Site  Spec  Total    By l ...
   1   Sr    1:28   1:1(s)   2:4(p)   5:9(d)   10:16(f) 17:17(s) 18:20(p) 21:25(d) 26:28(p)
   2   Ti   29:49   29:29(s) 30:32(p) 33:37(d) 38:38(s) 39:41(p) 42:46(d) 47:49(p)
   3   O    50:62   50:50(s) 51:53(p) 54:58(d) 59:59(s) 60:62(p)
   4   O    63:75   63:63(s) 64:66(p) 67:71(d) 72:72(s) 73:75(p)
   5   O    76:88   76:76(s) 77:79(p) 80:84(d) 85:85(s) 86:88(p)
 suham :  89 augmentation channels, 89 local potential channels  Maximum lmxa=4

Selecting the titanium d basis functions (the list syntax can be quite complex see here):

mpirun -np 24 lmf --band~col=33:37,42:46 sto| tee log_lmf-bands-col

A second set can be included as ~col2, if desired. Then plotting proceeds:

echo -8 10 8 6 | plbnds -fplot -ef=0 -scl=13.6 -lt=1,col=0,0,0,colw=1,0,0 -lbl=X,G,M,R,G bnds.sto

which is the same as before except that the default colour and the weighted colour have been defined (black and red, in rgb).

fplot -f plot.plbnds
mv fplot.ps ti-d-states.ps

Optical properties using lmf

lmf allows the calculation of dielectric function and joint densities of states. These are controlled by the OPTICS ctrl category, “lmf –input” reports:

 OPTICS_MODE            opt    i4       1,  1     default = 0
   0: do nothing;
   1: generate linear eps_2
   2: generate linear eps_2, spin 2
      Add 20 for Second Harmonic Generation
   8: eps_2 for nonequilibrium absorption.  KT, IMREF are required inputs.
   9: eps_2 for photoluminescence.  KT, IMREF are required inputs.
  -1: generate joint DOS
  -2: generate joint DOS, spin 2 (LTET=3)
  -3: generate up-down joint DOS (LTET=3)
  -4: generate down-up joint DOS (LTET=3)
  -5: generate spin-up single DOS (LTET=3)
  -6: generate spin-dn single DOS (LTET=3)
  -8: Joint DOS for nonequilibrium absorption.  KT, IMREF are required inputs.
  -9: Joint DOS for photoluminescence.  KT, IMREF are required inputs.

For example, Im $\epsilon$ evaluated between 0 and 1 Ryd, in 1001 bins:

OPTICS  MODE=1 NPTS=1001 WINDOW=0 1

lmf again will now print the following after the k-integration step:

 OPTINQ:  Im eps(w) using occ bands (1,88) and unocc (1,88)
          1001 points in energy window (0,1)  Efermi=-0.076969
 tetfbz  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
          Reduce bands to range (9,31)
 TETWTQ:  total number of (occ,unocc) pairs = 66560, 130 or fewer per qp
 OPTINT:  writing opt file ...

The range of bands included in the calculation can be reduced using the options OPTICS_FILBND and OPTICS_EMPBND. There are ways to decompose the output by $k$ or by band, see this tutorial for instructions.

The opt.sto file contains the energy and Im $\epsilon$ for polarisations along x,y,z (which are identical in this case). Copy the file to the local machine and plot using gnuplot, xmgrace or python.

Resolving the dielectric function precisely generally requires quite many k-points. Because of the LDA underestimate of the gap, this plot looks rather different to the experimental one.