# Superconducting magnesium diboride

Updated by: Hari Paudyal

Warning

The following example aims at providing physically meaningful results. Calculations can therefore take a significant amount of time. For quick calculations, look at the `EPW/test-suite/epw_super/`

folder.

## Running EPW to calculate the anisotropic superconducting gap of MgB_{2}

This example is based on the papers Phys. Rev. B 87, 024505 (2013) and Comp. Phys. Commun. 209, 116 (2016), which you should check for the physical interpretation of the results and a description of the main equations used in this tutorial.

```
The results in this example are obtained with Quantum_ESPRESSO v6.6 and the corresponding EPW v5.3.
```

The MgB_{2} example is located inside the `EPW/examples/mgb2/`

directory. Within this directory you will find three subdirectories: `pp/`

contains the Mg and B pseudopotentials; `phonons/`

contains the input files needed for calculating the phonon dispersions of MgB_{2}; `epw/`

contains the input files to run `epw.x`

on MgB_{2} .

Once `pw.x`

, `ph.x`

, and `epw.x`

has been compiled, we are ready to run the example.

## Calculating phonons

The phonon code of QE requires a ground-state self-consistent calculation:

```
mpirun -np 4 ../../../../bin/pw.x < scf.in > scf.out
```

The electronic band structure and DOS of MgB_{2} are shown in Figure 1 (a-b) below.

We will now proceed with a sequential run. Let us compute the phonon frequencies and eigenmodes using `ph.x`

:

```
mpirun -np 4 ../../../../bin/ph.x < ph.in > ph.out
```

This will provide us with the files containing the dynamical matrices and the variations of the self-consistent potentials at the 28 irreducibles **q**-points of the 6x6x6 uniform **q**-mesh. The next step is to copy `.dyn`

, `.dvscf`

files along with the `_ph0/MgB2.phsave`

folder inside the `save/`

folder. These files contain all the information that EPW needs from QE.

EPW expects these files to be labelled according to the number of the irreducible **q**-point in the list. In order to label the files correctly you can use the python script `pp.py`

. If you decide to proceed this way, just issue the statement:

```
python pp.py
```

from within the `phonons/`

folder. The script will ask you the prefix used for the QE calculations as well as the number of irreducible **q**-points computed. The script will rename all the required files and place them in the folded `save/`

.

Clean up unnecessary files and directories

```
rm -r _ph0/
rm -r MgB2.*
```

The phonon dispersion and PHDOS of MgB_{2} are shown in Figure 1 (c-d) below.

## Running EPW

We first have to perform one scf and one nscf calculation. To do this, we go inside the `epw`

directory and issue:

```
mpirun -np $NPROC ../../../../bin/pw.x < scf.in > scf.out
```

```
mpirun -np $NPROC ../../../../bin/pw.x -npool $NPOOL < nscf.in > nscf.out
```

(NOTE: $NPROC = $NPOOL = 128 is used)

The reason for the non-self consistent calculation is that EPW needs the wavefunctions on the full BZ on a grid between 0 and 1. For the nscf calculation the number of pools `$NPOOL`

has to be the same as the total number of cores `$NPROC`

You can then launch the EPW calculation:

```
mpirun -np $NPROC ../../../bin/epw.x -npool $NPOOL < epw.in > epw.out
```

The number of cores and pools have to be the same as for the `nscf.in`

run

The calculation of superconducting properties is triggered by the flag `eliashberg = .true.`

and will be accompanied by significant I/O.

The following EPW input variables are important in our case:

```
iverbosity = 2 ! 2 = verbose output for the SC part
ephwrite = .true. ! write ephmatXX (XX = pool dependent), egnv, freq, and ikmap files in ``MgB2.ephmat/`` directory
eliashberg = .true. ! solve the ME eqs.
laniso = .true. ! solve the anisotropic ME eqs.
limag = .true. ! solve the imag-axis ME eqs.
lpade = .true. ! solve ME eqs. on the real axis using Pade approximants
nsiter = 500 ! number of self-consistent iterations when solving the ME eqs.
conv_thr_iaxis = 1.0d-4 ! convergence threshold for solving ME eqs. on imag-axis
wscut = 1.0 ! upper limit over Matsubara frequency summation in ME eqs on imag-axis in [eV]
nstemp = 1 ! nr. of temperatures
temps = 10 ! list of temperetures OR (nstemp and temps = tempsmin tempsmax for even space mode) in [K]
muc = 0.16 ! effective Coulomb potential used in the ME eqs.
dvscf_dir = '../phonons/save' ! directory where the save folder generated by the ``pp.py`` python script is stored
nk1 = 6 ! dimensions of the coarse electronic grid
nk2 = 6
nk3 = 6
nq1 = 6 ! dimensions of the coarse phonon grid
nq2 = 6
nq3 = 6
mp_mesh_k = .true. ! use irreduciable electronic fine mesh
nkf1 = 60 ! dimensions of the fine electron grid
nkf2 = 60
nkf3 = 60
nqf1 = 30 ! dimensions of the fine phonon grid
nqf2 = 30
nqf3 = 30
```

We can take a look at the epw.out step by step

Fourier-transform the electron-phonon matrix elements from a coarse 6x6x6 to a dense 60x60x60

**k**-point grid and from a coarse 6x6x6 to a dense 30x30x30**q**-point grid

```
Using uniform q-mesh: 30 30 30
Size of q point mesh for interpolation: 27000
Using uniform MP k-mesh: 60 60 60
Size of k point mesh for interpolation: 20522
Max number of k points per pool: 162
```

We have a **q**-point grid with 27000 **q**-points (no symmetry) and a **k**-point grid with 20522 k-points (IBZ). Since we are doing a parallelization on **k**-point, we will have 162 **k**-points per core (this calculation is done on 128 cores).

Write on disk: the files described below in the

`MgB2.ephmat`

(are produced by setting`ephwrite=.true.`

). The files are formatted and required for solving the Migdal-Eliashberg equations. Because the electron-phonon matrix elements do not depend on the temperature at which the Migdal-Eliashberg equations are solved, the files can be reused in subsequent EPW calculations at different temperatures.the

`ephmatXX`

files (one per CPU) containing the electron-phonon matrix elements within the Fermi window (`fsthick`

) on the dense**k**and**q**grids,the

`freq`

file containing the phonon frequencies on the dense**q**grid,the

`egnv`

file containing the eigenvalues within the Fermi window on the dense**k**grid, andthe

`ikmap`

file containing the index of the**k**-points on the dense (IBZ) grid within the Fermi window.

```
Finish mapping k+sign*q onto the fine irreducibe k-mesh and writing .ikmap file
Nr irreducible k-points within the Fermi shell = 2570 out of 10261
Progression iq (fine) = 100/ 27000
Progression iq (fine) = 200/ 27000
....
....
Progression iq (fine) = 27000/ 27000
Fermi level (eV) = 0.747126479252863D+01
DOS(states/spin/eV/Unit Cell) = 0.340746408495801D+00
Electron smearing (eV) = 0.100000000000000D+00
Fermi window (eV) = 0.400000000000000D+00
Finish writing .ephmat files
```

Solve the anisotropic Migdal-Eliashberg equations on the imaginary frequency axis by setting the keywords

`eliashberg = .true.`

,`laiso = .true.`

, and`limag = .true.`

```
===================================================================
Solve anisotropic Eliashberg equations
===================================================================
....
Electron-phonon coupling strength = 0.7705746
Estimated Allen-Dynes Tc = 17.110246 K for muc = 0.16000
Estimated w_log in Allen-Dynes Tc = 55.592716 meV
Estimated BCS superconducting gap = 2.595027 meV
...
temp( 1) = 10.00000 K
Solve anisotropic Eliashberg equations on imaginary-axis
Total number of frequency points nsiw( 1) = 185
Cutoff frequency wscut = 1.0000
Size of allocated memory per pool: ~= 2.6247 Gb
iter ethr znormi deltai [meV]
1 2.659908E+00 1.747231E+00 5.142639E+00
2 5.652175E-01 1.747922E+00 5.468204E+00
....
10 5.493978E-05 1.743316E+00 5.761882E+00
Convergence was reached in nsiter = 10
```

Perform the analytic continuation of the solutions along the imaginary frequency axis to the real frequency axis by using Padé approximants (

`lpade = .true.`

). Note that the analytic continuation with the iterative procedure (`lacon = .true.`

) is not performed since this is very expensive computationally (hours to days).

```
Pade approximant of anisotropic Eliashberg equations from imaginary-axis to real-axis
Cutoff frequency wscut = 1.0000
pade Re[znorm] Re[delta] [meV]
166 1.754658E+00 4.997289E+00
Convergence was reached for N = 166 Pade approximants
```

In the following we describe the various output files, and we will use `XX`

to indicate the temperature at which the anisotropic gap equations are solved.

`MgB2.a2f`

This file contains the Eliashberg spectral function and cumulative electron-phonon coupling strength as a function of frequency \(\omega\) (meV), calculated using various smearing parameters (see Figure 2 below corresponding to the first smearing value).

`MgB2.a2f_iso`

This file contains the Eliashberg spectral function as a function of frequency \(\omega\) (meV) . The second column is the Eliashberg spectral function corresponding to the first smearing value in

`MgB2.a2f`

. The remaining columns in`MgB2.a2f_iso`

(3 x number of atoms) contain the mode-resolved Eliashberg spectral function (there is no specific information on which modes correspond to which atomic species).

`MgB2.phdos`

This file contains the total phonon density of states as a function of frequency \(\omega\) (meV) calculated using various smearing parameters. The smearing values are the same as in

`MgB2.a2f`

.

`MgB2.phdos_proj`

This file contains the phonon density of states as a function of frequency \(\omega\) (meV). The second column is the total phonon density of states corresponding to the first smearing value in

`MgB2.phdos`

. The remaining columns in`MgB2.phdos_proj`

(3 x number of atoms) contain the mode-resolved phonon density of states (there is no specific information on which modes correspond to which atomic species).

`MgB2.imag_aniso_XX`

This file contains 5 columns: the frequency along the imaginary-axis \(\omega\) (eV), the energy eigenvalue \(E_{n{\bf k}}\) for band \(n\) and wavevector \({\bf k}\) relative to the Fermi energy (eV), the quasiparticle renormalization function \(Z_{n{\bf k}}\), the superconducting gap \(\Delta_{n{\bf k}}\) (eV), and the quasiparticle renormalization function in the normal state \(Z_{n{\bf k}}\). The gap is shown in the left panel of Figure 3 below.

`MgB2.pade_aniso_XX`

This file contains 6 columns: the frequency along the imaginary-axis \(\omega\) (eV), the energy eigenvalue \(E_{n{\bf k}}\) for band \(n\) and wavevector \({\bf k}\) relative to the Fermi energy (eV), the real part of the quasiparticle renormalization \(ReZ_{n{\bf k}}\), the imaginary part of the quasiparticle renormalization \(ImZ_{n{\bf k}}\), the real part of the superconducting gap \(Re\Delta_{n{\bf k}}\) (eV), and the imaginary part of the superconducting gap \(Im\Delta_{n{\bf k}}\) (eV). The superconducting gap is calculated along the real axis and is obtained via Padé approximants. This file is written if

`iverbosity=2`

in the input file, and the result is shown in the right panel of Figure 3 below.

`MgB2.acon_aniso_XX`

This files are generated by setting

`lacon = .true.`

, which contain similar information as`MgB2.pade_aniso_XX`

.

`MgB2.imag_aniso_gap0_XX`

This file contains the distribution of the anisotropic superconducting gap \(\Delta_{n{\bf k}}(\omega=0)\) (eV) on the Fermi surface, as obtained from the imaginary-axis calculation. This is shown in Figure 4 below.

`MgB2.pade_aniso_gap0_XX`

This file contains similar information as in

`MgB2.imag_aniso_gap0_XX`

, but in this case the gap is calculated on the real axis using Padé approximants.

`MgB2.qdos_XX`

This file contains the quasiparticle density of states in the superconducting state relative to the DOS in the normal state \(N_{\rm S}(\omega)/N_{\rm N}(\omega)\) as a function of frequency \(\omega\) (eV). This is shown in Figure 5 below.

`MgB2.imag_aniso_gap_FS_XX`

This file contains the superconducting gap on the Fermi surface. The first three columns contain the 3 Cartesian coordinates of each

**k**-point, the 4th column is the band index within the chosen energy window, the 5th column is the energy difference with respect to the Fermi level. Only states within`fsthick`

of the Fermi energy are considered. The 6th column is the superconducting gap in meV.

`MgB2.imag_aniso_gap0_XX_YY.cube`

This file contains the same information as

`MgB2.imag_aniso_gap_FS_XX`

, in a format compatible with VESTA. This is suitable for mapping the gap on the Fermi surface. Here`YY`

represents the band number within the chosen energy window during the EPW calculation. A 3D rendering of superconducting gap on top of the Fermi surface sheets is shown in Figure 6 below. See section**Note for expert users**below if you want to reproduce this plot.

`MgB2.lambda_aniso`

This file contains 4 columns: KS eigenvalue \(E_{nk}\) relative to the Fermi level (eV), \(\lambda_{n{\bf k}}\), the wavevector \({\bf k}\), and the band index \(n\).

`MgB2.lambda_k_pairs`

This file contains the normalized distribution of the anisotropic electron-phonon coupling strength \(\lambda_{n{\bf k}}\) on the Fermi surface. This is shown in Figure 7 below.

`MgB2.lambda_pairs`

This file contains the normalized distribution of the anisotropic electron-phonon coupling strength \(\lambda_{n{\bf k},n'{\bf k}'}\) on the Fermi surface.

`MgB2.fe_XX`

This file contains the free energy in the superconducting state.

`MgB2.lambda_FS`

This file contains \(\lambda_{n{\bf k}}\) on the Fermi surface. The first three columns represent the Cartesian coordinates of the

**k**-point, the 4th column is the band index within the specified energy window, the 5th column is the energy relative to the Fermi level. Only states within the`fsthick`

of the Fermi energy are considered. The 6th column is \(\lambda_{n{\bf k}}\).

`MgB2.lambda_YY.cube`

This file contains the same information as

`MgB2.lambda_FS`

but can be read by VESTA for direct visualization. Note that`YY`

is the band index within the specified energy window.

## Important Notes

**ephwrite=.true.**does not work with random**k**or**q**grids and requires**nkf1**,**nkf2**,**nkf3**to be multiple of**nqf1**,**nqf2**,**nqf3**.**mp_mesh k = .true.**specifies that only the irreducible points for the dense**k**grid are used. This significantly decreases the computational cost when solving the Migdal-Eliashberg equations.If the Migdal-Eliashberg equations are solved in a separate run from the one in which the

`ephmatXX`

,`freq`

,`egnv`

, and`ikmap`

files were generated and saved in`prefix.ephmat`

directory, the code requires to use the same number of CPUs as the number of`ephmatXX`

files. If you forget this the code will stop, asking to use`npool`

equal to the number of`ephmatXX`

files.**lpade = .true.**requires**limag = .true.****lacon = .true.**requires both**limag = .true.**and**lpade = .true.****wscut**gives the upper limit (in eV) of the summation over the frequencies on the imaginary axis in the Migdal-Eliashberg equations (**limag = .true.**). Note that the input variable**wscut**is ignored if the number of frequency points is given using the input variable**nswi**. In this case, the number of frequency points in the summation is the same irrespective of the temperature.**temps(1)**,**temps(2)**, … define the temperatures at which the Migdal-Eliashberg equations are evaluated. Note that the temperatures can also be defined using**nstemp**,**temps(1) = min_temp**,**temps(2)=max_temp**input variables.If temperatures larger than the critical temperature T

_{c}are specified in the input file, the code will print a WARNING for such a temperature since the Migdal-Eliashberg equations may not have solution at that point.**imag_read**works if**limag = .true.**and**laniso = .true.**and it allows the code to read from file the superconducting gap and renormalization function on the imaginary axis at specific temperature`XX`

from file**MgB2.imag_aniso_XX**. The temperature is specified as**temps(1) = XX**in the EPW input file.**imag_read**can be used to:solve the anisotropic Migdal-Eliashberg equations on the imaginary axis at temperatures greater than

`XX`

using as a starting point the superconducting gap estimated at temperature`XX`

.obtain the solutions of the Migdal-Eliashberg equations on the real axis with

**lpade = .true.**or**lacon = .true.**starting from the imaginary axis solutions at temperature`XX`

.write to file the superconducting gap on the Fermi surface in cube format at temperature

`XX`

for**iverbosity = 2**. The generated output files are**MgB2.imag_aniso_gap0_XX_YY.cube**, where`YY`

is the band number within the chosen energy window during the EPW calculation.

## Restart tricks (this requires to use the same number of cores used as in the original run)

- Restart from an interrupted
**q**-point while writing`ephmatXX`

files Required files:

`prefix.epmatwp`

,`prefix.ukk`

,`crystal.fmt`

,`epwdata.fmt`

,`dmedata.fmt/vmedata.fmt`

,`restart.fmt`

, and`selecq.fmt`

(`selecq.fmt`

is only reqired if`selecqread = .true.`

otherwise it will be re-created).Input setup:

ep_coupling = .true. elph = .true. kmaps = .true. epbwrite = .false. epbread = .false. epwwrite = .false. epwread = .true. ! reads .epmatwp, .fmt, .ukk files wannierize = .false. ephwrite = .true.

- Restart from an interrupted
- Restart from
`ephmatXX`

files Required files:

`prefix.ephmat`

directory (which contains`egnv, freq, ikmap, ephmatXX`

files),`selecq.fmt`

, and`crystal.fmt`

Input setup:

ep_coupling = .false. elph = .false. kmaps = .true. epbwrite = .false. epbread = .false. epwwrite = .false. epwread = .true. wannierize = .false. ephwrite = .false.

- Restart from

## Note for expert users

To reproduce Figure 6 above (showing the map of the superconducting gap on the Fermi surface), you can follow one of the two ways described below;

- Fermi surface calculation with EPW in
`.cube`

format (fast): If you want to calculate the Fermi surfaces only, you can restart your calculation by reading

`MgB2.epmatwp`

,`MgB2.ukk`

, and`*.fmt`

files using following input flags;

- Fermi surface calculation with EPW in

```
fermi_plot =.true
mp_mesh_k = .true. ! works for both .true. and .false.
nkf1 = 60 ! dimensions of the fine electron grid, should be the same used in SC calculations
nqf2 = 60
nqf3 = 60
nqf1 = 1 ! dimensions of the fine phonon grid, we only need at \Gamma for the Fermi surface calculation
nqf2 = 1 ! works for any dimesions
nqf3 = 1
```

Trick: If you use “fermi_plot = .true.” in your input file for the superconducting calculation, both Fermi surface and superconducting calculations can be done in a single run.

This will generate **MgB2.fs_YY.cube** files that can be plotted with VESTA. To plot the Fermi surface, import **MgB2.fs_YY.cube** one by one in ‘Volumetric Data’ and choose ‘Multiply to current data’(if you have more than one files). Further, to color the Fermi surface with the supercondicting gap, import **MgB2_imag_aniso_gap0_XX.YY.cube** (obtained from the anisotropic SC calculation with `iverbosity = 2`

) in ‘Volumetric Data’ –> ‘Surface coloring’ one by one and choose ‘Add to current data’.

**OR**

Fermi surface calculation (in

`.xsf`

format) using the python script from the following tar file (an`nscf`

calculation without symmetry on the full uniform grid is expensive):

The following video explains how to use the scripts in conjunction with VESTA to obtain Figure 6.

This video tutorial can also be found here.