Back to the main Index

Carbon in diamond structure

In this lesson, we discuss how to compute the electron-phonon (EPH) self-energy and the renormalization of the electronic states due to EPH interaction. We start with a very brief discussion of the basic equations and some technical details about the Abinit implementation. 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.

Please note that the code for the E-PH self-energy is still under active development and lacks a systematic validation.

This new implementation differs from the one described in the official [Abinit tutorial](https://docs.abinit.org/tutorial/tdepes/) in several points. First of all, the new code computes the electron-phonon matrix elements directly using the input wavefunctions and the DFPT potentials while the official version employs a python script to post-process the matrix elements produced by the DFPT code and uses the first order change of the KS wavefunctions to reduce the number of bands.

The new approach, on the other hand, is based on the sum over empty states approach and is implemented in Fortran for better efficiency and improved parallel scalability. The new algorithm can take advantage of symmetries to reduce the number of k-points/q-points as well as the number of atomic perturbations that must be computed explicitly.

- 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} \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 EPH 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 representation are given by:

where $n$ ($f$) are the temperature-dependent occupation factors for phonons (electrons)
and $\eta$ is a small positive number (*zcut*).
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_{nb}^{\nu}(\kk, \qq=\Gamma)$ elements (see references for a more detailed discussion).

The QP energies are usually obtained with:

\begin{equation} \ee^{QP}_{b\kk} \approx \ee^{KS}_{b\kk} + Z_{b\kk}\,\Sigma^{ep-h}(\ee^{KS}_{b\kk}) \end{equation}where the renormalization factor $Z$ is given by:

\begin{equation} Z_{b\kk} = \Bigl ( 1 - \Re\dfrac{\partial\Sigma^{e-ph}}{\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, also because in this example we will focus on 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 electron-phonon matrix elements and build the self-energy matrix elements. A schematic representation of a typical E-PH workflow 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 [1]:

```
# 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 [2]:

```
from lesson_sigeph import build_flow
```

Please read the code carefully, in particular the comments.
Don't worry if the meaning of some input variables is not immediately clear
as we will try to clarify the most technical parts in the rest of this notebook.

In [3]:

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

Out[3]:

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 EPH 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 [4]:

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

We now print the input of the last `EphTask`

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

, `varrf`

and `vareph`

sections.

In [5]:

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

Out[5]:

The input does not contain any `irdwfk` or `getwfk` variable. AbiPy will add the required `ird*`
variables at runtime and create symbolic links to connect this task to its parents.

As usual, it is much easier to understand what is happening if we plot a graph
with the individual tasks and their connections.
Since there are several nodes in our graph, we first focus on the EPH part
and then we "zoom-out" to unveil the entire workflow.
Let's have a look at the parents of the last `EphTask`

with:

In [6]:

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

Out[6]:

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 [7]:

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

Out[7]:

As you can see, this section of the `flow`

consists of three `EphTasks`

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

.
We can easily check this with:

In [8]:

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

or, alternatvely, we select the last work and call `get_vars_dataframe`

with the list of variable names:

In [9]:

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

Out[9]:

There is another `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 [10]:

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

In [11]:

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

Out[11]:

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 [12]:

```
flow.get_graphviz()
```

Out[12]:

Remember that the $q$-point mesh cannot be chosen arbitrarily
since all $q$ wavevectors should connect two $k$ points of the grid used for the electrons.
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 [13]:

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

The task `w0/t0`

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

In [14]:

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

It's not surprising to see that diamond is a indirect gap semiconductor.
The KS band structure gives a direct gap at $\Gamma$ of ~5.6 eV while the fundamental gap
is ~4.2 eV (VBM at $\Gamma$, CBM at ~[+0.357, +0.000, +0.357] in reduced coordinates).
These values have been computer by ApiPy on the basis of the k-path employed to compute the band structure
in the `NscfTask`

:

In [15]:

```
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 [16]:

```
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 [17]:

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

This is the input file for `mrgddb`

produced automatically by python:

In [18]:

```
!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 [19]:

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

In [20]:

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

Let's open the `DDB`

file with:

In [21]:

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

and then invoke `anaddb`

to compute phonon bands and DOS:

In [22]:

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

Finally we plot the results with:

In [23]:

```
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 [24]:

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

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 [25]:

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

The output reveals that the direct QP gaps are smaller
than the corresponding KS values, even at $T=0$.
Still it is difficult to understand what is happening without a *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 [26]:

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

In diamond, the EPH interaction at T=0 decreases the direct gap at $\Gamma$ with respect
to the KS value. Moreover the QP gap decreases with T.
This behaviour, however, cannot be generalized in the sense that there are systems
in which the direct gap *increases* with T.

So far we have been focusing on the QP direct gaps but we can also look at the QP results for the CBM and the VBM:

In [27]:

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

In [28]:

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

We can also plot the real and the imaginary part of $\Sigma_{nk}^{eph}(\omega)$ and the spectral function $A_{nk}(\omega)$ obtained for the different temperatures with:

In [29]:

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

In [30]:

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

Converging these quantities is not an easy task since they are very sensitive to the q-point sampling. Still, we can see that the spectral function has a dominant peak that is shifted with respect to the initial KS energy. The spectral function of the CBM has more structures than the VBM...

In some cases, we need to access the raw data.
For example we may want to build a table with the QP results for all bands at fixed k-point, spin.
The code below creates a pandas `DataFrame`

and selects only the entries at 0 K.

In [31]:

```
df = sigeph.get_dataframe_sk(spin=0, kpoint=[0, 0, 0], ignore_imag=True)
df[df["tmesh"] == 0.0]
```

Out[31]:

In [32]:

```
# To plot the QP results as a function of the KS energy:
#sigeph.plot_qps_vs_e0();
```

In [33]:

```
# To plot KS bands and QP results:
#sigeph.plot_qpbands_ibzt();
```

Our initial goal was to analyze the convergence of the QP corrections as a function of `nbsum`

.
This means that we should now collect multiple `SIGEPH.nc`

files, extract the data,
order the results and visualize them.
Fortunately, we can use the `SigEphRobot`

(or the `abicomp.py`

script)
to automate the most boring and repetitive operations.

The code below, for instance, tells our robot to open all the `SIGEPH`

files
located within `flow_diamond/w2`

(fixed q-mesh, different `nband`

).
The syntax used to specify multiple directories should be familiar to Linux users:

In [34]:

```
robot_bsum = abilab.SigEPhRobot.from_dir_glob("flow_diamond/w2/t*/outdata/")
robot_bsum
```

Out[34]:

Let's change the labels of the files by replacing file paths with strings constructed
from `nbsum`

:

In [35]:

```
robot_bsum.remap_labels(lambda sigeph: "nbsum: %d" % sigeph.nbsum)
```

Out[35]:

and check that everything is OK by printing the convergence parameters:

In [36]:

```
robot_bsum.get_params_dataframe()
```

Out[36]:

As expected, the only parameter that is changing in this set of netcdf files is `nbsum`

.
Now all the files are in memory and we can start our convergence study.

There are several quantities that depend on T thus complicating a bit the analysis.
So we use the `itemp` index to select the first temperature.

To plot the convergence of QP gaps vs the `nbsum`

parameter, use:

In [37]:

```
robot_bsum.plot_qpgaps_convergence(sortby="nbsum", itemp=0);
```

The convergence of the QP gaps versus `nbsum`

is not variational
and very slow: the same run with 20 bands (not shown here) would give completely different results.
Note that QP gaps (energy differences) are usually easier to converge than absolute QP energies.
A similar behaviour is observed in $GW$ calculations as well.

We can also analyze the convergence of the individual QP terms for given $(\sigma, {\bf{k}}, n)$ with:

In [38]:

```
robot_bsum.plot_qpdata_conv_skb(spin=0, kpoint=(0, 0, 0), band=3, sortby="nbsum");
```

In [39]:

```
#robot_bsum.plot_qpdata_conv_skb(spin=0, kpoint=(0, 0, 0), band=4, sortby="nbsum");
```

and perform a similar analysis for $\Sigma^{eph}_{n{\bf k}}(\omega)$ and $A_{n{\bf k}}(\omega)$:

In [40]:

```
robot_bsum.plot_selfenergy_conv(spin=0, kpoint=(0 , 0, 0), band=4, itemp=0);
```

Note how the largest variations are observed in the real part of the self-energy while the imaginary part is rather insensitive to the number of empty states. Why?

In the previous section, we have analyzed the convergence of the QP results as a function of `nbsum`

for a fixed q-point mesh.
The main conclusion was that the convergence with respect to the number of bands is slow
and that even 200 bands are not enough to converge the QP gap.

Here we use *all* the `SIGEPH.nc`

files produced by the `flow`

to study
how the final results depend on `nbsum`

and `nqibz`

.

The main goal is to understand whether these two parameters are correlated or not.
In other words we want to understand whether we really need to perform a convergence study in which
both parameters are increased or whether we can assume some sort of *weak correlation* that allows
us to *decouple* the two convergence studies.

Let's start by loading all the `SIGPEH.nc`

files located inside `flow_diamond`

with:

In [41]:

```
robot_qb = abilab.SigEPhRobot.from_dir("flow_diamond")
robot_qb
```

Out[41]:

and replace the labels with:

In [42]:

```
robot_qb.remap_labels(lambda sigeph: "nbsum: %d, nqibz: %d" % (sigeph.nbsum, sigeph.nqibz))
```

Out[42]:

All the calculations in `w2`

are done with a 4x4x4 q-mesh, while in `w3`

we have EPH calculations
done with a 8x8x8 q-mesh.
Let's print a table with the parameters of the run, just to be sure...

In [43]:

```
robot_qb.get_params_dataframe()
```

Out[43]:

We can start to analyze the data but there is a small problem because, strictly speaking,
we are dealing with a function of two variables. (`nbsum`

and `nqibz`

).
Let's analyze the convergence of the QP direct gap as a function of `nbsum`

and use the `hue`

argument
to group results obtained with different q-meshes:

In [44]:

```
robot_qb.plot_qpgaps_convergence(sortby="nbsum", hue="nqibz", itemp=0);
```

The figure reveals that the largest variation is due to the q-point sampling,
nonetheless the two curves as function of `nbsum`

have very similar trend
(they are just converging to different values)
This suggests that one could simplify considerably the convergence study by decoupling the two parameters.
In particular, one could perform an initial convergence study with respect to `nbsum`

with a relatively coarse q-mesh
to find a reasonable value then *fix* `nbsum`

and move to the convergence study with respect to the q-point sampling
and the zcut broadening.

A similar approach has been recently used by van Setten to automate $GW$ calculation.

In [45]:

```
#robot_qb.plot_qpfield_vs_e0("qpeme0", itemp=0, sortby="nbsum", hue="nqibz");
```

In [46]:

```
#robot_qb.plot_qpdata_conv_skb(spin=0, kpoint=(0, 0, 0), band=3, sortby="nbsum", hue="nqibz");
```

In [47]:

```
#robot_qb.plot_selfenergy_conv(spin=0, kpoint=(0.5, 0, 0), band=4, itemp=0, sortby="nbsum", hue="nqibz");
```

In [48]:

```
#robot_qb.plot_qpgaps_t(); #sortby="nbsum", plot_qpmks=True, hue="nqibz");
```

In [49]:

```
#robot_qb.plot_qpgaps_t(sortby="nbsum", plot_qpmks=True, hue="nqibz");
```

In [50]:

```
#robot_qb.plot_qpgaps_convergence(itemp=0, sortby="nbsum", hue="nqibz");
```

In [51]:

```
#robot_qb.plot_qpfield_vs_e0("ze0", itemp=0, sortby="nbsum", hue="nqibz");
```

In the previous paragraphs, we discussed how to compute the QP corrections as function of $T$
for a set of k-points in the IBZ.
In many cases, however, we would like to visualize the effect of the self-energy
on the *electronic dispersion* using a high-symmetry k-path.
There are several routes to address this problem.

The most accurate method would be performing multiple self-energy calculations
with different shifted k-meshes so as to sample the high-symmetry k-path
and then merging the results.
Alternatively, one could try to *interpolate* the QP energies obtained in the IBZ.

It is worth noting that the QP energies must fulfill the symmetry properties of the point group of the crystal:

\begin{equation} \epsilon({\bf k + G}) = \epsilon({\bf k}) \end{equation}and

\begin{equation} \epsilon(\mathcal O {\bf k }) = \epsilon({\bf k}) \end{equation}where $\bf{G}$ is a reciprocal lattice vector and $\mathcal{O}$ is an operation of the point group.

Therefore it is possible to employ the star-function interpolation by Shankland, Koelling and Wood in the improved version proposed by Pickett to fit the ab-initio results. This interpolation technique, by construction, passes through the initial points and satisfies the basic symmetry property of the band energies.

It should be stressed, however, that this Fourier-based method can have problems in the presence of band crossings that may cause unphysical oscillations between the ab-initio points. To reduce this spurious effect, we prefer to interpolate the QP corrections instead of the QP energies. The corrections, indeed, are usually smoother in k-space and the resulting fit is more stable.

To interpolate the QP corrections we have to provide a SIGEPH file with the QP energies computed
in the entire irreducible zone.
A netcdf file obtained with `nband == 100`

is already provided in the github repository.
Note that the QP corrections in this case have been obtained by setting $Z = 1$

In [52]:

```
sigeph_ibz = abilab.abiopen("ibz_SIGEPH.nc")
```

These are the QP energies and the KS bands in the irreducible wedge:

In [53]:

```
sigeph_ibz.plot_qpbands_ibzt();
```

Now we extract the KS bands along the k-path from the GSR file produced by the second NscfTask

In [54]:

```
with abilab.abiopen("flow_diamond/w0/t1/outdata/out_GSR.nc") as gsr:
ks_ebands_kpath = gsr.ebands
```

Finally, we **interpolate the QP corrections** for the first two temperatures by calling:

In [55]:

```
tdep = sigeph_ibz.interpolate(ks_ebands_kpath=ks_ebands_kpath, itemp_list=[0, 1])
```

The small value of the MAE indicates that the fit passes through the ab-initio points (as expected).
This, however, does not exclude the possibility of unphysical oscillations between the input data points.

Finally, one can plot the interpolated QP bands and the ab-initio KS dispersion with:

In [56]:

```
tdep.plot();
```

Let's zoom in on the region around the gap:

In [57]:

```
tdep.plot(ylims=(-1, 6));
```

Back to the main Index

In [ ]:

```
```