# 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_{3}Si_{6}; 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:

- The Green’s function packages
**lmgf**and**lmpg**compute partial DOS. - The DOS the 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. - The optics modes in the
**lmf**and**lm**codes also enable you to resolve DOS in other forms

- Total DOS in Co

Set up input file

```
blm --mag --ctrl=ctrl --wsitex --noshorten co
lmfa --basfile=basp co
blm --mag --ctrl=ctrl --wsitex --noshorten --gmax=8.5 --nk=-2000 co
```

Self-consistent density

```
lmfa --basfile=basp co
lmf ctrl.co
```

Total Co dos

```
lmf ctrl.co --quit=dos --dos@npts=2001@window=-1,1
echo 5 8 -10 10 | pldos -esclxy=13.6 -ef=0 -fplot -lst=1 -lst2 dos.co
fplot -pr10 -f plot.dos
open fplot.ps
```

- Partial DOS in Cr
_{3}Si_{6}

Set up input file and self-consistent density:

```
blm --loc=0 --mto=1 --ctrl=ctrl --wsitex cr3si6
lmfa --basfile=basp cr3si6
blm --loc=0 --mto=1 --ctrl=ctrl --wsitex --gmax=7.1 --nk=-200 cr3si6
lmfa --basfile=basp cr3si6
lmf cr3si6
```

Partial DOS resolved by l

```
lmf cr3si6 -vnkabc=-5000 --quit=rho --pdos~mode=1~sites=1,4~lcut=2,1
lmdos cr3si6 -vnkabc=-5000 --quit=rho --pdos~mode=1~sites=1,4~lcut=2,1 --dos:npts=1001:window=-1,1
echo 1.4 5 -8 6 | pldos -fplot~sh~long~open~tmy=.25~dmin=0.40~xl=E -esclxy=13.6 -ef=0 -lst="1;2;3;4;5;" dos.cr3si6
open fplot.ps
```

Partial DOS resolved by l and m

```
rm mixm.cr3si6
lmf cr3si6 --quit=rho --pdos~nl=3
lmdos cr3si6 --dos:npts=1001:window=-1,1 --pdos~nl=3
echo .5 3 -6 6 | pldos -fplot~sh~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
open fplot.ps
```

- DOS and Core-level Spectroscopy in Fe

Set up input file and self-consistent density

```
blm --mag --ctrl=ctrl --wsitex --noshorten fe
lmfa --basfile=basp fe
blm --mag --ctrl=ctrl --wsitex --noshorten --gmax=8.3 --nk=16 --nit=20 fe
lmfa --basfile=basp fe
lmf ctrl.fe
```

Mulliken DOS

```
lmf fe --quit=rho -vso=f --mull~mode=2~nl=3
lmdos fe -vso=f --mull~mode=2~nl=3 --dos:npts=1001:window=-.7,.8
cp dos.fe dos-mull.fe
echo .5 5 -6 5 | pldos -esclxy=13.6 -ef=0 -fplot -lst="9;11;13;15;17" -lst2 dos-mull.fe
fplot -pr10 -f plot.dos
open fplot.ps
```

Partial DOS, and comparison to Mulliken DOS

```
lmf fe --quit=rho -vso=f --pdos~mode=2~nl=3
lmdos fe -vso=f --pdos~mode=2~nl=3 --dos:npts=1001:window=-.7,.8
cp dos.fe dos-pdos.fe
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
open fplot.ps
```

Core-level spectroscopy

```
lmf fe --quit=rho --cls:1,1,2 -vmet=2 --dos:npts=1001:window=-.7,.8
mv dos.fe tdos.fe
lmdos --cls --dos:wtfn=cls:npts=1001:window=-.7,.8 fe
mv dos.fe dos-cls.fe
catdos dos-cls.fe -s1/10 tdos.fe
echo .5 10 -6 6 | pldos -esclxy=13.6 -fplot -lst="1,3,5;7" -lst2 dos.dat
fplot -pr10 -f plot.dos
open fplot.ps
```

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

Set up self-consistent density

```
cp ~/lm/fp/test/ctrl.crn .
lmfa gas
mpirun -n 8 lmf crn -vnit=50
```

Core-level spectroscopy

```
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
echo .25 8 0 1 | pldos -fplot -lst="1;3;5" -lst2="2;4;6" dos.crn
fplot -disp -pr10 -f plot.dos
open fplot.ps
```

### Table of Contents

- Preliminaries
- Introduction
- Making total DOS using band programs lmf, lm, or tbe
- 1. Total DOS in elemental Co
- 2. Partial DOS in Cr3Si6
- 3. DOS and Core-level Spectroscopy in Fe
- 4. Core-level spectroscopy in the presence of a core hole in CrN
- Further exercises
- References

### Preliminaries

This tutorial assumes you have cloned and built the Questaal repository (located here). For the purpose of demonstration, *~/lm* will refer to top-level (source) directory of the cloned repository. In practice, this directory can be named differently. Questaal executables such **lmf**, **lmdos**, **pldos**, and **catdos** are required assumed to be in your path.

You will also need a postscript viewer. This document assumes you are using the apple-style **open** command to view postscript files.

### Introduction

The density-of-states is given by a sum over states *i* as [2]

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 into component parts. Note that an eigenstate is normalized: . 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 . The full energy-dependence of 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 *l* that contribute to the DOS: and the energy derivative . The (,) 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 *l* and *m* quantum numbers by a compound index *L*, and labeling and respectively as and , the eigenfunction can be projected onto an augmentation sphere centered at **R** as

denotes a projection onto augmentation sphere **R**, ranges from 0 to 1 (and if local orbitals are present, encompasses them), up to the augmentation cutoff **lmxa**.

Coefficients (which is determined from a solution of the secular matrix) then represent a particular kind of decomposition of . Assuming the (, ) 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

where is a linear combination of , (and possibly local orbitals) normalized so that .

The make up partial contributions to , Eq. (1). The contribution to from a particular partial wave 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:*

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.

**tbe**is not an augmented wave method; this kind of decomposition is not possible.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

The are augmented smooth Hankel functions.

We can decompose through coefficients . In this case the are eigenvectors of the **lmf** hamiltonian: they diagonalize both the hamiltonian and overlap. In matrix form

If the overlap matrix were diagonal, it is evident from Eq. (5) that the eigenvectors would satisfy . The overlap is not diagonal; however there is a generalization

is a vector; there is one eigenvector for each of the components. Thus is a square matrix that can be inverted.

The sum over components in Eq. (5) evaluates to 1 for a particular state , so we can decompose or resolve the unit norm into separate elements, resolving by , and . Decomposition by is not so meaningful (**lmf** contracts over ; while **lm** and **tbe** only have a single ), but resolving by or 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

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.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.

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.

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.

**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 l but not by m).**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** will generate DOS in a particular energy window. Tag **BZ_DOS** specifies the energy window, **BZ_NPTS** the number of points.

If **BZ_DOS** *is* specified in the input file, **lmf** will use the specified window. Otherwise **lmf** will select the window as follows. It makes a rough estimate of the Fermi level from the first *k* point, subtracts 0.5 Ry from the first eigenvalue, and adds 0.5 Ry to the estimate for the Fermi level.

However, if you further use the command-line switch `--no-fixef0`

, default values are used. You can find them with

```
lmf --input | grep BZ_DOS
```

If **BZ_NPTS** is specified, it uses the specified value for the number of points. Otherwise it uses a default, which you can find by invoking

```
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 --mag --ctrl=ctrl --wsitex --noshorten co
lmfa --basfile=basp co
blm --mag --ctrl=ctrl --wsitex --noshorten --gmax=8.5 --nk=-2000 co
```

See this page for command-line arguments to **blm** and this page for arguments to **lmfa**.

The mesh density cutoff set by `--gmax=8.5`

can not be determined in advance, but a good value can be obtained from the output of **lmfa**, as explained in the introductory tutorials. **blm** is run twice, once to set up **lmfa**, and again to set the final version of the input file, *ctrl.co*.

`--nk=-2000`

sets the number of *k* points. The negative value is a flag telling **lmf** to find a set of (*n*_{1},*n*_{2},*n*_{3}) divisions of the reciprocal lattice such that the number of microcells *n*_{1}×*n*_{2}×*n*_{3} is approximately 2000, while rendering the microcells as close to equidimension as possible (**G*** _{i}*/

*n*as uniform as possible). 2000 points makes a reasonably fine mesh, good enough to generate a reasonably smooth DOS with the tetrahedron method. 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.

_{i}#### Self-consistent Co density

Make the density self-consistent:

```
$ lmfa --basfile=basp co
$ lmf ctrl.co
```

You should get a reasonably self-consistent density in 10 iterations. The last line of the file *save.ext* should read

```
c mmom=3.1596775 ehf=-5564.3100279 ehk=-5564.310028
```

The magnetic moment/atom is then calculated to be 3.1596775/2, close to the experimental moment (1.6 *μ _{B}*).

#### Total DOS in Co

Use the command line argument **--dos** to generate the total DOS with **lmf**:

```
$ 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 -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.

**pldos** takes four arguments from standard input:

**dmx**DOS upper bound in the figure**ht**approximate height of figure, in cm**emin**minimum energy to draw (left point of abscissa)**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^{−1}to eV^{−1}.**−ef=0**Shift the abscissa, putting the Fermi energy at 0.**−fplot**Set up input for**fplot**. This entails the following:- Create file
*dosp.dat*for the spin-1 dos, written in the standard Questaal format for 2D arrays. - Create a corresponding file
*dosp2.dat*for the spin-2 dos (applicable only to spin polarized cases).

*Note:*The spin-2 dos are scaled by −1 to make it convenient for drawing the figure (see below). - Create an
**fplot**script*plot.dos*

- Create file
**−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; but when the dos is resolved into components there will be many.**−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 **fplot**. Create a handsome picture with the following

```
$ fplot -pr10 -f plot.dos
$ open fplot.ps
```

The majority spin dos is shown above *y*=0, the minority spin below it. (Energy is on the *x* axis 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*and the minority

_{F}*d*straddling it.

Add this line to *ctrl.co* :

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

Copy the original *dos.co* to a backup and invoke **lmf** without a command-line argument a

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

The last line compares the two dos. There should be no difference.

However, remove tag **NEVMX=999** and remake the dos. Now a small difference appears near **emax**. **lmf** does not necessarily compute all the eigenvalues and eigenvectors. When **NEVMX** is present it specifies how many eigenfunctions to make. Now all eigenvalues are found since the rank of the hamiltonian is 36, much less than 999, and there is no loss. (Command line switch **--dos** also causes **lmf** to generate all the eigenvalues.)

### 2. Partial DOS in Cr3Si6

#### Input file and self-consistency in Cr3Si6

Cr_{3}Si_{6} (AKA CrSi_{2}) 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 steps up the input file and generate a self-consistent density. The purpose for the first invocation of the (**blm**,**lmfa**) pair is solely to determine **gmax**:

```
$ blm --loc=0 --mto=1 --ctrl=ctrl --wsitex cr3si6
$ lmfa --basfile=basp cr3si6
$ blm --loc=0 --mto=1 --ctrl=ctrl --wsitex --gmax=7.1 --nk=-200 cr3si6
$ lmfa --basfile=basp cr3si6
$ lmf cr3si6
```

**--loc=0**suppresses**lmfa**’s search for deep lying states to be treated as local orbitals.**--mto=1**specifies a single-kappa LMTO basis (minimal basis).**--nk=-200**specifies a somewhat coarse*k*mesh of 7×7×4 divisions (24 inequivalent points)**--ctrl=ctrl**tells**blm**to write the input file directly*ctrl.cr3si6*.**--wsitex**tells**blm**to write site positions in crystal coordinates. as they are in the*init*file.**–basfile=basp**tells**lmfa**to write the*basp*file directly to*basp.cr3si6*.

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

```
x ehf=-9761.745587 ehk=-9761.7455858
```

#### 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.

DOS can be decomposed by site **R**, by **R** and *l* within **R**, and by **R** and both *l* and *m* within **R**.

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 *l* 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: *l*>2 for Cr and *l*>1 for Si is of limited interest. Also, most often, it is not necessary to resolve each *l* into individual *m* components. Here we resolve DOS by *l*, leaving the *m* resolution to Further Exercises. Since the three Cr and six Si atoms should be equivalent by symmetry, we will generate *l* resolved DOS for only the first Cr (site 1) and the first Si atom (site 4). (Symmetry is explored in more detail in Further Exercises.) Run this command:

```
$ lmf cr3si6 -vnkabc=-5000 --quit=rho --pdos~mode=1~sites=1,4~lcut=2,1
```

**-vnkabc=-5000**makes a much finer*k*mesh so that 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*l*=0,1,2 and Si DOS for*l*=0,1 ; see command line options.

Data is written into file *moms.cr3si6*.

Next run **lmdos**. You *must* include the **--pdos** switch in exactly the same way you ran **lmf** to make *moms.cr3si6*.

```
$ lmdos cr3si6 -vnkabc=-5000 --quit=rho --pdos~mode=1~sites=1,4~lcut=2,1 --dos:npts=1001:window=-1,1
```

If you leave off the **--dos** switch you will be prompted for three numbers that define the energy mesh (number of points, minimum and maximum energies).

**lmdos** should generate the following output:

```
ASADOS: reading weights from file MOMS
expecting file to be resolved by l
file has 5 channel(s)
Using npts=1001 emin=-1 emax=1
IOMOMQ: read 308 qp efermi=0.163599 vmtz=0.000000
```

Don’t worry that the output refers to “ASADOS.” **lmdos** handles both ASA and FP cases. **lmdos** exists 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
```

Generate a postscript file using **pldos**:

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

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

The standard input tells **pldos** to:

- limit the maximum DOS to 1.5 eV
^{−1} - A panel should not exceed 5 cm
- DOS spans (−8,6) eV

The **-fplot** switch tells **pldos** to:

- create files and an input for the
**fplot**(/docs/misc/fplot) utility - Run
**-fplot**in a subshell to generate the*fplot.ps*directly - The other
**pldos**command line switches modify the style of the figure.

The remaining switches do the following

**-esclxy=13.6**converts abscissa to eV and ordinate to eV^{−1}**-ef=0**puts the Fermi level at 0**-lst=”1;2;3;4;5”**creates 5 panels

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 has a moment of 2.2

*μ*in the ground state).

_{B}```
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
```

Construct the input file. The first two lines are present to determine a value for **gmax**.

```
$ blm --mag --ctrl=ctrl --wsitex --noshorten fe
$ lmfa --basfile=basp fe
$ blm --mag --ctrl=ctrl --wsitex --noshorten --gmax=8.3 --nk=16 --nit=20 fe
```

Make the Fe density self-consistent:

```
$ lmfa --basfile=basp fe
$ lmf ctrl.fe
```

You should find that the at the last (11^{th}) iteration the density is self-consistent with a magnetic moment of 2.23 *μ _{B}*.

*Note:* the magnetic moment is not variational in the density, so the value is more sensitive to *k* convergence than the total energy. With 13 divisions you should find a value of 2.168 *μ _{B}*; with 24 divisions you should find a value of 2.218

*μ*.

_{B}#### Mulliken analysis in Fe

Here we perform Mulliken analysis. Use `--mull~mode=0`

to resolve by **R** only – not meaningful here since there is only a single site; Use `--mull~mode=1`

to resolve by **R** and *l*. We use `--mull~mode=2~nl=3`

to resolve by both *l* and *m*, and limit the number of *l* to *spd*.

```
$ lmf fe --quit=rho -vso=f --mull~mode=2~nl=3
$ lmdos fe -vso=f --mull~mode=2~nl=3 --dos:npts=1001:window=-.7,.8
$ cp dos.fe dos-mull.fe
```

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
$ fplot -pr10 -f plot.dos
$ open fplot.ps
```

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

Mulliken DOS are confined to a window of about (−4,3) eV. Note that the three and two states 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 3*z*^{2}−1 and *x*^{2}−*y*^{2} states correspond to orbitals 7,9. With the interleaving of spins, the spin 1 states correspond to channels (9,11,15) and to channels (13,17) in *dos-mull.fe*. You can see that panels 1,2,4 appear identical, as do panels 3,5 [5].

**lmf** symmetrizes the Mulliken DOS by rotating each irreducible *k* point to all points in the star, so that the full Brillouin zone is used. The DOS should be the same whether symmetry operations are used or not.

To test to what extent this is actually the case, remake the DOS just calculated, but this time writing to an **mcx**-readable format:

```
lmdos fe -vso=f --mull~mode=2~nl=3 --dos:npts=1001:window=-.7,.8:rdm
cp dos.fe dos-mull-sym.fe
```

Remake the DOS without symmetry operations:

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

Find the globally maximum *change* in DOS in any channel:

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

This shows that there are some numerical errors even in something that should be formally exact.

For a visual comparison, try the following, which compare the minority-spin *d* channels with and without symmetry

```
fplot -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
open fplot.ps
```

The list of columns (**1+10:1+18:2**) follows standard Questaal syntax for integer lists, and correspond to columns 11,13,15,17,19 (in this format columns are shifted by one to accommodate the energy in the first column.)

The red and black curves should be nearly identical. Also the three and two states should be identical. Visually inspect the difference in the states as follows:

```
$ 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
$ open fplot.ps
$ fplot -y 0,10 -colsy 15 dos-mull-sym.fe -lt 2,col=1,0,0 -colsy 19 dos-mull-sym.fe
$ open fplot.ps
```

The wave function rotator has not yet been implemented in the noncollinear case [5]. To see the effect of SO coupling, the original calculation (without SO coupling) must be run without symmetry operations to compare against the result with SO coupling:

The following steps accomplish this:

```
lmf fe --quit=rho -vso=f --nosym --mull~mode=2~nl=3
lmdos fe -vso=f --nosym --mull~mode=2~nl=3 --dos:npts=1001:window=-.7,.8:rdm
cp dos.fe dos-mull-noso.fe
lmf fe --quit=rho -vso=t --nosym --mull~mode=2~nl=3
lmdos fe -vso=t --nosym --mull~mode=2~nl=3 --dos:npts=1001:window=-.7,.8:rdm
cp dos.fe dos-mull-so.fe
```

Without spin orbit coupling, the *xy* and *xz* DOS should be identical. But SO coupling reduces the symmetry and this is no longer the case. Make the following figure:

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

It draws the *xy* without SO in black, the *xy* with SO in red, and the *xz* with SO as a dotted green line. You can see how SO modifies the *xy* orbital. SO modifies the *xz* orbital in a similar, but slightly different manner, owing to the reduction in symmetry.

#### Partial DOS in Fe

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

```
lmf fe --quit=rho -vso=f --pdos~mode=2~nl=3
lmdos fe -vso=f --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
$ open fplot.ps
```

The first command makes four panels consisting of , , and 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 and channels, while the channels are nearly identical.

#### Core-level spectroscopy in Fe

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 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 2*p* level, which can be specified as **--cls:1,1,2** for atom 1, *l*=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 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.008083 -0.037621 0.018932
1 -0.000471 0.000136 -0.057507 0.014659
...
```

Make the CLS and rename the file:

```
$ lmdos --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
$ open fplot.ps
```

The two panels look very similar.

*Notes:*- 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 this tutorial.
- 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 Further 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 state, in CrN [5].

The self-consistent calculation proceeds with an electron missing from the N 1 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*. 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 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 gas
$ mpirun -n 8 lmf crn -vnit=50
```

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 *l* 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
```

#### Core-hole spectroscopy

You must specify a particular core level, which can be done in several ways. For CLS of the 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. 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
```

It generates *dosp.dat*, easily read by standard graphics packages, and a script *plot.dos* file readable by the **fplot** graphics package in particular. Make and view a postscript figure of the CLS DOS:

```
$ fplot -disp -pr10 -f plot.dos
$ open fplot.ps
```

### Further exercises

- This exercise picks up on Partial DOS in Cr3Si6 resolved by
*l*, but now resolving DOS by both*l*and*m*.

Run **lmf** again, this time globally limiting the maximum number of *l*’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 *l* transform into one another for the different sites. The sum over all *m* for a particular *l* 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
$ open 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 orbital (the *d* orbital) also should not depend on the Cr atom. This orbital is the middle one (7 for Cr_{1}, 16 for Cr_{2}, 25 for Cr_{3}). 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
open 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)'`

For Cr

_{3}Si_{6}, try improving the basis and see the effect on the DOS. Remove switches**--loc=0 --mto=1**from the command line when running**blm**.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

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

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.

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.