Inverse Participation Ratio

This tutorial demonstrates how to calculate the inverse participation ratio (IPR) in Cr. The IPR is a measure of the localization of an eigenfunction and is defined as $\Omega / \int d\mathbf{r} \lvert\psi^{n\mathbf{k}}(\mathbf{r})\rvert^4$, for normalized eigenstate $\psi^{n\mathbf{k}}$, i.e. $\int d\mathbf{r} \lvert\psi^{n\mathbf{k}}(\mathbf{r})\rvert^2 = 1$. Here $\Omega$ is the unit cell volume.

Preliminaries

This tutorial uses several Questaal executables, e.g. blm, lmfa, lmf, pfit, plbnds, and fplot. They are assumed to be in your path.

This tutorial is aimed at showing how to calculate the IPR. We will not bother with self-consistency, as it doesn’t affect the explanation. For a potential we will just use a Mattheis construction (overlapped atomic densities).

Tutorial

Cr is an antiferromagnet, with an incommensurate spin-density wave superposed on it. We will ignore the SDW part and treat Cr as a simple antiferromagnet.

Copy the box below to init.cr.

# Init file for Cr (Wein2K equilibrium volume 11.773A, http://molmod.ugent.be/sites/default/files/deltadftcodes/supplmat/SupplMat-WIEN2k.pdf)
LATTICE
    ALAT=5.41632 PLAT= 1 0 0   0 1 0   0 0 1
SPEC
    ATOM=Crup  MMOM=0,0,1
    ATOM=Crdn  MMOM=0,0,-1
SITE
    ATOM=Crup  X=0 0 0
    ATOM=Crdn  X=.5 .5 .5

Create a ctrl file and a symmetry lines file

blm init.cr --upcase  --simple --wsite --nk,met --ctrl=ctrl
lmchk ctrl.cr --syml~n=71~mq~lblq:G=0,0,0,M=1/2,1/2,0,X=0,1/2,0~lbl=GXM

Get the free atom density and Fermi level

lmfa ctrl.cr
lmf cr --quit=rho --rs=0

Generate the energy bands (bnds.h5)

lmf cr '--band~h5~scol:(z==24):l=2~scol2:(z==24):l=0~fn=syml' --rs=0
cp bnds.h5 bnds0.h5

Note that the instruction specifies two Mulliken projections of the bands: the first selects out Cr d character, the second Cr s character. These we will use as color weights when drawing the bands later. The second instruction preserves bnds.h5 because we want this information, but the file itself will be (unintentionally) overwritten when making the IPR.

Make the inverse participation ratio on for the k points in syml.cr.
Note: lmf can be replaced by an MPI instruction, e.g. mpirun -n 16 lmf. In this case the code execute quickly and you won’t have to wait long even with a scalar calculation.

lmf cr '--band~h5~fn=syml' --pwz~all~rps=2/3~norm~pprt~evrnge=5,20~ipar > out.ipr

writes the inverse participation ratio to file evec.h5.

Note: By default --pwz works with envelope function part of the eigenfunctions only, omitting the difference between the augmented basis functions and the envelope function inside any augmentation sphere. That the envelope function by itself is not normalized can be mitigated by normalizing the quasi-eigenfunction (i.e. contraction of the eigenvector with the unaugmented envelope functions, rather than with the true basis functions). However the shape of the quasi-eigenfunction inside the augmentation spheres is different from the true eigenfunction. This is especially problematic if the true partial waves are sharply peaked, as they are for d orbitals in transition metals. Since the inverse participation ratio is calculated for the purpose of measuring “peakiness” of the eigenfunction, the result is unreliable.

To make a better estimate, --pwz has an option to augment the bare envelope functions with pseudo partial waves. True partial waves are replaced by waves constructed from a simple pseudopotential: a constant inside a radius r_ps, where r_ps is selected to make the value and slope of the pseudo-wave match the true one at r_ps. (The pseudo partial wave inside r_ps is a spherical Bessel function.) One consequence of matching value and slope is that the norm of the pseudo wave should equal that of the true one. Pseudo partial waves are much smoother than true ones, and they are much more easily representable in terms of plane waves. With this form, an (approximate, but reasonable) plane wave representation of the eigenfunction is possible, which is useful in several contexts, including the calculation of $\int d\mathbf{r} \lvert\psi^{n\mathbf{k}}(\mathbf{r})\rvert^4$ needed for the inverse participation ratio of this tutorial, and also Compton scattering experiments, as well as plotting of charge densities.

To cause the --pwz switch to include (pseudo)augmentation part, add an option rps=#. It will construct a pseudopotential that creates a nodeless partial wave for radii between 0 and #. (The true partial wave is kept for radii larger than #).

For the interpretation of the other modifiers to --pwz~all~rps=2/3~norm~pprt~evrnge=5,20~ipar, see the command line options page.

See what evec.h5 contains (h5ls evec.h5). Some of the key elements should be the eigenvalues, inverse participation ratio, and a list of 72 k points

eval                     Dataset {72, 1, 56}
ipar                     Dataset {72, 1, 56}
qirr                     Dataset {72, 3}

Inspect out.ipr. You should see tables like this:

 ib    eval     IPR
  5  -2.7276  14.4311  |   8  -0.0344   2.1115  |  11   0.2724   6.8014  |  14   0.3261   6.8460  |  17   0.5006  10.2828  |  20   1.4727   4.6536
  6  -2.7276  14.4311  |   9   0.0988   4.2044  |  12   0.2724   6.8014  |  15   0.3631   9.1850  |  18   0.5561   6.2882  |
  7  -0.1621   1.5977  |  10   0.2261   5.6411  |  13   0.3261   6.8460  |  16   0.4396  10.5194  |  19   1.4727   4.6536  |

The inverse participation ratio is unity for a constant eigenfunction, and becomes larger the more sharply peaked (and therefore more localized) it is. The IPR is very large for states 5 and 6 (bands 1-6 are high-lying core 3p orbitals).

You can use the inverse participation ratio as you see fit. Here we will merge it with the Mulliken weights generated with the bands file, bnds.h5.

To do this, use a feature of the plbnds editor that will merge ipar residing in evec.h5, into the bnds file an additional color weight. Restore the bnds file from before, and merge the IPR as a color weight:

cp bnds0.h5 bnds.h5
plbnds --h5 --edit~ipr@neg@scl=2~quit bnds.h5

The raw IPR isn’t suitable as a color weight, as it ranges between 1 and ∞. For that reason The inverse of IPR is merged, scaled by a factor scl. Without further modification, the result would range between 1 for a constant eigenstate to 0 for an infinitely localized one. For plotting purposes the picture is clearer if we assign 0 for a constant eigenstate and 1 for an infinitely localized one. The @neg option makes this change, so with the switches given, the third color weight becomes  1−2/IPR, since scl=2.

In the above, plbnds appends a third color weight writes the output to bnds.ipr.h5. The command below displays the weights

mcx -r:h5:id=colwt:c=0,0,1,0:o=0,0,2,0 bnds.ipr.h5 -px

Weights are nonzero only for bands 5 through 20, and the value varies between 0 and 1.

Make a picture of the bands with

plbnds --h5 -range=-10,6 -fplot~scl=.3,.5~ts=.5 -ef=0 -scl=13.6 -lt=1,bold=3,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=1,1,1 bnds.ipr.h5

This makes fplot.ps, which you can view with your favorite reader. It should look like the figure shown.

Green depicts s character, black p character and red d character. The more bleached a band is, the larger the IPR. You can see that the states around 2 eV have the highest IPR, as expected since d orbitals are atomic-like and confined, while s and p orbitals are extended.