Making the dynamical GW self energy
Table of Contents
 Table of Contents
 Preliminaries
 Command summary
 Introduction
 Theory
 Make the GW selfenergy
 The selfenergy postprocessor
 Dynamical selfenergy editor
 Compare interacting and independentparticle densityofstates in Fe
 Spectral Function of Fe near the H point
 Simulation of photoemission midway between the Γ and H points
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 GW executables are in ~/bin/code2.
Command summary
Repeat the steps for LDA selfconsistency and QSGW selfconsistency 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 selfenergy Σ^{0} (file sigm), make sure your charge density is selfconsistent.
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 selfenergy step:
lmgwsc wt code2 sym metal tol=1e5 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)
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 lmfgwsreadable 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 selfenergy was made with the following command:
$ lmgwsc wt code2 sym metal tol=1e5 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 function at q=0 directly from the GW output files SEComg.UP and SEComg.DN. (For nonmagnetic calculations, only SEComg.UP is made).
 Use lmfgws to generate the interacting densityofstates (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
Theory
Z factor renormalization
Begin with a noninteracting Green’s function G_{0}, defined through an hermitian, energyindependent exchangecorrelation 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_{xc}^{j} 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 \Sigma (k,\omega )/\partial \omega } \right_{\omega ^j } .$defines the Z factor. The dependence of ω^{j} and Z^{j} on k is suppressed.
Define the QP peak as the value of ω 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ω/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 + 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 Z^{j} on k. Write dω^{j}/dk as $v_0^j$, and use the definition of Z 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 selfenergy, Σ(ω).
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\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(\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}}$Make the GW selfenergy
The 1shot GW selfenergy 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 selfenergy 
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 selfenergy is made. Do:
$ lmgwsc wt code2 sym metal tol=1e5 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_{n}, 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 selfenergy.
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.
The selfenergy postprocessor
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 which lmfgws can read.
It also has some limited ability to make spectral functions directly, as described next.
Generate spectral functions for q=0
SEComg.UP and SEComg.DN contain the diagonal matrix element Σ_{jj}(k,ω) 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, Σ_{ii}(k,ω)=0 so the spectral function would be proportional to δ(ω−ω^{*}), where ω* is the QP level (see Theory section).
Interactions give Σ_{ii}(k,ω) an imaginary part which broadens out the level, and in general, ReΣ_{ii}(k,ω) shifts and renormalizes the quasiparticle weight by Z. As noted in the Theory section, there is no shift if V_{xc}^{j} is the QSGW selfenergy; there remains, however, a reduction in the quasiparticle weight. This will be apparent when comparing the interacting and noninteracting DOS.
The spectral tool has a limited ability to convert raw files SEComg.{UP,DN} into spectral functions, which this section demonstrates.
Do the following:
$ spectral eps=.005 domg=0.003 'cnst:iq==1&eqp>10&eqp<30'
Commandline arguments have the following meaning:

–eps=.005 : 0.005 eV is added to the imaginary part of the selfenergy. 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 (kpoint 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 k_{n}. A sec files takes the following format:
# ib= 5 iq= 1 Eqp= 0.797925 q= 0.000000 0.000000 0.000000
# omega omegaEqp Re sigmvxc Im sigmvxc int A(w) int A0(w) A(w) A0(w)
0.2721160D+02 0.2641368D+02 0.6629516D+01 0.1519810D+02 0.2350291D04 0.6897219D08 0.7774444D02 0.2281456D05
0.2720858D+02 0.2641065D+02 0.6629812D+01 0.1520157D+02 0.4701215D04 0.1379602D07 0.7776496D02 0.2281979D05
...
spectral also makes the kintegrated 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 qintegrated 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 qintegrated dos to file dos.dn ...
Setup for lmfgws
Starting from SEComg.UP (and SEComg.DN in the magnetic case) generated by the GW code, use spectral to generate se. This is described in the next section.
Dynamical selfenergy editor
lmfgws is the dynamical selfenergy editor, which peforms a variety of manipulations of the GW selfenergy Σ(k_{n},ω) for different purposes. It 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 selfenergy se.ext in special a format written by the spectral tool.
For definiteness this section assumes that ext is fe. Starting from SEComg.UP (and SEComg.DN in the magnetic case) generated by the GW code, use spectral to generate se.fe:
$ spectral ws nw=1
$ ln s se se.fe
 ws tells spectral to write the selfenergy 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 selfenergy on the frequency mesh it was generated; no interpolation takes place.
Try starting the dynamical selfenergy editor:
$ lmfgws ctrl.fe `cat switchesforlm` 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 Options prompt waiting for another instruction. The editor will print a short summary of instructions if you type ? <RET>.
Editor instructions
The following summarizes the instruction set of the dynamical selfenergy editor.
 readsek [fn]
reads the selfenergy from an ASCII file. In the absence fn, the file name defaults to se.ext. 
readsekb [fn]
reads the selfenergy from a binary file. In the absence fn, the file name defaults to seb.ext.  units Ry
select Rydberg units for which data is to be output 
units eV
select electron vole units for which data is to be output 
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}(k_{n},ω) 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=nq=#1,#2,#3 ib=list [getev[=#1,#2,#3]] [nw=ndomg=#] [isp=#] [range=#1,#2]
Make Σ(ω) and A(ω) for given q and range of bands.
Required arguments are: iq=n index to q_{n}, from list in QIBZ. Alternatively specify q by:
 q=#1,#2,#3 qpoint in units of 2π/alat. lmfgws will interpolate Σ(q_{n}) to any q.
 ib=list Sum together A^{j}(ω) 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.
 range=#1,#2 Generate spectral function in a specified energy window (#1,#2)
 pepeqp iq=nq=#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).
Required arguments are: iq=n index to q_{n}, from list in QIBZ. Alternatively specify q by:
 q=#1,#2,#3 qpoint 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_{0}
 range=#1,#2 Generate spectral function in a specified energy window (#1,#2)

savesea [fn]
saves spectrum DOS or selfenergy + spectral function, in ASCII format. In the absence fn, the file name defaults to seia.ext or seia2.ext when writing band and kresolved spectral functions (se or pe) and to sdos.ext or sdos2.ext when writing spectrum dos (dos). 
savese [fn]
saves qinterpolated selfenergy + 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 a batch mode by stringing instructions together:
$ lmfgws ctrl.fe `cat switchesforlm` 'sfuned~first command~second command~...'
~ is the delimiter separating instructions (you can use another character). lmfgws will parse through all the commands given, and return to the Option prompt, unless the editor encounters “quit” instruction a
or q
. See the next section for an example.
Compare interacting and independentparticle densityofstates in Fe
This section uses the selfenergy editor, lmfgws, to interpolate Σ(k_{n},ω) to a fine k and ω mesh to obtain a reasonably well converged densityofstates.
$ lmfgws fe `cat switchesforlm` '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 asgiven 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.
Note: for now, copy the lmfgenerated DOS to your working directory.
cp ~/lm/gwd/test/fe/dosp.fe dosp.fe
The following uses the fplot tool to draw a figure comparing the DOS generated the three different ways. 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,{ymaxdy},{ymax} p0 x {emin},{emax} y {dmin},{dmax} tmy 1 1p
colsy 3 lt 1,bold=3,col=.5,.5,.5 qr sdos.fe
colsy 2 lt {ltdos} ord y qr dosp.fe
colsy 2 lt 1,bold=3,col=1,0,0 qr 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_{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.
Spectral Function of Fe near the H point
This example computes the selfenergy 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 switchesforlm` '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 switchesforlm` '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 6 lt 1,col=0,0,0 seia.fe colsy 7 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 can see a weak plasmon peak near −8 eV.
Simulation of photoemission midway between the Γ and H points
This test simulates an ARPES measurement for a point approximately midway between Γ and H, at q≅0,0,0.45 near where bands cross the Fermi level. It is done in two ways:
 with the interacting spectral function, but without SO coupling
 with The QP spectral function broadened by 0.01 eV, and including SO coupling.
$ lmfgws fe `cat switchesforlm` 'sfuned~units=eV~eps .01~readsek~evsync~pe q=0,0,0.450980392 ib=2,3 nw=10 getev nqf=220 ke0=1395+14.8 isp=2~savesea~q'
$ lmfgws fe vso=1 `cat switchesforlm` 'sfuned~units=eV~eps .01~qpse~evsync~peqp q=0,0,0.450980392 ib=1:8 nw=10 getev nqf=220 ke0=1395+14.8 isp=1~savesea~q'
The first command writes a file pes2.fe, the second pesqp.fe.
… need to complete
{::comment}
se ib=2
\# ib=2 qp = 0.050000 0.910000 0.010000 eps=0.01000 eqp=2.508988
\# omega Re sigmvxc Im sigmvxc int A(w) int A0(w) A(w) A0(w)
se ib=2,3
\# ib=2,3 qp = 0.050000 0.910000 0.010000 eps=0.01000
\# omega A(w) A0(w)
Edit This Page