# ASA calculation for CsPbI3 in cubic and pseudocubic structures

In this tutorial,

1. The energy band structure of CsPbI3 is computed in both an ideal cubic structure, and also a pseudocubic structure, using the LDA code lmf.
2. Corresponding calculations are carried out using the ASA code lm, and results compared.
3. A QSGW calculation is carried out, which modifies the band structure and widens the bandgap.

This tutorial is self-contained. It is required to prepare necessary inputs to the Levenberg-Marquardt tutorial which makes small adjustments to the ASA potential parameters to fit the ASA bands to QSGW bands.

You may wish to go through the introductory tutorial for the ASA, if you haven’t already done so.

** Cubic structure **

** Pseudocubic structure **

** QSGW **

** Levenberg-Marquardt **

### Introduction

CsPbI3 is a close cousin of NH3CH3PbI3 (usually called MAPI). They are in the family of perovskites that have attracted a great deal of attention recently because of their considerable potential as efficient, low-cost solar cells.

CsPbI3 has, on average, a cubic structure. However, the PbI3 cage flexes rather dramatically in time, roughly on the scale of a phonon period, distorting the underlying cubic structure and changing the instantaneous band structure. MAPI flexes in a similar way, but CsPbI3 is simpler so this tutorial is written for it. The effect of distortions on the band structure is significant, and it is the subject of this tutorial. It is based on a paper analyzing the effects of nuclear motion in CH$_3$NH$_3$PbI$_3$ and CsPbI$_3$.

This tutorial develops the energy bands for CsPbI3 in its ideal cubic structure, and a pseudocubic structure which approximates the flexing of the system at room temperature.

### Cubic structure

Copy the following box to file site0.cspbi

% site-data vn=3.0 xpos fast io=15 nbas=5 alat=11.82540057 plat= 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
#                        pos
Pb        0.0000000   0.0000000   0.0000000
Cs        0.5000000   0.5000000   0.5000000
I         0.5000000   0.0000000   0.0000000
I         0.0000000   0.5000000   0.0000000
I         0.0000000   0.0000000   0.5000000


and the box below to site0.dcspbi.

% site-data vn=3.0 xpos fast io=15 nbas=5 alat=1.88972688 plat= 6.2887979 0.0 0.0 0.0015427 6.2287686 0.0 0.1397332 -0.0003656 6.3725603
#                        pos
Pb        0.4448840  -0.5000200  -0.0223770
Cs        0.0094815   0.0000333  -0.4695814
I         0.4094870  -0.5000180   0.4751450
I         0.3992710  -0.0000420   0.0236640
I        -0.0676570   0.4999950  -0.0775210


and the box below to syml.cspbi.

#Template file generated by Syml_template.py
31    0.500000   0.500000   0.000000    0.500000   0.500000   0.500000  M to R
41    0.500000   0.500000   0.500000    0.000000   0.000000   0.000000  R to G


#### LDA energy bands, calculated with lmf

Use site0.cspbi to create the input file

blm --pr55 --nfile --rdsite=site0 --nosort --shorten=no --nk=4 --gw --wsitex@fn=site1 --gmax=8 cspbi
cp actrl.cspbi ctrl.fp
cp site1.cspbi site.fp


Command line switches to blm are explained here.

Make the system self-consistent and draw energy bands

lmfa ctrl.fp
cp basp0.fp basp.fp
lmfa ctrl.fp
lmf ctrl.fp >out.fp
cp syml.cspbi syml.fp
lmf ctrl.fp --band~mq~fn=syml


#### LDA-ASA energy bands, calculated with lm

The ASA doesn’t have a Mattheis construction to supply a reasonable starting guess the charge density, as lmf does. On the other hand, because it is the ASA, it is able to store information about the density in a highly compact way1.

Extract l resolved sphere charges generated by lmf
The ASA contracts all the information about the density into energy moments of the sphere charges 1, for each l. This step is not necessary, but it makes use of moments generated by lmf as input to the ASA. Thus it supplies a reasonable starting guess for the density.
blm --pr55 --nfile --rdsite=site0 --nosort --shorten=no --nk=4 --asa~pfree~rmaxs=10 --wsitex@fn=site cspbi
cp actrl.cspbi ctrl.asa
cp site.cspbi site.asa
lmf ctrl.fp --asars2
cp rsta.fp rsta.asa
lm ctrl.asa -vnit=0 --rs=1


Note that there is missing charge. From the end of the output of the last step you should see

 Sum Q=-6.262774  Emad=-0.724614(-3.110233)  Vmtz=-0.279207


We will address the deviation from charge neutrality by adding empty spheres, and adding charges to them to satisfy charge neutrality.

Autogenerate empty spheres
Use blm to remake the ctrl.asa, but this time locating empty spheres.
blm --pr55 --nfile --rdsite=site0 --nosort --shorten=no --nk=4 --asa~pfree~rmaxs=10 --wsitex@fn=site cspbi --findes~rmin=1.5~omaxi=0~resize=20~nspmx=3


Note the --asa switch.

• pfree inserts a line into ctrl.cspbi to keep the continuous principal quantum number Pl from falling too low
• rmaxs=10 sets the range of the tight binding structure constants. The entire family of command line switches is explained here.

Compare site.cspbi and site.asa. They should be the same except for the addition of empty spheres.

Use the –zerq command-line switch to tell lm to add charges to empty spheres in a manner that charge neutrality is satisfied.

cp site.cspbi site.asa
cp actrl.cspbi ctrl.asa


You should see these lines near the beginning of stdout

 asazerq : system charge is now -6.262774
adding charge 0.956974 to s channel, class 1
adding charge 0.423982 to s channel, class 2


lm added electrons to the s channel in the empty spheres, in proportion to their volumes, so that the total extra charge makes the system neutral. You should see this table at the end of stdout

 Class        Qtot       Qbak       Vmad     Vh(Rmax)    V(Rmax)
Pb        -1.324686   0.000000   0.516047  -0.286793  -0.807596
Cs        -1.361714   0.000000   0.991409   0.166127  -0.285513
I         -1.192124   0.000000   0.697869  -0.024630  -0.551172
E          0.956974   0.000000  -0.421284   0.176218  -0.299758
E1         0.423982   0.000000  -0.312501   0.034746  -0.434927


The E sphere has twice the volume of the E1 sphere, so it gets twice the charge.

We now have a good starting point for a self-consistent ASA calculation.

Self-consistent ASA, first attempt
Drive the density to self-consistency, and make the energy bands.
lmstr ctrl.asa
rm -f mixm.asa sv.asa log.asa
lm ctrl.asa -vnit=100 > out.asa


The density converged smoothly, as can be seen by inspecting the first column of sv.asa.

See this page for a description of the program flow in lm, how it manages the charge density and iterates towards self-consistency.

Create the energy bands:

cp syml.cspbi syml.asa
lm ctrl.asa -vnit=100 --band~mq~col=1:9~col2=10:18~col3=19:45~fn=syml


LDA (blue) and first attempt at LDA-ASA energy bands of CsPbI3. Note the spurious band at 3 eV in the ASA.

The energy bands compare favorably to the full LDA calculation, except for a spurious flat state near 3 eV.

This is a Pb 6d state. You can establish this out by using using color weights to distinguish the Cs, Pb, and I characters. Invoke

lm ctrl.asa --pr55 --quit=ham


The table tells you to choose color weights 1:9 for Pb, 10:18 for Cs and 19:45 for I.

lm ctrl.asa --band~mq~col=1:9~col2=10:18~col3=19:45~fn=syml
echo -8,6,5,10 | plbnds -fplot~sh -ef=0 -scl=13.6 -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 -lbl=M,R,G bnds.asa


The flat state at 3 eV is bright red, showing that is of Pb character. You can make further use of the color weights to establish that it is a Pb d state (we won’t do that here).

Self-consistent ASA

The spurious state or “ghost band,” 2 must be removed before proceeding. In this particular instance the problem can be fixed simply by reducing the overlap between the Pb and empty spheres. But this is a special case. For pedagogical reasons, we point out three other ways to solve the problem.

1. Downfold the Pb 6d.
2. Freeze the Pb 6d continuous principal quantum number Pd at a higher value than the one reached through the self-consistency cycle 2 (PMIN=-1 is designed to guard against this problem but it wasn’t sufficient.)
3. Replace the Pb 6d with Pb 5d.

Option 1 : Downfold the Pb 6d.
Set  IDXDN  to 2 in the Pb d channel and iterate to self-consistency
cat actrl.cspbi | sed 's/$$ATOM=Pb.*IDXDN=0,0$$,0/\1,2/'   > ctrl.asa
rm -f mixm.asa sv.asa log.asa
lm ctrl.asa -vnit=100 > out


Draw the energy bands of the downfolded hamiltonian

lm ctrl.asa -vnit=100 --band~mq~col=1:4~col2=5:13~col3=14:40~fn=syml
echo -8,6,5,10 | plbnds -fplot~sh -ef=0 -scl=13.6 -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 -lbl=M,R,G bnds.asa


The spurious state is gone; now the bands compare favorably with the full LDA result.

Note: this method has one drawback: namely that the layer Green’s function code lmpg does not have downfolding capability.

Option 2. Freeze the Pb 6d Pnu at a higher linearization energy
Set  IDMOD  to 1 in the Pb d channel; replace P the Pb d (now at about 6.15) with 6.2, and iterate to self-consistency
sed -i 's/$$2 6.[0-9]*$$ / 2   6.2       /' pb.asa
cat actrl.cspbi | sed 's/$$ATOM=Pb.*$$/\1 IDMOD=0,0,1/'   > ctrl.asa
rm -f mixm.asa sv.asa log.asa
lm ctrl.asa -vnit=100 > out
lm ctrl.asa -vnit=100 --band~mq~col=1:9~col2=10:18~col3=19:27~fn=syml
echo -8,6,5,10 | plbnds -fplot~sh -ef=0 -scl=13.6 -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 -lbl=M,R,G bnds.asa


These bands are practically indistinguishable from Option 1.

Option 3 : Replace the Pb 6d with Pb 5d.
Replace P on the Pb with 5d. So far we have treated the Pb 5d as a core state. This will improve the treateent of the the Pb 5d, but in the region of the Fermi level, Pb has more 6d than 5d character near the Fermi level, so we might expect this method to be a little less accurate.

In setting up the starting conditions the Pb moments are now completely out of whack, since Pb 6d and Pb 5d are very different. Rather than try to guess Pb moments, instead we can remove the Pb atom file, let lm choose some moments taken from the free atom, and scale the charge on that atom to enforce charge neutrality.

cat actrl.cspbi | sed 's/$$ATOM=Pb.*$$/\1 P=0,0,5.93/'   > ctrl.asa
rm pb.asa
lm ctrl.asa -vnit=0 --zerq~qin~ic=1
rm -f mixm.asa sv.asa log.asa
lm ctrl.asa -vnit=100 > out
lm ctrl.asa -vnit=100 --band~mq~col=1:9~col2=10:18~col3=19:45~fn=syml


Agreement now is bit worse, as the Figure shows. Pb d partial waves do have some effect near the Fermi level.

Here we adopt option 1, and also adjust the ES positions to reduce overlaps. This is automatically done with the lmchk utility.

cp site.asa site1.asa
lmchk -vfile=1 ctrl.asa --mino~z --wsite@short@fn=site
cat actrl.cspbi | sed 's/$$ATOM=Pb.*IDXDN=0,0$$,0/\1,2/'   > ctrl.asa


The largest atom-ES overlap, initially 16%, drops to 12.8%. Make the system self-consistent and draw the energy bands.

rm -f mixm.asa sv.asa log.asa
lmstr ctrl.asa
lm ctrl.asa -vnit=100
lm ctrl.asa -vnit=100 --band~mq~col=1:4~col2=5:13~col3=14:40~fn=syml
echo -8,6,5,10 | plbnds -fplot~sh -ef=0 -scl=13.6 -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 -lbl=M,R,G bnds.asa


LDA (blue) and LDA-ASA energy bands of CsPbI3, with the Pb 6d folded down.

Bands no longer have a spurious state at 3 eV, and they are in reasonable agreement with the full LDA result.

### Pseudocubic structure

Compare site0.dcspbi and site0.cspbi. They are basically the same structure, with relatively modest distortions, but:

1. alat and plat are distributed differently: plat(i,i) is about 6.2, corresponding to the (slightly distorted) cube edge in Angstroms. This data was extracted from a VASP POSCAR file.
2. Some sites are translated by approximate integer multiples of lattice vectors (approximate because of distortions)
3. The order of the Iodine sites is permuted.

It is convenient to render the site file of the distorted lattice as near as possible to the undistorted one. Then comparisons are simpler to make. A combination of switches in the blm and lmscell make this easy work.

#### LDA energy bands, pseudocubic structure

blm can eliminate difference (1) with the --scala switch, and most of the differences in lattice translations elminated with the --xshft switch.

blm --pr55 --scala=6.25772999 --nfile --rdsite=site0 --nk=4 --gw --wsitex@fn=site1 --gmax=8 dcspbi --xshftx=-.5,.5,0
cp actrl.dcspbi ctrl.dfp
cp site1.dcspbi site1.dfp


site1.dfp and site.fp are much closer to each other. The superlattice editor lmscell --stack can undo the permutation of Iodine sites with the ~sort option, and it can clean the up the one remaining difference in the translation vector (Cs) with the ~addpos option.

lmscell -vfile=1 ctrl.dfp --stack~sort@targ=1,2,5,4,3~addpos@dp=0,1,1@targ=2~wsitex@short@fn=site


site.dfp and site.fp are as similar as they can be: the remaining differences are all connected with distortions, which is the property of interest.

Use lmf to make self-consistent and draw the energy bands. Here we can re-use the basis set from the undistorted lattice.

cp basp.fp basp.dfp
lmfa ctrl.dfp
lmf ctrl.dfp >out.dfp
cp syml.cspbi syml.dfp
lmf ctrl.dfp --band~mq~fn=syml


Notice that the gap in the pseudocubic structure is 1.58 eV, almost 0.3 eV larger than the gap of the cubic structure. This change is important for the hybrid perovskite devices.

To find the gap, do

grep gap out.fp out.dfp


#### LDA-ASA energy bands, pseudocubic structure

The procedure follows that of the cubic structure, but we can reuse site and ctrl files from prior calculations.

Extract l resolved sphere charges from the lmf calcualtion
cp ctrl.asa ctrl.dasa
cp site.dfp site.dasa
lmf ctrl.dfp --asars2
cp rsta.dfp rsta.dasa
lm ctrl.dasa -vnit=0 --rs=1

Take the empty spheres from the cubic lattice, and use the overlap minimizer to optimize their positions:
  cp ctrl.asa ctrl.dasa
cp site.dfp site1.dasa
tail -11 site.asa >> site1.dasa
sed -i 's/nbas=5/nbas=16/' site1.dasa
lmchk -vfile=1 ctrl.dasa --mino~z --wsite@short@fn=site


This creates a working site file. The volume of the pseudocubic lattice is slightly larger than the cubic one. So sphere volumes no long quite fill the cell volumes, a requirement for the ASA.

This is remedied with sphere resizer, invoked by the command-line argument --sfill

  lmchk dasa --sfill~sclwsr=11~lock=1,1,1~omax=.18,.26,.26~wsrmx=3.3 --pr20


You should see this table

   spec  name        old rmax    new rmax     ratio
1   Pb          3.300000    3.300000    1.000000
2   Cs          3.300000    3.300000    1.000000
3   I           3.300000    3.300000    1.000000
4   E           3.203246    3.239422    1.011294
5   E1          2.441964    2.469543    1.011294


Use your text editor to cut and paste the modified E and E1 sphere radii into ctrl.dasa. Alternatively do it with sed:

  sed -i 's/$$ATOM=E .*R= *$$$$[.0-9]*$$/\13.239422/'  ctrl.dasa
sed -i 's/$$ATOM=E1 .*R= *$$$$[.0-9]*$$/\12.469543/' ctrl.dasa

Self-consistent ASA
Drive the density to self-consistency, and make the energy bands.
lm ctrl.dasa -vnit=0 --zerq~qin~ic=1
lmstr ctrl.dasa
rm -f mixm.dasa sv.dasa log.dasa
lm ctrl.dasa -vscr=0 -vnit=100 --pr31,20 > out.dasa
cp syml.cspbi syml.dasa
lm ctrl.dasa -vnit=100 --band~mq~col=1:4~col2=5:13~col3=14:40~fn=syml


LDA (blue) and LDA-ASA energy bands of CsPbI3, with the Pb 6d folded down. Left: cubic structure. Right: pseudocubic structure.

Distortions widen the gap. The figure shows that the ASA tracks shifts in the bands with displacement fairly well.

### Footnotes and references

1 The ASA does not keep the density itself, but contracts all the information about the density into energy moments of the sphere charges, $Q_{0l}$, $Q_{1l}$, and $Q_{2l}$, of the density of states for each $l$, and the continuous principal quantum number. $P_l$. It can do this because of the special structure of the ASA. When the density is needed to make the potential, it is constructed from the $Q_{0{\ldots}2l}$ and the $P_l$.

If no data is available for a particular atom, the ASA codes lm, lmgf, and lmpg select $Q_{0l}$ from charges in the atom, and sets $Q_{1{\ldots}2l}$ to zero. It reads preset values for the $Q_{0l}$ and $P_l$ from a lookup table. This is will make a crude estimate of the density, but it is usually sufficient to make a starting guess, which can be iterated to the self-consistent values.

Another alternative is to extract the moments from the same species of a different self-consistent calculation, as a starting guess for the self-consistent moments.

A third alternative is to extract the moments $Q_{0l}$, and boundary conditions $P_{l}$, from a full-potential calculation run by lmf. The present tutorial follows this strategy.

2 Normally the continuous principal quantum numbers. $P_l$ are allowed to float to band center of gravity (the $P_l$ at which $Q_{1l}$ vanishes), but when the partial wave is far removed from the Fermi level, this can cause “ghost bands” to appear. One guard against this is to restrict the $P_l$, and not let it fall below the free-electron value. Tag HAM_PMIN is designed for this purpose. Another guard is to freeze the $P_l$ to a fixed value, using SPEC_ATOM_IDMOD. Another way is to downfold the $P_l$. You can tell the ASA codes to downfold a particular state with SPEC_ATOM_IDXDN.