Back to the main Index

Aluminium

In this lesson, we discuss how to compute the phonon linewidths and the Eliashberg function in metals. 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 these calculations.

In the second part, we use the AbiPy objects to analyze the results. We start from the simplest case of a single calculation whose most important results are stored in the A2F.nc netcdf file. Then we show how to use the A2fRobot to handle multiple netcdf files and perform convergence studies.

- A bit of theory
- Technical details
- Building the Flow
- Electronic properties
- Vibrational properties
- Linewidths and Eliashberg function
- Convergence studies with A2FRobot

[back to 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}$ $\newcommand{\dd}{{\,\text{d}}}$

The phonon linewidths are proportional to the electron phonon coupling, and depend on the phonon wavevector $q$ and the branch index $\nu$:

\begin{equation} \gamma_{\qq\nu} = 2\pi \omega_{\qq\nu} \sum_{mn\kk} |g_{mn\nu}(\kk, \qq)|^2 \delta(\ee_{\kpq m}) \delta(\ee_{\kk n}) \end{equation}Throughout these notes we shall use Hartree atomic units ($e = \hbar = m_e = 1$), the Fermi level is set to zero.

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^\KS |\psi_{n\kk}\> \end{equation}For further details about $\Delta_{\qq\nu} V^\KS$ and their Fourier interpolation see the previous lesson on the EPH self-energy and Phys. Rev. B 78, 045124.

The Eliashberg function, $\alpha^2F(\omega)$, is similar to the density of states of the phonons, $F(\omega)$, but is weighted according to the coupling of the phonons to the electrons:

\begin{equation} \alpha^2F(\omega) = -\dfrac{1}{N_F} \sum_{\qq\nu} \dfrac{\gamma_{\qq\nu}}{\omega_{\qq\nu}} \delta(\ww - \ww_{\qq \nu}) \end{equation}The first inverse moment of $\alpha^2F(\omega)$ gives the total coupling strength, or mass renormalization factor, $\lambda$.

\begin{equation} \lambda = \int \dfrac{\alpha^2F(\omega)}{\omega}\dd\omega = \sum_{\qq\nu} \lambda_{\qq\nu} \end{equation}\begin{equation} \lambda_{\qq\nu} = \dfrac{\gamma_{\qq\nu}}{\pi N_F \ww_{\qq\nu}^2} \end{equation}From $\lambda$, using the McMillan formula
as modified by Allen and Dynes in Phys. Rev. B 12 905,
one can estimate the superconducting critical temperature $T_c$ in the **isotropic** case.

where the logarithmic moment, $\omega_{log}$, is defined by:

\begin{equation} \omega_{log} = \exp \Biggl [ \dfrac{2}{\lambda} \int \dfrac{\alpha^2F(\omega)}{\omega}\log(\omega)\dd\omega \Biggr ] \end{equation}The formula contains the $\mu^*$ parameter which approximates the effect of the Coulomb interaction.
It is common to treat $\mu^*$ as an *adjustable parameter* to reproduce (explain) the experimental $T_c$
from the ab-initio computed $\lambda$.

It is worth noting that, if we assume constant e-ph matrix elements in the BZ, the phonon linewidths
become proportional to the so-called *nesting* function defined by:

Roughly speaking, there are two factors entering into play in the equation for the phonon linewidths:
the behaviour in (k, q) space of the matrix elements and the *geometrical* contribution related
to the shape of the Fermi surface (described by the nesting term).

The input variables *optdriver = 7* and *eph_task = 1* activate the computation of:

for all q-points in the IBZ as determined by the `ddb_ngqpt`

input variable.
In the above equation, $p$ is a short-hand notation for atom index and direction,
$\vec d$ is the phonon displacement and $\tilde\gamma_{p p'}(\qq)$ is given by:

The $\tilde\gamma_{p p'}(\qq)$ matrix has the same symmetries as the dynamical matrix
and the elements in the *full* BZ can be obtained by *rotating* the initial set of q-points in the IBZ.

Once $\tilde\gamma(\qq)$ is known in the *full* BZ, one can Fourier transform to real space with:

and use Fourier interpolation to obtain $\gamma_{\qq\nu}$ everywhere in the BZ at low cost (a similar approach is used for the dynamical matrix and phonon frequencies).

Thanks to the relatively inexpensive Fourier interpolation, one can obtain the phonon linewidths along
an arbitrary q-path and evaluate the Eliashberg function on a q-mesh (specified by *ph_ngqpt*)
that can be made **much denser** than the initial ab-initio DDB sampling (specified by *ddb_ngqpt*) and
even denser than the `eph_ngqpt_fine`

q-mesh used to interpolate the DFPT potentials.

Note that the calculation of phonon linewidths in metals is made difficult
by the slow convergence of the double-delta integral over the Fermi surface.
It is therefore convenient to use a coarse k-mesh to calculate phonons with DFPT on a suitable q-grid
and then use a denser k-mesh to perform the integration over the Fermi surface.
The resolution in q-space can be improved by interpolating the DFPT potentials via the
`eph_ngqpt_fine`

input variable as discussed in the previous EPH lesson.

The general theory of electron-phonon coupling and Eliashberg superconductivity is reviewed by P.B. Allen and B. Mitrovic in Theory of Superconducting Tc. The first implementations similar to that in Abinit are those in Savrasov and Liu.

From the previous discussion, it should be clear that a typical calculation of phonon linewidths requires:

- An initial set of KS wavefunctions and eigenvalues.
- The knowledge of the dynamical matrix $D(\qq)$ from which one can obtain frequencies and displacements everywhere in the BZ via Fourier interpolation.
- A set of DFPT potentials given in the IBZ

Thanks to these three ingredients, the code can compute the EPH matrix elements and perform the integration over the Fermi surface. A schematic representation of a typical workflow with Abinit is given in the figure below:

The `mrgddb`

and `mrgdv`

are Fortran executables whose main goal is to produce the final files
required in the EPH code.

Note, in particular, how the computation of the WFK file and the DFPT part of the graph are now decoupled. This is the approach we are going to implement with AbiPy to converge the phonon linewidths in Al:

Compute DFPT phonons on a 4x4x4 q-mesh with a coarse 8x8x8 k-sampling

Generate 3 WFK files on a much denser k-mesh (x16, x24, x32)

Run the EPH code with

- one of the WFK files generated in point 2.
- interpolated DFPT potentials (from the initial 4x4x4 to a 8x8x8 q-mesh)

Compute the Eliashberg function on the

`ph_ngqpt`

mesh via Fourier interpolation.- Analyze the convergence of the results wrt
`nkpt`

.

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_eph_al`

module to build our AbiPy flow.

In [2]:

```
from lesson_eph_al 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 [5]:

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

Out[5]:

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 *phonon linewidths* with different k-point samplings.

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

for phonon calculations is given in
lesson_dfpt.

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 then print the input of the last `EphTask`

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

and `eph`

sections.

In a nutshell: we read a WKF on a 32x32x32 k-mesh and a DDB on a 4x4x4 q-mesh,
activate the computation of phonon linewidths with `optdriver`

and `eph_task`

,
interpolate the potentials onto a 8x8x8 q-mesh with `eph_ngqpt_fine`

and set other variables
for the computation of the phonon linewidths around the Fermi surface.

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 mainly focus on the EPH part.

Let's have a look at the parents of the last `EphTask`

with:

In [6]:

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

Out[6]:

We do not repeat here the detailed description of the Flow because it is very similar in spirit to what has been already done in the previous EPH lesson. Hopefully, you are already familiar with the AbiPy philosophy.

In [7]:

```
print("ngkpt in EPH tasks:", [task.input["ngkpt"] for task in flow[2]])
```

In [8]:

```
flow[2].get_vars_dataframe("ngkpt")
```

Out[8]:

For completeness, we show the entire flow with the different connections:

In [9]:

```
flow.get_graphviz()
```

Out[9]:

You may wonder why we have so many tasks especially if you are used to the multi-dataset philosophy of Abinit.
Indeed if you are a `ndtset` master, it should not be so difficult for you
to write a *single input file* that performs the same kind of convergence study.
Still you have to consider that the multi-dataset approach does not scale well:
this kind of calculation has several independent steps
that could be executed in parallel whereas abinit datasets are executed sequentially.
This makes a huge difference when running on clusters with hundreds or thousands of CPUs because
the *time to solution* can be considerably reduced with this kind of parallelism.

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

with:

```
flow.build_and_pickle_dump()
```

and then use the `abirun.py`

script to launch the entire calculation:

```
abirun.py flow_eph_al scheduler
```

You may want to run this example in the terminal if you've already installed and configured AbiPy and Abinit on your machine. The calculation requires ~12 minutes on a poor 1.7 GHz Intel Core i5 (50% of the time is spent in the last task to compute the phonon linewidths with the 32x32x32 k-mesh)

If you prefer to skip this part, jump to next section in which we focus on the post-processing of the results. Note that most of the output files are already available in the github repository so it is possible to try the AbiPy post-processing tools without having to run the flow. Some DVDB, DDB, PHDOS, and PBSTS files are missing, but their absence will not prevent running the present tutorial to the end. In particular, one can use the command line and the commands:

```
abiopen.py FILE
```

to open the file inside ipython,

```
abiopen.py ut_A2F.nc --expose
```

to visualize the EPH results and finally,

```
abicomp.py a2f flow_ep_al/
```

to compare multiple `A2F.nc`

files with the robot and ipython.

In [10]:

```
!find flow_eph_al/ -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 [11]:

```
with abilab.abiopen("flow_eph_al/w0/t1/outdata/out_GSR.nc") as gsr:
ebands_kpath = gsr.ebands
# NB: Fermi level set to zero
ebands_kpath.plot();
```

In [12]:

```
#print(ebands_kpath)
```

To better understand whatâ€™s happening at $\ee_F$, we can plot the bandwidths with:

In [13]:

```
ebands_kpath.boxplot();
```

There is only one band **index** crossing the Fermi level if we exclude a tiny portion with band index 2 (see the region close to the W point).
As phonon lifetimes are very sensitive to the Fermi surface sampling,
it is a good idea to analyze the convergence of the electronic DOS,
in particular the behavior in the region around $\ee_F$.

Let's use the `GsrRobot`

to load all the `GSR`

files produced by the `w0`

work:

In [14]:

```
gsr_robot = abilab.GsrRobot.from_dir_glob("flow_eph_al/w0/t*/")
gsr_robot
```

Out[14]:

In total, we have 5 `GSR`

files but `w0/t1`

computed the energies along the k-path and the
DOS **requires** a homogeneous sampling.
Let's remove the file for which this is not the case from the robot with:

In [15]:

```
gsr_robot.pop_label("flow_eph_al/w0/t1/outdata/out_GSR.nc")
gsr_robot
```

Out[15]:

and use a small lambda function to change the labels associated to the files so that we have the number of k-points in the IBZ:

In [16]:

```
gsr_robot.remap_labels(lambda gsr: "nkpt: %d" % gsr.nkpt)
```

Out[16]:

Now we can finally compare the electronic DOS obtained with the different k-meshes:

In [17]:

```
gsr_robot.combiplot_edos(xlims=[-15, +5]);
```

Clearly, the 8x8x8 k-mesh used to compute the density is not enough to converge the DOS at $\ee_F$. Remember, however, that we have decided to use a minimalistic sampling in the GS/DFPT run to keep the calculation manageable. In real life, one should use a much denser k-sampling for the GS/DFPT and this is particularly true if we are trying to relax the structure. Let's forget about this technical point and focus on the DOS obtained with the other two k-meshes.

As you can see, even if 145 k-points in the IBZ are not enough, the DOS is becoming *smoother*
and starts to resemble the one of the free-electron gas (as a matter of fact the band dispersion of Al is not
that different from the ideal free-electron model provided that BZ folding is taken into account).

Visual inspection suggests that the k-sampling becomes *acceptable* at and beyond 24x24x24 (413 nkpt).

The convergence of the DOS at the Fermi level does not necessarly imply convergence in the final EPH results.
This is something that should be checked explicitly by looking at the behaviour of the final observables
as a function of the input parameters.

In the introduction, we mentioned that there are two factors governing the strengh of the E-PH coupling in metals:
the behaviour in (k, q)-space of the matrix elements and the *geometrical* contribution due to the Fermi surface.

To understand the *geometrical* contribution, it is useful to visualize the Fermi surface.
Unfortunately, graphical applications usually require KS eigenvalues on a
$\Gamma-$centered k-mesh in the *full* BZ whereas ab-initio codes usually work with KS states
in the *irreducible* wedge.

Fortunately, we can use AbiPy to reconstruct the KS eigenvalues in the full BZ:

In [18]:

```
# Select the last GSR file in thr robot.
eb3d = gsr_robot.abifiles[-1].ebands.get_ebands3d()
```

and then use matplotlib to visualize the Fermi energy isosurfaces:

In [19]:

```
eb3d.plot_isosurfaces();
```

Note that, at present, the matplotlib version is only able to display isosurfaces in the unit cell of the reciprocal lattice. To visualize isosurfaces in the first BZ, one can export the data into BXSF format and then call xcrysden with:

In [20]:

```
#eb3d.xcrysden_view()
```

This is the image you should get with xcrysden:

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

:

In [21]:

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

In [22]:

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

In the same directory, we have the `DVDB`

file containing the independent DFPT potentials

In [23]:

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

In [24]:

```
#!cat flow_eph_al//w1/outdata/mrgdvdb.stdin
```

Let's open the `DDB`

file computed on the 4x4x4 q-mesh with:

In [25]:

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

and then invoke `anaddb`

to compute phonon bands and DOS:

In [26]:

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

Finally we plot the results with:

In [27]:

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

The vibrational spectrum seems OK but remember that we are *enforcing* the acoustic sum rule
with `asr`

.
Since our input parameters are underconverged, its a good idea to compare the spectrum with/without `asr`

:

In [28]:

```
ph_plotter = ddb.anacompare_asr()
```

In [29]:

```
ph_plotter.combiplot(units="cm-1");
```

Not so good! The breaking of the ASR is quite large.
This is mainly due to the use of a too small value for `ecut`.
In real life, we should increase `ecut` and rerun the DFPT part but since this is a tutorial
aiming at showing how to perform EPH calculations, we ignore this convergence issue
keeping in mind that we should redo all our calculations with larger ecut before submitting the final
version of the paper!

We can now finally turn our attention to the phonon linewidths and the Eliashberg function.

We have generated a pair of `DDB`

-`DVDB`

files on a 4x4x4 q-mesh
and **three** `WFK`

files with a much denser k-sampling (x16, x24, x32).
In total we have three EPH calculations done with different k-meshes to analyze.

Let's focus on the output files produced with the 16x16x16 k-mesh by the first `EphTask`

in `w2/t0`

:

In [30]:

```
!ls flow_eph_al/w2/t0/outdata
```

The most important results are stored in:

*out_A2F.nc*: main results in netcdf format*out_PHDOS.nc*: phonon DOS and projections over atoms and directions*out_PBSTS.nc*: phonon band structure along the q-path

There is also a bunch of text files with the same results in text format if you are a gnuplot/xmgrace aficionado...

Let's get a quick summary of the most important results with:

In [31]:

```
a2fnc = abilab.abiopen("flow_eph_al/w2/t0/outdata/out_A2F.nc")
print(a2fnc)
```

In the `E-PH calculation`

we have the value of $\lambda$ and $\omega_{log}$ computed
with a 16x16x16 k-mesh for electrons and two q-meshes.
...

The value of $\lambda$ is smaller (almost a factor 2) with respect to other values reported in the literature, likely due to the coarse 12x12x12 k-sampling. We will investigate this problem in the next section. For the time being, we prefer to focus on the visualization of the results with AbiPy.

Let's use matplotlib to plot the Eliashberg function obtained with the two q-meshes:

In [32]:

```
a2fnc.plot_a2f_interpol();
```

This Eliashberg function obtained on the [24, 24, 24] q-mesh looks nicer, in particular
we see the appearance of Van Hove singularities.
As expected, the integral $\lambda(\omega)$ is less sensitive to the interpolation in q-space.
We conclude that the fact that our $\lambda$ is too small when compared with other ab-initio calculations ($\lambda \approx 0.4$)
is not related to the q-sampling but to the *quality* of our phonon linewidths that in turn is related
to the description of the FS.

We can also visualize the phonon linewidths with a scatter plot:

In [33]:

```
a2fnc.plot_with_a2f(units="meV", what="gamma");
#a2fnc.plot(units="meV", what="gamma");
```

To plot the band energies used in the calculation:

In [34]:

```
#a2fnc.plot_ebands();
```

The notation $\alpha^2 F(\omega)$ was introduced because the Eliashberg function is usually proportional to the phonon DOS $F(\omega)$. There are, however, exceptions so it is useful to plot $\alpha^2F(\omega)$, $F(\omega)$ and their ratio $\alpha^2$.

It is just a matter of passing the phonon DOS computed from the `DDB`

file in the previous section
to the `plot_a2`

method:

In [35]:

```
a2fnc.a2f_qcoarse.plot_a2(phdos=phdos);
```

To conclude this section: our results for $\lambda$ and $\alpha^2F(\omega)$ look reasonable but we are still far from the results reported in previous works. Perhaps we are not completely converged and we should analyze in more detail what happens if we increase the k-point sampling. Fortunately we have already computed these results in our Flow so it is just a matter of using AbiPy to compare multiple calculations.

In this section, we use the `A2fRobot`

to analyze the convergence behaviour of our results
with respect to the k-point sampling in the double-delta integral.

Let's ask our robot to open all the `A2F`

files located within the `flow_eph_al/`

directory.

In [36]:

```
robot = abilab.A2fRobot.from_dir("flow_eph_al/")
robot
```

Out[36]:

We know that all these calculations have been done with different values of `nkpt`

.
Let's change the labels of the files by replacing file paths with more informative strings:

In [37]:

```
robot.remap_labels(lambda ncfile: "nkpt: %d" % ncfile.nkpt)
```

Out[37]:

and print a pandas `DataFrame`

with the parameters of the calculation can help to understand the data:

In [38]:

```
robot.get_params_dataframe()
```

Out[38]:

As usual, it is much easier to analyze the convergence of scalar quantities. Since we are mainly interested in $T_c$-related properties, it makes sense to print a table with the value of $\lambda$ and $\omega_{log}$ extracted from the different calculations:

In [39]:

```
robot.get_dataframe(with_params=False)
```

Out[39]:

This table gives the values integrated on the `eph_ngqpt`

mesh as well as the values obtained
by Fourier interpolating the results on the `ph_ngqpt`

mesh.
Note the **big jump** in $\lambda$ when we go from the 18x18x18 k-mesh to the 36x36x36 k-mesh (~0.2 --> ~0.4).
This clearly shows that our initial estimate for $\lambda$ obtained with a 18x18x18 k-mesh was really bad!
(Did I tell you that these calculations are very sensitive to the k-point sampling?)

If you prefer figures instead of tables with numbers, just use:

In [40]:

```
robot.plot_a2fdata_convergence(sortby="nkpt");
```

When converging with respect to the number of k-points, it is common to plot the physical results as function of $\frac{1}{N_{kpt}}$. Let's define a small function that tells our robot how to sort the results:

In [41]:

```
def inv_nkpt(a2f_file):
"""$\dfrac{1}{N_k}$"""
return 1/ a2f_file.nkpt
robot.plot_a2fdata_convergence(sortby=inv_nkpt); #, hue="eph_fsmear");
```

At this point, our estimate for $\lambda$ should be somewhere in [0.39, 0.42]
that compares much better with the value of 0.44 reported by Savrasov.
**Most importantly**, our results started to converge (although slowly).
Now we know that a serious calculation of the phonon linewidths of Al
would require something around 32x32x32 k-points (this is indeed the mesh used by Savrasov in their paper).

So far, we have been focusing on $\lambda$ but what about the convergence of $\alpha^2F(\omega)$?

In [42]:

```
robot.plot_a2f_convergence();
```

In [43]:

```
print(robot.abifiles[-1])
```

Other post-processing tools:

In [44]:

```
#robot.gridplot_a2f(tight_layout=True);
```

In [45]:

```
#robot.plot_a2f_convergence(sortby="nkpt", hue="eph_fsmear", tight_layout=True);
```

In [46]:

```
#robot.plot_lambda_convergence(what="gamma", sortby="nkpt", hue="eph_fsmear", tight_layout=True);
```

In [47]:

```
#robot.abifiles[-1].plot_a2f_interpol();
```

In [48]:

```
#robot.abifiles[-1].a2f_qcoarse.plot_tc_vs_mustar(start=0.1, stop=0.2);
```

In [49]:

```
#robot.gridplot_phbands_with_a2f(units="meV");
```

In [50]:

```
#robot.plot_a2fdata_convergence(sortby="nkpt", hue=None);
```

Back to the main Index

In [ ]:

```
```