# QSGW Tutorial for magnetic bcc Fe

This tutorial describes how to make a QSGW self-energy, and draw energy bands for elemental Fe, a bcc metal.

Several tutorials build on this one, in particular a tutorial for the transverse spin susceptibility and one for dynamical self-energy.

### Preliminaries

Executables blm, lmfa, and lmf are required and are assumed to be in your path; similarly for the QSGW script lmgwsc; and the binaries it requires should be in subdirectory code2.

When a text editor is required, the tutorial uses the nano text editor.

The band-plotting part uses the plbnds tool and to make the figure, the fplot graphics tool.

To view the postscript file, this document assumes you are using the apple-style open command.

It proceeds in a manner similar to the basic QSGW tutorial so you may want to be familiar with that tutorial. It also forms a starting point for other tutorials, for example the calculation of the dynamical self-energy, another for k-resolved DOS, and as well as combined of k-resolved and Mulliken-resolved DOS - going over these tutorials after or during this tutorial may help you better follow both.

### Command summary

1. LDA self-consistency (starting from init.fe)
nano init.fe
blm --nit=20 --nk=16 --gmax=7.9 --mag --nkgw=8 --gw --ctrl=ctrl fe
lmfa fe
cp basp0.fe basp.fe
lmf fe > out.lmf


Make the energy bands and save in bnds.lda:

nano syml.fe
lmf fe --quit=band
lmf fe --band:fn=syml
cp bnds.fe bnds.lda

1. QSGW self-consistency
lmfgwd --jobgw=-1 --sigw --ib=1:9 ctrl.fe
nano GWinput
rm mixm.fe
lmgwsc --mpi=6,6 --sym --metal --tol=1e-5 --getsigp fe > out.lmgwsc
grep more out.lmgwsc


Make the energy bands and save in bnds.gw:

lmf fe --quit=band
lmf fe --band:fn=syml
cp bnds.fe bnds.gw


Create and view a postscript the energy bands using plbnds and fplot.

echo -2,2 /  | plbnds -wscl=1,.8 -fplot -ef=0 -scl=13.6 --nocol -dat=gw -spin2 bnds.gw
echo -2,2 /  | plbnds -wscl=1,.8 -fplot -ef=0 -scl=13.6 --nocol -lbl=H,N,G,P,H,G -spin2 -dat=lda bnds.lda

rm -f plot.2bands
echo "% char0 colr=3,bold=5,clip=1,col=1,.2,.3" >>plot.2bands
echo "% char0 colb=2,bold=4,clip=1,col=.2,.3,1" >>plot.2bands
$cp basp0.fe basp.fe  In this case no local orbitals were generated (inspect basp.fe for any PZ present). lmfa does not need to be run a second time. #### The Fe 4d state It turns out the high-lying Fe 4d state affects the GW potential slightly, enough to affect states near the Fermi level by 0.05-0.1 eV. It matters little in the LDA, unless very high accuracy is demanded. In contrast to the LDA, both occupied and unoccupied states contribute to the GW self-energy. lmfa detects that a GW calculation is intended because of the token autobas_gw, and adjusts basp0.fe accordingly. Inspect this file. You should see the following  Fe RSMH= 1.561 1.561 1.017 1.561 EH= -0.3 -0.3 -0.3 -0.3 RSMH2= 1.561 1.561 1.017 EH2= -1.1 -1.1 -1.1 P= 4.738 4.538 3.892 4.148 5.102 PZ= 0 0 4.4  PZ=0 0 4.5 adds a 4d local orbital to the basis. #### Self-consistency in the LDA Make the LDA density self-consistent: $ lmf fe > out.lmf


This step overlaps the atomic density taken from file atm.fe generated by lmfa, generates an output density, and iterates the until the input and output densities are the same (self-consistency).

lmf should converge in 12 iterations, with the self-consistency density in rst.fe. Inspect save.fe:

h mmom=2.0619219 ehf=-2541.0404312 ehk=-2541.0267069
i mmom=2.1115096 ehf=-2541.0378418 ehk=-2541.0395735
...
c mmom=2.2003954 ehf=-2541.040546 ehk=-2541.0405493


The first line prints what is obtained from the trial Mattheis construction that overlaps free atomic densities. The moment gradually increases from 2.0 (the guessed moment) to 2.20.

At self-consistency, the Harris-Foulkes and Kohn-Sham total energies are nearly identical. (If they are not, something is wrong with the calculation!) Note that the magnetic moment (2.20 μB), is very close to the experimental value.

#### LDA Energy bands

In this section we compute the energy bands, which we will compare later with the QSGW result. For this tutorial use symmetry lines given in the box below. Copy its contents into syml.fe.

51   0   0   0      1   0   0             Gamma to H   (Delta)
51   1   0   0     .5  .5   0             H to N       (G)
51   0  .5  .5     .5  .5  .5             N to P       (D)
51  .5  .5  .5      0   0   0             P to Gamma   (Lambda)
0    0 0 0  0 0 0


To get the Fermi level corresponding to the density rst.fe without overwriting it, do:

$lmf fe --quit=band  Note that the Fermi level is -0.000584, very close to the last iteration in the self-consistency cycle. You can find this doing $ grep Fermi out.lmf


The Fermi level is saved in file wkp.fe. Make the energy bands with

$lmf fe --band:fn=syml$ cp bnds.fe bnds.lda


The latter command renames bnds.fe for future use.

### 2. QSGW calculation for Fe

#### Setup: the GWinput file

The GW codes require a separate GWinput file. Create a template with:

$lmfgwd --jobgw=-1 --sigw --ib=1:9 ctrl.fe  Switch --sigw makes small changes to file.ext in preparation to make the dynamical self-energy. It isn’t necessary for the QSGW calculation of this tutorial, but it is needed for the follow-on tutorial that makes dynamical self-energy and interacting spectral functions. Switch --ib=1:9 is used in 1-shot GW mode. It isn’t relevant for this tutorial, but sets list of QP states for which the dynamical self-energy is made in the dynamical self-energy tutorial. We are particularly interested in Fermi liquid properties, involving states near the Fermi surface. The raw GWinput template will generate a reasonable self-energy, but the 3s and 3p core levels affect the Fermi surface enough that we need to treat them at a little better level of approximation than the template gives you. Edit GWinput: $ nano GWinput


In the core part of the product basis you should see these lines:

  atom   l    n  occ unocc   ForX0 ForSxc :CoreState(1=yes, 0=no)
1    0    1    0    0      0    0    ! 1S *
1    0    2    0    0      0    0    ! 2S
1    0    3    0    0      0    0    ! 3S
1    1    1    0    0      0    0    ! 2P
1    1    2    1    0      0    0    ! 3P


With switches as given, no core level participates in the product basis, polarizability or self energy, except that the 3p is included in the product basis (occ=1). For accurate description of the Fermi surface the 3s also needs to be included; moreover, both 3s and 3p need to be included in the polarization (ForX0 and self-energy (ForSxc).

Caution: core levels are calculated indpendently of the valence levels, and there is a slight residual nonorthogonality that can cause problems if the core levels are too shallow. This can be a serious issue and a safer approach is to include these levels in the valence through local orbitals, though in this case the levels are deep enough that the present treatment is adequate.

Change the 3s and 3p lines as follows:

  atom   l    n  occ unocc   ForX0 ForSxc :CoreState(1=yes, 0=no)
...
1    0    3    1    0      1    1    ! 3S
...
1    1    2    1    0      1    1    ! 3P


#### QSGW self-consistency

Run the QSGW script as follows:

$rm mixm.fe$ lmgwsc --sym --metal --tol=1e-5 --getsigp fe > out.lmgwsc


or faster

$rm mixm.fe$ lmgwsc --mpi=6,6 --sym --metal --tol=1e-5 --getsigp fe > out.lmgwsc


The QSGW cycle should complete in a couple of hours, after 9 iterations. When it is finished, do

$grep more out.lmgwsc  You should see the following:  lmgwsc : completed iteration 0 of 999 more=T Mon Oct 31 09:05:33 GMT 2016 elapsed wall time 15.4m (0.3h) phpdl1 lmgwsc : iter 1 of 999 RMS change in sigma = 5.57E-03 Tolerance = 1e-5 more=T Mon Oct 31 09:22:09 GMT 2016 elapsed wall time 32.0m (0.5h) phpdl1 lmgwsc : iter 2 of 999 RMS change in sigma = 2.32E-03 Tolerance = 1e-5 more=T Mon Oct 31 09:38:12 GMT 2016 elapsed wall time 48.0m (0.8h) phpdl1 lmgwsc : iter 3 of 999 RMS change in sigma = 4.91E-04 Tolerance = 1e-5 more=T Mon Oct 31 09:54:15 GMT 2016 elapsed wall time 64.0m (1.1h) phpdl1 lmgwsc : iter 4 of 999 RMS change in sigma = 9.56E-04 Tolerance = 1e-5 more=T Mon Oct 31 10:11:04 GMT 2016 elapsed wall time 80.9m (1.3h) phpdl1 lmgwsc : iter 5 of 999 RMS change in sigma = 6.64E-04 Tolerance = 1e-5 more=T Mon Oct 31 10:28:16 GMT 2016 elapsed wall time 98.1m (1.6h) phpdl1 lmgwsc : iter 6 of 999 RMS change in sigma = 8.96E-05 Tolerance = 1e-5 more=T Mon Oct 31 10:45:00 GMT 2016 elapsed wall time 114.8m (1.9h) phpdl1 lmgwsc : iter 7 of 999 RMS change in sigma = 7.11E-05 Tolerance = 1e-5 more=T Mon Oct 31 11:01:18 GMT 2016 elapsed wall time 131.1m (2.2h) phpdl1 lmgwsc : iter 8 of 999 RMS change in sigma = 9.82E-06 Tolerance = 1e-5 more=F Mon Oct 31 11:18:00 GMT 2016 elapsed wall time 147.8m (2.5h) phpdl1  The self-consistent cycle ends when the RMS change in Σ0 falls below the specified tolerance (–tol=1e-5). #### QSGW Energy bands Σ0 is an effective non interacting potential with one-particle levels, similar to Hartree Fock or the LDA. The band structure can be drawn in the same way as in the LDA. First, find the Fermi level: $ lmf fe --quit=band


You should get a table like this one:

 BZINTS: Fermi energy:      0.056498;   8.000000 electrons;  D(Ef):   16.413
Sum occ. bands:   -0.6890437  incl. Bloechl correction:   -0.000979
Mag. moment:       2.269378


The Fermi level is higher than in the LDA value (-0.000599); it suggests that the work function would be somewhat smaller. The magnetic moment (2.27 μB) comes out slightly larger as well. A better converged QSGW calculation (calculating Σ0 with more k divisions) reduces this value to about (2.2 μB), very similar to what the LDA gets.

Generate the band structure, and copy the file to bnds.gw

$lmf fe --band:fn=syml$ cp bnds.fe bnds.gw


### Compare QSGW and LDA energy bands

At this point the LDA (bnds.lda) and QSGW (bnds.gw) energy bands should in your working directory, containing bands along four symmetry lines (Γ-H,  H-N,  N-P,  and P-Γ). For a brief description of the structure of these files, see here.

Use the plbnds tool render both data sets into files (bnd[1234].lda) for the four panels of LDA bands, and (bnd[1234].gw) for the four panels of QSGW bands. This tutorial will be concerned with the bands near the Fermi level. Restrict the range to EF±2 eV, and focus on the minority spin bands:

echo -2,2 /  | plbnds -wscl=1,.8 -fplot -ef=0 -scl=13.6 --nocol -dat=gw -spin2 bnds.gw
echo -2,2 /  | plbnds -wscl=1,.8 -fplot -ef=0 -scl=13.6 --nocol -lbl=H,N,G,P,H,G -spin2 -dat=lda bnds.lda


Each command generates a script (plot.plbnds), which fplot will turn into postscript files. The script draws four frames, one for each symmety line.

Here we adapt the second plot.plbnds, and create a new script plot.2bands that combines the bands from the two calculations in one figure. To distinguish bands the GW and LDA bands will be drawn respectively in red dots and blue dashed lines. fplot selects line type and colors with the -lt instruction.

To make script use the same color in each frame, it is convenient to make use of the file preprocessor’s ability to assign and use character variables. Execture the instructions in the box below. It creates a new file plot.2bands with character variables colr and colb. They contain strings that modify line types, notably the color (see here for a quick reference).

$rm -f plot.2bands$ echo "% char0 colr=3,bold=5,clip=1,col=1,.2,.3" >>plot.2bands
$echo "% char0 colb=2,bold=4,clip=1,col=.2,.3,1" >>plot.2bands  Next, append plot.plbnds to plot.2bands, adapt it to draw two kinds of bands in each frame, each with its own color. This could be accomplished with a text editor, but it is convenient to accomplish it with an awk command: $ awk '{if ($1 == "-colsy") {sub("-qr","-lt {colg} -qr");sub("dat","green");sub("green","gw");sub("colg","colr");print;sub("gw","lda");sub("colr","colb");print} else {print}}' plot.plbnds >> plot.2bands  Compare plot.plbnds and plot.2bands. $ diff plot.plbnds plot.2bands


Character variables are declared at the top. Each line in the script that draws bands has been split into two lines, one for bndn.lda with line type modifier {colb} and another for bndn.gw with line type modifier {colr}.

After preprocessing the script will contain instructions explained in the fplot manual, e.g. -colsy 2:6,   lt 3,…,   and -qr. To see how the script appears after preprocessing, do

$rdfile plot.2bands  $ fplot -f plot.2bands
\$ open fplot.ps


QSGW makes substantial shifts relative to the LDA. The QSGW bands generated in this tutorial disagree somewhat with experiment, because the QSGW potential isn’t quite converged. When well converged, agreement with the available experimental data in the Fermi liquid regime is excellent, though a considerable discrepancy with LDA remains.