Making the dynamical GW self energy

Table of Contents


Preliminaries


This tutorial assumes you have completed a QSGW calculation for Fe, following this tutorial, which requires that the GW script lmgwsc is in your path, along with the executables it requires. In addition it requires spectral and lmfgws.

This tutorial assumes the your build dirrectory is ~/lm, and the executables are in ~/bin and ~/bin/code2.

Command summary


Repeat the steps for LDA self-consistency and QSGW self-consistency in the [Fe tutorial[/tutorial/gw/qsgw_fe/]; see Command summary.

If you have already done so without removing any files, you can skip those steps.

If on the other hand you have retained the QSGW self-energy Σ0 (file sigm), make sure your charge density is self-consistent.

lmf fe > out.lmf

If you also retained the density restart file rst.fe this step is not necessary. Make all the inputs (e.g. screened coulomb interaction) up to the self-energy step:

lmgwsc --wt --code2 --sym --metal --tol=1e-5 --getsigp --stop=sig fe

This will give you everything you need to make Σ(k,ω)\Sigma(\mathbf{k},\omega).

Make the spectral function *files SEComg.UP (and SEComg.DN)

export OMP_NUM_THREADS=8
export MPI_NUM_THREADS=8
~/bin/code2/hsfp0_om --job=4 > out.hsfp0

Postprocessing step translating SEComg.{UP,DN} into lmfgws-readable form

spectral --ws --nw=1
ln -s se se.fe
... to be finished

Introduction

This tutorial starts after a QSGW calculation for Fe has been completed, in this tutorial.

The QSGW static self-energy was made with the following command:

$ lmgwsc --wt --code2 --sym --metal --tol=1e-5 --getsigp fe

Note: until that tutorial is written, perform the setup as follows (where ~/lm is your Questaal source directory)

~/lm/gwd/test/test.gwd --mpi=#,# fe 4

This tutorial will do the following:

Theory

Z factor renormalization

Begin with a noninteracting Green’s function G0, defined through an hermitian, energy-independent exchange-correlation potential Vjxc(k).   j refers to a particular QP state (pole of G0). There is also an interacting Green’s function, G.

The contribution to G0 from QP state j is

G0j(k,ω)=1ωωj(k)G_0^j(k,\omega ) = \frac{1}{\omega - \omega^j(k)}

where ωj(k)\omega^j(k) is the pole of G0G_0.

Write the contribution to G from QP state j as

Gj(k,ω)=1ωωjΣ(k,ω)+Vxcj(k)G^j(k,\omega) = \frac{1}{\omega - \omega^j - \Sigma (k,\omega ) + V^j_{xc}(k)}

Note that this equation is only true if Σ\Sigma is diagonal in the basis of noninteracting eigenstates. We will ignore the nondiagonal elements of Σ(k,ω)\Sigma(k,\omega). Note that if Vxcj is constructed by QSGW, this is a very good approximation, since ReΣ(k,ω)=Vxcj(k){\mathrm{Re}\Sigma (k,\omega ){=}V^j_{xc}(k)} at ω=ωj(k)\omega{=}\omega^j(k). Approximate G by its coherent part:

Gj,coh(k,ω)=1ωωjReΣ(k,ωj)+Vxcj(k)(ωωj)(11/Zj)iImΣ(k,ω)G^{j,\mathrm{coh}}(k,\omega) = \frac{1}{\omega - \omega^j - \mathrm{Re} \Sigma (k,\omega^j) + V^j_{xc}(k) - (\omega - \omega^j)(1 - 1/Z^j) - i\mathrm{Im} \Sigma (k,\omega )}

where

11/Zj=Σ(k,ω)/ωωj.1 - 1/{Z^j} = \left. {\partial \Sigma (k,\omega )/\partial \omega } \right|_{\omega ^j } .

defines the Z factor. The dependence of ωj and Zj on k is suppressed.

Define the QP peak as the value of ω where the real part of the denominator vanishes.

(ωωj)/Zj=ReΣ(k,ωj)Vxcj(k)({\omega^*} - \omega^j)/Z^j = \mathrm{Re} \Sigma(k,\omega^j) - V^j_{xc}(k)

and so

ω=ωj+Zj(ReΣ(k,ωj)Vxcj(k)){\omega^*} = \omega^j + Z^j\left( {\mathrm{Re} \Sigma (k,\omega^j) - V^j_{xc}(k)} \right)

Note that in the QSGW case, the second term on the r.h.s. vanishes by construction: the noninteracting QP peak corresponds to the (broadened) pole of G.

The group velocity is /dk. For the interacting case it reads

dωdk=dωjdk+ddkZj(ReΣ(k,ωj)Vxcj(k))\frac{d\omega^*}{dk} = \frac{d\omega ^j }{dk} + \frac{d}{dk}Z^j \left( {\text{Re}\Sigma(k,\omega^j) - V_{xc}^j (k)} \right)

Use the ratio of noninteracting and interacting group velocities as a definition of the ratio of inverse masses. From the chain rule

m0mdωdk/dωjdk=1+Zj(ωReΣ(k,ωj)ωjdωjdk+kReΣ(k,ωj)ωjkVxcj(k))+(dZjdk)(ReΣ(k,ωj)Vxcj(k))\frac{m_0}{m^*} \equiv \frac{d\omega^*}{dk}/\frac{d\omega^j}{dk} = 1 + Z^j \left( \frac{\partial}{\partial\omega} \left. \text{Re}\Sigma (k,\omega ^j ) \right|_{\omega^j} \frac{d\omega ^j }{dk} + \frac{\partial}{\partial k}\left. \text{Re}\Sigma(k,\omega^j) \right|_{\omega ^j } - \frac{\partial }{\partial k}V_{xc}^j (k) \right) + \left(\frac{dZ^j}{dk}\right) \left(\text{Re}\Sigma(k,\omega^j) - V_{xc}^j (k) \right)

Ignore the dependence of Zj on k. Write j/dk as v0jv_0^j, and use the definition of Z to get

m0m=1+1v0jZj((11/Zj)v0j+kReΣ(k,ωj)ωjkVxcj(k))\frac{m_0}{m^*} = 1 + \frac{1}{v_0^j}Z^j \left( {\left( {1 - 1/{Z^j}} \right) v_0^j + \frac{\partial}{\partial k} \left.\text{Re}\Sigma(k,\omega^j) \right|_{\omega^j} - \frac{\partial }{\partial k}V_{xc}^j(k)} \right)

So

m0m=Zj+Zjv0j(kReΣ(k,ωj)ωjkVxcj(k))\frac{m_0}{m^*} = Z^j + \frac{Z^j}{v_0^j }\left( {\frac{\partial}{\partial k}\left. \text{Re}\Sigma (k,\omega^j) \right|_{\omega ^j } - \frac{\partial }{\partial k}V_{xc}^j (k)} \right)

In the QSGW case the quantity in parenthesis vanishes. Thus QSGW there is no “mass renormalization” from the ω-dependent self-energy, Σ(ω).

Coherent part of the spectral function

Write Gj,coh(k,ω)G^{j,\mathrm{coh}}(k,\omega) as

Gj,coh(k,ω)=[(ωωj)Zj1ReΣ(k,ωj)+Vxc(k)iImΣ(k,ω)]1G^{j,\mathrm{coh}}(k,\omega) = \left[(\omega - \omega^j){Z^j}^{-1} - \mathrm{Re} \Sigma (k,\omega^j) + {V_{xc}}(k) - i\mathrm{Im} \Sigma (k,\omega )\right]^{-1}

Rewrite as

Gj,coh(k,ω)=ZjωωiZImΣ(k,ω)=Zjωω+iZImΣ(k,ω)(ωω)2+(ZjImΣ(k,ω))2G^{j,\mathrm{coh}}(k,\omega) = \frac{Z^j}{\omega - \omega^* - iZ\mathrm{Im} \Sigma (k,\omega )} = Z^j\frac{\omega - \omega^* + iZ\mathrm{Im} \Sigma (k,\omega )}{(\omega - {\omega^*})^2 + (Z^j\mathrm{Im} \Sigma (k,\omega ))^2}

Using the standard definition of the spectral function, e.g. Hedin 10.9:

A(ω)=π1ImG(ω)A(\omega ) = {\pi ^{ - 1}}\left| {\mathrm{Im} G(\omega )} \right|

the approximate spectral function is

which shows that the spectral weight of the coherent part is reduced by Z.

Akj,coh(ω)=ZjπZjImΣ(k,ω)(ωω)2+(ZjImΣ(k,ω))2A_k^{j,\mathrm{coh}}(\omega ) = \frac{Z^j}{\pi}\frac{Z^j\mathrm{Im} \Sigma (k,\omega )}{(\omega - {\omega^*})^2 + (Z^j\mathrm{Im} \Sigma (k,\omega ))^2}

Simulation of Photoemission

(needs cleaning up)

Energy conservation : requires (see Marder, p735, Eq. 23.58)

ω=Ekin+φsEb\hbar\omega=E_{kin}+{\varphi_s}-E_b

where Eb is the binding energy and Ekin+φsE_{kin}+{\varphi_s} is the energy of the electron after being ejected. (Marder defines EbE_{b} with the opposite sign, making it positive).

Momentum conservation : The final wave vector kf of the ejected electron must be equal to its initial wave vector, apart from shortening by a reciprocal lattice vector to keep kf in the first Brillouin zone.

Let EkinE_{kin} be the energy on exiting the crystal, φs\varphi_s the work function and EbE_b and V0V_0 are called the electron binding energy and “inner potential.”

Then

22m(k2+k2)=Ekin+V0, where Ekin=ωφs+Eb\frac{\hbar^2}{2m}(k_\parallel^2 + k_\bot^2) = E_{kin} + V_0, \text{ where } E_{kin} = \hbar \omega - \varphi_s + E_b

The total momentum inside the crystal, k+k\mathbf{k}_\parallel{+}\mathbf{k}_\bot, is linked to the kinetic energy measured outside the crystal through Eq.(1). The kinetic energy is linked to the binding energy through the equation Ekin=ωEbφa{E_{kin}}=\hbar\omega-{E_b}-{\varphi_a} where φa{\varphi_a} is the work function of the analyzer. Usually φa=φs{\varphi_a}{=}{\varphi_s}. The Fermi level is defined such that Eb=0E_b{=}0. The inner potential is defined by scanning the range of photon energy under the constraint of normal emission: then the Γ\Gamma-point can be identified and by using Eq.(1), and the inner potential experimentally determined.

The momentum of the particle in free space is

2k022m=Ekin\frac{\hbar ^2 k_0^2 }{2m} = E_{kin}

Resolve kf\mathbf{k}_f into components parallel and perpendicular to the surface

kf=k+k\mathbf{k}_f = \mathbf{k}_\parallel + \mathbf{k}_\bot

After passing through the surface, kf\mathbf{k}_f is modified to kf\overline{\mathbf{k}}_f; this is what is actually measured.

The conservation condition requires

k02=k¯2+k¯2k_0^2 = \bar k_\parallel^2 + \bar k_\bot^2

k\mathbf{k}_\parallel is conserved on passing through the surface; thus k¯=k\bar k_\parallel{=} k_\parallel. k\mathbf{k}_\bot is not conserved; therefore

k¯=k02k2\bar k_\bot = \sqrt{k_0^2{-}k_\parallel^2}

The wave number shift is then

Δk=(kk)k^\Delta{\mathbf{k}} = (\overline{k}_\bot-{k}_\bot)\hat{\mathbf{k}}_\bot

and the crystal momentum actually being probed by the experiment is

kf=kfΔk{\mathbf{k}}_f = \overline{\mathbf{k}}_f - \Delta{\mathbf{k}}

Make the GW dynamical self-energy

The 1-shot GW self-energy maker, hsfp0, has a mode (--job=4) make the dynamical Σ(k,ω). Some changes to GWinput are needed. lmfgwd will automatically make these changes if you used switch --sigw in the QSGW tutorial.

With your text editor, modify GWinput. Change these two lines:

 --- Specify qp and band indices at which to evaluate Sigma

into these four lines:

***** ---Specify the q and band indices for which we evaluate the omega dependence of self-energy ---
   0.01 2   (Ry) ! dwplot omegamaxin(optional)  : dwplot is mesh for plotting.
                   : this omegamaxin is range of plotting -omegamaxin to omegamaxin.
                   : If omegamaxin is too large or not exist, the omegarange of W by hx0fp0 is used.

Also change these lines

*** Sigma at all q -->1; to specify q -->0.  Second arg : up only -->1, otherwise 0
  0  0

to

*** Sigma at all q -->1; to specify q -->0.  Second arg : up only -->1, otherwise 0
  1  0

If you have removed intermediate files, you must remake them up to the point where the self-energy is made. Do:

$ lmgwsc --wt --code2 --sym --metal --tol=1e-5 --getsigp --stop=sig fe

This step is not necessary if you have completed the QSGW Fe tutorial without removing any files.

The next step will make Σ(kn,ω) on a uniform energy mesh −2 Ry < ω < 2 Ry, spaced by 0.01 Ry at irreducible points kn, for QP levels specified in GWinput. This is a fairly fine spacing so the calculation is somewhat expensive.

  • Run hsfp0 (or better hsfp0_om) in a special mode --job=4 to make the dynamical self-energy.
export OMP_NUM_THREADS=8
export MPI_NUM_THREADS=8
~/bin/code2/hsfp0_om --job=4 > out.hsfp0

This step should create SEComg.UP and SEComg.DN. These files contain Σ(k,ω), albeit in a not particularly readable format.

spectral, the self-energy translator

spectral is a postprocessor that reads SEComg.UP (and SEComg.DN in the spin polarized case).

Its main purpose is to translate these files into file se.ext which lmfgws can read.

SEComg.UP and SEComg.DN contain the diagonal matrix element Σjj(k,ω)\Sigma_{jj}(\mathbf{k},\omega) for each QP level j, for each irreducible point kn in the Brillouin zone, on a uniform mesh of points ω as specified in the GWinput file of the last section. If the absence of interactions, Σjj(k,ω)=0\Sigma_{jj}(\mathbf{k},\omega)=0 so the spectral function would be proportional to δ(ωω*), where ω* is the QP level (see Theory section).

Interactions give Σjj(k,ω)\Sigma_{jj}(\mathbf{k},\omega) an imaginary part which broadens out the level, and in general, ReΣjj(k,ω)\mathrm{Re}\Sigma_{jj}(\mathbf{k},\omega) shifts and renormalizes the quasiparticle weight by Z. As noted in the Theory section, there is no shift if VxcjV_\mathrm{xc}^j is the QSGW self-energy Σjj0(k,ω)\Sigma_{jj}^0(\mathbf{k},\omega); there remains, however, a reduction in the quasiparticle weight. This will be apparent when comparing the interacting and noninteracting DOS.

1. Setup for lmfgws

Starting from SEComg.UP (and SEComg.DN in the magnetic case) generated by hsfp0, use spectral to generate the  se  file as described in the lmfgws tutorial below.

2. Use spectral to directly generate spectral functions for q=0

spectral also has a limited ability to directly generate spectral functions from raw output SEComg.{UP,DN} which this section demonstrates.

Do the following:

$ spectral --eps=.005 --domg=0.003 '--cnst:iq==1&eqp>-10&eqp<30'

Command-line arguments are described here. In this context they have the following meaning:

  • --eps=.005 :   0.005 eV is added to the imaginary part of the self-energy. This is needed because as ω→0,  ImΣ→0. Peaks in A(k,ω) become infinitely sharp for QP levels near the Fermi level.

  • --domg=.003 :   interpolates Σ(kn,ω) to a finer frequency mesh. ω is spaced by 0.003 eV. The finer mesh is necessary because Σ varies smoothly with ω, while A will be sharply peaked around QP levels.

  • --cnst:expr :   acts as a constraint to exclude entries in SEComg.{UP,DN} for which expr is zero.
    expr is an integer expression using standard Questaal syntax for algebraic expressions. It can that can include the following variables:

    • ib (band index)
    • iq (k-point index)
    • qx,qy,qz,q (components of q, and amplitude)
    • eqp (quasiparticle energy, in eV)
    • spin (1 or 2)

    The expression in this example, iq==1&eqp>-10&eqp<30, does the following:
       generates spectral functions only for the first k point (the first k point is the Γ point)
       eliminates states below the bottom of the Fe s band (i.e. shallow core levels included in the valence through local orbital)
       eliminates states 30 eV or more above the Fermi level.

spectral writes files sec_ibj_iqn.up and sec_ibj_iqn.dn which contain information about the G for band j and the k point kn. A sec files takes the following format:

# ib=   5  iq=   1  Eqp=   -0.797925  q=   0.000000   0.000000   0.000000
#     omega         omega-Eqp     Re sigm-vxc    Im sigm-vxc      int A(w)      int A0(w)       A(w)           A0(w)
  -0.2721160D+02 -0.2641368D+02 -0.6629516D+01  0.1519810D+02  0.2350291D-04  0.6897219D-08  0.7774444D-02  0.2281456D-05
  -0.2720858D+02 -0.2641065D+02 -0.6629812D+01  0.1520157D+02  0.4701215D-04  0.1379602D-07  0.7776496D-02  0.2281979D-05
  ...

spectral also makes the k-integrated DOS. However, the k mesh is rather coarse and a better DOS can be made with lmfgws. See below.

 spectral: read 29 qp from QIBZ
 Dimensions from file(s) SEComg.(UP,DN):
 nq=1  nband=9  nsp=2  omega interval (-27.2116,27.2116) eV with (-200,200) points
 Energy mesh spacing = 136.1 meV ... interpolate to target spacing 3 meV.  Broadening = 5 meV
 Spectral functions starting from band 1, spin 1, for 9 QPs

          file            Eqp      int A(G)   int A(G0) rat[G] rat[G0]
      sec_ib1_iq1.up   -8.743948     0.8473     0.9999     T     T
      sec_ib2_iq1.up   -1.674888     0.8251     0.9999     T     T
      sec_ib3_iq1.up   -1.674819     0.8251     0.9999     T     T
      sec_ib4_iq1.up   -1.674753     0.8251     0.9999     T     T
      sec_ib5_iq1.up   -0.795683     0.8278     0.9999     T     T
      sec_ib6_iq1.up   -0.795556     0.8278     0.9999     T     T
      sec_ib7_iq1.up   24.572881     0.7374     0.9994     T     T
      sec_ib8_iq1.up   24.572882     0.7374     0.9994     T     T
      sec_ib9_iq1.up   24.572884     0.7374     0.9994     T     T

 writing q-integrated dos to file dos.up ...
 Spectral functions starting from band 1, spin 2, for 9 QPs

          file            Eqp      int A(G)   int A(G0) rat[G] rat[G0]
      sec_ib1_iq1.dn   -8.458229     0.8447     0.9998     T     T
      sec_ib2_iq1.dn    0.015703     0.8718     0.9999     T     T
      sec_ib3_iq1.dn    0.016072     0.8700     0.9999     T     T
      sec_ib4_iq1.dn    0.016437     0.8688     0.9998     T     T
      sec_ib5_iq1.dn    1.552938     0.8363     0.9998     T     T
      sec_ib6_iq1.dn    1.553722     0.8364     0.9999     T     T
      sec_ib7_iq1.dn   24.695801     0.7317     0.9994     T     T
      sec_ib8_iq1.dn   24.695832     0.7317     0.9994     T     T
      sec_ib9_iq1.dn   24.695865     0.7317     0.9994     T     T

 writing q-integrated dos to file dos.dn ...

Dynamical self-energy editor lmfgws

lmfgws is the dynamical self-energy editor, which peforms a variety of postprocessing of the GW or DMFT self-energy Σ(kn,ω)\Sigma(\mathbf{k}_n,\omega) for different purposes. The collection of points kn are typically supplied for a regular mesh. This need not be the case, but when the points do not correspond to a regular mesh the parts of the editor that require k interpolation are not allowed. You must tell the editor that you are not using a uniform mesh (see irrmesh in the instructions below).

lmfgws requires the same files lmf needs to compute the QSGW band structure, e.g. ctrl.ext and sigm.ext. In addition it requires the dynamical self-energy se.ext in special a format written by the spectral utility.

Setting up lmfgws with the spectral utility

For definiteness this section assumes that ext is fe. Starting from SEComg.UP (and SEComg.DN in the magnetic case) generated by hsfp0, use spectral to generate se.fe:

$ spectral --ws --nw=1
$ ln -s se se.fe
  • --ws tells spectral to write the self-energy to file se for all k points, in a special format designed for lmfgws. Individual files are not written. It must be renamed to se.ext for use by lmfgws.
  • --nw=1 tells spectral to write the self-energy on the frequency mesh it was generated; no interpolation takes place.

Try starting the dynamical self-energy editor:

$ lmfgws ctrl.fe `cat switches-for-lm` --sfuned

You should see:

 Welcome to the spectral function file editor.  Enter '?' to see options.

 Option :

The editor operates interactively. It reads a command from standard input, executes the command, and returns to the  Option prompt waiting for another instruction. The editor will print a short summary of instructions if you type  ? <RET>.

Editor instructions

This sections documents the instruction set of the dynamical self-energy editor. Codes that can generate the input for this editor are spectral and lmfdmft.

  • readsek[flags]  |  readsekb[flags]   [fn]
    reads the self-energy from an ASCII (or binary) file. In the absence fn, the file name defaults to the ASCII file se.ext (readsek), or the binary seb.ext (readsekb).
    The structure of the file is documented here. Data is read in the basis of 1-particle eigenfunctions for whatever states are supplied in the file. Some points of note:
    1. Data is stored for a collection of k points; the list of points is written in the file. These points may, or may not constitute a uniform mesh of points.
    2. QP levels are stored relative to the chemical potential (which may, but need not, be written in the header).
    3. Only the diagonal elements of the potentials are read. The full complement of static potentials consist of the static QSGW self-energy Σ0\Sigma^0, the Fock exchange VxxV_\mathrm{xx}, and VxcLDAV_\mathrm{xc}^\mathrm{LDA}.
    4. The se file may, but need not, contain these potentials. For example, none are supplied by lmfdmft.
    5. Optional flags are strung together and separated by a delimiter, taken is the first character, e.g.  @ .
      Note: if you are operating the editor in batch mode, be sure to distinguish this delimiter from the batch mode delimiter.
      • @irrmesh  points are not on a regular  k  mesh : no  k  interpolations allowed
      • @ib=list    after reading data from file, pare bands read from file to those in list
      • @minmax    print minimum and maximum QP levels for each band
      • @readef   file chemical potential becomes Fermi level
      • @makeqpse  Not documented yet
  • units Ry  |  units eV
    select Rydberg units or electron volt units for which data is to be output (default=Ry).
    Note: the se file can store data in either eV or Ry units; lmfgws internally converts it to whatever units you select.

  • evsync
    replace quasiparticle levels read from se.ext by recalculating them with the same algorithm lmf uses.

  • eps val
    add a constant val to Im Σ, needed to broaden spectral functions so that integrations are tractable.

  • ef ef
    Use ef for the Fermi level, overriding the internally calculated value.

  • dos   [nq=#1,#2,#3]   [nw=#|domg=#]   [range=#1,#2]   [isp=#]
    Integrate spectral function to make both the QP and spectrum DOS. Options are:

    • nq=#1,#2,#3     Interpolate Σj(kn,ω) to a new uniform mesh of k points, defined by (#1,#2,#3) divisions.
    • nw=n                    Refine the given energy mesh by interpolating Σ to an n multiple of the given energy mesh. n must be an integer.
    • range=#1,#2     Generate DOS in a specified energy window (#1,#2), in eV.
    • isp=i                      Generate DOS for spin i (1 or 2). Default value is 1.
  • se   iq=n|q=#1,#2,#3|allq|band[~args]   ib=list|ibx=list   [getev[=#1,#2,#3]]   [nw=n|domg=#]   [isp=#]   [range=#1,#2]
    Make Σ(ω) and A(ω) for given q and range of bands.
         Required arguments are:
    • iq=n                            index to qn, from list in QIBZ. Alternatively specify q by:
    • q=#1,#2,#3           q-point in units of 2π/alat. lmfgws will interpolate Σ(qn) to any q.
    • allq        Σ(ω) is made for all q in the irreducible BZ and written to disk
    • band       A(ω), Σ(ω) are made for qp along symmetry lines and written to disk.
              Use this mode to draw interacting energy bands, in conjunction with plbnds −sp
              Optional ~args are parsed like options of the --band switch
    • ib=list       Sum together Aj(ω) derived from QP states j in list. See here for the syntax of integer lists.
              ibx=list is similar to ib=list, but Aj(ω) is resolved by band, writing each Aj(ω) in succession.
      Options are:
    • getev                         Do not interpolate energy but calculate it at q.
    • getev=#1,#2,#3 Generate evals on independent mesh with #1,#2,#3 divisions of uniformly spaced points.
    • nw=n                         Refine the given energy mesh by interpolating Σ to an n multiple of the given energy mesh. n must be an integer.
    • range=#1,#2         Generate spectral function in a specified energy window (#1,#2)
  • pe|peqp   iq=n|q=#1,#2,#3   ib=#   [getev[=#1,#2,#3]]   [nw=#|domg=#]   [nqf=#]   [ke0=#]   [isp=i]   [range=#1,#2]
    Model ARPES for given q and band(s).
    pe uses the spectrum self-energy, while peqp uses just the quasiparticle hamiltonian. Final-state effects are folded into both. Only the latter works with SO coupling now.
    Required arguments are:
    • iq=n                            index to qn, from list in QIBZ. Alternatively specify q by:
    • q=#1,#2,#3           q-point in units of 2π/alat. lmfgws will interpolate Σ(qn) to any q.
    • ib=list                       Sum together PE spectrum derived from QP states j in list. See here for the syntax of integer lists.
      Options are:
    • getev                        Do not interpolate energy but calculate it at q.
    • getev=#1,#2,#3 Generate evals on independent mesh with #1,#2,#3 divisions of uniformly spaced points.
    • nw=n                         Refine the given energy mesh by interpolating Σ to an n multiple of the given energy mesh. n must be an integer.
    • isp=i                          Generate spectra for spin i (1 or 2). Default value is 1.
    • nqf=n                       number of mesh points for final state integration. Default is 200.
    • ke0=#                       kinetic energy of emitted electron. KE+V0=ℏω−φs+V0
    • range=#1,#2   Generate spectral function in a specified energy window (#1,#2)
  • qpse
    Generates the Quasiparticle “self-energy” (in practice the QP levels relative to the Fermi level).

  • savesea [fn]
    saves spectrum DOS or self-energy + spectral function, in ASCII format. In the absence fn, the file name defaults to seia.ext or seia2.ext when writing band and k-resolved spectral functions (se or pe) and to sdos.ext or sdos2.ext when writing spectrum dos (dos).

  • savese [fn]
    saves q-interpolated self-energy + spectral function in binary format. In the absence fn, the file name defaults to seib.ext.

  • q
    quits the editor unless information has generated that has not been saved. Program terminates.
  • a
    (abort) unconditionally quits the editor. Program terminates.

You can also run the editor in batch mode by stringing instructions together separated by a delimiter:

$ lmfgws ctrl.fe `cat switches-for-lm` '--sfuned~first command~second command~...'

The delimiter ( ~  in this case), is the first character following --sfuned. lmfgws will parse through all the commands sequentially until it encounters “quit” instruction ( ~a  or  ~q ) which causes it to exit. If no exit instructions are encountered, lmfgws returns to the interactive mode and prompts you with  Option : .

Compare interacting and independent-particle density-of-states in Fe

This section uses the self-energy editor, lmfgws, to interpolate Σ(kn,ω) to a fine k- and ω- mesh to obtain a reasonably well converged density-of-states.

$ lmfgws fe `cat switches-for-lm` '--sfuned~units eV~readsek~eps .030~dos isp=1 range=-10,10 nq=32 nw=30~savesea~q'

This invocation runs lmfgws in batch mode, and writes the spectral and noninteracting DOS to file sdos.fe. The editor’s instructions do the following (documented here):

  • units eV
    Set units to eV; spectrum DOS will be written in eV.
  • readsek
    Read se.fe
  • eps .030
    Add 30 meV smearing to Im Σ
  • dos isp=1 range=-10,10 nq=32 nw=30
    Make the DOS for spin 1, in the energy range (-10,10) eV, interpolating Σ to a k mesh 32×32×32 divisions, and refining the energy mesh by a factor of 30. The as-given k mesh is 8×8×8 divisions.
  • savesea
    Write the DOS.
  • q
    Exit the editor.

Notes:

  • The mesh is very fine, so the interpolation takes a little while (2 to 3 minutes). The frequency and k meshes are both pretty fine and the DOS is rather well converged, as the figure below demonstrates.
  • The spectrum DOS is written to file sdos.fe. Columns 1,2,3 are ω, A(ω), and A0(ω), respectively.
  • A0(ω) should compare directly to the DOS calculated as a byproduct of lmf.

You can make the QP DOS yourself, but to speed things up just copy it from the build directory to your working directory.

cp ~/lm/gwd/test/fe/dosp.fe dosp.fe

The following script draws a figure comparing the DOS generated the three different ways, using the fplot utility. Cut and paste the contents of the box below into script file plot.dos.

% char0 ltdos="1,bold=3,col=0,0,0"
% var ymax=1.4 dy=0.4 dw=.00 ymax+=dy emin=-10 emax=5 ef=0
fplot

% var ymax-=dy+dw dy=0.4 dmin=0 dmax=3
 -frme 0,1,{ymax-dy},{ymax} -p0 -x {emin},{emax} -y {dmin},{dmax} -tmy 1 -1p
 -colsy 3 -lt 1,bold=3,col=.5,.5,.5 sdos.fe
 -colsy 2 -lt {ltdos} -ord y -qr dosp.fe
 -colsy 2 -lt 1,bold=3,col=1,0,0 sdos.fe
 -lt 2,bold=3,col=0,0,0,2,.5,.05,.5 -tp 2~{ef},{dmin},{ef},{dmax}

k-integrated spectral function in Fe

$ fplot -f plot.dos
$ open fplot.ps   [choose your postscript file viewer]

Notes on the figure:

  • The black line (col=0,0,0) is the noninteracting DOS generated by lmf.
  • The grey line (col=.5,.5,.5) is the noninteracting DOS A0(ω), generated by lmfgws
  • The red line (col=1,0,0) is the interacting DOS A(ω), generated by lmfgws
  • Grey and black lines nearly coincide, as they should if the DOS is well converged. Note that the black line was generated from energy bands with the tetrahedron method, the other effectively by integrating G0(k,ω) by sampling with a smearing of 30 meV.
  • The noninteracting DOS at the Fermi level is D(EF)≅1/eV (one spin). The Stoner criterion for the onset of ferromagnetism is I×D(EF)>1, where I is the Stoner parameter, which DFT predicts to be approximately 1 eV for 3d transition metals. Combining DOS for the two spins would indicate that the Stoner criterion is well satisfied.
  • The interacting DOS is smoothed out, and is roughly half the amplitude of the noninteracting DOS. This is also expected: the Z factor for the d states is about 0.5.

Spectral Function of Fe near the H point

This example computes the self-energy for a q point near the H point. It is calculated from band 2 for the majority spin and bands 2,3 for the minority spin. These bands were chosen because of their proximity to the Fermi level.

$ lmfgws fe `cat switches-for-lm` '--sfuned~units=eV~eps .01~readsek~evsync~se q=1.05,2.91,1.01 ib=2 nw=10 getev=12 isp=1~savesea~q'
$ lmfgws fe `cat switches-for-lm` '--sfuned~units=eV~eps .01~readsek~evsync~se q=1.05,2.91,1.01 ib=2,3 nw=10 getev=12 isp=2~savesea~q'

The first command writes a file seia.fe, the second seia2.fe

The following makes a picture comparing A (solid lines) and A0 (dashed lines), majority spin (black) and minority spin (red)

$ fplot -x -9,5 -y 0,1 -colsy 2 -lt 1,col=0,0,0 seia.fe -colsy 3 -lt 2,col=0,0,0 seia.fe -colsy 2 -lt 1,col=1,0,0 seia2.fe -colsy 3 -lt 2,col=1,0,0 seia2.fe
$ open fplot.ps   [choose your postscript file viewer]

k-integrated spectral function in Fe

You should see a weak plasmon peak in the majority spin band near −8 eV.

Interacting band structure

This block uses lmfgws to generates the band structure of the interacting Green’s function, i.e. the k-resolved spectral function along symmetry lines similar to a band plot for a noninteracting G0G_0. Peaks in the spectral function correspond to the band structure; the plot can be compared directly to the bands of the noninteracting G0. Use syml.fe from that tutorial, or use file syml2.fe, which contain the symmetry lines as appear in Figure 1 of this Phys. Rev. B paper. Invoke lmfgws in batch mode as follows:

$ cp ~/lm/gwd/test/fe/syml2.fe .
$ lmfgws fe `cat switches-for-lm` '--sfuned~units=eV~readsek~eps .01~evsync=6~se band:fn=syml2 ib=1:10 nw=10 getev=12 isp=1 range=-10,10'

The self-energy editor carries out the following

  • units eV
    Set units to eV
  • readsek
    Read se.fe
  • eps .01
    Add 10 meV smearing to Im Σ
  • evsync
    refresh quasiparticle levels read from se.fe by recalculating them.
  • se   band:fn=syml2   ib=1:10   nw=10   getev=12   isp=1   range=-10,10
    Generate the self-energy and spectral function A(k,ω)A(\mathbf{k},\omega) along symmetry lines given in file syml2.fe. Include bands 1-10, and generate A(k,ω)A(\mathbf{k},\omega) on a frequency mesh 10× finer than the one in se.fe. getev refines the k-mesh to a 12×12×12 mesh, and using that mesh to interpolate bands along symmetry lines in syml2.fe. Genearte bands in an energy window [−10,10] eV.

lmfgws writes a file, spq.fe.

Invoke plbnds in “spectral function mode:”

$ plbnds -sp~atop=10~window=-4,4 spq.fe

It will generate a gnuplot script file gnu.plt together with a data file spf.fe.

Run gnuplot

$ gnuplot gnu.plt
$ open spf.ps   [choose your postscript file viewer]

to generate and view postscript file spf.ps.


Edit This Page