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 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:
68 0.5 0.0 0.0 0.0 0.0 0.0 # X--G 96 0.0 0.0 0.0 0.5 0.5 0.0 # G--M 68 0.5 0.5 0.0 0.5 0.5 0.5 # M--R 118 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.0 # terminates
Each line of the file describes a section of a path in . 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 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 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
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 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 or by band, see this tutorial for instructions.
The opt.sto file contains the energy and Im 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.