QSGW Tutorial for magnetic bcc Fe


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.

Table of Contents


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
awk '{if ($1 == "-colsy") {sub("-qr","-lt {colg} -qr");sub("lda","green");sub("green","gw");sub("colg","colr");print;sub("gw","lda");sub("colr","colb");print} else {print}}' plot.plbnds >> plot.2bands

fplot -f plot.2bands
open fplot.ps


This tutorial carries out a QSGW calculation for Fe, a ferromagnet with a fairly large moment of 2.2μB.

Quasiparticle Self-Consistent GW (QSGW) is a special form of self-consistency within GW. Its advantages are briefly described in the overview, and also in the summary of the basic QSGW tutorial. Self-consistency is rather necessary in magnetic systems, because the magnetic moment is not reliably described by 1-shot GW, and is somewhat ill-defined. For the most part magnetic moments predicted by QSGW are significantly better than the LDA. Moments in itinerant magnets, however, tend to be overestimated because diagrams with spin fluctuations (absent in GW) tend to reduce the magnetic moment. This is true for both QS_GW_ and LDA.

This tutorial proceeds in a manner similar to the basic QSGW tutorial.

1. Self-consistent LDA calculation for Fe

QSGW requires a starting point. The LDA is a reasonably accurate, and convenient starting point, indeed it is good enough in its own right for many magnetic systems.

Following other LDA tutorials, the input file is built with the blm tool, which requires file init.fe.

Input file setup

Fe is a bcc metal with a lattice constant 5.408 a0 and a magnetic moment of 2.2 μB. The LDA needs some breaking of the symmetry to stabilize a nonzero magnetic moment, so we supply a trial moment of 2μB.

Cut and paste the contents in the box below into init.fe.

# Init file for Fe
  SPCGRP=Im-3m  A=5.408
# ALAT=5.408 PLAT= -0.5 0.5 0.5   0.5 -0.5 0.5   0.5 0.5 -0.5

   ATOM=Fe MMOM=0,0,2     # Trial spin moment on the d channel
   ATOM=Fe X=0 0 0

Supply either the space group (Im-3m) or the bcc lattice vectors.

The input file setup and self-consistent cycle are very similar to the PbTe tutorial; review it first if you are not familiar with the structure of the input file, ctrl.fe, or how to set up an LDA calculation from structural information.

Run blm this way:

$ blm --nit=20 --mag --gw --nk=16 --nkgw=8 --gmax=7.9 fe
$ cp actrl.fe ctrl.fe

The switches do the following:

  • --nit=20    : sets the maximum number of iterations in lmf self-consistency cycle.
    Not required, the default (10 iterations) is about how many are needed are needed to make it self-consistent
  • --mag     : tells blm that you want to do a spin-polarized calculation.
  • --gw      : tells blm to prepare for a GW calculation. The effect is to make the basis set a bit larger than usual and sets the basis Hankel energies deeper than is needed for LDA. This is to make the basis short enough ranged so that the quasiparticlized self-energy Σ0 can be smoothly interpolated between k points.
  • --nk=16     : sets the k mesh. You must supply the mesh.
    We use a rather fine mesh here, because Fe is a transition metal with a high density-of-states near the Fermi level.
  • --nkgw=8     : the k mesh for GW. If you do not supply this mesh, it will use the lmf mesh. The self-energy varies much more smoothly with k than does the kinetic energy; 16 divisions is expensive and overkill.
  • --gmax=7.9 : Plane-wave cutoff. You must supply this number. It is difficult to determine in advance; however you can leave it out at first and run lmfa. You sometimes need to run lmfa twice anyway, since it may find some local orbitals and change the valence-core partitioning. (See the LDA tutorial for PbTe).

Note: blm writes the input file template to actrl.fe, to avoid overwriting a the ctrl.fe, which you may want to preserve.

Free atom density

Generate the free atom density and copy the basis set information it generates to basp.fe. See the PbTe tutorial for further explanation.

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

Additional exercises

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

  2. Test the PBE functional

Edit This Page