Partial Densities of States and Core-Level Spectroscopy

This tutorial explains how to generate total and partial densities-of-states (DOS) with the Questaal package, using the band codes lmf, lm, and tbe. Total DOS is simplest and can be generated in a normal band cycle without postprocessing; it is demonstrated for Co below. Resolving DOS into components can provide much useful information and it can be decomposed in multiple ways. This tutorial demonstrates two common forms in the Questaal suite: projection of the eigenfunction onto partial waves in augmentation spheres, which is demonstrated for Cr₃Si₆; and Mulliken analysis. Partial wave and Mulliken decomposition are compared in the Fe portion of this tutorial.

Core-level spectroscopy [1] describes the excitation of an core electron, calculated by Fermi’s golden rule which involves the matrix element of the dipole operator with the core and valence wave functions. It is closely related to the partial DOS and is computed in much the same way, as described in the Fe tutorial below. The CLS should be calculated in the presence of a core hole, where the density redistributes. This is done in the CrN tutorial on this page.

Other forms DOS that Questaal can calculate, but not described in this tutorial are:

1. The Green’s function packages lmgf and lmpg compute partial DOS.
2. The DOS is the independent-particle approximation to the spectral function. The spectral function for the interacting case can be calculated in the GW framework, and also the DMFT framework. The noninteracting and interacting forms are compared at the GW level in this tutorial.
3. The optics modes in the lmf and lm codes also enable you to resolve DOS in other forms.

Table of Contents

• Preliminaries
• Introduction
  • Partial Waves
  • Mulliken Analysis
  • How information is assembled for analysis
  • Notes on Partial wave and Mulliken analysis, and their relative merits
• Making total DOS using band programs lmf, lm, or tbe
• 1. Total DOS in elemental Co
  • Building a Co input file
  • Self-consistent Co density
  • Total DOS in Co
• 2. Partial DOS in Cr3Si6
  • Input file and self-consistency in Cr3Si6
  • Partial DOS in Cr3Si6 resolved by l
• 3. DOS and Core-level Spectroscopy in Fe
  • Input file and Self-consistent density in Fe
  • Mulliken analysis in Fe
  • Partial DOS in Fe
  • Core-level spectroscopy in Fe
• 4. Core-level spectroscopy in the presence of a core hole in CrN
  • CrN input file
  • Self-consistent density in CrN
  • Core-hole spectroscopy in CrN
• Additional exercises
• References

Preliminaries

This tutorial assumes you have cloned and built the Questaal repository (located here). Questaal executables such lmf, lmdos, pldos, and catdos are required assumed to be in your path.

Note: most jobs in this tutorial fast enough that MPI is not needed. However lmf can be run with MPI, e.g. invoke lmf as mpirun -n 8 lmf.

Introduction

The density-of-states $D(E)$ is given by a sum over states i as [2] $$D(E) = \sum\nolimits_i\delta(E{-}E_i) \qquad (1)$$

Band methods lmf, lm, and tbe work in a different manner than the Green’s function methods, lmgf and lmpg. They can evaluate Eq. (1) directly by approximating the δ-function with a Gaussian function. This method (sometimes called gaussian sampling) is simple and safe but is slow to converge with k. Convergence can be greatly accelerated with Methfessel and Paxton’s polynomial generalization of the Gaussian, but it is more cumbersome than the tetrahedron method, which is also implemented in the band programs. We use the tetrahedron method here.

This tutorial lmf to generate DOS, but lm and tbe perform similar functions.

Programs lm, lmf, and tbe have a facility to resolve, or decompose, the eigenfunction of a particular eigenfunction $\psi_i(E_i,\mathbf{r})$ into component parts. Note that an eigenstate is normalized: $\smallint|\psi_i(E_i,\mathbf{r})|^2d^3\mathbf{r}{=}1$. Decomposition amounts to resolving the unit norm of the wave function in different ways. A myriad of ways are possible [3]: Questaal offers two kinds, “partial waves” and “Mulliken analysis.” Core-level spectroscopy (rather closely related to the partial wave analysis), is also explained here.

Partial Waves

The eigenfunction inside an augmentation sphere is given by solutions to the radial wave equation $\phi_l(E,r)$. The full energy-dependence of $\phi_l(E,r)$ is approximated by to linear order, expanded to first order in a Taylor series in energy. Thus, inside an augmentation sphere there are two partial waves for a particular site R and angular momentum $\ell$ that contribute to the DOS: $\phi_{\mathbf{R}l}(\mathbf{r})$ and the energy derivative $\dot\phi_{\mathbf{R}l}(\mathbf{r})$. The ($\phi_{\mathbf{R}l}$,$\dot\phi_{\mathbf{R}l}$) pair are assumed to completely span the hilbert space inside augmentation sphere R (unless there is an additional wave from a local orbital). See this page for the definition of the lmf basis set.

Denoting the $\ell$ and m quantum numbers by a compound index $\ell$, and labeling $\phi_{\mathbf{R}l}$ and $\dot\phi_{\mathbf{R}l}$ respectively as $\phi_{0\mathbf{R}l}$ and $\phi_{1\mathbf{R}l}$, the eigenfunction can be projected onto an augmentation sphere centered at R as $$P_\mathbf{R}\, \psi_i(E_i,\mathbf{r}) = \sum\limits_{\alpha L} {A_{\alpha\mathbf{R}L}^{(i)}} \, \phi_{\alpha\mathbf{R}l}(r) Y_L(\mathbf{\hat r}) \qquad (2)$$

$P_\mathbf{R}$ denotes a projection onto augmentation sphere R, $\alpha$ ranges from 0 to 1 (and if local orbitals are present, encompasses them), $L$ up to the augmentation cutoff lmxa.

Coefficients ${A_{\alpha\mathbf{R}L}^{(i)}}$ (which is determined from a solution of the secular matrix) then represent a particular kind of decomposition of $\psi _i(E,\mathbf{r})$. Assuming the ($\phi$, $\dot\phi$) basis is complete, this decomposition is independent of basis set. However, it does depend on the augmentation radius. In sum Eq. (2) can be expressed in terms of the energy-dependent partial wave as $$P_\mathbf{R}\, \psi _i(E,\mathbf{r}) = \sum\limits_{\alpha L} {C_{\mathbf{R}L}^{(i)}} \, \phi_{\mathbf{R}l}(E,r) Y_L(\mathbf{\hat r}) \qquad (3)$$

where $\phi_{\mathbf{R}l}(E,r)$ is a linear combination of $\phi$, $\dot\phi$ (and possibly local orbitals) normalized so that $\smallint |\phi_i(E,\mathbf{r})Y_L(\mathbf{\hat r})|^2d^3\mathbf{r}{=}1$.

The ${C_{\mathbf{R}L}^{(i)}}$ make up partial contributions to $D(E)$, Eq. (1). The contribution to $D(E)$ from a particular partial wave $\phi_{\mathbf{R}l}$ is well defined, positive and less than 1, since contributions from all partial waves at most sum to one and the interstitial also adds a positive a contribution.

Notes:

1. In the ASA, with space-filling spheres, the sum of partial waves comprises the total wave function, and the separate contributions to sum to 1.

2. tbe is not an augmented wave method; this kind of decomposition is not possible.

3. The projection yields a physically measurable quantity (at least in principle). This stands in contrast to the Mulliken projection, which depends on chice of basis set.

Mulliken Analysis

Mulliken analysis is a decomposition of an eigenfunction into the separate orbital contributions. The eigenfunction is written as a linear combination of lmf basis functions $\chi_{\alpha\mathbf{R}L}(\mathbf{r})$ $$\psi_i(E_i,\mathbf{r}) = \sum\limits_{\alpha\mathbf{R}L} {z_{\alpha\mathbf{R}L}^{(i)}\,\chi_{\alpha\mathbf{R}L}(\mathbf{r})} \qquad (4)$$

The $\chi_{\alpha\mathbf{R}L}$ are augmented smooth Hankel functions.

We can decompose $D(E)$ through coefficients $z_{\alpha\mathbf{R}L}^{(i)}$. In this case the $z_{\alpha\mathbf{R}L}^{(i)}$ are eigenvectors of the lmf hamiltonian: they diagonalize both the hamiltonian and overlap. In matrix form $$(z^{(i)})^{\dagger} \, H \, z^{(i)} = E_i \quad \mathrm{and} \quad (z^{(i)})^{\dagger} \, O \, z^{(i)} = 1 \qquad (5)$$

If the overlap matrix were diagonal, it is evident from Eq. (5) that the eigenvectors would satisfy $\sum\nolimits_{\alpha {\mathbf{R}}L} [z_{\alpha \mathbf{R}L}^{(i)}]^\dagger z_{\alpha {\mathbf{R}}L} = 1$. The overlap is not diagonal; however there is a generalization $$\sum\nolimits_{\alpha {\mathbf{R}}L} {(z_{\alpha \mathbf{R}L}^{(i)}})^{-1} z_{\alpha {\mathbf{R}}L}^{} = 1 \qquad (6)$$

$z_{\alpha \mathbf{R}L}^{(i)}$ is a vector; there is one eigenvector $z^{(i)}$ for each of the ${\alpha\mathbf{R}L}$ components. Thus $z$ is a square matrix that can be inverted.

The sum over components ${\alpha\mathbf{R}L}$ in Eq. (5) evaluates to 1 for a particular state $i$, which is a statement of normalization of the eigenfunction. We can decompose or resolve the unit norm into separate elements (“domposition of the norm”), resolving by $\alpha$, $\mathbf{R}$ and $L$. Decomposition by $\alpha$ is not so meaningful (lmf contracts over $\alpha$; while lm and tbe only have a single $\alpha$), but resolving by $\mathbf{R}$ or $L$ can offer a great deal of physical insight. This decomposition is used to assign color weights in band plots. Examples can be found in the plbnds manual and in the lmf band plotting tutorial. How to do it in the partial dos context will be shown below.

How information is assembled for analysis

The following outlines the general procedure for making partial wave analysis. Steps are explained in more detail later, in the examples.

Both Mulliken and partial pave analysis enable decomposition of the unit norm into partial contributions associated with a particular site R and L=lm character. The band programs (lm, lmf, and tbe) will accumulate weights for a partial wave or Mulliken analysis in the course of a usual band pass, by adding a command-line switch --pdos or --mull. Both switches have numerous options that can select or group a subset of all states, to contract over m (leaving the resolution by R and l) or over both l and m, resolving by R only. This is described in more detail below. The band program will write a file moms.ext to disk with information about the partial decomposition.

With moms.ext in hand, run lmdos with exactly the same switch --pdos or --mull, including any modifiers. This should generate a file dos.ext with the requisite information. By default this file will take traditional standard format for DOS files; but you can change the format.

The pldos utility is designed to read DOS files, and select out or combined particular DOS, and either make a postscript file directly (for quick and dirty results) or format the data in easily read formats.

Notes on Partial wave and Mulliken analysis, and their relative merits

1. Mulliken analysis has somewhat imprecise meaning, as the results are dependent on the user’s choice of basis. However, to the basis functions do resemble atomic orbitals, especially for d and f electrons, they are a useful tool.

2. As noted above, partial wave analysis is approximately independent of basis, except for the choice of augmentation radius. As such, it is often preferable to Mulliken analysis.

3. The Jigsaw Puzzle Orbital basis is very short ranged, so association with a particular atomic orbital is more clear. The distinction between partial wave and Mulliken analysis will be much smaller.

4. In the ASA, with space-filling spheres, the sum of partial waves comprises the total wave function, and the separate contributions to sum to 1.

5. lm typically generate the moments file moms.ext automatically as part of the band pass. However the moments are generated for inequivalent classes only; the weights are ordered by class instead of by site. You can run lmpg after a band pass (without argument --pdos) to generate partial DOS for each class (typically resolved by $\ell$ but not by $m$).

6. tbe is not an augmented wave method; partial wave decomposition is not possible.

Making total DOS using band programs lmf, lm, or tbe

The simplest DOS is the total DOS/cell (not resolved into any components). This is automatically generated when you turn on  SAVDOS=T  in category  BZ. A band program (lmf, lm, or tbe) will generate DOS in a particular energy window, on a uniform mesh of points.

Note: you can also cause lmf to generate dos using the command-line argument  --dos. Modifiers to this switch allow you to control the energy mesh (and format) of the dos file.

lmf --input | grep BZ_DOS

lmf --input | grep BZ_NPTS

lmf writes the DOS to dos.ext, normally in the traditional standard format for dos files. You can reformat it yourself to use your favorite graphics package, or use pldos utility to format the dos into standard Questaal format for two-dimensional arrays, which are more easily read by other graphics packages.

At this stage, you can use your favorite graphics package to draw a figure from files dosp.dat and dosp2.dat. Alternatively pldos will have written an fplot script; you can immediately create a postscript file using fplot.

1. Total DOS in elemental Co

Building a Co input file

Copy the contents of the box below into file init.co.

LATTICE
% const a=4.730 c=7.690
  ALAT={a} PLAT= 1.0 0.0 0.0 0.5 0.8660254 0.0 0.0 0.0 {c/a}
SPEC
   ATOM=Co   MMOM=0,0,1.6
SITE
   ATOM=Co X=0 0 0
   ATOM=Co X=1/3 1/3 1/2

Construct the input file in the usual manner, see for example the Si tutorial or the PbTe tutorial:

blm --brief --mag --ctrl=ctrl --wsitex --noshorten --nk~met co

See this page for command-line arguments to blm. In particular, --nk~met sets the number of k points, telling blm to use an internal default blm appropriate for typical metals. Inspect the ctrl file, and you can find these lines

% const met=5 nk1=10 nk2=nk1 nk3=6 gmax=8.9  # BZ-related variables
...
      NKABC=  {nk1},{nk2},{nk3}       # 1 to 3 values

blm selects the k-mesh to be n₁×n₂×n₃ divisions, with n₁=n₂=10 and n₃=6, for 600 points in the full Brillouin zone (with symmetry it reduces to 36 irreducible k points).

This makes a reasonably fine mesh, good enough to generate a reasonably smooth DOS with the tetrahedron method (which lmf normally uses). Coarser meshes will cause the DOS to be much less smooth and this is especially severe if the integration is performed by simple sampling integration. For a really smooth DOS, use a finer mesh.

Self-consistent Co density

Make the density self-consistent:

lmfa co
lmf ctrl.co > out.lmfsc

You should get a reasonably self-consistent density in 10 iterations (confirm with grep DQ out.lmfsc). The last line of the file save.ext should read

c mmom=3.2369729 ehf=-5564.3100904 ehk=-5564.3100902

The magnetic moment/atom is then calculated to be 3.24/2, close to the experimental moment (1.6 μ_B).

Total DOS in Co

Run lmf with the command line argument  --dos  to generate the total DOS

lmf ctrl.co --quit=dos --dos@npts=2001@window=-1,1

--quit=dos tells lmf to stop after the DOS is generated.

Near the end of the standard output the following line should appear:

 ... Generating total DOS

The pldos utility will extract and reconfigure the contents of dos.co, saving the data in a more palatable format :

echo 5  8  -10  10 | pldos -esclxy=13.6 -ef=0 -fplot~tmy=1 -lst=1 -lst2 dos.co
     ↑  ↑   ↑   ↑
   dmx ht emin emax

pldos reads data in the traditional dos file format the band codes normally use. It can make a postscript figure directly, or be used as a preparatory step for fplot or another graphics package. We do the latter here. The figure will put energy on the abscissa and DOS on the ordinate.

pldos takes four arguments from standard input:

• dmx   DOS upper bound in the figure. On account of --esclxy, units will be in (eV-cell)⁻¹
• ht    approximate height of figure, in cm
• emin   minimum energy to draw (left point of abscissa). On account of --esclxy, units will be in eV.
• emax   maximum energy to draw (right point of abscissa)

Note: no interactive input is required from the command as written above. However, you can run pldos in an interactive mode by entering simply pldos dos.co.

Switches to pldos have the following effect:

• −esclxy=13.6   scales the abscissa (energy) to convert it from Ry to eV, and the ordinate (dos) converting it from  Ry⁻¹  to  eV⁻¹.
• −ef=0   Shift the abscissa, putting the Fermi energy at 0.
• −fplot   Set up input for fplot. This entails the following:
  1. Create file dosp.dat for the spin-1 dos, written in the standard Questaal format for 2D arrays.
  2. Create a corresponding file dosp2.dat for the spin-2 dos (applicable only to spin polarized cases).
     Note: −fplot can accept numerous arguments, in this case tmy=1, which sets the spacing for tic marks on the ordinate.
  3. Create an fplot script plot.dos
• −lst=1   Select which channels in the dos file (dos.co) to combine for the majority spin.
       In this case there is only a single channel per spin (there will be many DOS is resolved in parts, as shown later)
• −lst2    Select which channels in the dos file (dos.co) to combine for the second spin.
       −lst2 without arguments tells pldos to copy the list from -lst, incrementing each element by 1.
       Since (majority,minority) dos are interleaved, it simply generates the spin-2 channels counterpart to spin-1.
• dos.co    causes pldos to read DOS from dos.co, formatted in the way lmf usually writes dos files.

The order of switches is not important, but the file name specifying DOS must come last.

File dosp.dat contains the spin-1 dos (majority in this case), and dosp2.dat the negative of the spin-2 (minority) dos, written in Questaal’s standard 2D array format.

Use your favorite graphics package to draw a figure from files dosp.dat and dosp2.dat. Alternatively, use fplot : pldos has already created a script plot.dos for it. Create a handsome picture with the following

fplot -pr10 -f plot.dos

fplot generates a postscript file fplot.ps. Open this file with a postscript viewer of your choice, e.g. gs fplot.ps.

Click here to view the total DOS for Co

Note: The spin-2 dos are scaled by −1 to make it convenient for drawing the figure with majority DOS above zero and minority DOS below it. Energy is on the abscissa with the Fermi level at 0. Co d bands dominate near the Fermi level E_F: they form two broad peaks with the majority d falling completely below E_F and the minority d straddling it.

The DOS at the Fermi energy is essentially the Pauli susceptibility $\chi_P$. The Stoner criterion for a ferromagnetic instability is $I{\times}\chi_P{>}1$. It happens that the Stoner $I$ is 1 eV for the 3d transition metals, to a good approximation. Inspecting the DOS at the Fermi level, it is apparent that the Stoner criterion is well satisfied, and Co should spontaneously magnetize. Indeed the self-consistent calculation we just ran shows that Co is ferromagnetic in DFT.

BZ  SAVDOS=t NPTS=2001 DOS=-1,1 NEVMX=999

mv dos.co dos.bk
lmf ctrl.co --quit=dos
diff dos.co dos.bk

2. Partial DOS in Cr3Si6

Input file and self-consistency in Cr3Si6

Cr₃Si₆ (AKA CrSi₂) is a transition metal silicide with a small band gap [4], measured to be between 0.27 and 0.67 eV. The three Cr and six Si atoms are symmetry-equivalent.

To set up the computational conditions, copy the following box into init.cr3si6.

# Init file for Cr3Si6
LATTICE
    ALAT=8.37 PLAT= sqrt(3/4) -0.5 0.0 0.0 1.0 0.0 0.0 0.0 1.43369176
SITE
    ATOM=Cr   X= 1/2   1/2  -1/6
    ATOM=Cr   X= 0     1/2   1/6
    ATOM=Cr   X= 1/2   0     1/2
    ATOM=Si   X= 1/3   1/6   1/6
    ATOM=Si   X=-1/3  -1/6   1/6
    ATOM=Si   X=-1/6   1/6  -1/6
    ATOM=Si   X= 1/6  -1/6  -1/6
    ATOM=Si   X= 1/6   1/3   1/2
    ATOM=Si   X=-1/6  -1/3   1/2

In this tutorial we will use a small, single-kappa basis. It is not necessary, but it speeds up the calculation with minimal effect on the accuracy since the system is fairly close-packed.

The following sets up the input file and generate a self-consistent density.

blm --brief --ctrl=ctrl --wsitex --nk~met --autobas~loc=0~mto=1 cr3si6
lmfa cr3si6
lmf cr3si6 > out.lmfsc

• --autobas~loc=0 suppresses lmfa’s search for deep lying states to be treated as local orbitals.
• --autobas~mto=1 specifies a single-kappa LMTO basis (minimal basis).
• --nk~met causes blm to selects a default k-mesh of 6×6×4 divisions (21 inequivalent points) This is sufficient to make a good self-consistent density, but it is a little on the coarse for a smooth DOS. (This is a reflection of a general observation that the total energy is variational in the density but the separate contributions kinetic energy and potential energy are not.) When making partial DOS we will make the mesh finer.
• --ctrl=ctrl tells blm to write the input file directly to ctrl.cr3si6.
• --wsitex tells blm to write site positions in crystal coordinates, as they are written in the init file.

When lmf’s completes execution you should find that the 10^th and final iteration is (nearly) self-consistent with and RMS DQ=1.45e-5 (confirm with grep DQ out.lmfsc). You should find that the last line of file save.cr3si6 is

x ehf=-9761.7455954 ehk=-9761.745594

The ‘x’ indicates that lmf did not reach convergence according to the tolerances set in the ctrl file. However, DQ=1.45e-5 is quite acceptable and that ehf and ehk are essentially identical is another good sign that self-consistency has been reached.

Partial DOS in Cr3Si6 resolved by l

Note: A figure for partial DOS similar to the one generated here is shown in the pldos manual. Other features of the --pdos switch and switches to pldos are used there.

You have a choice to resolve DOS can site R, or more finely resolved by both R and angular momentum $\ell$ within R, or even at a finer level by R, $\ell$ and the azimuthal quantum number $m$.

Note: lmf uses real versions of spherical harmonics by default, and the $m$ projections reflect this.

The partial DOS requires a three-step process. lmdos, the utility that generates the partial DOS, requires energy bands and weights generated by a band program (lmf in this case), as explained in the Introduction. These weights are generated when lmf is run with an additional command-line argument --pdos. As a final step, once lmdos is invoked to create the partial DOS, the information must be rendered into a figure. In this tutorial we use the Questaal utility pldos, though you may want to use a different utility.

First, try using --pdos without any modifiers

rm mixm.cr3si6
lmf cr3si6 --quit=rho --pdos

At the beginning of the band pass you should see this line

 sumlst:  Partial DOS mode 2 (all sites lm-projected)  9 sites 171 channels

Each of the 9 sites will decomposed into DOS, resolved by both $\ell$ and $m$. There is a grand total of 171 channels, because DOS are expanded to lmxa=4 for Cr (25 channels) and lmxa=3 for Si (16 channels), as the input file specifies. This is overkill: $\ell$>2 for Cr and $\ell$>1 for Si is of limited interest. Also, most often, it is not necessary to resolve each $\ell$ into individual $m$ components. Here we resolve DOS by $\ell$, leaving the $m$ resolution to Additional Exercises. Since the three Cr and six Si atoms should be equivalent by symmetry, we will generate $\ell$ resolved DOS for only the first Cr (site 1) and the first Si atom (site 4). (Symmetry is explored in more detail in Additional Exercises.) Run this command:

lmf cr3si6 -vnk1=12 -vnk3=8 --quit=rho --pdos~mode=1~sites=1,4~lcut=2,1

• -vnk1=12 -vnk3=8 assigns values to preprocessor variables, which are later used when parsing the token NKABC= {nk1},{nk2},{nk3} to rendering the mesh 12×18×8 divisions (95 inequivalent points). This is double the mesh used for self-consistency, and it is enough to make the DOS look smooth.
• --quit=rho stops lmf after the DOS had been generated
• --pdos~mode=1~sites=1,4~lcut=2,1 makes Cr DOS for $\ell$=0,1,2 and Si DOS for $\ell$=0,1 ; see command line options. DOS are generated only for the first Cr and Si sites.

Information is written into a file in hdf5 format, moms.h5. To see what is in it, try h5ls moms.h5. dosw contains partial DOS weights; it is dimensioned to hold 95 k points, one spin, 81 bands and 5 channels.

Next run lmdos. Use the same command line switches used for lmf.

lmdos cr3si6 -vnk1=12 -vnk3=8 --quit=rho --pdos~mode=1~sites=1,4~lcut=2,1 --dos:npts=1001:window=-1,1

Note there is an extra --dos switch. If you leave it off the you will be prompted at the standard input for three numbers that define the energy mesh (number of points, minimum and maximum energies). Putting it on the command-line with --dos avoids the need to deal with lmdos interactively.

lmdos should generate the following output:

 ASADOS: reading weights from file moms.h5
         expecting file to be resolved by l
         file has 5 channel(s)
         Using npts=1001  emin=-1  emax=1

Don’t worry that the output refers to “ASADOS.” lmdos doesn’t care where the DOS weights and energy bands originated from. lmdos exits with a table showing how the panels are broken out:

 Channels in dos file generated by LMDOS:
 site class label   spin-1
    1    1   Cr     1:3
    4    2   Si     4,5

It is a handy table telling you which DOS are associated with which channel.

To render the $\ell$-resolved DOS into a postscript, use the pldos utility

echo 1.4 5 -8    6 | pldos -fplot~long~open~tmy=.25~dmin=0.40~xl=E -esclxy=13.6 -ef=0 -lst="1;2;3;4;5;" dos.cr3si6
      ↑  ↑   ↑   ↑
    dmx  ht emin emax

fplot generates a postscript file fplot.ps. Open this file with a postscript viewer of your choice, e.g. gs fplot.ps.

This instruction makes 5 panels, corresponding to Cr s, p, d, and Si s and p.

As explained previously, pldos takes four arguments from standard input:

• dmx   DOS upper bound in the figure. On account of --esclxy, units will be in (eV-cell)⁻¹
• ht    approximate height of figure, in cm
• emin   minimum energy to draw (left point of abscissa). On account of --esclxy, units will be in eV.
• emax   maximum energy to draw (right point of abscissa)

The -fplot switch tells pldos to:

• create files and an input for the fplot(/docs/misc/fplot) utility
• The arguments to -fplot modify the style of the figure as explained here.

The remaining switches do the following

• -esclxy=13.6 converts abscissa to eV and ordinate to eV⁻¹
• -ef=0 shifts the energy axis to put the Fermi level at 0
• -lst=”1;2;3;4;5” tells lmdos which DOS to combine into a single data set, and how many data sets to make. It uses  “;“  as a separator to tell lmdos to start a new data set. Each data set gets its own panel. By contrast, -lst=”1,2,3;4,5” would make 2 panels, combine channels 1,2,3 in the first and 4,5 in the second.

Click here to view the l-resolved DOS for Cr3si6

The figure shows that there is a small gap, and that the DOS at the valence and conduction band edges are dominated by Cr d, with some Si p mixed in.

3. DOS and Core-level Spectroscopy in Fe

In this tutorial for Fe, partial wave analysis, Mulliken analysis, and core-level spectroscopy are compared.

Input file and Self-consistent density in Fe

Copy the following into init.fe. As an initial guess, we use a trial moment of 2 μ_B. (Fe is known to have a magnetic moment of 2.2 μ_B/atom in the ground state).

LATTICE
% const a=5.4235
  ALAT={a} PLAT= -0.5 0.5 0.5   0.5 -0.5 0.5   0.5 0.5 -0.5
SPEC
   ATOM=Fe   MMOM=0,0,2
SITE
   ATOM=Fe X=0 0 0

The following sets up the input file and generate a self-consistent density.

blm --brief --mag --ctrl=ctrl --wsitex --noshorten --nk~met fe
lmfa fe
lmf ctrl.fe > out.lmfsc

Confirm that self-consistency is reached by looking at the RMS DQ in the last iteration (confirm with grep DQ out.lmfsc). It should be 1.7e-5. Also look at the last line of save.fe (tail save.fe). You should see a result similar to the following

c mmom=2.2472806 ehf=-2541.0382763 ehk=-2541.038279

That ehf and ehk are essentially identical is another good sign that self-consistency has been reached.

Note: the magnetic moment is not variational in the density, so the value is more sensitive to k convergence than the total energy. Note that blm selected a default mesh of 11×11×11 divisions. Try invoking blm with a finer k mesh, e.g. replace --nk~met with --nk~16.
With 16 divisions you should find mmom is 2.225 μ_B;
with 24 divisions you should find mmom is 2.218 μ_B.
The total energy also changes a little (about 7e-5 eV).

This confirms that for purposes of making the self-consistent density, blm_{:.exec}’s default k mesh is adequate.

Mulliken analysis in Fe

Mulliken analysis and partial DOS are very similar. Perform Mulliken analysis by adding --mull to the the lmf command line, analogous to --pdos for partial DOS (section 2).

As with partial DOS, you have a choice to resolve DOS can site R, or more finely resolved by both R and angular momentum $\ell$ within R, or even at a finer level by R, $\ell$ and the azimuthal quantum number $m$. We use --mull~mode=2~nl=3 here to resolve by both $\ell$ and $m$, and truncate the $\ell$ to spd.

As in the partial DOS case, you (1) make the DOS weights with lmf, (2) use lmdos to combine the energy levels and weights into the

lmf -vnk1=20 fe --quit=rho --mull~mode=2~nl=3
lmdos -vnk1=20 fe --quit=rho --mull~mode=2~nl=3 --dos:npts=1001:window=-.7,.8
cp dos.fe dos-mull.fe

(The last step is there to preserve this dos for later in the tutorial.)

You should see this table just before lmdos exits:

 Channels in dos file generated by LMDOS:
 site class label   spin-1                       spin-2
    1    1   Fe     1:17:2                       2:18:2

Spin channels are interleaved, so odd channels hold spin-1 (majority spin) and even channels minority spin Mulliken dos.

Draw a figure showing each of the 5 d channels

echo .5  5  -6 5 | pldos -esclxy=13.6 -ef=0 -fplot -lst="9;11;13;15;17" -lst2 dos-mull.fe
      ↑  ↑   ↑   ↑
    dmx  ht emin emax
fplot -f plot.dos

View the figure with your favorite postscript viewer, e.g. gs fplot.ps.

The figure below is drawn with the convention that spin-1 DOS is shown above y=0, spin-2 DOS below it.

Click here to view the Mulliken DOS for Fe

DOS is large where d states are present, in the window about (−5,2) eV. Fe has cubic symmetry, which means that the three $t_{2g}$ DOS (panels 1,2,4) and the two $e_g$ DOS (panels 3,5) should be identical. The former are xy, yz, and xz states, which correspond to orbitals 5,6,8 in Questaal’s ordering, and the latter 3z²−1 and x²−y² states correspond to orbitals 7,9. With the interleaving of spins, the spin 1 $t_{2g}$ states correspond to channels (9,11,15) and $e_g$ to channels (13,17) in dos-mull.fe. You can see that panels 1,2,4 appear identical, as do panels 3,5 [5].

lmdos -vnk1=20 fe --quit=rho --mull~mode=2~nl=3 --dos:npts=1001:window=-.7,.8:rdm
cp dos.fe dos-mull-sym.fe

lmf -vnk1=20 fe --quit=rho --mull~mode=2~nl=3 --nosym
lmdos -vnk1=20 fe --quit=rho --mull~mode=2~nl=3 --nosym --dos:npts=1001:window=-.7,.8:rdm
cp dos.fe dos-mull-nosym.fe

mcx dos-mull-sym.fe dos-mull-nosym.fe -- -abs -max:g

plot -colsy 1+10:1+18:2 dos-mull-sym.fe -colsy 1+10:1+18:2 -lt 2,col=1,0,0 dos-mull-nosym.fe
gs fplot.ps

fplot -colsy 11 dos-mull-sym.fe -lt 2,col=1,0,0 -colsy 13 dos-mull-sym.fe -lt 3,col=0,1,0 -colsy 17 dos-mull-sym.fe
gs fplot.ps
fplot -y 0,10 -colsy 15 dos-mull-sym.fe -lt 2,col=1,0,0 -colsy 19 dos-mull-sym.fe
gs fplot.ps

echo 'SYMGRP find SOC=t' >> ctrl.fe

lmf -vnk1=20 fe --quit=rho --mull~mode=2~nl=3 -vso=t
lmdos -vnk1=20 fe --quit=rho --mull~mode=2~nl=3 -vso=t --dos:npts=1001:window=-.7,.8:rdm
cp dos.fe dos-mull-so-sym.fe

fplot -frme 0,1,0,.7 -frmt th=3,1,1 -colsy 11 dos-mull-sym.fe -colsy 11 -lt 2,col=1,0,0 dos-mull-so-sym.fe -colsy 13 -lt 3,col=0,1,0 dos-mull-so-sym.fe

Partial DOS in Fe

The partial DOS proceeds in much the same way as Mulliken analysis presented in the preceding section.

lmf -vnk1=20 fe --quit=rho --pdos~mode=2~nl=3
lmdos -vnk1=20 fe --quit=rho --pdos~mode=2~nl=3 --dos:npts=1001:window=-.7,.8
cp dos.fe dos-pdos.fe

You can draw the partial DOS in the same manner as the Mulliken DOS. In the steps below we combine the two so the similarities and differences can be compared.

echo .999  6  -6 5 | pldos -esclxy=13.6 -ef=0 -fplot~ext=mull~dmin=.4~tmy=.25 -lst="1;3,5,7;9,11,15;13,17" -lst2 dos-mull.fe
echo .999  6  -6 5 | pldos -esclxy=13.6 -ef=0 -fplot~ext=pdos~dmin=.4~tmy=.25 -lst="1;3,5,7;9,11,15;13,17" -lst2 dos-pdos.fe
awk '{if ($NF == "dosp.pdos") {print; sub("pdos","mull");sub("{ltdos}","2,bold=3,col=1,0,0");print} else if ($NF == "dosp2.pdos") {print; sub("pdos","mull");sub("{ltdos}","2,bold=3,col=1,0,0");print} else {print}}' plot.dos  > plot.dos2
fplot -pr10 -f plot.dos2

View the figure with your favorite postscript viewer, e.g. gs fplot.ps.

The first command makes four panels consisting of $s$, $p_x{+}p_y{+}p_z$, $t_{2g}$ and $e_{g}$ Mulliken channels; the second does the same for partial wave decomposition. awk combines the two kinds of data to make a single fplot script, plot.dos2. fplot creates the postscript figure.

Mulliken and partial wave analysis show modest differences in the $s$ and $p$ channels, while the $d$ channels are nearly identical.

Core-level spectroscopy in Fe

You can tell lmf to compute EELS, (also called core-level spectroscopy) by adding --cls to the command-line. This switch acts in a manner roughly similar to the --mull and --pdos switches. In the CLS case, matrix element $\langle\psi_c|\mathbf{r}|\psi_v\rangle$ of the core level with the valence states are calculated, and the number of channels is the number of core-level states. The weights file is written to cls.fe, and includes the matrix element.

You must specify a particular core level. There are several ways to do it. Here we consider the excitation from the 2p level, which can be specified as  --cls:1,1,2  for atom 1, $\ell$=1, n=2.

Run the following to simultaneously generate cls.fe (--cls) for CLS and the total dos (--dos). The two will be compared later.

lmf -vnk1=20 fe --quit=rho --cls:1,1,2 -vmet=2 --dos:npts=1001:window=-.7,.8
mv dos.fe tdos.fe

The --cls tag requires BZ_METAL=2 or BZ_METAL=3.

The second command renames the file containing the total DOS, since that file will be overwritten CLS.

You should see a table of matrix elements printed just before lmf exits

 CLS atom 1 (Fe) n=2 l=1
 Spin 1 ..
 vcdmel: ecor0=-50.757934 ecore=-50.757934
 (not including electrostatic potential shift)
   l   <u|core>     <s|core>    <u|r|core>  <s|r|core>
   0    0.013933   -0.008084   -0.037621    0.018937
   1   -0.000471    0.000136   -0.057505    0.014664
 ...

Make the CLS and rename the file:

lmdos -vnk1=20 --cls --dos:wtfn=cls:npts=1001:window=-.7,.8 fe
mv dos.fe dos-cls.fe

lmdos converts the cls weights into k-integrated spectra in much the same way it converts the moms into DOS. You must supply lmdos with a switch telling it that this is a CLS calculation and also tell it to read data from file cls.

cls.fe should contain three channels per spin for x, y, and z. Concatenate it and the DOS into a single file. They have different units, so the latter is scaled by 1/10 to make them about the same size:

catdos dos-cls.fe -s1/10 tdos.fe

catdos concatenates the two DOS, writing the output to file dos.dat with 8 channels. The first 6 channels are the CLS, the last two are the total DOS.

Draw two panels, one with the three CLS combined and the other with DOS

echo  .5 10 -6 6 |  pldos -esclxy=13.6 -fplot -lst="1,3,5;7" -lst2 dos.dat
fplot -pr10 -f plot.dos

View the figure with your favorite postscript viewer, e.g. gs fplot.ps.

The two panels look similar.

Notes:

1. the final state is better described in terms the self-consistent density in the presence of a core hole (sudden approximation). In other words, you should compute the matrix elements with the core partially occupied. We ignore that step here, but see the tutorial below.
2. the $m$ resolved CLS should be calculated without symmetry operations, since the code does not rotate irreducible k points to the star of k [4] for CLS. See Additional exercises.

4. Core-level spectroscopy in the presence of a core hole in CrN

EELS, also known as core-level spectroscopy, involves the excitation of a core electron to an excited state. In this tutorial we demonstrate core-level spectroscopy for the N 1$s$ state, in CrN [5].

The self-consistent calculation proceeds with an electron missing from the N 1$s$ core, which corresponds to the ‘sudden approximation’ (system relaxes instantaneously from electron ejected out of a core hole).

CrN input file

The electron should be ejected from a single, isolated N atom. We must use periodic boundary conditions, so we simulate it with a supercell of 4 N atoms, with one (Nh) differentiated as having a core hole.

This tutorial uses an already-built file from the source directory ~/lm Note: your source directory might be in a different place. Substitute ~/lm for the top-level directory where you installed Questaal.

This file is essentially similar to the input used in Ref [5]. A tutorial detailing the steps required to generate a basic input file can be found here.

cp ~/lm/fp/test/ctrl.crn .

Inspect ctrl.crn, and note in particular this line in species Nh:

  C-HOLE=1s C-HQ=-1,-1

These tags tell lmf to fractionally occupy the Nh $1s$ state, with one missing electron and a magnetic moment of −1.

Self-consistent density in CrN

This test runs a little slowly, so we use lmf in the MPI mode. If you only have the serial mode installed, just remove mpirun 8.

lmfa crn
mpirun -n 8 lmf crn -vnit=50 > out.lmfsc

In the output of lmfa, the part concerning Nh, you should see a table like this:

 Species Nh:  Z=7  Qc=1  R=1.800000  Q=0  mom=-1
 mesh:   rmt=1.800000  rmax=23.783012  a=0.03  nr=223  nr(rmax)=309
         Add core hole:  kcor=1  lcor=0  qcor=-1  amom=-1
  Pl=  2.5     2.5     3.5     4.5     spn 2   2.5     2.5     3.5     4.5
  Ql=  1.0     2.0     0.0     0.0     spn 2   1.0     2.0     0.0     0.0

The core charge has only one electron; the net magnetic moment of the system is −1 because the core hole has a magnetic moment.

lmf also indicates that a core hole is present. These lines of its output:

 site  5  z=  7.0  rmt= 1.80000  nr=223   a=0.030  nlml=16  rg=0.450  Vfloat=T
 core hole:  kcor=1  lcor=0  qcor=-1  amom=-1

indicate that the core hole is present on atom 5 with principal and $\ell$ quantum numbers 1 and 0, respectively.

lmf should converge to self-consistency in 32 iterations and terminate with the following:

c nit=50 mmom=-.775742 ehf=-8797.2800063 ehk=-8797.2800322

Thus the core hole’s magnetic moment (−1) is partially screened by the valence electrons.

Core-hole spectroscopy in CrN

You must specify a particular core level, which can be done in several ways. For CLS of the $1s$ state on Nh (the fifth atom), do:

mpirun -n 8 lmf --rs=1,0 --cls:5,0,1 -vnit=1 -vmetal=2 -vnk=8 crn
lmdos --dos:cls:window=0,1:npts=1001 --cls crn -vnk=8

lmf makes weights and stores dos decorated by matrix elements in cls.crn. Some matrix elements are printed just before lmf exits:

 CLS atom 5 (Nh) n=1 l=0
 Spin 1 ..
 vcdmel: ecor0=-29.419700 ecore=-29.419695
 (not including electrostatic potential shift)
   l   <u|core>     <s|core>    <u|r|core>  <s|r|core>
   0   -0.000090   -0.000273    0.086293   -0.081045
   ...

lmdos converts cls.crn weights into k-integrated spectra. You must supply lmdos with a switch telling it that this is a CLS calculation and also tell it to read data from file cls cls.crn. This file should contain three channels per spin, for x, y, and z . Convert the result into a more convenient form using the pldos utility:

echo .25 8 0 1 | pldos -fplot -lst="1;3;5" -lst2="2;4;6" dos.crn
fplot -pr10 -f plot.dos

View the figure with your favorite postscript viewer, e.g. gs fplot.ps.

Click here to view the core-level spectroscopy for CrN

Note the difference in x, y, and z spectroscopies.

Additional exercises

1) Try making the m-resolved partial DOS in Co. The moment is 1.6 μ_B, and since all the majority d states are filled, slightly less than two minority states are empty. It is a worthwhile exercise to find out which of the 5 d states are the unoccupied ones.

2) This exercise picks up on Partial DOS in Cr3Si6 resolved by $\ell$, but now resolving DOS by both $\ell$ and $m$.

Run lmf again, this time globally limiting the maximum number of $\ell$’s to 3 at any site. This time we compute DOS for all 9 sites. The number of channels should be 81:

rm mixm.cr3si6
lmf cr3si6 --quit=rho --pdos~nl=3
lmdos cr3si6 --dos:npts=1001:window=-1,1 --pdos~nl=3

Here is the breakdown of channels:

 Channels in dos file generated by LMDOS:
 site class label   spin-1
    1    1   Cr     1:9
    2    1   Cr     10:18
    3    1   Cr     19:27
    4    2   Si     28:36
    5    2   Si     37:45
    6    2   Si     46:54
    7    2   Si     55:63
    8    2   Si     64:72
    9    2   Si     73:81

While the three Cr and six Si atoms should be equivalent by symmetry, the orbitals of a given $\ell$ transform into one another for the different sites. The sum over all m for a particular $\ell$ should be the same for symmetry equivalent atoms, but any particular m may look different.

The following steps combine the 5 Cr d orbitals on each atom, into three panels (panel 1 for the first Cr, panel 2 for the second Cr, panel 3 for the third Cr). Then we can check to see how well this invariance is kept.

echo 40 5 -.5 .5 | pldos -ef=0 -fplot -lst="5:9;14:18;23:27" dos.cr3si6

This sets up DOS in three panels.  --lst tells lmdos which DOS to combine into a single data set, and how many data sets to make. It uses  “;“  as a separator to tell lmdos to start a new data set. Each data set gets its own panel. DOS in channels 5:9, 14:18, and 23:27 are each added together to create DOS respectively for the first, second, and third panels. Note from the table above that these channels correspond to d orbitals for the first, second and third Cr atoms.

DOS is drawn on a Ry energy scale in this case. pldos creates a data file dosp.dat in the standard Questaal format for two-dimensional arrays. pldos also generates a script file plot.dos readable by fplot. To see a picture do:

fplot -pr10 -f plot.dos
gs fplot.ps

You can compare directly the last three columns in dosp.dat, to check how similar they are. This is easily accomplished with the mcx calculator:

mcx dosp.dat -e2 x2-x3 x2-x4

The differences are much smaller than the dos itself; but evidently there is some difference. This is largely an artifact of incomplete k convergence. Repeat the calculation with more k points

lmf cr3si6 --quit=rho --pdos~nl=3 -vnkabc=-1000
lmdos cr3si6 --dos:npts=1001:window=-1,1 --pdos~nl=3 -vnkabc=-1000
echo 40 5 -.5 .5 | pldos -ef=0 -fplot -lst="5:9;14:18;23:27" dos.cr3si6
mcx dosp.dat -e2 x2-x3 x2-x4

and the error becomes much smaller.

The Co $3z^2{-}1$ orbital (the $m{=}0$ d orbital) also should not depend on the Cr atom. This orbital is the middle one (7 for Cr₁, 16 for Cr₂, 25 for Cr₃). Do the following:

echo 20 5 -.5 .5 | pldos -ef=0 -fplot -lst="7;16;25" dos.cr3si6
fplot -pr10 -f plot.dos
open fplot.ps

The three panels should look nearly identical.

The following creates a figure with the following panels:

• The s Co orbital (first atom)
• The sum of Co p orbitals (first atom)
• The 5 Co d orbitals, each given its own panel (first atom)
• The Si s orbital (fourth atom)
• The sum of Si p orbitals (fourth atom)

This makes a grand total of 9 panels.

The command below sets up figure in eV units.

echo .5 3 -6 6 | pldos -fplot~long~open~tmy=.125~dmin=0.20~xl=E -esclxy=13.6 -ef=0 -lst="1;2,3,4;5;6;7;8;9;28;29:31" dos.cr3si6

Modifiers to the −fplot switch alter plot.dos, to “prettify” the figure when fplot generates it. To see what they do, run pldos with no arguments.

Make and display a postscript figure:

fplot -pr10 -f plot.dos
gs fplot.ps

You should see a figure like the one shown below.

The first and second panels (Cr s and p) show very little DOS near the Fermi level. Panels 3 through 7 show the five Cr d DOS; they dominate the electronic structure near the Fermi level (shown by the blue dot-dashed line). The Si p also makes a significant contribution.

From an aesthetic perspective, the autogenerated script plot.dos makes a reasonable figure but some tweaking is needed.

• number labelling for adjacent panels collide with each other. This can easily be rectified by editing plot.dos and making a global change dmax=0.5 → dmax=0.499.

• Labels are needed and the energy axis label (-xl  E) needs some improvement. Try one of:

  -xl  '&\{E} (eV)'
  -xl  "~\{w} (eV)"
  -lbl 3.8,-.25  rd '~\{w} (eV)'
  

3) For Cr₃Si₆, try improving the basis and see the effect on the DOS. Remove switches --loc=0 --mto=1 from the command line when running blm.

4) For Cr₃Si₆, try improving the basis and see the effect on the DOS.

1. The x, y and z components of the CLS in Fe should be equivalent, but they are not, because wave functions are not rotated [5].

   Rerun the calculation suppressing symmetry operations (lmf --nosym and lmdos --nosym) and confirm that they are essentially identical. Try also the CrN case. There the symmetry operations have no effect.

References

1. A. T. Paxton, M. van Schilfgaarde, M. MacKenzie and A. J. Craven, ``The Near-edge Structure in Energy-Loss Spectroscopy: Many-Electron and Magnetic Effects in Transition Metal Nitrides and Carbides,’’ J. Phys. Cond. Mat.12, 729 (2000).

2. Equation (1) only applies to noninteracting effective hamiltonians. In the interacting case the δ-function gets broadened, as described in this tutorial. It is similarly broadened in the Coherent Potential Approximation.

3. One decomposition of the DOS is to resolve the charge density by energy, or just evaluate the charge density for a single state. We do not consider such a decomposition in this tutorial.

4. H. Lange H, Phys. Status Solidi b201, 3 (1997)

5. That the Fe xy, xz, and yz DOS come out the same is non-trivial. To resolve DOS by m the full Brillouin zone must be used. lmf rotates the eigenstate at an irreducible point in the Brillouin to each point in the star of k, to accumulate partial DOS.