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.

../_images/MgB2_BS_PH_qe6.6.png

Figure 1: (a) Electronic band structure, (b) DOS, (c) phonon dispersion relations, and (d) PHDOS of bulk MgB2. The calculations are executed using the experimental lattice parameters. The inelastic X-ray scattering data (black dots in (d)) at 300~K are from A. Shukla et al..

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.

    1. the ephmatXX files (one per CPU) containing the electron-phonon matrix elements within the Fermi window (fsthick) on the dense k and q grids,

    2. the freq file containing the phonon frequencies on the dense q grid,

    3. the egnv file containing the eigenvalues within the Fermi window on the dense k grid, and

    4. the 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).

../_images/MgB2_a2F_qe6.6.png

Figure 2: Calculated isotropic Eliashberg spectral function \(\alpha^2 F(\omega)\) of MgB2 (blue solid line), and cumulative electron-phonon coupling strength \(\lambda\) (black dashed line). The results are integrated on a homogeneous fine k-point grid containing 60x60x60 points, and a fine homogeneous q-mesh with 30x30x30 points.

  • 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.

../_images/MgB2_aniso2_qe6.6.png

Figure 3: Calculated energy-dependent superconducting gap of MgB2 from the anisotropic Migdal-Eliashberg theory at 10~K (blue), 35~K (green) and 45~K (red). The left panel shows the gap along the imaginary energy axis and the right panel is along the real energy axis (approximate analytic continuation via Padé functions). The Coulomb potential is set to \(\mu_c^*=0.16\). Only the KS states n**k** with an energy within 0.1~eV from the Fermi level are shown.

  • 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.

../_images/MgB2_gaps_qe6.6.png

Figure 4: Calculated anisotropic superconducting gap \(\Delta_{n\mathbf{k}}(\omega=0)\) of MgB2 on the Fermi surface as a function of temperature. The Coulomb potential is set to \(\mu_c^*=0.16\). For each temperature the histograms indicate the number of states on the Fermi surface with that superconducting gap energy. The dashed lines are guides to the eye. The two superconducting gaps vanish at the calculated critical superconducting temperature \(T_{\rm c}\) of 53~K [note this value is about 36% higher than in the experiment].

  • 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.

../_images/MgB2_tunnel_qe6.6.png

Figure 5: Quasiparticle DOS in the superconducting state for various temperatures [obtained using \(\Delta(\omega=0)\)]. The dashed black line is the normal state DOS taken from the top of Figure 1 and normalized to 1 at the Fermi level. The superconducting DOS \(N_{\rm S}(\omega)\) was scaled so that its high energy tail becomes 1 beyond the highest phonon energy. The two-gap nature of MgB2 is clearly seen, and can be associated with the \(\sigma\) and \(\pi\) sheets on its Fermi surface, as shown in Figure 6 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.

../_images/MgB2_FS_band1and2.png

Figure 6: The superconducting energy gap of MgB2 calculated at T=10 K mapped on the Fermi surface. The Fermi surface consists of two \(\sigma\) sheets along the \(\Gamma-\Gamma\) lines, and two \(\pi\) sheets along the K-M and the H-L lines.

  • 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.

../_images/MgB2_k_pairs_qe6.6.png

Figure 7: Distribution of the electron-phonon coupling strength \(\lambda_{\bf k}\) for MgB2.

  • 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 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:
    1. 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.

    2. 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.

    3. 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)

  1. 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.
    
  2. 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.
    

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;

  1. 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_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

  1. 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):

fermi.tar

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.