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 MgB2
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 MgB2 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 MgB2; epw/
contains the input files to run epw.x
on MgB2 .
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 MgB2 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 MgB2 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 settingephwrite=.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.
, andlimag = .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 inMgB2.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 inMgB2.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 asMgB2.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. HereYY
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 thatYY
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
, andikmap
files were generated and saved inprefix.ephmat
directory, the code requires to use the same number of CPUs as the number ofephmatXX
files. If you forget this the code will stop, asking to usenpool
equal to the number ofephmatXX
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 Tc 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 temperatureXX
.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, whereYY
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
, andselecq.fmt
(selecq.fmt
is only reqired ifselecqread = .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 q -point while writing
- Restart from
ephmatXX
files Required files:
prefix.ephmat
directory (which containsegnv, freq, ikmap, ephmatXX
files),selecq.fmt
, andcrystal.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 (annscf
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.