Back to the main Index

Carbon in diamond structure

In this lesson, we discuss how to compute the electron-phonon (e-ph) self-energy and the renormalization of the electronic states due to the e-ph interaction. We start with a very brief discussion of the basic equations and some technical details about the Abinit implementation of the EPH code. Then we show how to build an AbiPy flow to automate the multiple steps required by this calculation.

In the second part, we use the AbiPy objects to analyze the results. We start from the simplest case of a single EPH calculation whose most important results are stored in the SIGEPH.nc file. Then we show how to use SigEPhRobot to handle multiple netcdf files and perform convergence studies. Finally, we present a possible approach to interpolate the QP corrections to obtain temperature-dependent band structures along an arbitrary high-symmetry k-path.

- A bit of Theory
- Building a Flow for EPH self-energy calculations
- Electronic properties
- Vibrational properties
- E-PH self-energy
- Using SigEphRobot for convergence studies
- Convergence study with respect to nbsum and nqibz
- Interpolating QP corrections with star functions

$\newcommand{\kk}{{\mathbf k}}$ $\newcommand{\qq}{{\mathbf q}}$ $\newcommand{\kpq}{{\kk+\qq}}$ $\newcommand{\RR}{{\mathbf R}}$ $\newcommand{\rr}{{\mathbf r}}$ $\newcommand{\<}{\langle}$ $\renewcommand{\>}{\rangle}$ $\newcommand{\KS}{{\text{KS}}}$ $\newcommand{\ww}{\omega}$ $\newcommand{\ee}{\epsilon}$

The electron-phonon matrix elements are defined by:

\begin{equation} %\label{eq:eph_matrix_elements} g_{mn}^{\nu}(\kk,\qq) = \dfrac{1}{\sqrt{2\omega_{\qq\nu}}} \<\psi_{m \kpq} | \Delta_{\qq\nu} V |\psi_{n\kk}\> \end{equation}where the perturbation potential, $\Delta_{\qq\nu} V$, is related to the phonon displacement and the first order derivative of the KS potential:

\begin{equation} \partial_{\kappa\alpha,\qq} v^{KS} = \Bigl ( \partial_{\kappa\alpha,\qq}{v^H} + \partial_{\kappa\alpha,\qq}{v^{xc}} + \partial_{\kappa\alpha,\qq}{v^{loc}} \Bigr ) + \partial_{\kappa\alpha,\qq}{V^{nl}} \ = \partial_{\kappa\alpha,\qq}{v^{scf}} + \partial_{\kappa\alpha,\qq}{V^{nl}} \end{equation}The first term in parentheses depends on the first order change of the density (output of the DFPT run).
The second term derives from the non-local part of the atom-centered (pseudo)potentials
and can be easily computed without having to perform a full DFPT calculation.
Hereafter, we use the term *DFPT potential* to refer to the above expression and
the term *self-consistent DFPT potential* to indicate the contribution that depends on the DFPT density
(the term in parentheses).

Throughout these notes, we shall use Hartree atomic units ($e = \hbar = m_e = 1$).

The e-ph self-energy consists of two terms (Fan-Migdal + Debye-Waller):

\begin{equation} \Sigma^{e-ph}(\ww) = \Sigma^{FM}(\ww) + \Sigma^{DW} \end{equation}where only the first part (FM) is frequency dependent.

The *diagonal* matrix elements of the Fan-Migdal self-energy in the KS basis set are given by:

where $n_{\qq\nu}$ ($T$) are the temperature-dependent occupation factors for phonons (electrons)
and $\eta$ is a small positive number (*zcut* input variable).
These occupation factors introduce a dependence on T in the self-energy and therefore in the QP results.
For the sake of simplicity, the dependence on T will not be explicitly mentioned in the other equations.

The Debye-Waller term is given by:

\begin{equation} \Sigma_{n\kk}^{DW} = \sum_{\qq\nu m} (2 n_{\qq\nu} + 1) \dfrac{g_{mn\nu}^{2,DW}(\kk, \qq)}{\ee_{n\kk} - \ee_{m\kk}} \end{equation}where $g_{mn\nu}^{2,DW}(\kk, \qq)$ is an effective matrix element that, under the rigid-ion approximation, can be expressed in terms of the "standard" $g_{nm}^{\nu}(\kk, \qq=\Gamma)$ elements (see references for a more detailed discussion).

The QP energies are usually obtained with the linearized QP equation:

\begin{equation} \ee^{QP}_{n\kk} \approx \ee_{n\kk} + Z_{n\kk}\,\Sigma^{eph}(\ee_{n\kk}) \end{equation}where the renormalization factor $Z$ is given by:

\begin{equation} Z_{n\kk} = \Bigl ( 1 - \Re\dfrac{\partial\Sigma^{eph}}{\partial{\ee_{KS}}} \Bigr )^{-1} \end{equation}The above expression assumes that the KS states and eigenvalues are relatively close the QP results. A more rigourous approach would requires solving the non-linear QP equation in the complex plane with possibly multiple solutions.

Finally, the spectral function is given by:

\begin{equation} A_{n\kk}(\ww) = \dfrac{1}{\pi} \dfrac {|\Im \Sigma^{eph}_{n\kk}(\ww)|} {\bigl[ \ww - \ee_{n\kk} - \Re{\Sigma^{eph}_{n\kk}(\ww)} \bigr ] ^2+ \Im{\Sigma^{eph}_{n\kk}(\ww)}^2} \end{equation}EPH calculations require very dense BZ meshes to converge and the calculation of the DFPT potentials represent a significant fraction of the overall computing time. To reduce the computational cost, Abinit employs the Fourier-based interpolation technique proposed by Eiguren and Ambrosch-Draxl in Phys. Rev. B 78, 045124.

The DFPT potentials can be interpolated by introducing:

\begin{equation} W(\rr,\RR) = \dfrac{1}{N_\qq} \sum_\qq e^{-i\qq\cdot\RR} \partial v^{scf}_\qq(\rr) \end{equation}This quantity is usually short-ranged and the *interpolated* potential at an arbitrary point $\tilde q$
is obtained via:

In polar materials, particular care must be taken due to the long-range behavior of $W$. We will not elaborate on this point further but dynamical quadrupoles are important in diamond.

For the previous discussion, it should be clear that the evaluation of the e-ph self-energy requires:

An initial set of KS wavefunctions and eigenvalues including

**several empty**statesThe knowledge of the dynamical matrix $D(\qq)$ from which one can obtain phonon frequencies and displacements everywhere in the BZ via FT interpolation.

A set of DFPT potentials in the IBZ.

Thanks to these three ingredients, the EPH code can compute the e-ph matrix elements and evaluate the self-energy matrix elements. A schematic representation of a typical EPH workflow for ZPR computations is given in the figure below:

The `mrgddb`

and `mrgdv`

are Fortran executables whose main goal is to *merge* the partial files
obtained at the end of a single DFPT run for given q-point and atomic perturbation
and produce the final database required by the EPH code.

Note that, for each quasi-momentum $\kk$, the self-energy $\Sigma_{n\kk}$ is given by an integral over the BZ. This means that:

The WFK file must contain the $\kk$-points and the bands where the QP corrections are wanted as well as enough empty state to perform the summation over bands.

The code can (Fourier) interpolate the phonon frequencies and the DFPT potentials over an arbitrarily dense q-mesh. Still, for fixed $k$-point, the code requires the ab-initio KS states at $\kk$ and $\kpq$.

**Providing a WKF file on a k-mesh that fulfills this requirement is responsability of the user**.

Since we are computing the contribution to the electron self-energy due to the interaction with phonons,
it is not surprising to encounter variables that are also used in the $GW$ code.
More specifically, one can use the `nkptgw`

, `kptgw`

, `bdgw`

variables to select
the electronic states for which the QP corrections are wanted (see also `gw_qprange`

)
while `zcut`

defines for the complex shift.
If `symsigma`

is set to 1, the code takes advantage of symmetries to reduce the BZ integral to an appropriate
irreducible wedge defined by the little group of the k-point (calculations at $\Gamma$ are therefore
much faster than for other low-symmetry k-points).

You might find additional material, related to the present section, in the following references:

Before starting, we need to import the python modules used in this notebook:

In [35]:

```
# Use this at the beginning of your script so that your code will be compatible with python3
from __future__ import print_function, division, unicode_literals
import numpy as np
import warnings
warnings.filterwarnings("ignore") # to get rid of deprecation warnings
from abipy import abilab
abilab.enable_notebook() # This line tells AbiPy we are running inside a notebook
import abipy.flowtk as flowtk
# This line configures matplotlib to show figures embedded in the notebook.
# Replace `inline` with `notebook` in classic notebook
%matplotlib inline
# Option available in jupyterlab. See https://github.com/matplotlib/jupyter-matplotlib
#%matplotlib widget
```

and a function from the `lesson_sigeph`

module that will build our AbiPy flow.

In [36]:

```
from lesson_sigeph import build_flow
```

In [37]:

```
abilab.print_source(build_flow)
```

Out[37]:

```
def build_flow(options):
"""
C in diamond structure. Very rough q-point mesh, low ecut, completely unconverged.
The flow computes the ground state density and a WFK file on a 8x8x8 k-mesh including
empty states needed for the self-energy. Then all the independent atomic perturbations
for the irreducible qpoints in a 4x4x4 grid are obtained with DFPT.
Finally, we enter the EPH driver to compute the EPH self-energy.
"""
workdir = options.workdir if (options and options.workdir) else "flow_diamond"
# Define structure explicitly.
structure = abilab.Structure.from_abivars(
acell=3*[6.70346805],
rprim=[0.0, 0.5, 0.5,
0.5, 0.0, 0.5,
0.5, 0.5, 0.0],
typat=[1, 1],
xred=[0.0, 0.0, 0.0, 0.25, 0.25, 0.25],
ntypat=1,
znucl=6,
)
# Initialize input from structure and norm-conserving pseudo provided by AbiPy.
gs_inp = abilab.AbinitInput(structure, pseudos="6c.pspnc")
# Set basic variables for GS part.
gs_inp.set_vars(
ecut=25.0, # Too low, shout be ~30
nband=4,
tolvrs=1e-8,
)
# The kpoint grid is minimalistic to keep the calculation manageable.
# The q-mesh for phonons must be a submesh of this one.
gs_inp.set_kmesh(
ngkpt=[8, 8, 8],
shiftk=[0.0, 0.0, 0.0],
)
# Build new input for NSCF calculation with k-path (automatically selected by AbiPy)
# Used to plot the KS band structure and interpolate the QP corrections.
nscf_kpath_inp = gs_inp.new_with_vars(
nband=8,
tolwfr=1e-16,
iscf=-2,
)
nscf_kpath_inp.set_kpath(ndivsm=10)
# Build another NSCF input with k-mesh and empty states.
# This step generates the WFK file used to build the EPH self-energy.
nscf_empty_kmesh_inp = gs_inp.new_with_vars(
nband=210, # Too low. ~300
nbdbuf=10, # Reduces considerably the time needed to converge empty states!
tolwfr=1e-16,
iscf=-2,
)
# Create empty flow.
flow = flowtk.Flow(workdir=workdir)
# Register GS + band structure parts in the first work
work0 = flowtk.BandStructureWork(gs_inp, nscf_kpath_inp, dos_inputs=[nscf_empty_kmesh_inp])
flow.register_work(work0)
# Generate Phonon work with 4x4x4 q-mesh
# Reuse variables from GS input and let AbiPy handle the generation of the input files
# Note that the q-point grid is a sub-grid of the k-point grid (here 8x8x8)
ddb_ngqpt = [4, 4, 4]
ph_work = flowtk.PhononWork.from_scf_task(work0[0], ddb_ngqpt, is_ngqpt=True,
tolerance={"tolvrs": 1e-6}) # This to speedup DFPT
flow.register_work(ph_work)
# Build template for self-energy calculation. See also v8/Input/t44.in
# The k-points must be in the WFK file
#
eph_inp = gs_inp.new_with_vars(
optdriver=7, # Enter EPH driver.
eph_task=4, # Activate computation of EPH self-energy.
ddb_ngqpt=ddb_ngqpt, # q-mesh used to produce the DDB file (must be consistent with DDB data)
nkptgw=1,
kptgw=[0, 0, 0],
bdgw=[1, 8],
# For more k-points...
#nkptgw=2,
#kptgw=[0, 0, 0, 0.5, 0, 0],
#bdgw=[1, 8, 1, 8],
tmesh=[0, 200, 5], # (start, step, num)
zcut="0.2 eV", # Too large. needed to get reasonable results due to coarse q-mesh.
nfreqsp=301, # Compute A(w)
)
# Set q-path for Fourier interpolation of phonons.
eph_inp.set_qpath(10)
# Set q-mesh for phonons DOS.
eph_inp.set_phdos_qmesh(nqsmall=16, method="tetra")
# EPH part requires the GS WFK, the DDB file with all perturbations
# and the database of DFPT potentials (already merged by PhononWork)
deps = {work0[2]: "WFK", ph_work: ["DDB", "DVDB"]}
# Now we use the EPH template to perform a convergence study in which
# we change the q-mesh used to integrate the self-energy and the number of bands.
# The code will activate the Fourier interpolation of the DFPT potentials if eph_ngqpt_fine != ddb_ngqpt
for eph_ngqpt_fine in [[4, 4, 4], [8, 8, 8]]:
# Create empty work to contain EPH tasks with this value of eph_ngqpt_fine
eph_work = flow.register_work(flowtk.Work())
for nband in [100, 150, 200]:
new_inp = eph_inp.new_with_vars(eph_ngqpt_fine=eph_ngqpt_fine, nband=nband)
eph_work.register_eph_task(new_inp, deps=deps)
flow.allocate()
return flow
```

OK, the function is a little bit long but it is normal as we are computing
in a single workflow the *electronic* properties, the *vibrational* spectrum
and the *electron-phonon* self-energy by varying the number of q-points and the
number of empty states.

Note that we have already encountered similar flows in the previous AbiPy lessons.
The calculation of electronic band structures is
discussed in
lesson_base3,
an example of `Flow`

for phonon calculations is given in
lesson_dfpt.
An AbiPy Flow for the computation of $G_0W_0$ corrections is discussed in
lesson_g0w0.

The novelty is represented by the generation of the `EphTasks`

in which we have to specify several variables related to phonons and the e-ph self-energy.
For your convenience, we have grouped the variables used in the `EphTask`

in sub-groups:

Hopefully things will become clearer if we build the flow and start to play with it:

In [38]:

```
flow = build_flow(options=None)
flow.show_info()
```

`EphTask`

.
Please read carefully the documentation of the variables, in particular
those in the `vargw`

, `varrf`

and `vareph`

sections.

In [39]:

```
print(flow[-1][-1])
flow[-1][-1].input
```

<EphTask, node_id=242751, workdir=flow_diamond/w3/t2>

Out[39]:

##############################################

#### SECTION: basic

##############################################

ecut 25.0

nband 200

tolvrs 1e-08

ngkpt 8 8 8

kptopt 1

nshiftk 1

shiftk 0.0 0.0 0.0

##############################################

#### SECTION: eph

##############################################

eph_task 4

ddb_ngqpt 4 4 4

tmesh 0 200 5

ph_ndivsm 10

ph_nqpath 12

ph_qpath

0.0 0.0 0.0

0.5 0.0 0.5

0.5 0.25 0.75

0.375 0.375 0.75

0.0 0.0 0.0

0.5 0.5 0.5

0.625 0.25 0.625

0.5 0.25 0.75

0.5 0.5 0.5

0.375 0.375 0.75

0.625 0.25 0.625

0.5 0.0 0.5

ph_intmeth 2

ph_wstep 0.0001 eV

ph_ngqpt 16 16 16

ph_qshift 0 0 0

ph_nqshift 1

eph_ngqpt_fine 8 8 8

##############################################

#### SECTION: gstate

##############################################

optdriver 7

##############################################

#### SECTION: gw

##############################################

nkptgw 1

kptgw 0.0000000000 0.0000000000 0.0000000000

bdgw 1 8

zcut 0.2 eV

nomegasf 301

##############################################

#### STRUCTURE

##############################################

natom 2

ntypat 1

typat 1 1

znucl 6

xred

0.0000000000 0.0000000000 0.0000000000

0.2500000000 0.2500000000 0.2500000000

acell 1.0 1.0 1.0

rprim

0.0000000000 3.3517340250 3.3517340250

3.3517340250 0.0000000000 3.3517340250

3.3517340250 3.3517340250 0.0000000000

#### SECTION: basic

##############################################

ecut 25.0

nband 200

tolvrs 1e-08

ngkpt 8 8 8

kptopt 1

nshiftk 1

shiftk 0.0 0.0 0.0

##############################################

#### SECTION: eph

##############################################

eph_task 4

ddb_ngqpt 4 4 4

tmesh 0 200 5

ph_ndivsm 10

ph_nqpath 12

ph_qpath

0.0 0.0 0.0

0.5 0.0 0.5

0.5 0.25 0.75

0.375 0.375 0.75

0.0 0.0 0.0

0.5 0.5 0.5

0.625 0.25 0.625

0.5 0.25 0.75

0.5 0.5 0.5

0.375 0.375 0.75

0.625 0.25 0.625

0.5 0.0 0.5

ph_intmeth 2

ph_wstep 0.0001 eV

ph_ngqpt 16 16 16

ph_qshift 0 0 0

ph_nqshift 1

eph_ngqpt_fine 8 8 8

##############################################

#### SECTION: gstate

##############################################

optdriver 7

##############################################

#### SECTION: gw

##############################################

nkptgw 1

kptgw 0.0000000000 0.0000000000 0.0000000000

bdgw 1 8

zcut 0.2 eV

nomegasf 301

##############################################

#### STRUCTURE

##############################################

natom 2

ntypat 1

typat 1 1

znucl 6

xred

0.0000000000 0.0000000000 0.0000000000

0.2500000000 0.2500000000 0.2500000000

acell 1.0 1.0 1.0

rprim

0.0000000000 3.3517340250 3.3517340250

3.3517340250 0.0000000000 3.3517340250

3.3517340250 3.3517340250 0.0000000000

`EphTask`

with:

In [40]:

```
flow[-1][-1].get_graphviz()
```

Out[40]:

The graph shows that the `EphTask`

uses the `WFK`

produced by a NSCF calculation
as well as two other files generated by the `PhononWork`

.
If this is not the first time you run DFPT calculations with Abinit, you already know that the `DDB`

is essentially a text file with the independent entries of the dynamical matrix.
This file is usually produced by merging partial `DDB`

files with the `mrgddb`

utility.
For further information about the `DDB`

see this notebook.

The `DVDB`

, on the other hand, is a Fortran binary file containing an independent set
of DFPT potentials (actually the first-order variation of the self-consistent KS potential
due to an atomic perturbation characterized by a q-point and the [`idir, ipert`

] pair).
The `DVDB`

is generated by merging the `POT`

files produced at the end of the DFPT run
with the `mrgdv`

tool.

Now let's turn our attention to the last `Work`

:

In [41]:

```
flow[-1].get_graphviz()
```

Out[41]:

`flow`

consists of three `EphTasks`

, all with the same parents.
This is the part that computes the e-ph self-energy with a fixed q-mesh and different values of `nband`

.
We can easily check this with:

In [42]:

```
print("nband:", [task.input["nband"] for task in flow[-1]])
print("q-mesh for self-energy:", [task.input["eph_ngqpt_fine"] for task in flow[-1]])
```

nband: [100, 150, 200] q-mesh for self-energy: [[8, 8, 8], [8, 8, 8], [8, 8, 8]]

`get_vars_dataframe`

with the list of variable names:

In [43]:

```
flow[-1].get_vars_dataframe("nband", "eph_ngqpt_fine")
```

Out[43]:

nband | eph_ngqpt_fine | class | |
---|---|---|---|

w3_t0 | 100 | [8, 8, 8] | EphTask |

w3_t1 | 150 | [8, 8, 8] | EphTask |

w3_t2 | 200 | [8, 8, 8] | EphTask |

`work`

made of `EphTasks`

with similar organization, the only difference being
that we use the ab-initio q-mesh for the DFPT potentials (no Fourier interpolation here).

In [44]:

```
print("nband:", [task.input["nband"] for task in flow[-2]])
print("q-mesh for self-energy:", [task.input["eph_ngqpt_fine"] for task in flow[-2]])
```

nband: [100, 150, 200] q-mesh for self-energy: [[4, 4, 4], [4, 4, 4], [4, 4, 4]]

In [45]:

```
flow[-2].get_vars_dataframe("nband", "eph_ngqpt_fine")
```

Out[45]:

nband | eph_ngqpt_fine | class | |
---|---|---|---|

w2_t0 | 100 | [4, 4, 4] | EphTask |

w2_t1 | 150 | [4, 4, 4] | EphTask |

w2_t2 | 200 | [4, 4, 4] | EphTask |

At this point, we should have a good understanding of the input files required by the `EphTasks`

.
There is only one point missing! How do we get the `WFK`

, `DDB`

, `DVDB`

files?

Well, in the old days we used to write several input files one for each q-point in the IBZ, link files and merge the intermediate results manually but now we can use python to automate the entire procedure:

In [46]:

```
flow.get_graphviz()
```

Out[46]:

More specifically, the q-mesh associated to the
`DDB`

(`DVDB`

) must be a submesh of the k-mesh for electrons.
Moreover the EPH code can compute self-energy corrections only for the states that are available in the `WFK`

file.

Now we can generate the directories and the input files of the `Flow`

by executing the
`lesson_sigeph.py`

script that will generate the flow_diamond directory with all the input
files required by the calculation.

Then use the `abirun.py`

script to launch the calculation with the scheduler:

```
abirun.py flow_diamond scheduler
```

You may want to run this example in the terminal if you have already installed and configured AbiPy and Abinit on your machine. The calculation requires ~8 minutes on a poor 1.7 GHz Intel Core i5 (2 minutes for all the EPH tasks, the rest for the DFPT section).

If you prefer to skip this part, jump to next section in which we focus on the post-processing of the results. Note that the output files are already available in the github repository so it is possible to play with the AbiPy post-processing tools without having to run the flow. In particular, one can use the command line and the commands:

```
abiopen.py FILE
```

to open the file inside ipython,

```
abiopen.py out_SIGEPH.nc --expose
```

to visualize the self-energy results and finally,

```
abicomp.py sigeph flow_diamond/w3/
```

to compare multiple `SIGEPH.nc`

files with the robot inside ipython.

In [47]:

```
!find flow_diamond/ -name "out_GSR.nc"
```

`w0/t0`

computed the electronic band structure on a high-symmetry k-path.
Let's plot the bands with:

In [48]:

```
with abilab.abiopen("flow_diamond/w0/t1/outdata/out_GSR.nc") as gsr:
ebands_kpath = gsr.ebands
ebands_kpath.plot(with_gaps=True);
```

`NscfTask`

:

In [49]:

```
print(ebands_kpath)
```

In the `NscfTask`

`w0/t2`

we have computed the KS eigenstates including several empty states.
This is the `WFK`

file we are going to use to build the `EPH`

self-energy.

Let's plot the bands in the IBZ and the associated DOS with:

In [50]:

```
with abilab.abiopen("flow_diamond/w0/t2/outdata/out_GSR.nc") as gsr:
ebands_4sigma = gsr.ebands
# 8x8x8 is too coarse to get nice DOS so we increase the gaussian smearing
# (we are mainly interested in the possible presence of ghost states --> sharp peaks)
ebands_4sigma.plot_with_edos(ebands_4sigma.get_edos(width=0.6));
```

The figure shows that our value of `nband`

correspond to ~400 eV above the valence band maximum (VBM).

Note that it is always a good idea to look at the behaviour of the high-energy states when running self-energy calculations, especially the first time you use a new pseudopotential.

**Ghost states**, indeed, can appear somewhere in the high-energy region causing sharp peaks in the DOS.
The presence of these high-energy ghosts is not clearly seen in standard GS calculations
but they can have drammatic consequences for many-body calculations based on sum-over-states approaches.

This is one of the reasons why the pseudopotentials of the pseudo-dojo project have been explicitly tested for the presence of ghost states in the empty region (see here for further details)

Now we turn our attention to the vibrational properties.
AbiPy has already merged all the independent atomic perturbations in `flow_diamond/w1/outdata/out_DDB`

:

In [51]:

```
!find flow_diamond/ -name "out_DDB"
```

This is the input file for `mrgddb`

produced automatically by python:

In [52]:

```
!cat flow_diamond//w1/outdata/mrgddb.stdin
```

The first line gives the name of the output DDB file followed by a title, the number of partial DDB files we want to merge and their path.

A similar input file must be provided to `mrgdv`

, the only difference is that DDB files are not replaced
by `POT*`

files with the SCF part of the DFPT potential (Warning: Unlike the out_DDB files,
the out_DVDB file and the associated mrgdvdb.stdin files are available
in the present tutorial only if you have actually issued `abirun.py flow_diamond scheduler`

above.
Still, the files needed beyond the two next commands are available by default, so that you can continue to examine this tutorial.)

In [53]:

```
!find flow_diamond/ -name "out_DVDB"
```

flow_diamond//w1/outdata/out_DVDB

In [54]:

```
!cat flow_diamond//w1/outdata/mrgdvdb.stdin
```

Let's open the `DDB`

file with:

In [55]:

```
ddb = abilab.abiopen("flow_diamond/w1/outdata/out_DDB")
print(ddb)
```

and then invoke `anaddb`

to compute phonon bands and DOS:

In [56]:

```
phbst, phdos = ddb.anaget_phbst_and_phdos_files()
```

Finally we plot the results with:

In [57]:

```
phbst.phbands.plot_with_phdos(phdos);
```

The maximum phonon frequency is ~0.17 eV that is much smaller that the KS fundamental gap. This means that at T=0 we cannot have scattering processes between two electronic states in which only one phonon is involved. This however does not mean that the electronic properties will not be renormalized by the EPH interaction, even at T=0 (zero-point motion effect).

So far, we managed to generate a `WFK`

with empty states on a 8x8x8 k-mesh,
and a pair of `DDB`

-`DVDB`

files on a 4x4x4 q-mesh.
We can now finally turn our attention to the EPH self-energy.

Let's focus on the output files produced by the first `EphTask`

in `w2/t0`

:

In [58]:

```
!ls flow_diamond/w2/t0/outdata
```

out_EBANDS.agr out_PHBANDS.agr out_PHDOS.nc out_OUT.nc out_PHBST.nc out_SIGEPH.nc

where:

*out_PHBST.nc*: phonon band structure along the q-path*out_PHDOS.nc*: phonon DOS and projections over atoms and directions*out_SIGEPH.nc*: $\Sigma^{eph}_{nk\sigma}(T)$ matrix elements for different temperatures $T$

There are several physical properties stored in the `SIGEPH`

file.
As usual, we can use `abiopen`

to open the file and `print(abifile)`

to get a quick summary
of the most important results.

Let's open the file obtained with the 8x8x8 q-mesh and 200 bands (our best calculation):

In [59]:

```
#sigeph = abilab.abiopen("flow_diamond/w2/t0/outdata/out_SIGEPH.nc")
sigeph = abilab.abiopen("flow_diamond/w3/t2/outdata/out_SIGEPH.nc")
print(sigeph)
```

*graphical* representation of the results.
Let's use matplotlib to plot the **difference** between the QP and the KS direct gap at $\Gamma$ as a function or T:

In [60]:

```
sigeph.plot_qpgaps_t();
```

*increases* with T.

In [61]:

```
# Valence band max
sigeph.plot_qpdata_t(spin=0, kpoint=0, band_list=[3]);
```

In [62]:

```
# Conduction band min
sigeph.plot_qpdata_t(spin=0, kpoint=0, band_list=[4]);
```

Note how the real part of the self-energy at the VBM goes up with T while the CBM goes down with T, thus explaining the behavior of the QP gap vs T observed previously (well one should also include $Z$ but the real part gives the most important contribution...)

We do not discuss the results for the imaginary part because these values requires a very dense sampling of the BZ to give reasonable results.

In [63]:

```
sigma = sigeph.get_sigeph_skb(spin=0, kpoint=[0, 0, 0], band=3)
sigma.plot_tdep(xlims=[-1, 2]);
```

In [64]:

```
sigeph.get_sigeph_skb(spin=0, kpoint=[0, 0, 0], band=4).plot_tdep(xlims=(-1, 2));
```