Tutorial on Partial Densities of States

This tutorial explains how to generate partial densities-of-states with the Questaal package. The total density-of-states (DOS) can be decomposed in multiple ways. Questaal has two forms, projecting the DOS into partial waves onto augmentation spheres, and by Mulliken analysis. Programs lmf, lm, lmgf, lmpg, tbe can generate DOS. This tutorial will focus mainly on lmf, using elemental Co as a demonstration.

Table of Contents

Construct the Co 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

Make the Co density self-consistent

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

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 directory of the cloned repository. In practice, this directory can be named differently. Questaal executables such lmf, lmdos, and fplot 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 [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 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 [2]: Questaal offers two kinds, “partial waves” and “Mulliken analysis.” Core-level spectroscopy (rather closely relted 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 . is expanded to linear order in a Taylor series in energy, so 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 completely span the hilbert space inside augmention 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 repectively as and , the eigenfunction can be projected onto an augmentation sphere centered at R as

denotes a projection onto augmenation 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:

  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.

Mulliken Analysis

Mulliken analysis is a decomposition of an eigenfunction into the separate orbital contributions. The eigenfunction is writtten 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 may be treated 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 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

  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 l but not by m).

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

How make 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 immediatly create a postscript file using fplot.

Building the input file for Co

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 (n1,n2,n3) divisions of the reciprocal lattice such that the number of microcells n1×n2×n3 is approximately 2000, while rendering the microcells as close to equidimension as possible (Gi/ni 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.

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:
    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: The spin-2 dos are scaled by −1 to make it convenient for drawing the figure (see below).
    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; 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 fplot can read to create a handsome picture, as shown below.

Total DOS for Co

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 EF: they form two broad peaks with the majority d falling competely below EF and the minority 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.)

Input file and self-consistent density for Cr3Si6

Cr3Si6 (AKA CrSi2) is a transition metal silicide with a small band gap [3], 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 the coordinates in the site as multiples of the lattice vectors, 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 10th 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

DOS can be decomposed by site R, by R and l wihin 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. Run lmf again, this time limiting the maximum number of l’s to 3 at any site. Now the number of channels should be 81

rm mixm.cr3si6
lmf cr3si6 --quit=rho --pdos~nl=3

Data is written into file moms.cr3si6. Next run lmdos:

lmdos cr3si6 --dos:npts=1001:window=-1,1 --pdos~nl=3

You must include the --pdos switch in exactly the same way you used it to make moms.cr3si6.

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 and m
         file has 81 channel(s)
         Using npts=1001  emin=-1  emax=1
 IOMOMQ: read 24 qp  efermi=0.160300  vmtz=0.000000

Don’t worry that the output refers to “ASADOS.” lmdos handles both ASA and FP cases. Next comes a list of channels indicating how the grand total of 171 channels is divided up.

 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

Also the three Cr and six Si atoms should be equivalent by symmetry; however 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.

The steps following 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 Cr1, 16 for Cr2, 25 for Cr3). 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.

Partial DOS for Cr3Si6

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.5dmax=0.499.
  • Labels are needed and the energy axis label (-xl  E) needs some improvment. Try one of:

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

Core-level spectroscopy (EELS), Mulliken analysis, partial DOS in Fe

Self-consistent density

Copy the following into init.fe. We use a trial moment of 2 μB. Fe has a moment of 2.2 μB 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

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 (11th) iteration the density is self-consistent with a magnetic moment of 2.23 μB. Note: the magnetic moment is not variational, so the value is more sensitive to k convergence than the total energy. With nkabc=13 you should find a value of 2.168 μB; with 24 divisions you should find a value of 2.218 μB.

Partial DOS

… to be completed

Mulliken analysis

… to be completed

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

Further exercises

  1. For Cr3Si6, try improving the basis and see the effect on the DOS. Remove switches --loc=0 --mto=1 from the commmand line when running blm.

References

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

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

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


Edit This Page