Superconducting properties¶
Note
Hands-on based on QE-v7.5 and EPW-v6.0
In this tutorial we compute superconducting properties of three metals by solving the Migdal-Eliashberg equations with EPW: the isotropic case in fcc Pb (Exercise 1), the anisotropic case in MgB2 (Exercise 2), and the anisotropic full-bandwidth case in bcc Nb with sparse-ir Matsubara sampling (Exercise 3). Theory background: Phys. Rev. B 87, 024505 (2013) and Phys. Rev. B 110, 064505 (2024).
Prerequisites¶
Quantum ESPRESSO (pw.x, ph.x) and EPW
(epw.x); the helper scripts pp.py,
kmesh.pl and nscf2supercond.x; norm-conserving
pseudopotentials for Pb, Mg, B, Nb in a shared pseudo/
directory. Exercise 3 additionally needs a sparse-ir object file
generated with the Python sparse-ir package (see
sparse-ir-fortran).
Exercise 1: Isotropic superconductivity in fcc Pb¶
We compute the superconducting gap and critical temperature of fcc Pb by solving the isotropic Migdal-Eliashberg equations, and then extract the nesting function and the mode-resolved electron-phonon coupling along a q-path.
Step 1¶
Run a self-consistent calculation on a homogeneous 14x14x14 k-point grid and a phonon calculation on a homogeneous 3x3x3 q-point grid. Move to the phonon directory and submit the job:
cd exercise1/phonon
mpirun -np 8 pw.x -nk 4 -in scf.in > scf.out
mpirun -np 8 ph.x -nk 4 -in ph.in > ph.out
Note. The input degauss used here is large, reduce it for your system.
scf.in
&control
calculation = 'scf'
restart_mode = 'from_scratch'
prefix = 'pb'
pseudo_dir = '../../pseudo/'
outdir = './'
/
&system
ibrav = 2
celldm(1) = 9.222558
nat = 1
ntyp = 1
ecutwfc = 90
ecutrho = 360
occupations = 'smearing'
smearing = 'mp'
degauss = 0.025
/
&electrons
diagonalization = 'david'
mixing_beta = 0.7
conv_thr = 1.0d-12
/
ATOMIC_SPECIES
Pb 207.2 pb_s.UPF
ATOMIC_POSITIONS crystal
Pb 0.00 0.00 0.00
K_POINTS automatic
14 14 14 0 0 0
ph.in
&inputph
prefix = 'pb'
fildyn = 'pb.dyn'
fildvscf = 'dvscf'
tr2_ph = 1.0d-17
ldisp = .true.
nq1 = 3
nq2 = 3
nq3 = 3
/
During the run, check the irreducible q-point grid written in
ph.out:
Dynamical matrices for ( 3, 3, 3) uniform grid of q-points
( 4 q-points):
N xq(1) xq(2) xq(3)
1 0.000000000 0.000000000 0.000000000
2 -0.333333333 0.333333333 -0.333333333
3 0.000000000 0.666666667 0.000000000
4 0.666666667 -0.000000000 0.666666667
Step 2¶
Collect the .dyn, .dvscf, and patterns files
into a new save/ directory using the helper script
pp.py distributed with EPW:
python3 /path/to/EPW/bin/pp.py
The script will ask for the prefix of the calculation (here
"pb").
Step 3¶
Run a non-self-consistent calculation on a uniform 3x3x3 Gamma-centered k-grid written in crystal coordinates, then perform the EPW calculation that Wannierises the electronic structure and solves the isotropic Migdal-Eliashberg equations.
cd ../epw1-2
mpirun -np 8 pw.x -nk 8 -in scf.in > scf.out
mpirun -np 8 pw.x -nk 8 -in nscf.in > nscf.out
mpirun -np 8 epw.x -nk 8 -in epw1.in > epw1.out
nscf.in
&control
calculation = 'bands'
restart_mode = 'from_scratch'
prefix = 'pb'
pseudo_dir = '../../pseudo/'
outdir = './'
verbosity = 'high'
/
&system
ibrav = 2
celldm(1) = 9.222558
nat = 1
ntyp = 1
ecutwfc = 90
ecutrho = 360
occupations = 'smearing'
smearing = 'mp'
degauss = 0.025
nbnd = 10
/
&electrons
diagonalization = 'david'
mixing_beta = 0.7
conv_thr = 1.0d-12
/
ATOMIC_SPECIES
Pb 207.2 pb_s.UPF
ATOMIC_POSITIONS crystal
Pb 0.000000000 0.000000000 0.000000000
K_POINTS crystal
27
0.00000000 0.00000000 0.00000000 3.703704e-02
0.00000000 0.00000000 0.33333333 3.703704e-02
...
The full 27-point k-list for the uniform 3x3x3 crystal grid can be
generated with the kmesh.pl utility distributed with
Wannier90:
/path/to/wannier90/utility/kmesh.pl 3 3 3
epw1.in
&inputepw
prefix = 'pb'
outdir = './'
dvscf_dir = '../phonon/save' ! directory with .dyn, .dvscf and
! prefix.phsave/patterns.xx.yy files
ep_coupling = .true.
elph = .true.
epwwrite = .true.
epwread = .false.
etf_mem = 0
wannierize = .true.
nbndsub = 4
bands_skipped = 'exclude_bands = 1-5'
num_iter = 300
dis_froz_min = -3
dis_froz_max = 13.5
dis_win_max = 50
proj(1) = 'Pb:sp3'
wdata(1) = 'bands_plot = .true.'
wdata(2) = 'begin kpoint_path'
wdata(3) = 'G 0.00 0.00 0.00 X 0.00 0.50 0.50'
wdata(4) = 'X 0.00 0.50 0.50 W 0.25 0.50 0.75'
wdata(5) = 'W 0.25 0.50 0.75 L 0.50 0.50 0.50'
wdata(6) = 'L 0.50 0.50 0.50 K 0.375 0.375 0.75'
wdata(7) = 'K 0.375 0.375 0.75 G 0.00 0.00 0.00'
wdata(8) = 'G 0.00 0.00 0.00 L 0.50 0.50 0.50'
wdata(9) = 'end kpoint_path'
wdata(10)= 'bands_plot_format = gnuplot'
fsthick = 0.4 ! Fermi window thickness [eV]
degaussw = 0.1 ! energy-conserving delta smearing [eV]
degaussq = 0.5 ! phonon-sum smearing [meV]
ephwrite = .true.
eliashberg = .true.
liso = .true. ! isotropic Migdal-Eliashberg
limag = .true. ! imaginary-axis equations
lpade = .true. ! Pade continuation to real axis
lacon = .true. ! iterative real-axis continuation
nsiter = 500
npade = 20
conv_thr_iaxis = 1.0d-4
conv_thr_racon = 1.0d-3
wscut = 0.1 ! Matsubara cutoff [eV]
muc = 0.1 ! effective Coulomb parameter
temps = 0.3 0.9 1.5 2.1 2.7 3.3 3.9 4.0 4.2 4.4 4.5 4.6
nk1 = 3
nk2 = 3
nk3 = 3
nq1 = 3
nq2 = 3
nq3 = 3
mp_mesh_k = .true.
nkf1 = 18
nkf2 = 18
nkf3 = 18
nqf1 = 18
nqf2 = 18
nqf3 = 18
/
Note 1. EPW calculations with ephwrite = .true. require that the fine k and q grids be commensurate, i.e.
nkf1,nkf2,nkf3 must be integer multiples of
nqf1,nqf2,nqf3.
Note 2. If pw.x silently adds extra k-points to the non-self-consistent list, replace calculation='nscf' with
calculation='bands' as shown above.
Note 3. The non-self-consistent run requires the charge density produced by the preceding self-consistent run. Either rerun the SCF step as shown above, or create pb.save/ and copy
charge-density.dat and data-file-schema.xml from the
phonon/pb.save/ directory.
With this input, EPW performs the following steps:
Fourier-transforms the electron-phonon matrix elements from the coarse 3x3x3 grids to a dense 18x18x18 k-grid and 18x18x18 q-grid. In
epw1.outthis appears asUsing uniform q-mesh: 18 18 18 Size of q point mesh for interpolation: 5832 Using uniform MP k-mesh: 18 18 18 Size of k point mesh for interpolation: 390
Pre-computes the q-points that lie inside the Fermi window set by
fsthick(a q-point is selected if at least one band eigenvalue at k+q falls within the window for some k).Writes the files
pb.ephmat/ephmatXX,freq,egnv, andikmapthat are required to solve the Migdal-Eliashberg equations. This step is driven byephwrite = .true..Solves the isotropic Migdal-Eliashberg equations on the imaginary frequency axis self-consistently at each temperature listed in
temps. The self-consistent renormalisation function Z and gap function Delta enter a coupled integral equation on the Matsubara frequencies; the iteration ends when eitherconv_thr_iaxisornsiteris reached. Expect to see output like=================================================================== Solve isotropic Eliashberg equations =================================================================== Electron-phonon coupling strength = 1.4242792 Estimated Allen-Dynes Tc = 3.168499 K for muc = 0.10000 Estimated w_log in Allen-Dynes Tc = 2.519536 meV Estimated BCS superconducting gap = 0.480551 meV temp( 1) = 0.30000 K Solve isotropic Eliashberg equations on imaginary-axis Total number of frequency points nsiw( 1) = 616 iter ethr znormi deltai [meV] 1 3.130497E+00 2.278915E+00 6.137465E-01 ... Convergence was reached in nsiter = 8Performs the analytic continuation of the imaginary-axis solution to the real frequency axis with Pade approximants (
lpade = .true.) and with the iterative procedure (lacon = .true.). The iterative procedure is more expensive but less sensitive to the number of Matsubara frequencies kept.
At the end of the calculation, EPW writes one set of output files per
temperature. In the following XX stands for the temperature tag.
Step 4¶
The following files are produced:
pb.a2f– isotropic Eliashberg spectral function (alpha2`F) and cumulative electron-phonon coupling strength as a function of frequency, for several phonon smearing values. Written when ``eliashberg = .true.`.pb.a2f_proj– the same spectral function resolved by atomic contributions (3 x number of atoms extra columns).pb.imag_iso_XX– three columns: the Matsubara frequency (eV), the renormalisation Z, and the gap Delta (eV). Written whenliso = limag = .true..pb.pade_iso_XX– five columns giving the real frequency, real and imaginary parts of Z, real and imaginary parts of Delta (eV). Written whenlpade = .true..pb.acon_iso_XX– same layout aspb.pade_iso_XX, obtained by the iterative continuation (lacon = .true.).
Plot the gap on the imaginary and real axes at T = 0.3 K. A sample
gnuplot script fig1.plt is included in the tutorial
directory:
gnuplot fig1.plt
evince fig1.pdf
Superconducting gap of fcc Pb at 0.3 K. Left: imaginary-axis data (columns 1:3 of pb.imag_iso_000.30). Right: real-axis data from pb.pade_iso_000.30 (Pade) and pb.acon_iso_000.30 (iterative).
Step 5¶
Extract the leading-edge value of the gap at each temperature by
reading the second line (the smallest Matsubara frequency) of every
pb.imag_iso_XX file, pulling column 4 (the gap in eV),
converting to meV, and pairing the result with the temperature tag
XX from the filename. A small shell or Python helper does
the job; any script that produces a two-column
pb.imag_iso_gap0 file (temperature in K, gap in meV) will
work with the gnuplot templates shipped with the tutorial.
./script_gap0_imag
gnuplot fig2.plt
evince fig2.pdf
Isotropic superconducting gap of Pb on the imaginary axis as a function of temperature. Left: result of the tutorial settings. Right: converged reference from Phys. Rev. B 87, 024505 (2013).
The same leading-edge extraction applied to pb.pade_iso_XX
and pb.acon_iso_XX produces the real-axis counterparts; the
three curves should agree within the numerical precision of the
continuation.
Step 6¶
Close to Tc the superconducting gap becomes small and the
Migdal-Eliashberg equations reduce to a linear eigenvalue problem for
the gap function. The critical temperature is the value at which the
leading eigenvalue of the resulting kernel equals unity. EPW solves
this problem directly, using only the Eliashberg spectral function
pb.a2f, when tc_linear = .true..
mpirun -np 1 epw.x -nk 1 -in epw2.in > epw2.out
epw2.in
&inputepw
prefix = 'pb'
outdir = './'
ep_coupling = .false.
elph = .false.
epwwrite = .false.
epwread = .true.
wannierize = .false.
ephwrite = .false.
fsthick = 0.4
degaussw = 0.1
degaussq = 0.5
eliashberg = .true.
liso = .true.
limag = .true.
lpade = .false.
lacon = .false.
fila2f = 'pb.a2f'
tc_linear = .true.
tc_linear_solver = 'power' ! 'power' or 'lapack'
nstemp = 21
temps = 0.25 5.25 ! evenly spaced points when two values
! are given with nstemp > 2
nsiter = 500
conv_thr_iaxis = 1.0d-4
wscut = 0.1
muc = 0.1
nk1 = 3
nk2 = 3
nk3 = 3
nq1 = 3
nq2 = 3
nq3 = 3
nkf1 = 18
nkf2 = 18
nkf3 = 18
nqf1 = 18
nqf2 = 18
nqf3 = 18
mp_mesh_k = .true.
/
Note. When the linearised equations are solved from the Eliashberg spectral function the ephmatXX, freq,
egnv, and ikmap files are not used. A single MPI rank
suffices. This shortcut is only available for the isotropic equations.
Extract the temperature-dependence of the leading eigenvalue from the
summary table in epw2.out (the block printed under the header
Max. eigenvalue): the first column is the temperature in K
and the second column is the maximum eigenvalue. Save the pairs to
data_max_eigenvalue.dat and plot.
./script_max_eigenvalue
gnuplot fig4.plt
evince fig4.pdf
Tc is the temperature at which the extracted maximum eigenvalue crosses unity.
Step 7¶
The mode-resolved electron-phonon coupling along a q-path and the
Fermi-surface nesting function are useful to understand which phonon
branches dominate the coupling. Reuse the interpolated matrix elements
from Step 3 by symlinking pb.epmatwp and the *.fmt
files into a new directory:
cd ../epw3-4
ln -s ../epw1-2/pb.epmatwp .
ln -s ../epw1-2/wigner.fmt .
ln -s ../epw1-2/vmedata.fmt .
ln -s ../epw1-2/epwdata.fmt .
ln -s ../epw1-2/crystal.fmt .
ln -s ../epw1-2/pb.ukk .
mpirun -np 8 epw.x -nk 8 -in epw3.in > epw3.out
mpirun -np 8 epw.x -nk 8 -in epw4.in > epw4.out
epw3.in
&inputepw
prefix = 'pb'
outdir = './'
dvscf_dir = '../phonon/save'
ep_coupling = .true.
elph = .true.
epwwrite = .false.
epwread = .true.
etf_mem = 0
wannierize = .false.
nbndsub = 4
fsthick = 0.4
degaussw = 0.1
degaussq = 0.5
phonselfen = .true. ! compute the phonon self-energy
delta_approx = .true. ! double-delta approximation
nk1 = 3
nk2 = 3
nk3 = 3
nq1 = 3
nq2 = 3
nq3 = 3
filqf = 'pb_band.kpt'
nkf1 = 18
nkf2 = 18
nkf3 = 18
/
The file epw4.in is identical except that
phonselfen is replaced by nest_fn = .true.:
epw4.in (differences only)
nest_fn = .true. ! calculate the nesting function
Note. When any of elecselfen, phonselfen,
specfun_el, or specfun_ph is true the full
(non-reduced) fine k-mesh is used, so mp_mesh_k must be
false (or left at its default).
Phonon self-energy and mode-resolved coupling. With
phonselfen = .true. and delta_approx = .true., EPW
evaluates, for every q-point on the path specified by filqf,
the phonon linewidth and the branch-resolved electron-phonon coupling
strength lambdaq,nu in the double-delta approximation:
the coupling for each phonon branch nu is obtained as a Brillouin-zone
integral of the squared matrix element divided by the phonon frequency,
with two delta functions restricting the initial and final electronic
states to the Fermi surface.
A typical excerpt from epw3.out looks like:
===================================================================
Phonon (Imaginary) Self-Energy in the Migdal Approximation
===================================================================
Fermi Surface thickness = 0.400000 eV
Golden Rule strictly enforced with T = 0.025852 eV
Gaussian Broadening: 0.100000 eV, ngauss= 1
DOS = 0.296296 states/spin/eV/Unit Cell at Ef= 11.757717 eV
ismear = 1 iq = 43 coord.: 0.37500 0.37500 0.75000 ...
-------------------------------------------------------------------
lambda___( 1 )= 0.245003 gamma___= 0.006201 meV omega= 5.2143 meV
lambda_tr( 1 )= 0.215219 gamma_tr= 0.005447 meV omega= 5.2143 meV
lambda___( 2 )= 0.223415 gamma___= 0.008930 meV omega= 6.5530 meV
lambda_tr( 2 )= 0.111656 gamma_tr= 0.004463 meV omega= 6.5530 meV
lambda___( 3 )= 0.265418 gamma___= 0.016608 meV omega= 8.1988 meV
lambda_tr( 3 )= 0.241549 gamma_tr= 0.015114 meV omega= 8.1988 meV
lambda___( tot )= 0.733836
lambda_tr( tot )= 0.568425
-------------------------------------------------------------------
The default temperature for the printout is 300 K, but
lambdaq,nu within the double-delta approximation is
independent of temperature. The full set of branch- and q-point-resolved
values is written to lambda.phself.300.000K.
Nesting function. The nesting function
fnest`(q) measures how strongly pairs of states separated
by q overlap on the Fermi surface. With ``nest_fn = .true.`,
EPW evaluates, for every q-point on the path, the Brillouin-zone
integral that counts such k and k+q pairs restricted to the Fermi level
by two delta functions. The corresponding excerpt of epw4.out
is:
===================================================================
Nesting Function in the double delta approx
===================================================================
Fermi Surface thickness = 0.400000 eV
Gaussian Broadening: 0.100 eV, ngauss= 1
DOS = 0.296296 states/spin/eV/Unit Cell at Ef= 11.757717 eV
iq = 1 coord.: 0.00000 0.00000 0.00000 wt: 1.00000
-------------------------------------------------------------------
Nesting function (q)= 0.417115E+03 [Adimensional]
-------------------------------------------------------------------
Extract the q-path-resolved nesting function from epw4.out into
a two-column file pb.nesting_fn (q-point index, nesting
function) with a small grep/awk pipeline; the gnuplot
template fig5.plt then plots fnest`(q) together
with the branch-resolved lambda:sub:`q,nu read from
lambda.phself.300.000K.
./script_nesting_fn
gnuplot fig5.plt
evince fig5.pdf
q-point dependence along a high-symmetry path of the Fermi-surface nesting function fnest`(q) and of the branch-resolved electron-phonon coupling lambda:sub:`q,nu for the three acoustic and optical modes of Pb.
Note. Comparing fnest`(q) with lambda:sub:`q,nu is a useful diagnostic for decomposing the total isotropic coupling lambda into the branches and q-regions that contribute the most.
Exercise 2: Anisotropic superconductivity in MgB2¶
MgB2 is the standard test case for anisotropic Migdal-Eliashberg calculations because of the two distinct superconducting gaps that open on the sigma- and pi-derived Fermi sheets. This exercise solves the anisotropic gap equations in both the Fermi-surface restricted (FSR) and full-bandwidth (FBW) formulations, computes the quasiparticle density of states, and optionally lets the chemical potential shift self-consistently with temperature.
Move to the Exercise 2 directory:
cd ../../exercise2
Step 1¶
Run SCF on a 12x12x12 grid and phonons on a 3x3x3 q-grid:
cd phonon
mpirun -np 24 pw.x -nk 6 -in scf.in > scf.out
mpirun -np 24 ph.x -nk 6 -in ph.in > ph.out
Note. The smearing used in this tutorial is larger than in published work to keep the run time short. Converge it before production runs.
scf.in
&control
calculation = 'scf'
restart_mode = 'from_scratch'
prefix = 'mgb2'
pseudo_dir = '../../pseudo/'
outdir = './'
/
&system
ibrav = 4
celldm(1) = 5.8260252227888
celldm(3) = 1.1420694129095
nat = 3
ntyp = 2
ecutwfc = 40
smearing = 'mp'
occupations = 'smearing'
degauss = 0.05
/
&electrons
diagonalization = 'david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-9
/
ATOMIC_SPECIES
Mg 24.305 Mg.pz-n-vbc.UPF
B 10.811 B.pz-vbc.UPF
ATOMIC_POSITIONS crystal
Mg 0.000000000 0.000000000 0.000000000
B 0.333333333 0.666666667 0.500000000
B 0.666666667 0.333333333 0.500000000
K_POINTS automatic
8 8 8 0 0 0
ph.in
&inputph
prefix = 'mgb2'
fildyn = 'mgb2.dyn.xml'
fildvscf = 'dvscf'
tr2_ph = 1.0d-16
ldisp = .true.
nq1 = 3
nq2 = 3
nq3 = 3
/
Step 2¶
python3 /path/to/EPW/bin/pp.py
Enter the prefix "mgb2" when prompted.
Step 3¶
Run the non-self-consistent calculation on a 6x6x6 uniform Gamma-centered grid and launch the EPW solver for the anisotropic gap in the FSR approximation.
cd ../epw1-FSR
mpirun -np 48 pw.x -nk 48 -in scf.in > scf.out
mpirun -np 48 pw.x -nk 48 -in nscf.in > nscf.out
mpirun -np 48 epw.x -nk 48 -in epw1.in > epw1.out
The 6x6x6 k-list for nscf.in can be generated with
kmesh.pl:
/path/to/wannier90/utility/kmesh.pl 6 6 6
nscf.in
&control
calculation = 'nscf'
prefix = 'mgb2'
pseudo_dir = '../../pseudo/'
outdir = './'
verbosity = 'high'
/
&system
ibrav = 4
celldm(1) = 5.8260252227888
celldm(3) = 1.1420694129095
nat = 3
ntyp = 2
ecutwfc = 40
smearing = 'mp'
occupations = 'smearing'
degauss = 0.05
nbnd = 16
/
&electrons
diagonalization = 'david'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-9
/
ATOMIC_SPECIES
Mg 24.305 Mg.pz-n-vbc.UPF
B 10.811 B.pz-vbc.UPF
ATOMIC_POSITIONS crystal
Mg 0.000000000 0.000000000 0.000000000
B 0.333333333 0.666666667 0.500000000
B 0.666666667 0.333333333 0.500000000
K_POINTS crystal
216
0.00000000 0.00000000 0.00000000 4.629630e-03
0.00000000 0.00000000 0.16666667 4.629630e-03
...
epw1.in
&inputepw
prefix = 'mgb2'
outdir = './'
dvscf_dir = '../phonon/save'
ep_coupling = .true.
elph = .true.
epwwrite = .true.
epwread = .false.
etf_mem = 0
wannierize = .true.
nbndsub = 5
num_iter = 500
dis_froz_max = 8.8
proj(1) = 'B:pz'
proj(2) = 'f=0.5,1.0,0.5:s'
proj(3) = 'f=0.0,0.5,0.5:s'
proj(4) = 'f=0.5,0.5,0.5:s'
iverbosity = 2
fsthick = 0.4
degaussw = 0.1
degaussq = 0.5
fermi_plot = .true.
ephwrite = .true.
eliashberg = .true.
laniso = .true.
limag = .true.
lpade = .true.
lacon = .false.
nsiter = 500
conv_thr_iaxis = 1.0d-4
wscut = 0.5
muc = 0.05
nstemp = 10
temps = 10 55 ! evenly spaced nstemp temperature points
nk1 = 6
nk2 = 6
nk3 = 6
nq1 = 3
nq2 = 3
nq3 = 3
mp_mesh_k = .true.
nkf1 = 40
nkf2 = 40
nkf3 = 40
nqf1 = 20
nqf2 = 20
nqf3 = 20
/
Note. By default the Migdal-Eliashberg equations are solved in the FSR approximation. The FBW treatment is activated in Step 8 below with fbw = .true..
EPW first interpolates the electron-phonon matrix elements onto the
dense 40x40x40 k-grid and 20x20x20 q-grid, restricts to q-points that
fall inside the Fermi window, writes the mgb2.ephmat/
directory, and then solves the anisotropic FSR Migdal-Eliashberg
equations at the temperatures listed in temps. The Pade
continuation to the real axis follows immediately. The iterative
continuation (lacon) is disabled for the anisotropic case
because it is expensive.
The Fermi surface files mgb2.fs_YY.cube (one per band in the
Fermi window) and mgb2.fs.frmsf are also written, since
fermi_plot = .true.. The .cube files can be opened
with VESTA and the
.frmsf file with FermiSurfer.
Step 4¶
Three coupling-related files are produced:
mgb2.lambda_pairs– anisotropic coupling lambda:sub:`n k, m k+q`(0) on the Fermi surface.mgb2.lambda_k_pairs– band- and wave-vector-resolved coupling lambda:sub:`n k`(0) on the Fermi surface.mgb2.a2f– isotropic alpha:sup:`2`F and cumulative lambda as a function of frequency for several phonon smearing values.
Note. Before plotting mgb2.a2f with gnuplot, comment out (with #) the first header line and the last seven lines that give smearing metadata.
gnuplot fig6.plt
evince fig6.pdf
Left: distribution of the anisotropic coupling lambdan k, m k+q`(0). Middle: lambda:sub:`n k`(0) on the Fermi surface. Right: Eliashberg spectral function and cumulative lambda from ``mgb2.a2f`.
Step 5¶
mgb2.imag_aniso_XX– imaginary-axis Matsubara frequency, Kohn-Sham eigenvalue relative to the Fermi level, renormalisation Zn k, gap Deltan k (eV). Written whenlaniso = limag = .true..mgb2.pade_aniso_XX– real-axis frequency, eigenvalue relative to the Fermi level, real/imaginary parts of Z, real/imaginary parts of Delta (eV). Written whenlpade = .true..
The file fig7.pdf can reach tens of megabytes, so convert it
to PNG before opening:
gnuplot fig7.plt
pdftoppm -png -r 100 fig7.pdf fig7
display fig7-1.png
Anisotropic superconducting gap of MgB2 at 10 K on the imaginary and real frequency axes.
Step 6¶
Collect the leading edge of the gap at every temperature from the
mgb2.imag_aniso_gap0_XX files and plot:
gnuplot fig8.plt
evince fig8.pdf
Anisotropic gap of MgB2 as a function of temperature. Left: present tutorial. Right: converged reference.
Step 7¶
EPW writes the superconducting quasiparticle density of states relative
to the normal-state density of states in files
mgb2.qdos_XX. Divide the second column of these files by the
Fermi-level density of states reported in epw1.out before
plotting.
gnuplot fig9.plt
evince fig9.pdf
Normalised superconducting quasiparticle density of states of MgB2 at 10 K.
Step 8¶
Drop the Fermi-window restriction by reusing the same electron-phonon
matrix elements with fbw = .true.. Reuse the
mgb2.ephmat/ directory, the spectral function mgb2.a2f,
and the *.fmt files from Step 3:
cd ../epw2-FBW
ln -s ../epw1-FSR/*fmt .
ln -s ../epw1-FSR/mgb2.ephmat .
ln -s ../epw1-FSR/mgb2.a2f .
mpirun -np 48 epw.x -nk 48 -in epw2.in > epw2.out
Only the flags that differ from ../epw1-FSR/epw1.in are shown:
epw2.in (differences only)
ep_coupling = .false.
elph = .false.
epwwrite = .false.
epwread = .true.
wannierize = .false.
ephwrite = .false.
fbw = .true.
In the FBW treatment the electronic states are no longer restricted to the vicinity of the Fermi surface: the Matsubara sums include a band eigenvalue dependence and an additional energy-shift function chi. The chemical potential is held fixed at the Fermi level in this step.
Note. The Eliashberg spectral function prefix.a2f is not strictly required; if it is absent, EPW will recompute it.
Extract the FBW temperature-dependent gap as in Step 5 and compare with the FSR result:
gnuplot fig11.plt
evince fig11.pdf
Anisotropic superconducting gap of MgB2 as a function of temperature from the FSR (dashed) and FBW (solid) calculations.
Note. Adding muchem = .true. lets the chemical potential adjust at each temperature. The files
mgb2.imag_aniso_gap0_XX.frmsf (one per temperature), written
when iverbosity = 2, can be opened with
FermiSurfer.
Exercise 3: Superconducting properties of bcc Nb by solving the anisotropic ME equations using sparse-ir sampling¶
Niobium has a large phonon bandwidth and a narrow spread of Matsubara frequencies of interest, which makes it a good candidate for sparse intermediate-representation (sparse-ir) sampling of the imaginary axis. In this exercise the anisotropic FBW Migdal-Eliashberg equations are solved twice: first with the semi-empirical Coulomb parameter muc, and then with the ab-initio Coulomb interaction computed on an outer-band grid.
The sparse-ir sampling is activated by gridsamp = 2; EPW reads
the precomputed basis functions and sampling points from a file pointed
at by filirobj. With sparse-ir sampling the number of
Matsubara points required for a given numerical accuracy is typically
much smaller than with uniform sampling, so runs that would be
infeasible on a dense Matsubara grid become tractable.
Move to the Exercise 3 directory:
cd ../../exercise3
Step 1¶
cd phonon
mpirun -np 24 pw.x -nk 6 -in scf.in > scf.out
mpirun -np 24 ph.x -nk 6 -in ph.in > ph.out
scf.in
&control
calculation = 'scf'
prefix = 'nb'
restart_mode = 'from_scratch'
pseudo_dir = '../../pseudo/'
outdir = './'
/
&system
ibrav = 3
celldm(1) = 6.252854867061436
nat = 1
ntyp = 1
nbnd = 24
ecutwfc = 40
occupations = 'smearing'
smearing = 'm-p'
degauss = 0.01
/
&electrons
diagonalization = 'cg'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-12
/
ATOMIC_SPECIES
Nb 92.91 Nb_ONCV_PBE-1.2.upf
ATOMIC_POSITIONS crystal
Nb 0.0000000000 0.0000000000 0.0000000000 0 0 0
K_POINTS automatic
12 12 12 0 0 0
ph.in
&inputph
prefix = 'nb'
fildvscf = 'dvscf'
outdir = './'
fildyn = 'Nb.dyn'
ldisp = .true.
nq1 = 4
nq2 = 4
nq3 = 4
tr2_ph = 1.0d-14
/
The q-point grid written to ph.out is:
Dynamical matrices for ( 4, 4, 4) uniform grid of q-points
( 8 q-points):
N xq(1) xq(2) xq(3)
1 0.000000000 0.000000000 0.000000000
2 0.000000000 -0.250000000 0.250000000
3 0.000000000 0.500000000 -0.500000000
4 -0.250000000 0.750000000 -0.500000000
5 0.000000000 0.000000000 0.500000000
6 0.000000000 0.750000000 -0.250000000
7 0.500000000 -0.500000000 0.500000000
8 0.000000000 0.000000000 -1.000000000
Step 2¶
python3 /path/to/EPW/bin/pp.py
Enter the prefix "nb" when prompted.
Step 3¶
Run the NSCF calculation on a 4x4x4 Gamma-centered grid and perform the EPW interpolation step.
cd ../epw1-lambda
mpirun -np 48 pw.x -nk 48 -in scf.in > scf.out
mpirun -np 48 pw.x -nk 48 -in nscf.in > nscf.out
mpirun -np 48 epw.x -nk 48 -in epw1.in > epw1.out
The 4x4x4 k-list for nscf.in is produced with:
/path/to/wannier90/utility/kmesh.pl 4 4 4
nscf.in
&control
calculation = 'nscf'
prefix = 'nb'
restart_mode = 'from_scratch'
pseudo_dir = '../../pseudo/'
outdir = './'
/
&system
ibrav = 3
celldm(1) = 6.252854867061436
nat = 1
ntyp = 1
nbnd = 24
ecutwfc = 40
occupations = 'smearing'
smearing = 'm-p'
degauss = 0.01
/
&electrons
diagonalization = 'cg'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-10
/
ATOMIC_SPECIES
Nb 92.91 Nb_ONCV_PBE-1.2.upf
ATOMIC_POSITIONS crystal
Nb 0.0000000000 0.0000000000 0.0000000000 0 0 0
K_POINTS crystal
64
0.00000000 0.00000000 0.00000000 1.562500e-02
0.00000000 0.00000000 0.25000000 1.562500e-02
0.00000000 0.00000000 0.50000000 1.562500e-02
...
epw1.in
&inputepw
prefix = 'nb'
outdir = './'
dvscf_dir = '../phonon/save'
ep_coupling = .true.
elph = .true.
epwwrite = .true.
epwread = .false.
etf_mem = 0
vme = 'dipole'
nbndsub = 9
bands_skipped = 'exclude_bands = 1:4'
wannierize = .true.
num_iter = 0
dis_win_min = 10
dis_froz_min = 11
dis_froz_max = 28
proj(1) = 'Nb:s'
proj(2) = 'Nb:p'
proj(3) = 'Nb:d'
iverbosity = 2
fsthick = 0.5
degaussw = 0.1
degaussq = 0.5
ephwrite = .true.
eliashberg = .true.
nk1 = 4
nk2 = 4
nk3 = 4
nq1 = 4
nq2 = 4
nq3 = 4
mp_mesh_k = .true.
nkf1 = 20
nkf2 = 20
nkf3 = 20
nqf1 = 20
nqf2 = 20
nqf3 = 20
/
After the interpolation step, EPW writes nb.ephmat/,
nb.a2f, the coupling pair files, and the *.fmt files
that will be reused in later steps.
Step 4¶
Plot the isotropic and band-resolved anisotropic coupling strengths and the Eliashberg spectral function:
gnuplot fig14.plt
evince fig14.pdf
Left: distribution of lambdan k, m k+q`(0) on the Fermi surface of Nb. Middle: lambda:sub:`n k`(0). Right: Eliashberg spectral function and cumulative lambda from ``nb.a2f`.
Step 5¶
Reuse the nb.ephmat/ directory and the *.fmt files
from Step 3 and run the FBW anisotropic solver with sparse-ir Matsubara
sampling.
cd ../epw2-mustar
ln -s ../epw1-lambda/*fmt .
ln -s ../epw1-lambda/nb.ephmat .
ln -s ../epw1-lambda/nb.a2f .
mpirun -np 96 epw.x -nk 96 -in epw2.in > epw2.out
Only the flags that differ from ../epw1-lambda/epw1.in are
shown:
epw2.in (differences only)
ep_coupling = .false.
elph = .false.
epwwrite = .false.
epwread = .true.
wannierize = .false.
iverbosity = 0 ! set to 2 to request .frmsf files
ephwrite = .false.
laniso = .true. ! anisotropic Migdal-Eliashberg
limag = .true. ! imaginary-axis equations
fbw = .true. ! full-bandwidth formulation
muchem = .true. ! update chemical potential
gridsamp = 2 ! sparse-ir Matsubara sampling
filirobj = '../../irobjs/ir_nlambda6_ndigit8.dat'
broyden_beta = -0.7 ! negative -> linear mixing
conv_thr_iaxis = 1.0d-4
nsiter = 500
wscut = 1.0
temps = 0.2 6.0 10.0 12.0 13.0 13.5
muc = 0.24
Note 1. With gridsamp = 2, the Matsubara frequencies are read from the sparse-ir object file, so the value of wscut is used only to set the upper limit of the real frequency axis during analytic continuation.
Note 2. EPW cannot generate the sparse-ir object file itself. Produce it beforehand with the Python sparse-ir package and point filirobj at it. The file name pattern
ir_nlambda<N>_ndigit<M>.dat encodes the two parameters
Lambda = 10N (maximum Matsubara frequency) and
epsilonIR = 10-M (basis accuracy).
Larger N or M yield more sampling points.
Note 3. From EPW v5.8 onwards, ephmatXX files can be reread with a different number of MPI pools than used when they were written. Older ephmat directories are not compatible with the current reader.
Note 4. For anisotropic FBW runs, linear mixing (broyden_beta < 0) can converge more reliably than Broyden mixing.
During the run EPW prints the actual number of Matsubara sampling points, which should be much smaller than the uniform-grid number at the same cutoff:
===================================================================
Solve full-bandwidth anisotropic Eliashberg equations
===================================================================
Start reading ir object file
Finish reading ir object file
Actual number of frequency points ( 1) = 48 for sparse-ir sampling
temp( 1) = 0.20000 K
Total number of frequency points nsiw( 1) = 48
Parameters for IR basis: Lambda = 1.00E+06, eps_IR = 1.00E-08
Maximum frequency = 131.8778 eV
iter ethr znormi deltai [meV] shifti [meV] mu [eV]
1 1.129460E+00 2.444918E+00 1.657449E+00 1.255267E+01 1.786864E+01
...
Convergence was reached in nsiter = 23
Compare this with the dense-grid behaviour reported in Exercise 2: at a comparable accuracy target, sparse-ir sampling reaches convergence with tens of frequencies, while the uniform-grid FBW run in Exercise 2 needs hundreds.
Step 6¶
gnuplot fig15.plt
evince fig15.pdf
Anisotropic superconducting gap of Nb on the Fermi surface as a function of temperature, from the FBW solver with sparse-ir sampling.
Step 7¶
To switch from the semi-empirical Coulomb pseudopotential to the
ab-initio Coulomb interaction, run an NSCF calculation on an outer-band
window and extract the eigenenergies with
nscf2supercond.x. The k-grid must be commensurate with the
nkf1-3 used in Step 3.
cd ../outerbands
mpirun -np 24 pw.x -nk 6 -in scf.in > scf.out
mpirun -np 24 pw.x -nk 6 -in nscf_outerbands.in > nscf_outerbands.out
/path/to/EPW/bin/nscf2supercond.x < nscf2supercond.in
nscf_outerbands.in
&control
calculation = 'bands'
prefix = 'nb'
restart_mode = 'from_scratch'
pseudo_dir = '../../pseudo/'
outdir = './'
verbosity = 'high'
/
&system
ibrav = 3
celldm(1) = 6.252854867061436
nat = 1
ntyp = 1
nbnd = 24
ecutwfc = 40
occupations = 'smearing'
smearing = 'm-p'
degauss = 0.01
/
&electrons
diagonalization = 'cg'
mixing_mode = 'plain'
mixing_beta = 0.7
conv_thr = 1.0d-10
/
ATOMIC_SPECIES
Nb 92.91 Nb_ONCV_PBE-1.2.upf
ATOMIC_POSITIONS crystal
Nb 0.0000000000 0.0000000000 0.0000000000 0 0 0
K_POINTS automatic
20 20 20 0 0 0
nscf2supercond.in
&bands
outdir = './'
prefix = 'nb'
filband = 'nb.bands.20x20x20.dat'
/
Step 8¶
Rerun the sparse-ir FBW solver with the Coulomb term accounting for outer bands. Reuse the interpolated matrix elements from Step 3:
cd ../epw3-outerbands
ln -s ../epw1-lambda/*fmt .
ln -s ../epw1-lambda/nb.ephmat .
ln -s ../epw1-lambda/nb.a2f .
mpirun -np 96 epw.x -nk 96 -in epw3.in > epw3.out
Only the flags that differ from ../epw2-mustar/epw2.in are
shown:
epw3.in (differences only)
broyden_beta = -0.5 ! linear mixing, smaller factor
icoulomb = 1 ! include the outside-fsthick states
filnscf_coul = '../outerbands/nb.bands.20x20x20.dat'
emax_coulomb = 15.0d0 ! outer window upper limit [eV]
emin_coulomb = -15.0d0 ! outer window lower limit [eV]
muc = 0.429 ! ab-initio Coulomb parameter
Note 1. icoulomb = 1 requires that the fine k-grid and the outer-band k-grid be commensurate. The 20x20x20 grid used here for nscf_outerbands.in matches nkf1-3 = 20 from Step 3; compatible alternatives are 10x10x10, 5x5x5, 4x4x4, etc.
Note 2. The value muc = 0.429 is the ab-initio Coulomb parameter for Nb taken from
Phys. Rev. B 101, 134511 (2020).
Note 3. Convergence with icoulomb = 1 is slower than the bare FBW run. A smaller mixing factor (broyden_beta = -0.5) often helps.
EPW reports which outer bands it selected from the NSCF file:
Start reading nscf file for Coulomb
Finish reading nscf file for Coulomb
k-grid read from ../outerbands/nb.bands.20x20x20.dat : 20 20 20
Nr irreducible k-points read from ../outerbands/nb.bands.20x20x20.dat : 256
Minimum eigenvalue of bands taken from the file (eV) = -3.6519037290E+01
Maximum eigenvalue of bands taken from the file (eV) = 6.5821232648E+01
emin_coulomb + "Fermi level" (eV) = 2.8752569141E+00
emax_coulomb + "Fermi level" (eV) = 3.2875256914E+01
Only states taken from nscf file between 2.875257 eV and 32.875257 eV
will be included in the Eliashberg calculations.
9 bands in the interval [emin_coulomb + "Fermi level", emax_coulomb + "Fermi level"]
Step 9¶
Plot the temperature-dependent gap from the ab-initio run and compare with the semi-empirical result of Step 5:
gnuplot fig16.plt
evince fig16.pdf
Anisotropic superconducting gap of Nb obtained with the semi-empirical muc = 0.24 (Step 5) and with the ab-initio muc = 0.429 plus outer-band Coulomb term (Step 8). Right panel: converged reference from Phys. Rev. B 110, 064505 (2024).
Appendix: Restart options¶
EPW supports three restart modes that cover the usual iteration patterns of a superconductivity workflow.
Restart from the Wannier-representation matrix elements (``prefix.epmatwp``, ``epwdata.fmt``). Useful when changing any of
nkf1-3,nqf1-3,fsthick, ordegaussw. Deleteprefix.ephmat,restart.fmt,selecq.fmt, andprefix.a2fbefore restarting. Required files:prefix.epmatwp,prefix.ukk,crystal.fmt,epwdata.fmt, andvmedata.fmt(ordmedata.fmt).ep_coupling = .true. elph = .true. epwwrite = .false. epwread = .true. wannierize = .false. ephwrite = .true.
The number of MPI pools may differ from the original run.
Restart an interrupted q-loop during ``ephmatXX`` writing. Required files:
prefix.epmatwp,prefix.ukk,crystal.fmt,epwdata.fmt,vmedata.fmt(ordmedata.fmt),restart.fmt, andselecq.fmt(selecq.fmtonly ifselecqread = .true.; otherwise it is regenerated).ep_coupling = .true. elph = .true. epwwrite = .false. epwread = .true. wannierize = .false. ephwrite = .true.
The number of MPI pools must match the original run.
Restart from the saved ``ephmatXX`` files. Required files: the
prefix.ephmat/directory (withegnv,freq,ikmap, and theephmatXXfiles),prefix.dos,selecq.fmt, andcrystal.fmt.ep_coupling = .false. elph = .false. epwwrite = .false. epwread = .true. wannierize = .false. ephwrite = .false.
From EPW v5.8 onwards,
ephmatXXfiles may be read back with a different number of MPI pools.