ASA calculation for CsPbI3 in cubic and pseudocubic structures

In this tutorial,

1. The energy band structure of CsPbI₃ 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.

Introduction

CsPbI₃ is a close cousin of NH₃CH₃PbI₃ (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.

CsPbI₃ has, on average, a cubic structure. However, the PbI₃ 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 CsPbI₃ 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 CsPbI₃ 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

LDA energy bands, calculated with lmf

Use site0.cspbi to create the input file

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

Command line switches to blm are explained here.

Set up symmetry lines for energy bands, file syml.fp. Command adpated from the symmetry lines page

lmchk ctrl.fp --syml~n=61~mq~wmq~lblq:G=0,0,0,R=1/2,1/2,1/2,X=0,0,1/2,M=0,1/2,1/2~lbl=MRG

Make the system self-consistent and compute energy bands

lmfa ctrl.fp
lmf ctrl.fp >out.fp
lmf ctrl.fp --band~mq~scol@atom=pb~scol2@atom=cs~scol3@atom=i~fn=syml

To create a postscript figure of the band structure, do

plbnds -range=-10,8 -fplot -ef=0 -scl=13.6 -dat=green -lt=1,col=0,0,0,colw=1,0,0,colw2=0,1,0,colw3=0,0,1 -lbl=M,R,G  bnds.fp
fplot -f plot.plbnds
gs fplot.ps

[Any postscript viewer can be used in place of the Ghostscript gs]

The states near the valence band maximum are iodine-like (blue), while the states at the conduction band minimum are Pb like (purple = red+blue).

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

Extract l resolved sphere charges generated by lmf

The ASA contracts all the information about the density into energy moments of the sphere charges ¹, 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 --brief --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=-5.836727  Emad=-0.559668(-2.648708)  Vmtz=-0.315457

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

Autogenerate empty spheres

Use blm to remake the ctrl.asa, but this time locating empty spheres.

blm --brief --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 P_l from falling too low (below the free electron value)
• 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. Next do a little cleanup:

cp site.cspbi site.asa
cp actrl.cspbi ctrl.asa
lmchk asa --mino~z --wsitex~short~fn=site1
cp site1.asa site.asa

The last two steps are not essential, but they move the E sites around to reduce sphere overlaps.

Now we are ready to create a starting potential Use the --zerq command-line switch to tell lm to add charges to empty spheres in a manner that charge neutrality is satisfied.

lm ctrl.asa -vnit=0 --zerq~es~add~qin --pr45

You should see these lines near the beginning of stdout

 asazerq : system charge is now -5.836727
           adding charge 0.962397 to s channel, class 1
           adding charge 0.368692 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.324811   0.000000   0.434942  -0.337686  -0.829351
 Cs        -0.935569   0.000000   0.835369   0.300759  -0.129149
 I         -1.192116   0.000000   0.658431  -0.036809  -0.527685
 E          0.962397   0.000000  -0.404977   0.204973  -0.278103
 E1         0.368692   0.000000  -0.252756   0.068981  -0.407160
 Sum Q=0.000000  Emad=0.432771(-3.399443)  Vmtz=-0.453235

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 in its iteration to self-consistency, as can be seen by looking for the rms change DQ in the moments from the out.asa (do grep DQ out.asa) or by inspecting the fourth 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.fp syml.asa
lm ctrl.asa --band~mq~scol@atom=pb~scol2@atom=cs~scol3@atom=i~fn=syml

LDA (blue) and first attempt at LDA-ASA energy bands of CsPbI₃.

The energy bands compare favorably to the full LDA calculation, except for a flat state in the FP case at -8 eV. It can be seen from color weights in the FP bands that this state is Cs derived, and further analysis shows it comes from the semicore Cs 5p orbital. This is included in the valence using local orbitals in lmf, but there is no facility for it in the ASA.

To preserve a complete setup, append the self-consistent moments to ctrl.asa

lmctl asa
cat log.asa >> ctrl.asa

To confirm that ctrl.asa and site.asa make a self-contained setup, copy these files to a clean directory and do:

lmstr asa
lm -vnit=0 asa
lm asa --quit=rho

You should see that the RMS DQ is small, less than 10⁻⁵.

Note: with a modest change in how the calculation is set up, you may find a spurious flat band near 3 eV; this is a spurious Pb 6d (“ghost”) state.

Self-consistent ASA

The spurious state or “ghost band,” ² 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 P_d at a higher value than the one reached through the self-consistency cycle ² (PMIN=-1 is designed to guard against this problem but it wasn’t sufficient.)
3. Replace the Pb 6d with Pb 5d.

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

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

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

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 --brief --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.fp syml.dfp
lmf ctrl.dfp --band~mq~scol@atom=pb~scol2@atom=cs~scol3@atom=i~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 calculation

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

Add empty spheres

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

You should see this table

   spec  name        old rmax    new rmax     ratio  lock
     1   Pb          3.429366    3.429366    1.000000 *
     2   Cs          3.500000    3.500000    1.000000 *
     3   I           3.429366    3.429366    1.000000 *
     4   E           3.155659    3.195825    1.012728
     5   E1          2.291887    2.321059    1.012728
   

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.bak 's/\(ATOM=E .*R= *\)\([.0-9]*\)/\13.195825/'  ctrl.dasa
  sed -i.bak 's/\(ATOM=E1 .*R= *\)\([.0-9]*\)/\12.321059/' 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

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

Next steps: compute the partial DOS in the ASA; see this tutorial

Footnotes and references

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

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