# 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 . 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 carries information about the inverse lifetime (the time for decay of a state to another state).

is in general nonlocal in both space and time. Nonlocality in time implies that depends on two time coordinates and ; however in equilibrium it depends only on the difference . Converting time to its Fourier (frequency) representation, depends only on a single frequency . Nonlocality in space implies that depends on two coordinates and . Translational invariance of a crystal implies that for translations between one unit cell and another, depends only on one (lattice) translation vector . With a Fourier from to representation, depends only on one vector. For points within a unit cell, depends on two coordinates and confined to a unit cell. Thus in general is written as .

The *GW* approximation indeed makes a fully nonlocal , and this tutorial described how to construct and analyze it. The main utility used to analyze is **lmfgws**. Dynamical Mean Field Theory (DMFT) is a single-site approximation; it thus makes 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 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, is nonlocal in time but local in space. It operates as a stand-alone package.

### Table of Contents

- Table of Contents
- Preliminaries
- Command summary
- Introduction
- Theory
- Make the GW dynamical self-energy
- spectral, the self-energy translator
- Dynamical self-energy editor lmfgws
- Compare interacting and independent-particle density-of-states in Fe
- Spectral Function of Fe near the H point
- Interacting joint Density-of-States and Optics
- Interacting band structure

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

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

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

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

- Generate spectral functions directly from files
*SEComg.UP*and*SEComg.DN*. (For nonmagnetic calculations, only*SEComg.UP*is made). - Use
**lmfgws**to generate the interacting density-of-states (DOS) from Im*G*, compare it to the noninteracting DOS from Im 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

### Theory

#### Z factor renormalization

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

The contribution to from QP state is

where is the pole of .

Write the contribution to from QP state as

Note that this equation is only true if is diagonal in the basis of noninteracting eigenstates. We will ignore the nondiagonal elements of . Note that if is constructed by QSGW, this is a very good approximation, since at . Approximate *G* by its coherent part:

where

defines the factor. The dependence of and on is suppressed.

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

and so

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 . For the interacting case it reads

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

Ignore the dependence of on . Write as , and use the definition of to get

So

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

#### Coherent part of the spectral function

Write as

Rewrite as

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

the approximate spectral function is

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

#### Simulation of Photoemission

(needs cleaning up)

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

where *E _{b}* is the binding energy and is the energy of the electron after being ejected. (Marder defines 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 be the energy on exiting the crystal, the work function and and are called the electron binding energy and “inner potential.”

Then

The total momentum inside the crystal, , 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 where is the work function of the analyzer. Usually . The Fermi level is defined such that . The inner potential is defined by scanning the range of photon energy under the constraint of normal emission: then the -point can be identified and by using Eq.(1), and the inner potential experimentally determined.

The momentum of the particle in free space is

Resolve into components parallel and perpendicular to the surface

After passing through the surface, is modified to ; this is what is actually measured.

The conservation condition requires

is conserved on passing through the surface; thus . is not conserved; therefore

The wave number shift is then

and the crystal momentum actually being probed by the experiment is

### 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 QS*GW* 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 --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 Σ(**k*** _{n}*,

*ω*) on a uniform energy mesh −2 Ry <

*ω*< 2 Ry, spaced by 0.01 Ry at irreducible points

**k**

*, for QP levels specified in*

_{n}*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
$(dirname $(which lmgw))/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 for each QP level *j*, for each irreducible point **k*** _{n}* 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, so the spectral function would be proportional to δ(

*ω*−

*ω*

^{*}), where

*ω** is the QP level (see Theory section).

Interactions give an imaginary part which broadens out the level, and in general, shifts and renormalizes the quasiparticle weight by *Z*. As noted in the Theory section, there is no shift if is the QS*GW* self-energy ; 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 Σ(**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:**: acts as a constraint to exclude entries in*expr**SEComg.{UP,DN}*for whichis zero.*expr*

is an integer expression using standard Questaal syntax for algebraic expressions. It can that can include the following variables:*expr*- 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_ib*j*_iq*n*.up* and *sec_ib*j*_iq*n*.dn* which contain information about the *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
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 performs a variety of postprocessing of the *GW* or DMFT self-energy for different purposes. The collection of points **k*** _{n}* 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 QS*GW* 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, the file name defaults to the ASCII file*fn**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:- 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. - QP levels are stored relative to the chemical potential (which may, but need not, be written in the header).
- Only the diagonal elements of the potentials are read. The full complement of static potentials consist of the static QS
*GW*self-energy , the Fock exchange , and . - The
*se*file may, but need not, contain these potentials. For example, none are supplied by**lmfdmft**. - 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.**@fn=**use*nam**nam*for self-energy file name**@useef**file chemical potential becomes Fermi level**@irrmesh**points are not on a regular**k**mesh : no**k**interpolations allowed**@ib=**after reading data from file, pare bands read from file to those in*list**list***@minmax**print minimum and maximum QP levels for each band**@dc=#**subtract double counting # from Re sigma(omega) after reading**@makeqpse**Not documented yet

- Data is stored for a collection of

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

replace quasiparticle levels read from*se.ext*by recalculating them with the same algorithm**lmf**uses.**eps***val*

add a constantto Im Σ, needed to broaden spectral functions so that integrations are tractable.*val***ef***Ef*| ef=*Ef*

Usefor the Fermi level, overriding the internally calculated value.*Ef*

*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 enterin Ry units.*Ef*

**dos|jdos|imeps [nq=#1,#2,#3] ib=***list*[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=**restrict contribution to spectra from QP states in*lst**list*.**nq=#1,#2,#3**Interpolate Σ(_{j}**k**,_{n}*ω*) to a new uniform mesh of**k**points, defined by (**#1,#2,#3**) divisions. Use between 1 and 3 numbers.**nw=**Refine the given energy mesh by interpolating Σ to an*n**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=**index to*n***q**, from list in_{n}*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**) to any_{n}**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 required arguments:**ib=**Sum together*list**A*(^{j}*ω*) derived from QP states*j*in*list*.

**ibx=**is similar to*list***ib=**, but*list**A*(^{j}*ω*) is resolved by band, writing each*A*(^{j}*ω*) in succession.

*Note:*with**ib**,**lmfgws**will … XXX Options are:**getev**Do not interpolate QP energy but calculate it at**q**.**a2qp**Extract the QP energy from the peak in and write it for the one-particle energy when writing*spq.ext*.**getev=#1,#2,#3**Generate evals on independent mesh with**#1,#2,#3**divisions of uniformly spaced points.**nw=**Refine the given energy mesh by interpolating Σ to an*n**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). 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=**index to*n***q**, from list in_{n}*QIBZ*. Alternatively specify**q**by:**q=#1,#2,#3****q**-point in units of 2*π*/alat.**lmfgws**will interpolate Σ(**q**) to any_{n}**q**.**ib=**Sum together PE spectrum derived from QP states*list**j*in. See here for the syntax of integer lists.*list*

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=**Refine the given energy mesh by interpolating Σ to an*n**n*multiple of the given energy mesh.*n*must be an integer.**isp=**Generate spectra for spin*i*(1 or 2). Default value is 1.*i***nqf=**number of mesh points for final state integration. Default is 200.*n***ke0=#**kinetic energy of emitted electron. KE+V0=ℏ*ω−φ*_{s}+V_{0}**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, the file name defaults to*fn**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, the file name defaults to*fn**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.

- Batch mode
- 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 Σ(**k*** _{n}*,

*ω*) 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*A*_{0}(*ω*), respectively. *A*_{0}(*ω*) 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}
```

```
$ 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*_{0}(*ω*), 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*_{0}(**k**,*ω*) by sampling with a smearing of 30 meV. - The noninteracting DOS at the Fermi level is
*D*(*E*)≅1/eV (one spin). The Stoner criterion for the onset of ferromagnetism is_{F}*I*×*D*(*E*)>1, where_{F}*I*is the Stoner parameter, which DFT predicts to be approximately 1 eV for 3*d*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 *A*_{0} (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.

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

Note that is a (weak) function of temperature since the Fermi function contains temperature.

In the limit of noninteracting particles and this expression reduces to the standard expression for joint density-of-states

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

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

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

```
fplot -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
```

In the absence of a vertex, is proportional to the joint DOS, decorated by the matrix elements of velocity operator, . The latter is usually calculated in terms of the momentum operator . In Ry units reads

The following makes using **lmf** with gaussian sampling integration.

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

The following computes for both static and interacting QS*GW* self-energies, using **lmfgws**.

```
lmf -vnk=32 fe `cat switches-for-lm` -vlteto=0 -voptmod=1 -vmefac=0 --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'
```

The **lmf** instructions generates 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, but for spectral functions from the static and interacting QS*GW* self-energies.

Draw a picture of the three independent calculations of :

```
fplot -frme 0,1,0,.7 -frmt th=3,1,1 -xl "~{w} (eV)" -x 0,10 -y 0,32 -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 , the abscissa from for *opt-lmf.fe* must be scaled (alternatively tell **lmfgws** to use Ry units).

You should see something similar to the figure shown below (the figure shown is smoother because **nq=64** divisions were used for the *k* mesh). For all three data, contributions are resolved into majority and minority parts. The physically relevant is the sum of the two.

Black ( computed by **lmf**) and Red ( computed by **lmfgws**, noninteracting case) are very similar. Dotted green is the corresponding computed in the RPA with the dynamical self-energy. There is a strong reduction of order because of loss of quasiparticle weight in the coherent part of

### 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 . Peaks in the spectral function correspond to the band structure; the plot can be compared directly to the bands of the noninteracting *G*_{0}. 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 along symmetry lines given in file*syml2.fe*. Include bands 1-10, and generate 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*.