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. For the interpretation of --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.