Spectral functions for LSCO

This tutorial explains how to generate densities-of-states and k-resolved spectral functions $A(\mathbf{k},\omega)$ from a local dynamical self-energies $\Sigma(\omega)$ generated by DFMT.

It has close similarities to generating spectral functions from dynamical self-energies $\Sigma(\mathbf{k},\omega)$ generated by the GW code; see for example this tutorial.

Table of Contents

• Table of Contents
• Command summary
• Preliminaries
• Interacting Energy bands
• Spectrum DOS

Command summary

touch ctrl.lscoq
rm -f *.lscoq
cp ~/lm/dmft/test/lscoq/{{ctrl,basp,sigm,rst,syml,syml2,indmfl&125;.lscoq,sig.inp%125; .

lmfdmft lscoq --rs=1,0 --ldadc=82.2 -job=1 --pade~nw=121~window=-4/13.606,4/13.606~icut=30,75
lmchk lscoq | grep -A 3 Qlat | tail -3 | awk '{print $4, $5, $6}' > qlat
~/lm/misc/symlinepoints -qlat=qlat -nk=101 syml.lscoq > syml2.lscoq
lmfdmft lscoq -vnkabc=8 --ldadc=82.2 --job=1 --gprt~band,fn=syml2~rdsigr=sig2~mode=19 --pr45
lmfgws  lscoq '--sfuned~units eV~readsek@useef@irrmesh@minmax@ib=30:38~eps 0.02~se band@fn=syml2 nw=2 isp=1 range=-4,4'
plbnds -sp~atop=10~window=-4,4 spq.lscoq ; gnuplot gnu.plt ; open spf.ps

lmf lscoq -vnkabc=12 --dos@npts=721@ef0@rdm@window=-6/13.606,6/13.606 --quit=dos
lmfdmft lscoq -vnkabc=8 --ldadc=82.2 --job=1 --gprt~rdsigr=sig2~mode=19
lmfgws  lscoq -vnkgw=8 '--sfuned~units eV~readsek@useef@ib=30:38~eps 0.02~dos isp=1 range=-4,4 getev=12 nw=4~savesea~q'

Preliminaries

This tutorial works with paramagnetic La₂CuO₄. QSGW code has no CPA-like capability (in contrast to DMFT, which may be thought of as a dynamical generalization of the Coherent Potential Approximation). At the GW level the material must be treated either as an ordered antiferromagnet, in which case it predicts La₂CuO₄ to be an insulator with a gap of ~4 eV, or nonmagnetically, in which case it predicts La₂CuO₄ to be a metal.

Experimentally, La₂CuO₄ is an insulator with an optical gap slightly less than 2 eV. Both LDA+DMFT and QSGW+DMFT yield an insulator in the paramagnetic state: the $x^2{-}y^2$ orbital (which in the nonmagnetic case forms a continuous band that crosses the Fermi level) gets split off.

This tutorial starts from a static nonlocal QSGW self-energy $\Sigma^0(\mathbf{k})$ and a local, dynamical $\Sigma(i\omega_n)$, projected onto the Cu d states, generated by DMFT for frequencies on the Matsubara axis. The materials system is La₂CuO₄; files are taken from a validation test taken from the Questaal installation.

It tutorial assumes your build directory is ~/lm, and the executables are your path. The test itself can be run by the script ~/lm/dmft/test/test.dmft.
It requires the following:

• lmfdmft^†  is used for Pade extrapolation of Σ to the real axis, and also to embed Σ into the QSGW bath.
• lmchk^†   is used to create file qlat, needed by symlinepoints
• ~/lm/misc/symlinepoints  maps k-point lines in syml.lscoq to syml2.lscoq.
• lmfgws^†   Makes $A(\mathbf{k},\omega)$ from the embedded $\Sigma(\mathbf{k},\omega)$. lmfgws uses the self-energy editor, described here
• plbnds^†    band plot utility, operating in the spectral function mode
• gnuplot^†
• open^†    [use any utility that views postscript files]

^†assumed to be in your path.

Interacting Energy bands

This section parallels the tutorial for the interacting band structure generated by GW. In the GW case, $\Sigma(\mathbf{k},\omega)$ is generated from QSGW on the real axis. Here $\Sigma(\mathbf{k},\omega)$ is generated from a local $\Sigma_d(\omega)$ generated by DMFT on the Matsubara axis.

To generate physical observables, we need to analytically continue $\Sigma(i\omega_n)$, given by DMFT for a set of discrete Matsubara frequencies, to the real axis. This is an important difference between Questaal’s QSGW and DMFT implementations.

In this tutorial we use lmfdmft’s built in Pade extrapolator to carry out the analytic continuation.
Note: analytic continuation is a tricky business. Pade works reliably only when the polynomial order is not too large, which means in can operate reliably only over a small energy window. Also the self-energy should be reasonably smooth.

Continue the local DMFT $\Sigma(i\omega)$ to the real axis in an energy window [−4,4] eV, using a Pade extrapolation:

$ lmfdmft lscoq --rs=1,0 --ldadc=82.2 -job=1 --pade~nw=121~window=-4/13.606,4/13.606~icut=30,75

The switches to --pade are delimited by  ~  which is the first character following --pade:

• ~nw=121
  Create an energy mesh of 121 points
• ~window=−4/13.606,4/13.606
  Use an energy window of [−4,4] eV (specify in Ry units)
• ~icut=30,75
  Specifies which Matsubara points to use in the Pade approximant. Use the first 30 points; thereafter increment the spacing by 1,2,4 … points.
  Continue until the counter passes point 75. In this case the Pade approximant includes Matsubara points 1:30,32,36,44,60

Next, we need a list of q-points where the spectral function $A(\mathbf{k},\omega)$ will be made (band plot). Here we use generic symmetry lines from syml.lscoq, which contain high-symmetry points in units of the reciprocal lattice vectors. The symlinepoints utility will convert these points to Cartesian coordinates and at the same time specify a number of points on each line to make the spacing between points approximately uniform. Do the following:

$ lmchk lscoq | grep -A 3 Qlat | tail -3 | awk '{print $4, $5, $6}' > qlat
$ ~/lm/misc/symlinepoints -qlat=qlat -nk=101 syml.lscoq > syml2.lscoq

symlinepoints writes points in Questaal’s symmetry line format to file syml2.lscoq.

Next, we need to embed the local DMFT self-energy $\Sigma_d(\omega)$ (now analytically continued to the real axis) in the bath defined by the one-particle QSGW hamiltonian. The next step will:

1. Determine the Fermi level, so the one-particle states can be referenced relative to it

2. Embed $\Sigma_d(\omega)$ into $\Sigma(\mathbf{k},\omega)$, in the basis of QSGW eigenfunctions, for each of the k points supplied through syml2.lscoq.

   $ lmfdmft lscoq -vnkabc=8 --ldadc=82.2 --job=1 --gprt~band,fn=syml2~rdsigr=sig2~mode=19 --pr45
   

−gprt tells lmfdmft to operate in a special purpose mode which creates and save some property of the Green’s function or self energy. Switches to –gprt are documented here.

lmfdmft first determines the chemical potential μ from $\Sigma_d$ on the Matsubara axis, avoiding errors in the Pade approximation. Next it makes another band pass for the points defined through syml2.lscoq. For each point it embeds $\Sigma_d(\omega)$ (now continued to the real axis) to construct the crystal $\Sigma(\mathbf{k},\omega)$, for the set of bands defined by the DMFT hybridization function. It writes $\Sigma(\mathbf{k},\omega)$ to se.lscoq.
Note: lmfdmft writes only the diagonal part of Σ to se.lscoq. It is assumed that the off-diagonal elements do not significantly modify the spectral function.

We are now in a position to make spectral functions. It is done by the dynamical self-energy editor lmfgws. $\Sigma(\mathbf{k},\omega)$ constructed either by DMFT and GW can use this utility to analyze spectral functions. For the GW case see this Fe tutorial.

lmfgws can work interactively by adding  --sfuned  unadorned by modifiers:

$ lmfgws  lscoq --sfuned

or it can be invoked in batch mode. We do the latter here. The following instruction will read se.lscoq, interpolate it to a finer frequency mesh, and write $A(\mathbf{k},\omega)$ to file spq.lscoq.

$ lmfgws  lscoq '--sfuned~units eV~readsek@useef@irrmesh@minmax@ib=30:38~eps 0.02~se band@fn=syml2 nw=2 isp=1 range=-4,4'

The editor’s instructions are delimited by the first character following --sfuned, which is  ~  in this case. They perform the following functions:

• units eV
  Select data is to be output in electron volt units

• readsek@flags
  Reads the self-energy from se.lscoq. Optional flags (separated by a delimiter) modify what is read.
  Note: you cannot use  ~  for the flags delimiter because it was already used to separate instructions to --sfuned (batch mode). The flags perform the following functions:

  • useef
    Tells the editor to use the chemical potential in se.lscoq. The QSGW and DMFT Fermi levels are different, and this switch can have some effect if lmfgws is used to interpolate between k points. However, in this context it is not because lmfgws already has the one-particle levels and self-energy relative to the Fermi level for every k point; there is no interpolation (in contradistinction to the DOS calculation
  • irrmesh
    Tells the editor that data is not being read on a regular mesh of points suitable for k integration. Operations that require k integration or k interpolation cannot be performed.
  • minmax
    Instructs the editor to print out the minimum and maximum values of QP levels among all the q points it reads. Values are printed for each band. This is for informational purposes only.
  • ib=30:38
    Restricts the bands entering into the spectral function to bands 30 through 38. These are the bands in the energy window [−4,4] eV, as the minmax printout will show.

• eps 0.02
  add a constant 0.02 eV to Im Σ. this broaden spectral functions that have small imaginary parts near their peak, rendering a cleaner display.

• se   band@fn=syml2   nw=2   range=−4,4
  Instructs the editor to generate spectral functions at the list of k points you specify.

  • band@fn=syml2
    generates a list of k points in lines defined through syml2.lscoq.  @  is a delimiter separating flags to band (In general the band instruction can have multiple flags). By default band operates in symmetry line mode, which read points symmetry lines from syml2.lscoq. Here we synchronize the k points with the same list for which lmfdmft made the self energy. That way $\Sigma$ will not need to be interpolated in q.
  • nw=n
    Refines the given energy mesh by interpolating Σ to an n multiple of the given energy mesh. n must be an integer.
  • range=#1,#2
    Specifies the energy window [#1,#2] where the spectral function is to be generated.

Finally, invoke plbnds in spectral function mode to generate an “interacting energy bands” figure.

$ plbnds -sp~atop=10~window=-4,4 spq.lscoq
$ gnuplot gnu.plt
$ open spf.ps   [choose your postscript file viewer]

plbnds will generate a gnuplot script file gnu.plt together with a data file spf.lscoq. gnuplot makes a postscript file, spf.ps.

The open merely displays spf.ps.

Spectrum DOS

The spectrum total DOS is a k-integrated version of the k resolved spectral function. Here we generate the noninteracting DOS from the QSGW hamiltonian, and compare it to the DMFT spectrum DOS.

Making the QSGW noninteracting DOS using lmf in the usual manner, using the --dos switch

$ lmf lscoq -vnkabc=12 --dos@npts=721@ef0@rdm@window=-6/13.606,6/13.606 --quit=dos

It creates file dos.lscoq. --quit=dos tells lmf to stop after the DOS has been written, and the modifies to --dos have the following meanings:

• npts=721  specifies number of energy points (the spacing between points is the same as for the spectrum dos)
• window=−6/13.606,6/13.606  energy window for DOS, in Ry. It is 3/2 wider than that generated for the spectrum DOS
• ef0     Energy window to be defined relative to Fermi level
• rdm    Write DOS in standard Questaal format for 2D arrays.

To make the spectrum DOS, $\Sigma(\omega)$ is needed on the real axis. The steps below are similar those used earlier in the tutorial for energy bands;

$ lmfdmft lscoq --rs=1,0 --ldadc=82.2 -job=1 --pade~nw=121~window=-4/13.606,4/13.606~icut=30,75
$ lmfdmft lscoq -vnkabc=8 --ldadc=82.2 --job=1 --gprt~rdsigr=sig2~mode=19
$ lmfgws  lscoq -vnkgw=8 '--sfuned~units eV~readsek@useef@ib=30:38~eps 0.02~dos isp=1 range=-4,4 getev=12 nw=4~savesea~q'

The main change to the lmfdfmt command is that --gprt has no band modifier, since we require the self-energy on a regular mesh to perform k-integration. lmfdfmt writes $\Sigma(\mathbf{k},\omega)$ to se.lscoq.

The main change to the lmfgws instruction is that dos is used instead of se in the self-energy editor. Note also the modifier getev=12. This interpolates the QP band structure to a finer mesh (12×12×12 divisions), by recalculating the one-particle levels for the finer mesh. The self-energy is linearly interpolated from the given 8×8×8 mesh to the 12×12×12 mesh. It writes file sdos.lscoq.

To compare the two DOS, you can use the fplot utility :

$ fplot -vef0=.39 sdos.lscoq -ab '(x1-ef0)*13.6' -ord x2/13.6 -lt 2,col=1,0,0 dos.lscoq 

The noninteracting DOS are written in Ry, on an absolute energy scale (i.e. not relative to the Fermi level).
-ab ‘(x1-ef0)*13.6’ shifts by ef0 and scales the abscissa from Ry to eV. -ord x2/13.6 scales the ordinate from Ry⁻¹ to eV⁻¹.

Note that we should use the Fermi level of the noninteracting DOS (0.264 Ry), but we used the Fermi level of the interacting one (0.392 Ry). The mostly noninteracting DOS near 4 eV line up properly that way.