#!/usr/bin/env python
# coding: utf-8
# Back to the main [Index](../index.ipynb)
#
#
Temperature-dependent QP band structures
#
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.
#
#
# ## Table of Contents
# [[back to top](#top)]
#
# * [A bit of Theory](#A-bit-of-theory-and-some-details-about-the-implementation)
# * [Building a Flow for EPH self-energy calculations](#Building-a-Flow-for-EPH-self-energy-calculations)
# * [Electronic properties](#Electronic-properties)
# * [Vibrational properties](#Vibrational-properties)
# * [E-PH self-energy](#E-PH-self-energy)
# * [Using SigEphRobot for convergence studies](#Using-SigEphRobot-for-convergence-studies)
# * [Convergence study with respect to nbsum and nqibz](#Convergence-study-with-respect-to-nbsum-and-nqibz)
# * [Interpolating QP corrections with star functions](#Interpolating-QP-corrections-with-star-functions)
#
#
#
# ## A bit of theory and some details about the implementation
# [[back to top](#top)]
#
# $\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$).
#
# ### EPH self-energy
#
# 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:
#
# \begin{equation}
# \Sigma^{FM}_{n\kk}(\ww) = \sum_{\qq\nu m} | g_{nm}^{\nu}(\kk, \qq) | ^ 2 \times
# \Biggl [
# \dfrac{n_{\qq\nu} + f_{m\kpq}}{\ww - \ee_{m\kpq} - \ww_{\qq\nu} + i\eta} +
# \dfrac{n_{\qq\nu} + 1 - f_{m\kpq}}{\ww + \ee_{m\kpq} + \ww_{\qq\nu} + i\eta}
# \Biggr ]
# \end{equation}
#
# 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}
#
#
# ### Interpolation of the DFPT potentials
#
# 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](http://dx.doi.org/10.1103/PhysRevB.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:
#
# \begin{equation}
# V^{scf}_{\tilde\qq}(\rr) = \sum_\RR e^{+i{\tilde \qq}\cdot\RR} W(\rr, \RR)
# \end{equation}
#
# 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.
#
#
# ## Technical details
#
# For the previous discussion, it should be clear that the evaluation of the e-ph self-energy requires:
#
# 1. An initial set of KS wavefunctions and eigenvalues including **several empty** states
#
# 2. The knowledge of the dynamical matrix $D(\qq)$ from which one
# can obtain phonon frequencies and displacements everywhere in the BZ via FT interpolation.
#
# 3. 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:
#
# 1. 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.
#
# 2. 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).
#
#
# ## Suggested references
# [[back to top](#top)]
#
# You might find additional material, related to the present section, in the following references:
#
# 1. [Temperature dependence of electronic eigenenergies in the adiabatic harmonic approximation](https://doi.org/10.1103/physrevb.90.214304)
# 2. [Temperature dependence of the electronic structure of semiconductors and insulators](https://aip.scitation.org/doi/abs/10.1063/1.4927081)
# 3. [Electron-phonon interactions from first principles](https://doi.org/10.1103/RevModPhys.89.015003)
# ## Building a Flow for EPH self-energy calculations
# [[back to top](#top)]
#
# 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
get_ipython().run_line_magic('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
#
# 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[37]:
abilab.print_source(build_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](https://nbviewer.jupyter.org/github/abinit/abitutorials/blob/master/abitutorials/base3/lesson_base3.ipynb),
# an example of `Flow` for phonon calculations is given in
# [lesson_dfpt](https://nbviewer.jupyter.org/github/abinit/abitutorials/blob/master/abitutorials/dfpt/lesson_dfpt.ipynb).
# An AbiPy Flow for the computation of $G_0W_0$ corrections is discussed in
# [lesson_g0w0](https://nbviewer.jupyter.org/github/abinit/abitutorials/blob/master/abitutorials/g0w0/lesson_g0w0.ipynb).
#
# 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()
# 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[39]:
print(flow[-1][-1])
flow[-1][-1].input
#
# 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[40]:
flow[-1][-1].get_graphviz()
# 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](../ddb.ipynb).
#
# 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()
# 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 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]])
# or, alternatvely, we select the last work and call `get_vars_dataframe` with the list of variable names:
# In[43]:
flow[-1].get_vars_dataframe("nband", "eph_ngqpt_fine")
# 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[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]])
# In[45]:
flow[-2].get_vars_dataframe("nband", "eph_ngqpt_fine")
# 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()
#
#
#
# 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](https://github.com/abinit/abitutorials)
# 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.
# ## Electronic properties
# [[back to top](#top)]
#
# Let's focus on the electronic properties first.
# In[47]:
get_ipython().system('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[48]:
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[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](http://www.pseudo-dojo.org/)
# have been explicitly tested for the presence of ghost states in the empty region
# (see [here](https://www.sciencedirect.com/science/article/pii/S0010465518300250?via%3Dihub) for further details)
# ## Vibrational properties
# [[back to top](#top)]
#
# 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]:
get_ipython().system('find flow_diamond/ -name "out_DDB"')
# This is the input file for `mrgddb` produced automatically by python:
# In[52]:
get_ipython().system('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]:
get_ipython().system('find flow_diamond/ -name "out_DVDB"')
# In[54]:
get_ipython().system('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.
# ## E-PH self-energy
# [[back to top](#top)]
#
# Let's focus on the output files produced by the first `EphTask` in `w2/t0`:
# In[58]:
get_ipython().system('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[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)
# 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[60]:
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[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.
# 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[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));
# 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[65]:
df = sigeph.get_dataframe_sk(spin=0, kpoint=[0, 0, 0], ignore_imag=True)
df[df["tmesh"] == 0.0]
# In[66]:
# To plot the QP results as a function of the KS energy:
#sigeph.plot_qps_vs_e0();
# In[67]:
# To plot KS bands and QP results:
#sigeph.plot_qpbands_ibzt();
# ## Using SigEphRobot for convergence studies
# [[back to top](#top)]
#
# 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[68]:
robot_bsum = abilab.SigEPhRobot.from_dir_glob("flow_diamond/w2/t*/outdata/")
robot_bsum
# Let's change the labels of the files by replacing file paths with strings constructed
# from `nbsum`:
# In[69]:
robot_bsum.remap_labels(lambda sigeph: "nbsum: %d" % sigeph.nbsum)
# and check that everything is OK by printing the convergence parameters:
# In[70]:
robot_bsum.get_params_dataframe()
# 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[71]:
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[72]:
robot_bsum.plot_qpdata_conv_skb(spin=0, kpoint=(0, 0, 0), band=3, sortby="nbsum");
# In[73]:
#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[74]:
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?
# ## Convergence study with respect to nbsum and nqibz
# [[back to top](#top)]
#
# 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[75]:
robot_qb = abilab.SigEPhRobot.from_dir("flow_diamond")
robot_qb
# and replace the labels with:
# In[76]:
robot_qb.remap_labels(lambda sigeph: "nbsum: %d, nqibz: %d" % (sigeph.nbsum, sigeph.nqibz))
# 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[77]:
robot_qb.get_params_dataframe()
# 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[78]:
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](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.96.155207)
# to automate $GW$ calculation.
# In[79]:
#robot_qb.plot_qpfield_vs_e0("qpeme0", itemp=0, sortby="nbsum", hue="nqibz");
# In[80]:
#robot_qb.plot_qpdata_conv_skb(spin=0, kpoint=(0, 0, 0), band=3, sortby="nbsum", hue="nqibz");
# In[81]:
#robot_qb.plot_selfenergy_conv(spin=0, kpoint=(0.5, 0, 0), band=4, itemp=0, sortby="nbsum", hue="nqibz");
# In[82]:
#robot_qb.plot_qpgaps_t(); #sortby="nbsum", plot_qpmks=True, hue="nqibz");
# In[83]:
#robot_qb.plot_qpgaps_t(sortby="nbsum", plot_qpmks=True, hue="nqibz");
# In[84]:
#robot_qb.plot_qpgaps_convergence(itemp=0, sortby="nbsum", hue="nqibz");
# In[85]:
#robot_qb.plot_qpfield_vs_e0("ze0", itemp=0, sortby="nbsum", hue="nqibz");
# ## Interpolating QP corrections with star-functions
# [[back to top](#top)]
# 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](http://dx.doi.org/10.1103/physrevb.38.2721)
# 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[86]:
sigeph_ibz = abilab.abiopen("ibz_SIGEPH.nc")
# These are the QP energies and the KS bands in the irreducible wedge:
# In[87]:
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[88]:
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[89]:
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[90]:
tdep.plot();
# Let's zoom in on the region around the gap:
# In[91]:
tdep.plot(ylims=(-1, 6));
# Back to the main [Index](../index.ipynb)