# QSGW Tutorial for magnetic bcc Fe

### Purpose

This tutorial begins with with a self-consistent LDA calculation for Fe, a bcc magnetic metal, and follows it with a QSGW calculation.

It proceeds in a manner similar to the basic QSGW 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.

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

### 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 fe
cp actrl.fe ctrl.fe
lmfa fe
cat basp0.fe | sed -e 's/$$Fe.*$$/\1 PZ=0 0 4.5'/ > 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
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 We could proceed directly to making a self-consistent LDA density with the setup as given. But when later proceeding to the GW part, it turns out the 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. (In contrast to the LDA, both occupied and unoccupied states contribute to the GW self-energy.) Anticipating that the GW code will need these states for a good description of the Fermi surface, edit basp.fe: $ nano basp.fe


You should see a line beginning with Fe. Append PZ=0 0 4.5 to the line so that it looks similar to

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


This adds a 4d state (local orbital) to the basis.

Note: To make this addition without a text editor, do

$cat basp0.fe | sed -e 's/$$Fe.*$$/\1 PZ=0 0 4.5'/ > basp.fe  #### 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.0619193 ehf=-2541.0404251 ehk=-2541.0266951
i mmom=2.173571 ehf=-2541.0406411 ehk=-2541.0395501
...
c mmom=2.2003534 ehf=-2541.0405422 ehk=-2541.0405455


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.000599, 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

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

1. Tetrahedron vs sampling vs fermi function, and METAL modes

2. Test the PBE functional