Making the dynamical GW self energy

An exact one-body description of the many-body Schrödinger equation requires a time-dependent, non-hermitian potential, called the self-energy $\Sigma$. This is because independent-particle states that solve a one-body Schrödinger equation are no longer eigenstates in the many-body case. Nevertheless many-body descriptions are usually characterized in terms of a one-particle basis. The imaginary part of $\Sigma$ carries information about the inverse lifetime (the time for decay of a state to another state).

$\Sigma$ is in general nonlocal in both space and time. Nonlocality in time implies that $\Sigma$ depends on two time coordinates $t$ and $t^\prime$; however in equilibrium it depends only on the difference $t-t^\prime$. Converting time to its Fourier (frequency) representation, $\Sigma$ depends only on a single frequency $\omega$. Nonlocality in space implies that $\Sigma$ depends on two coordinates $\mathbf{r}$ and $\mathbf{r}^\prime$. Translational invariance of a crystal implies that for translations between one unit cell and another, $\Sigma$ depends only on one (lattice) translation vector $\mathbf{T}$. With a Fourier from $\mathbf{T}$ to $\mathbf{k}$ representation, $\Sigma$ depends only on one $\mathbf{k}$ vector. For points within a unit cell, $\Sigma$ depends on two coordinates $\mathbf{r}$ and $\mathbf{r}^\prime$ confined to a unit cell. Thus in general $\Sigma$ is written as $\Sigma=\Sigma(\mathbf{k},\mathbf{r},\mathbf{r}^\prime,\omega)$.

The GW approximation indeed makes a fully nonlocal $\Sigma(\mathbf{k},\mathbf{r},\mathbf{r}^\prime,\omega)$, and this tutorial described how to construct and analyze it. The main utility used to analyze $\Sigma$ is lmfgws. Dynamical Mean Field Theory (DMFT) is a single-site approximation; it thus makes $\Sigma$ nonlocal in time but local in space. However its time-dependence is calculated to a higher level of theory than is done for the GW approximation. The DMFT package generates $\Sigma$ in a different way but, once created, can also use lmfgws to analyze it.

The Coherent Potential Approximation, or CPA replaces a true atom with some effective average atom (either alloys or atoms of one kind but different moment orientations) Even though it is a one-body approximation, the CPA construction implies that, like DMFT, $\Sigma$ is nonlocal in time but local in space. It operates as a stand-alone package.

Note: QSGW has another form of the self-energy, the “quasiparticlized” form, which has a nonlocal but static self-energy $\Sigma^0(\mathbf{k},\mathbf{r},\mathbf{r}^\prime)$ derived from the dynamical $\Sigma(\mathbf{k},\mathbf{r},\mathbf{r}^\prime,\omega)$. There is an editor for the quasiparticle form also, useful for other contexts. For applications of the static editor, see this page.

Table of Contents

• Table of Contents
• Preliminaries
• Command summary
• Theory
  • Z factor renormalization
  • Coherent part of the spectral function
  • Simulation of Photoemission
• The LQSGW Approximation
• Introduction
  • Setting up the input from distribution files
• Confirm the self-energy is self-consistent
  • Make sig.h5
  • One QSGW iteration
• Make the GW dynamical self-energy
  • Notes on the Fermi level and enforcement of the QSGW condition
  • Notes on interpolating the self-energy in k and ω
• spectral, the self-energy translator
  • Use spectral to directly generate spectral functions for q=0
• Dynamical self-energy editor lmfgws
  • Generate se.fe using the spectral utility
• Editor instructions
  • Invoking the editor: interactive vs batch mode
  • Instruction set
• Editor Application : Compare interacting and independent-particle density-of-states in Fe
• Editor Application : Spectral Function of Fe near the H point
• Editor Application : Interacting joint Density-of-States and Optics
• Editor Application : Interacting band structure
• Editor Application : Generating spectral functions from Brillouin zone unfolding of a supercell

Preliminaries

This tutorial assumes you have completed a QSGW calculation for Fe, following this tutorial. Alternatively you can start from a setup supplied the standard test suite. We will do the latter in this tutorial so you can run it without having to go through a prior one. The ctrl.ext files are a little different but the results are very similar.

This tutorial will do the following:

• Starting from an already-supplied QSGW self-energy Σ⁰, perform a single QSGW cycle to confirm it is self-consistent.
• Generate the dynamical self-energy Σ(k,ω), stored in files SEComg.UP and SEComg.DN.
• Generate spectral functions with the spectral utility.
• Use lmfgws dynamical self-energy postprocessor to generate the interacting density-of-states (DOS) from Im G, compare it to the noninteracting DOS from Im $G_0$ and to the noninteracting DOS generated as an output of an lmf band calculation.
• Use lmfgws to generate to calculate the spectral function A(k,ω) for k near the H point, and also simulate photoemission spectra

Also the instructions for dynamical self-energy editor are documented.

Command summary

Repeat the steps for LDA self-consistency and QSGW self-consistency in the Fe tutorial; 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 $\Sigma^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 --sym --metal --tol=1e-5 --getsigp --stop=sig fe

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

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

env OMP_NUM_THREADS=8 MKL_NUM_THREADS=8 $(dirname $(which lmgw))/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

Theory

Z factor renormalization

Begin with a noninteracting Green’s function $G_0$, defined through an hermitian, energy-independent exchange-correlation potential $V^j_{xc}(k)$. $j$ refers to a particular QP state (pole of $G_0$). There is also an interacting Green’s function, $G$.

The contribution to $G_0$ from QP state $j$ is $$G_0^j(k,\omega ) = \frac{1}{\omega - \omega^j(k)}$$

where $\omega^j(k)$ is the pole of $G_0$.

Write the contribution to $G$ from QP state $j$ as $$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 $\Sigma(k,\omega)$. Note that if $V^j_{xc}(k)$ is constructed by QSGW, this is a very good approximation, since ${\mathrm{Re}\Sigma (k,\omega ){=}V^j_{xc}(k)}$ at $\omega{=}\omega^j(k)$. Approximate G by its coherent part: $$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 $$1 - 1/{Z^j} = \left. {\partial \mathrm{Re} \Sigma (k,\omega )/\partial \omega } \right|_{\omega ^j } .$$

defines the $Z$ factor. The dependence of $\omega^j$ and $Z^j$ on $k$ is suppressed.

Define the QP peak as the value of $\omega$ where the real part of the denominator vanishes. $$({\omega^*} - \omega^j)/Z^j = \mathrm{Re} \Sigma(k,\omega^j) - V^j_{xc}(k)$$

and so $${\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 $d\omega/dk$. For the interacting case it reads $$\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 $$\frac{m_0}{m^*} \equiv \frac{d\omega^*}{dk}/\frac{d\omega^j}{dk} = 1 + \frac{Z^j}{v_0^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) \frac{1}{v_0^j}$$

Ignore the dependence of $Z^j$ on $k$. Write $\text{d}\omega^j/\text{d}k$ as $v_0^j$, and use the definition of $Z^j$ to get $$\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 $$\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 $G^{j,\mathrm{coh}}(k,\omega)$ as $$G^{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 $$G^{j,\mathrm{coh}}(k,\omega) = \frac{Z^j}{\omega - \omega^* - iZ^j \mathrm{Im} \Sigma (k,\omega )} = Z^j\frac{\omega - \omega^* + iZ^j \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(\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. $$A_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) $$\hbar\omega=E_{kin}+{\varphi_s}-E_b$$

where E_b is the binding energy and $E_{kin}+{\varphi_s}$ is the energy of the electron after being ejected. (Marder defines $E_{b}$ with the opposite sign, making it positive).

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

Let $E_{kin}$ be the energy on exiting the crystal, $\varphi_s$ the work function and $E_b$ and $V_0$ are called the electron binding energy and “inner potential.”

Then $$\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, $\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 ${E_{kin}}=\hbar\omega-{E_b}-{\varphi_a}$ where ${\varphi_a}$ is the work function of the analyzer. Usually ${\varphi_a}{=}{\varphi_s}$. The Fermi level is defined such that $E_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 $$\frac{\hbar ^2 k_0^2 }{2m} = E_{kin}$$

Resolve $\mathbf{k}_f$ into components parallel and perpendicular to the surface $$\mathbf{k}_f = \mathbf{k}_\parallel + \mathbf{k}_\bot$$

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

The conservation condition requires $$k_0^2 = \bar k_\parallel^2 + \bar k_\bot^2$$

$\mathbf{k}_\parallel$ is conserved on passing through the surface; thus $\bar k_\parallel{=} k_\parallel$. $\mathbf{k}_\bot$ is not conserved; therefore $$\bar k_\bot = \sqrt{k_0^2{-}k_\parallel^2}$$

The wave number shift is then $$\Delta{\mathbf{k}} = (\overline{k}_\bot-{k}_\bot)\hat{\mathbf{k}}_\bot$$

and the crystal momentum actually being probed by the experiment is $${\mathbf{k}}_f = \overline{\mathbf{k}}_f - \Delta{\mathbf{k}}$$

The LQSGW Approximation

Kutepov’s LQSGW theory of is a linearized form of QSGW. He constructs the quasiparticlized self-energy from a Taylor series around the origin. In his formalism (treating each band indpendently and suppressing band index, for simplicity of presentation) replaces the interacting $G$ $$G^{-1} (k,\omega ) = \omega + \mu - \epsilon - \Sigma(k,\omega)$$

by omitting the second order and higher terms of an expansion in $\Sigma$ $$\Sigma (k,\omega ) = \Sigma(k,0) + \omega\,\Sigma'(k,0) + \frac{1}{2}\omega^2\,\Sigma''(k,0) + ...$$

$G^{-1}$ simplifies to a linear function of $\omega$ $$G^{-1} (k,\omega ) = \bar{Z}^{-1}\, \omega + \mu - \epsilon - \Sigma(k,0)$$

We use a bar to denote the $Z$ factor since it is defined at $\omega{=}0$: $$1 - 1/{\bar{Z}^j} = \Sigma'(k,0)$$

Evidently $\epsilon - \mu + \Sigma(k,0)$ is the eigenvalue of a hamiltonian defined as the one-body part of $G^{-1}$, but adding the static part of $\Sigma$. The (linearized) energy-dependence of $\Sigma$ modifies this eigenvalue to read $$E-\mu = \bar{Z}[\epsilon-\mu + \Sigma(k,0)]$$

$E-\mu$ is identical to the QSGW $\omega^j$ if $\Sigma$ is a linear function of $\omega$.

Now let us retain the quadratic term and determine the shift in $E$ to estimate the difference between LQSGW and QSGW. Let us denote the LQSGW eigenvalue $E-\mu$ as $E_0$. Expanding $G^{-1}$ to second order we obtain, to lowest order in $\Sigma''(k,0)$ : $$G^{-1} = \omega - (E_0 + \frac{\bar{Z}}{2}E_0^2\,\Sigma''(k,0) )$$

Introduction

The GW calculations are performed using the main GW script, lmgw.sh, which is assumed to be in your path together with the other binaries compiled in your build directory.

This tutorial takes files from in the standard test suite. Here we refer to the top-level directory where Questaal is installed as $toplevel (e.g. ~/lmf).

Note: This tutorial follows the steps taken from the standard test case

$toplevel/gwd/test/test.gwd --mpi=#,# fe 4

Setting up the input from distribution files

As noted here, complications can arise if you do not start from a fresh directory. Either start from a fresh directory, or initialize your working directory with the following:

rm -f {wkp,bnds,rst,rst0,save,log,hssn,sigm,erange}.fe switches-for-lm QPU sigm.fe}
touch 1run meta sig.h5 mixsigma
rm -r *run meta *.h5 mix*

Copy the following setup files:

cp $toplevel/{gwd/test/fe/syml.fe,gwd/test/fe/syml2.fe,gwd/test/fetb/GWinput.gw,gwd/test/fetb/switches-for-lm,gwd/test/fetb/sigm.in,gwd/test/fetb/basp.fe,gwd/test/fetb/ctrl.fe,gwd/test/fetb/site.fe,gwd/test/fetb/rst.fe} .
cp GWinput.gw GWinput

Σ⁰ was generated using the tight-binding form of the basis set. Here we will use standard unscreened basis functions. This requires two modifications.

• In this setup, file  switches-for-lm  controls whether the basis is screened or not. (This file contains command-line switches   –tbeh –v8 -vrsham=2   which turns on the screened basis.) lmgw.sh, reads  switches-for-lm  and passes the contents as command-line switches to executables such as the one-body code lmf and the GW interface code lmfgwd.
• Σ⁰ (in file sigm.in) was generated in the screened basis. Thus we need to rotate the self-energy (sigm.in) to an unscreened form.

Make these two modifications as follows:

sed -i.bak 's/--tbeh --v8 -vrsham=2/--v8/' switches-for-lm
lmf --tbeh --v8 -vrsham=2 fe --scrsig:ifl=sigm.in:ofl=sigm.fe:in=s:out=u > out.unscr
cp sigm.fe sigm.unscr

Note the second line creates a file sigm.fe, in unscreened form.. The third line isn’t needed, but for better consistency in the generation of Σ(k,ω), you will want to save it; see below for an explanation

Confirm the self-energy is self-consistent

Make sig.h5

When we are starting with a QSGW self-energy (in file sigm.fe), lmgw.sh requires it only to make the eigenfunctions for the GW cyle. However the self-energy maker (lmsig, called within lmgw.sh) can use this information to check how much the self-energy changes relative to the starting point, and also to use the input-output Σ⁰ pair to accelerate convergence to self-consistency. lmsig does not read from binary sigm.fe, but from sig.h5. Eventually the binary file will be phased out, but for now sigm.fe should be consistent with sig.h5. As sig.h5 does not exist yet, make it as follows

ln -s sigm.fe sigm
s2s5

One QSGW iteration

This system is very simple, and you can run these calculations in serial mode, but you can save considerable time by running in parallel. We will do that here. For an explanation how to set up lmgw.sh to run in parallel for this system, see the QSGW tutorial for Fe.

The next step is not necessary, but you can speed up the self-energy generation by making a pqmap file.

lmfgwd ctrl.fe `cat switches-for-lm` --job=0 --lmqp~rdgwin --batch~np=4~pqmap@plot

Perform one iteration to confirm Σ⁰ is self-consistent

/home/ms4/lmf/s/lmgw.sh --incount --iter=1 --maxit=1 --sym --tol=1e-5 --split-w0 --mpirun-1 "env OMP_NUM_THREADS=16 MKL_NUM_THREADS=16 mpirun -n 1" --mpirun-n "env OMP_NUM_THREADS=4 MKL_NUM_THREADS=4 mpirun -n 4" ctrl.fe > lmgw.log

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

Confirm that the change in self-energy is small:

grep more lmgw.log

You should see a line similar to

lmgw: iter 1 of 1 completed in 149s, 149s from start. RMS change in sigma = 1.54e-5  Tol = 1e-5  more=F

This is a small change. You can also confirm that the input and output charge densities are very similar

grep DQ llmf

Make the GW dynamical self-energy

This section makes the dynamical self-energy Σ(k,ω). You must complete the prior step before doing this one. Note only the diagonal elements of Σ(k,ω) are calculated.

The instructions in this section should create SEComg.UP and SEComg.DN. These files contain Σ(k,ω), albeit in a not particularly readable format; they will be translated to another format when postprocessing steps are performed.

If you repeat this section, clean the directory by removing these files

touch se
rm -f SEComg.{UP,DN} se se.fe lsg out.spectral out2.spectral out.lmfgws

Note: If self-consistency were entire, the input Σ^0,in would match Σ^0,out generated by the GW cycle. This is almost, but not exactly the case. The difference is slight but not totally negligible very near E_F because lmsig makes Σ(k,ω), from the QP levels generated from Σ^0,in, and calculates E_F from them. (Note further E_F is computed by sampling since lmsig uses sampling to make Σ(k,ω).)

lmsig generates Σ(k,ω) from intermediate files generated in the previous step, which used Σ^0,in for their construction. However at the final stage of that step, rst and sigm have been updated, which means that E_F corresponding to the new Σ is a little different. It is not necessary to do this, but to be consistent, you should restore both rst and sigm before proceeding in the generation of Σ(k,ω).

cp rst0.fe rst.fe
cp sigm.unscr sigm

This is further eplained in a section below.

Make SEComg.UP and SEComg.DN either in parallel or serial mode. Do one of the following

lmsig --v8 --rdwq0 --dynsig~nw=401~wmin=-2.0~wmax=2.0~bmax=11 ctrl.fe >& lsg
env OMP_NUM_THREADS=4 MKL_NUM_THREADS=4 mpirun -n 4 lmsig --v8 --rdwq0 --dynsig~nw=401~wmin=-2.0~wmax=2.0~bmax=11 ctrl.fe >& lsg

In parallel mode it should execute in about 3 minutes on a standard Intel X86 machine.

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

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

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

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

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

export OMP_NUM_THREADS=8
export MPI_NUM_THREADS=8
$(dirname $(which lmgw))/code2/hsfp0_om --job=4 > out.hsfp0

Notes on the Fermi level and enforcement of the QSGW condition

The dynamical self-energy is computed as a function of frequency. To the extent the QSGW self-consistency condition is satisfied, that is Σ(k,ω=E^QP-E_F) = Σ⁽⁰⁾(k), the pole of the Green’s function is given solely by E^QP. For this condition to be properly satisfied.

Also, for reasons noted earlier it may not be exactly satisfed, for example if the self-consistency is not fully realized. In any case, the dynamical postprocessor lmfgws normally enforces this condition, unless you specify otherwise. Since the enforcement depends on knowledge of the quasiparticle E_F, and error is introduced if E_F is not precisely specified. Usually the uncertainty in E_F is unimportant, but there is structure near E_F on the scale of the uncertainty, this can be important. To work around this the lmfgws editor has instructions to allow you to specify it. The most automatic way is to tell it to find E_F from a finely spaced k mesh (see getef in the instructions below). You can also specify E_F by hand (see ef= in the instructions).

Once you specify E_F, lmfgws will enforce the QSGW condition, unless to tell it not to.

Notes on interpolating the self-energy in k and ω

For reasons of cost, only the diagonal element of Σ_jj of Σ_ij(k,ω) is calculated. Also for reasons of cost, Σ_jj is usually computed on a relatively coarse mesh, and for a limited number of ω. However, quantities often need to be determined on a finer k-mesh, and also a finer ω-mesh, which to accomplish Σ_jj must be known there as well. This is accomplished by an interpolation scheme. Normally when Σ_jj(k,ω) to a k where it is not known, lmfgws will shift either Σ_jj(k,ω) or Σ⁽⁰⁾(k) by a constant to enforce the QSGW condition. There is no other adjustment of the ω-dependence.

Also, in the absence of extra information, the QP levels E^QP must be interpolated. This error can be more severe than interpolating Σ. Fortunately lmfgws can generate E^QP for any k, and this is the recommended approach when you seek quantities on a k mesh other than the one used to generate Σ_jj(k,ω); see getef in the instructions below.

spectral, the self-energy translator

spectral is a utility that reads SEComg.UP (and SEComg.DN in the spin polarized case). SEComg.UP and SEComg.DN contain the diagonal matrix element $\Sigma_{jj}(\mathbf{k},\omega)$ for each QP level j, for each irreducible point k_n in the Brillouin zone, on a regular mesh of points ω as specified in the previous section. If the absence of interactions, $\Sigma_{jj}(\mathbf{k},\omega)=0$ so the spectral function would be proportional to δ(ω−ω^*), where ω* is the QP level (see Theory section).

Interactions give $\Sigma_{jj}(\mathbf{k},\omega)$ an imaginary part which broadens out the level, and in general, $\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 $V_\mathrm{xc}^j$ is the QSGW self-energy $\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.

Most often, spectral is use to translate SEComg.UP (SEComg.DN) into file se.ext which the dynamical self-energy postprocessor lmfgws can read. This is discussed in the below.

spectral also has a limited ability to create spectral functions in its own right, as described next.

Use spectral to directly generate spectral functions for q=0

spectral 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 Σ(k_n,ω) 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 Im G for band j and the k point k_n. 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
 ...
 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
...
 writing q-integrated dos to file dos.dn ...

Dynamical self-energy editor lmfgws

This section explains the features lmfgws has to construct various properties of the interacting Green’s function. Later sections will demonstrate some of these features for Fe.

lmfgws is the dynamical self-energy editor, which performs a variety of postprocessing of the GW or DMFT self-energy $\Sigma(\mathbf{k}_n,\omega)$ for different purposes. Typically $\Sigma$ is supplied on a uniform mesh of points $\mathbf{k}_n$. It can interpolate in both $\mathbf{k}$ and $\omega$ to provide information for any $\mathbf{k}$ and $\omega$.

lmf computes properties the noninteracting QSGW band structure, while lmfgws computes the corresponding interacting one. They require the same input files; in addition lmfgws requires the dynamical self-energy, se.ext. In the GW context, see above for its construction. You need use the file translater explained below to make se.ext.

lmfgws can:

• Generate the spectral function $A(\mathbf{k},\omega)$ at specified $k$ points
• compute the interacting density-of-states (DOS). Entails an integral of $A$ over $\mathbf{k}$.
• Simulate ARPES with a simple model for final-state scattering
• Compute the joint noninteracting and interacting DOS
• Compute the noninteracting and interacting imaginary part of the dielectric function
• Generate the interacting band structure along specified symmetry lines

Note 1: It is not required that the given k_n are supplied on a regular mesh, but in such a case some instructions of the editor (those that require k interpolation) are not allowed. Also in that a case you must tell the editor that you are not using a uniform mesh (see irrmesh in the instructions below).

Note 2: In Feb. 2022, lmfgws was parallelized to work with MPI. (If use the editor with MPI, typically you need to operate strictly in batch mode as well.) Most of the editor instructions have been parallelized, though not all of them. Run lmfgws as mpirun -np # lmfgws ...

These files will be made in later stages of this tutorial. If you have already run parts of this tutorial, delete these files:

rm -f sdos.fe seia.fe pes2.fe pesqp0.fe spq.fe spq-bnds.fe seq.fe jdos-lmf.fe jdosni.fe

Generate se.fe using the spectral utility

Use spectral to translate SEComg.UP (and SEComg.DN since this system is magnetic) into se.fe:

spectral --pr31 --ws --nw=1 > out2.spectral
cp 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 frequency interpolation takes place. This will be done by lmfgws.

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.

Invoking the editor: interactive vs batch mode

You can run the editor interactively, or in batch mode.

Interactive mode

To operate lmfgws in interactive mode, invoke

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

Batch mode

You can also run the editor in batch mode by stringing instructions together separated by a delimiter, e.g.

lmfgws fe '--sfuned~units=eV~eps@.01'

The delimiter ( ~  in this case), is the first character following --sfuned. lmfgws will parse through all the commands sequentially until it encounters “quit” instruction ( ~q ) which causes it to exit. If no such instruction is encountered, lmfgws returns to the interactive mode.

Instruction set

Note: Other delimiters may be used in place of  @  assumed below, but if you are operating the editor in batch mode, be sure to distinguish this delimiter from the batch mode delimiter. You can use a space as a delimiter here, but in batch mode, you need to enclose ‘--sfuned’ in single quotation marks. If you run lmfgws with mpirun, it may be safer to avoid using spaces, as mpirun may strip the quotes.

• readsek[flags]  |  readsekb[flags]  |  readeb[flags]  |  readebr[flags]
  reads the self-energy from an ASCII file (readsek), binary file (readsekb), or hdf5 file (readsekh). In the absence fn, the file name defaults to se.ext for the ASCII case, seb.ext for the binary case, and se.h5 for the hdf5 case. Data is read in the basis of 1-particle eigenfunctions for whatever states are supplied in the file. When reading a file generated by the QSGW cycle it typically reads header information, including number of spins in file; number of bands and k points; number of frequencies and the smallest and largest values; the list of bands stored in the file (a subset of all the bands may be stored); and the stored chemical potential if it is written into the file.

  The se file structure is documented here for both ASCII and hdf5 formats.
  
  readeb and readebr are instructions for the Brillouin zone unfolding mode. Its generation and usage proceed a bit differently, and information is read from a bands file, written in hdf5 format (see here for documentation). readeb and readebr are the same except the weight $w_{n\mathbf{k}}$ (see Eq. (5) in this tutorial) is read as written for readeb, while for readebr the weight is turned into $1{-}w_{n\mathbf{k}}$.
  
  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 $\Sigma^0$, the Fock exchange $V_\mathrm{x}$, and $V_\mathrm{xc}^\mathrm{LDA}$,

  4. The se file may, but need not, contain these potentials. se files generated lmfdmft or in the Brillouin zone unfolding mode of lmf contain only Σ(k,ω). On the other hand se files generated by the current GW implementation typically contain: (1) the k-points where Σ is stored; (2) the QP levels generated from the starting self-energy Σ^0,in; (3) the quasiparticlized self-energy Σ⁽⁰⁾; (4) the dynamical self-energy Σ(k,ω).

  5. For the GW-generated case, after reading the file it checks how well the QSGW condition is satisfied at self-consistency, namely that Σ(k,ω) should equal Σ⁽⁰⁾(k) at the QP level. The result is written to stdout in a line like this one:
     <(sig-vxc)@eQP>=-0.000004 RMS=0.000074 shift =-0.000948 RMS=0.029138 Adjust=F

  6. Optional flags are strung together and separated by a delimiter, taken as the first character. The sequences below assume the delimiter is  @ .\

     • @fn=nam   use nam for self-energy file name [for readeb and readebr the file should be a bands file in h5 format].
     • @useef    writes the file chemical potential into the Fermi level.
     • @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
     • @dc=#     subtract double counting from Re sigma(omega) after reading, mode # (for DMFT)
     • @scl=#     scale Σ(k,ω) and Σ⁽⁰⁾(k) by #.
     • @makeqpse  Not documented yet

• units@Ry  |  units@eV
  Select Rydberg units or electron volt units (default=Ry).
  Note: the se file can store data in either eV or Ry units; lmfgws internally converts it to whatever units you select.

• getef[@n=#1[#2,#3]@show]  |  getef=#1[,#2,#3]
  Finds the QSGW fermi level for a specified mesh  #1[#2,#3], and sets it as the default to be used in future operations. If no mesh is specified the mesh taken from se.ext is used.
  Optional show causes lmfgws to suppress assigning the calculated E_F to the default value.

• evsync[@frzsig]
  Replaces quasiparticle levels E^QP read from se.ext by recalculating them internally with lmf, shifting the energy zero to E_F.
  If E_F has already been supplied externally (typically by  readsek@useef,  ef=#, or  getef), that value is used. Otherwise it uses the internal E_F found when E^QP are generated.

  After synchronizing E^QP, the QSGW condition is imposed : Σ(k,ω)=Σ⁽⁰⁾(k) = 0 at E^QP. For QSGW, it is imposed by shifting Σ⁽⁰⁾(k) for QSGW generated se.ext. Optional argument frzsig causes lmfgws not to show how well the QSGW condition is imposed, but not to update the self-energy.

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

• ef=#  |  ef0=#
  Assign  #@useef  to the Fermi level, overriding the internally calculated value.
  Note: the order in which you use this switch is important. If you use the  ef  switch before  readsek, QP levels are shifted by  μ−Ef  when they are subsequently read (provided the chemical potential μ is supplied in the se file). If you use this switch after  readsek, no shifts are added. In such a case you likely want to realign the QP levels with evsync after  readsek. Always enter  #  in Ry units.
  Use ef0=# to set the QP Fermi level. In QSGW it is not distinct from ef, but it is in DMFT.

• dos|jdos|imeps   [getev=#1,#2,#3] | [nq=#1,#2,#3]   ib=list   wts=strn   [getev|getev[=#1,#2,#3]]   [nw=#|domg=#]   [range=#1,#2]   [isp=#]
  dos  integrates the spectral function to make both the QP and spectrum DOS (writes to sdos.ext or sdos2.ext).
  jdos  integrates either the QP or interacting spectral function to make the joint DOS (writes to jdos.ext).
  imeps  integrates either the QP or interacting spectral function to make Im ε (writes to opt.ext or optni.ext).
  Options are:
  • ib=lst        restrict contribution to spectra from QP states in list.
  • nq=#1,#2,#3    Evaluate the physical quantity on a specified mesh. This entails interpolating Σ_j(k_n,ω) to a new uniform mesh of k points, defined by (#1,#2,#3) divisions. Use between 1 and 3 numbers.
  • getev        Calculate the QP levels on the supplied mesh, and enforce the QSGW condition. The enforcement should have already been satisfied if evsync was invoked previously and the Fermi level is consistent with the mesh.
  • getev=#1,#2,#3  Calculate the QP levels on a specified mesh, and enforce the QSGW condition. Thus it is similar to nq=#1,#2,#3 except the QP part (i.e. QSGW part) of the hamiltonian need not be interpolated and the QSGW condition for Σ_j(k_n,ω) is enforced.
  • wts=strn     (Applies to dos only) Add extra columns to file sdos, writing projections of dos onto a subspace, defined through the orbitals of the basis set.
             The projection uses a Mulliken population analysis of the orbitals comprising the basis, to resolve a normalized eigenfunction into some portion of it.
             This is the same procedure used to make color weights for energy bands and it uses the same “color weights” syntax as the --band switch uses.
             Example:   wts#scol@atom=u#scol2@atom=co#scol3@atom=ge
             will generate a pair of columns each for three decompositions: orbitals belonging to U, Co, and Ge respectively.
             If the crystal has nothing but U, Co, and Ge, the sum of the three partial dos should add to unity.
             Hazard: wts requires two levels of delimiters ( @  and  #  in the example above).
             Take care that these delimiters are distinct from the delimiter separating arguments to dos,
             as well as the one separating batch commands if you are using batch mode.
  • 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.
  • kT         Temperature, units of omega (applies only to jdos and imeps).
  • a0         Spectra for noninteracting spectral function (only for jdos and imeps).
  • isp=#       Generate DOS for spin # (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]]   [a2qp]   [nw=n|domg=#]   [isp=#]   [range=#1,#2]
  Make Σ(ω) and A(ω) for given q and range of bands.
  Specify which q by:
  • iq=n        index to q_n, from list in QIBZ. Writes to seia.ext, or to seia2.ext for second spin.
  • q=#1,#2,#3    q-point in units of 2π/alat. lmfgws will interpolate Σ(q_n) to any q. Writes to seia.ext.
  • allq        Σ(ω) is made for all q in the irreducible bz and written to seq.ext.
  • band       A(ω), Σ(ω) are made for qp along symmetry lines and written to spq.ext.
             Use this mode to draw interacting energy bands, in conjunction with plbnds −sp
             Optional @args are parsed like options of the --band switch.
    Other arguments:
  • ib=list         Sum together A^j(ω) derived from QP states j in list.
  • ibx=list       ibx and ib perform the same function, but ibx writes out more information, e.g. A^j(ω) is resolved by band, writing each A^j(ω) in succession.
  • getev         Do not interpolate QP energy but calculate it at each q where Σ and A are to be evaluated.
  • getev=#1,#2,#3   Generate evals on independent mesh with #1,#2,#3 divisions of uniformly spaced points. E^QP may have to be interpolated if the mesh does not contain q.
  • a2qp         Extract the QP energy from the peak in $A(\omega)$ and write it for the one-particle energy when writing spq.ext.
  • 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).
  • wx=#         Interpolate Σ(ω) and A(ω) to frequency  #  and write to disk only data for this frequency (for Fermi surfaces).

• 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). Writes to pes.ext or pes2.ext.
  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 q_n, from list in QIBZ. Alternatively specify q by:
  • q=#1,#2,#3           q-point in units of 2π/alat. lmfgws will interpolate Σ(q_n) 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+V₀
  • range=#1,#2   Generate spectral function in a specified energy window (#1,#2)

• qsgwh
  Generates the Quasiparticle “self-energy” (in practice the QP levels relative to the Fermi level, treated as δ-functions in energy)

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

Editor Application : Compare interacting and independent-particle density-of-states in Fe

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

Once you have created file se.fe, do the following:

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 above):

• 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 (maybe a minute). You can use MPI, e.g. mpirun -n 12 lmfgws ... to speed up the calculation.
• Since the ω 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 A₀(ω), respectively.
• A₀(ω) 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 an already generated DOS from the build directory to your working directory.

cp $toplevel/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}

Draw the figure with

$ 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 A₀(ω), 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 G₀(k,ω) by sampling with a smearing of 30 meV.
• The noninteracting DOS at the Fermi level is D(E_F)≅1/eV (one spin). The Stoner criterion for the onset of ferromagnetism is I×D(E_F)>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.

Editor Application : Spectral Function of Fe near the H point

This example computes the self-energy for a q point near (1,1,0) — 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=5@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=5,6 nw=10 getev=12 isp=2~savesea~q'

The first command writes a file seia.fe, the second seia2.fe To interpret the command-line argument, refer to the editor manual.

The following makes a picture comparing A (solid lines) and A₀ (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]

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

Editor Application : Interacting joint Density-of-States and Optics

lmfgws can make the joint density-of-states (JDOS) and the macroscopic dielectric function. The joint DOS is given by $$D(\omega ) = \sum\nolimits_{\mathbf{k}} {\int {d\omega'[f(\omega) - f(\omega' + \omega )]\sum\nolimits_{ij} {\{ A_i ({\mathbf{k}},\omega' + \omega )A_j ({\mathbf{k}},\omega')\} } } }$$

Note that $D$ is a (weak) function of temperature since the Fermi function $f(\omega)$ contains temperature.

In the limit of noninteracting particles and $k_BT{\rightarrow}0$ this expression reduces to the standard expression for joint density-of-states $$D(\omega) \rightarrow \sum\limits_{\mathbf{k}} {\sum\limits_{i=\mathrm{occ}\atop j=\mathrm{unocc}} {\delta (\omega _i - \omega _j - \omega )} }$$

The following computes joint DOS (noninteracting case) using the lmf optics package. It renames the file for future comparison.

lmf -vnkabc=32 fe `cat switches-for-lm` -vlteto=0 -voptmod=-1 --quit=rho
cp jdos.fe jdos-lmf.fe

Note: you can use MPI for this step.

The following computes joint DOS for both static and interacting QS_GW_ self-energies, using lmfgws.

lmfgws fe `cat switches-for-lm` '--sfuned~units@eV~readsek~eps@.040~jdos@range=-10,10@nq=32@a0@nw=5~savesea~q'
lmfgws fe `cat switches-for-lm` '--sfuned~units@eV~readsek~eps@.040~jdos@range=-10,10@nq=32@nw=5~savesea~q'

The first command makes file jdosni.fe, the second jdos.fe.
Note: the jdos option can be run with MPI.

The following will make a postscript file, with frequency on the abscissa in eV:

fplot -frme 0,.7,0,.5 -x 0,10 -ab 'x1*13.6' -colsy 2,3 -ord y/13.6 jdos-lmf.fe -lt 2,col=1,0,0 -colsy 2,3 jdosni.fe -lt 3,bold=5,col=0,1,0 -colsy 2,3 jdos.fe

Black (joint DOS computed by lmf) and Red (noninteracting joint DOS by lmfgws) are very similar. Dotted green is the corresponding joint DOS with the dynamical self-energy. There is a strong reduction of order $Z^2$ because of loss of quasiparticle weight in the coherent part of $A(\mathbf{k},\omega)$

Optics are very similar the joint DOS. In the absence of a vertex, $\mathrm{Im} \epsilon(\omega)$ is proportional to the joint DOS, decorated by the matrix elements of velocity operator, $\lvert\langle {v_{ij}}\rangle\rvert^2$. The latter is usually calculated in terms of the momentum operator $\lvert\langle {p_{ij}}\rangle\rvert^2$. In Ry units $\mathrm{Im} \epsilon(\omega)$ reads $$\mathrm{Im} \epsilon(\omega ) = \frac{32\,\pi^2}{\omega^2 V} \sum\nolimits_{\mathbf{k}} {\sum\nolimits_{ij} \lvert\left< p_{ij}\right>\rvert^2 \int {d\omega'[f(\omega) - f(\omega' + \omega )] {\{ A_i ({\mathbf{k}},\omega' + \omega )A_j ({\mathbf{k}},\omega')\} } } }$$

The following makes $\mathrm{Im}\epsilon(\omega)$ using lmf with gaussian sampling integration.

lmf -vnkabc=32 fe `cat switches-for-lm` -vmefac=2 -vlteto=0 -voptmod=1 --quit=rho
cp opt.fe opt-lmf.fe

Note: lmf optics can be run with MPI.

The following computes $\mathrm{Im}\,\epsilon(\omega)$ for both static and interacting QSGW self-energies, using lmfgws.

lmf -vnkabc=32 fe `cat switches-for-lm` -vlteto=0 -voptmod=1 -vmefac=2 --quit=rho --opt:woptmc
lmfgws fe `cat switches-for-lm` '--sfuned~units@eV~readsek~eps@.040~imeps@range=-10,10@nq=32@a0@nw=5~savesea~q'
lmfgws fe `cat switches-for-lm` '--sfuned~units@eV~readsek~eps@.040~imeps@range=-10,10@nq=32@nw=5~savesea~q'

Note: these commands can be run with MPI.

lmf generates and stores matrix elements of the velocity operator in file optdatac.fe for lmfgws. The latter commands generate optni.fe and opt.fe. They are calculated in the same way as in the independent particle case, but for spectral functions from the static and interacting QSGW self-energies.

Draw a picture of the three independent calculations of Im ε:

fplot -x 0,10 -frme 0,.7,0,.5 -frmt th=3,1,1 -xl "~{w} (eV)" -x 0,10 -y 0,35 -ab 'x1*13.6' -colsy 2,5 opt-lmf.fe -lt 2,col=1,0,0 -colsy 2,5 optni.fe -lt 3,bold=5,col=0,1,0 -colsy 2,5 opt.fe

Note that lmf uses Ry units; we specified eV in the lmfgws instruction. Thus when comparing Im ε, the abscissa for opt-lmf.fe is scaled to eV.

You should see something similar to the figure shown below. For all three data, contributions are resolved into majority and minority parts. The physically relevant Im ε is the sum of the two.

Black (Im ε computed by lmf) and Red: (Im ε computed by lmfgws, noninteracting case) are very similar. Dotted green is the corresponding Im ε computed in the RPA with the dynamical self-energy. There is a strong reduction in intensity because of loss of quasiparticle weight in the coherent part of $A(\mathbf{k},\omega)$.

Editor Application : 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 $G_0$. Peaks in the spectral function correspond to the band structure; the plot can be compared directly to the bands of the noninteracting G₀. 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 $toplevel/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 (explained in more detail in editor instructions):

• 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(\mathbf{k},\omega)$ along symmetry lines given in file syml2.fe. Include bands 1-10, and generate $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~drawqp 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. The yellow line is the QP band structure from the noninteracting QSGW self-energy. By construction it falls at the peak of the interacting band structure. Blue shows the interacting band structure. Broadening goes to zero at the Fermi level (as it should, according to Fermi liquid theory). In the conduction band, broadening stays small, but in the valence band the broadening increases much more rapidly.

Editor Application : Generating spectral functions from Brillouin zone unfolding of a supercell

This application of lmfgws has its own tutorial. See the Additional excercises portion of this tutorial.