Interpolation of electron-phonon matrix elements¶

J. Lafuente-Bartolomé

Note

Hands-on based on QE-v7.5 and EPW-v6.0

In this tutorial we will learn to use the core capabilities of EPW. In particular, we will focus on how to check the quality of Wannier interpolation of physical quantities for a future production run.

Prerequisites¶

This tutorial assumes that Quantum ESPRESSO and EPW have already been downloaded and compiled. If you have not yet completed these steps, please follow the Installation and setup tutorial before continuing.

The examples below also assume that the environment variables used to run the executables (such as $QE, $RUN_PW, and $RUN_EPW) have been defined as described in the setup tutorial.

Note that variables defined with export only apply to the current terminal session. If you compiled the code previously but are starting a new session, you may need to redefine the run variables following the instructions in Defining run commands.

If the setup tutorial has been followed in the current session, the commands below should work directly in your terminal.

Exercise 1: Metallic Pb¶

In this exercise we will calculate physical properties of metallic Pb using Quantum ESPRESSO and repeat the same calculations through Wannier interpolation using EPW. We will then compare them to check the interpolation quality. After assessing interpolation quality, we will calculate the phonon linewidth for Pb.

Let us first move to the directory containing the input files for the coarse-grid phonon calculations:

cd phonons

Step 1: Run a self-consistent calculation on a homogeneous $8 \times 8 \times 8$ k-point grid and a phonon calculation on a homogeneous $3 \times 3 \times 3$ q-point grid using the following input files:

Note: The energy cutoff ecutwfc needed for convergence will be higher and should be checked.

scf.in

&CONTROL
  calculation='scf'
  prefix='pb',
  pseudo_dir = '../',
  outdir='./'
/
&SYSTEM
  ibrav=  2,
  celldm(1) = 9.2225583816,
  nat= 1,
  ntyp= 1,
  ecutwfc = 30.0
  occupations='smearing',
  smearing='marzari-vanderbilt',
  degauss=0.05
/
&ELECTRONS
  conv_thr =  1.0d-10
  mixing_beta = 0.7
/
ATOMIC_SPECIES
  Pb 207.2 pb_s.UPF
ATOMIC_POSITIONS
  Pb 0.00 0.00 0.00
K_POINTS {automatic}
  8 8 8 0 0 0

ph.in

&INPUTPH
  prefix = 'pb',
  fildyn = 'pb.dyn.xml',
  fildvscf = 'dvscf',
  ldisp = .true.,
  nq1 = 3,
  nq2 = 3,
  nq3 = 3,
  tr2_ph = 1.0d-16
/

The calculations can be launched with the following commands:

$RUN_PW < scf.in > scf.out
$RUN_PH < ph.in > ph.out

The keyword fildvscf in ph.in tells the code to write on file the change of the self-consistent potential due to phonon perturbations, $\partial_{\mathbf{q}\nu}V^{\rm scf}$, that is needed to compute the electron-phonon matrix elements. In the output file ph.out, locate the list of irreducible q points:

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

The list of irreducible $\mathbf{q}$ points is also written in the pb.dyn0 file. If you type ls, you can see pb.dynX.xml file containing the dynamical matrix has been produced for each irreducible $\mathbf{q}$ point. The dvscf files are all named pb.dvscf1 and are located inside the _ph0/pb.q_X/ folders, except for the onecorresponding to the first $\mathbf{q}$ point ($\Gamma$) that is located in _ph0/.

Note: The dynamical matrices files have been written in .xml format, which is advisable, because we specified it in the ph.in input with fildyn = ’pb.dyn.xml’.

Step 2: Calculate the phonon dispersion by using Fourier interpolation.

First, bring the interatomic force constants to real space using a Fourier transformation with q2r.x:

q2r.in

&INPUT
  zasr='simple',
  fildyn='pb.dyn.xml',
  flfrc='pb333.fc'
/

$RUN_Q2R < q2r.in > q2r.out

In the output, check that the FFT was successful:

fft-check success (sum of imaginary terms < 10^-12)

Then, calculate the phonon dispersion along the high-symmetry lines using matdyn.x:

matdyn.in

&INPUT
  asr              = 'simple',
  flfrc            = 'pb333.fc.xml',
  flfrq            = 'pb.freq',
  q_in_band_form   = .true.,
  q_in_cryst_coord = .true.
/
  6
  0.000 0.000 0.000  30 ! \Gamma
  0.500 0.000 0.500  30 ! X
  0.500 0.250 0.750  30 ! W
  0.500 0.500 0.500  30 ! L
  0.000 0.000 0.000  30 ! \Gamma
  0.375 0.375 0.750  30 ! K

$RUN_MATDYN < matdyn.in > matdyn.out

This produces the file named pb.freq.gp with the phonon frequencies along the path, expressed in cm$^{-1}$ units. This will be checked against the phonon frequencies along the same path from EPW later.

Step 3: Gather the .dyn and .dvscf files into a new save/ directory which EPW will read. The files in _ph0/pb.phsave/ containing the displacement patterns are also needed. This can be done using the pp.py python script which is included in the EPW/bin distribution:

python3 $QE/../EPW/bin/pp.py

The script will ask you to enter the prefix of your calculation (in this case, pb):

Enter the prefix used for PH calculations (e.g. diam)
pb

Note: If numpy is not already installed in your Python environment, the script may fail. In that case, install it with python3 -m pip install numpy and then run the script again.

Step 4: Run a non self-consistent calculation for electronic band structures using the charge density and other outputs from the previous self-consistent run:

cd ../band/
cp -rp ../phonon/pb.save/ .

nscf.in

&CONTROL
  calculation='bands'
  prefix='pb',
  pseudo_dir = '../',
  outdir='./'
/
  &SYSTEM
  ibrav=  2,
  celldm(1) = 9.2225583816,
  nat= 1,
  ntyp= 1,
  ecutwfc = 30.0
  occupations='smearing',
  smearing='marzari-vanderbilt',
  degauss=0.05
  nbnd=10
/
&ELECTRONS
  conv_thr =  1.0d-10
  mixing_beta = 0.7
/
ATOMIC_SPECIES
  Pb 207.2 pb_s.UPF
ATOMIC_POSITIONS
  Pb 0.00 0.00 0.00
K_POINTS crystal_b
  6
  0.000 0.000 0.000  30 ! \Gamma
  0.500 0.000 0.500  30 ! X
  0.500 0.250 0.750  30 ! W
  0.500 0.500 0.500  30 ! L
  0.000 0.000 0.000  30 ! \Gamma
  0.375 0.375 0.750  30 ! K

$RUN_PW < bands.in > bands.out

Then, run the program of bands.x to obtain the band structure data:

bands.in

&BANDS
  prefix = 'pb'
  lsym   = .false.
/

$RUN_BANDS < bandsx.in > bandsx.out

This will produce a file named bands.out.gnu, which will be later compared with the interpolated electronic band structure from EPW.

Step 5: Run a non self-consistent calculation on a homogeneous $6\times6\times6$ $\mathbf{k}$-point grid for generating the input of wave functions for EPW, and an EPW calculation for the Wannier interpolation of electronic band structure and phonon dispersion along the high-symmetry lines:

cd ../epw1
cp -rp ../phonon/pb.save/ .

nscf.in

&CONTROL
  calculation='bands'
  prefix='pb',
  pseudo_dir = '../',
  outdir='./'
/
&SYSTEM
  ibrav=  2,
  celldm(1) = 9.2225583816,
  nat= 1,
  ntyp= 1,
  ecutwfc = 30.0
  occupations='smearing',
  smearing='marzari-vanderbilt',
  degauss=0.05
  nbnd=10
/
&ELECTORNS
  conv_thr =  1.0d-10
  mixing_beta = 0.7
/
ATOMIC_SPECIES
  Pb 207.2 pb_s.UPF
ATOMIC_POSITIONS
  Pb 0.00 0.00 0.00
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      = 'pb',
  outdir      = './'
  dvscf_dir   = '../phonon/save'

  etf_mem     = 0

  elph        = .true.
  epbwrite    = .true.
  epbread     = .false
  epwwrite    = .true.
  epwread     = .false.

  wannierize  = .true.
  nbndsub     =  4
  bands_skipped = 'exclude_bands = 1-5'
  num_iter    = 300
  dis_win_max = 21
  dis_froz_max= 13.5
  proj(1)     = 'Pb:sp3'
  wannier_plot= .true.
  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 G 0.00 0.00 0.00'
  wdata(7)    = 'G 0.00 0.00 0.00 K 0.375 0.375 0.75'
  wdata(8)    = 'end kpoint_path'
  wdata(9)    = 'bands_num_points = 10'
  wdata(10)   = 'dis_num_iter      = 200'
  wdata(11)   = 'num_print_cycles  = 10'
  wdata(12)   = 'conv_tol = 1E-12'
  wdata(13)   = 'conv_window = 4'

  nkf1         = 1
  nkf2         = 1
  nkf3         = 1
  nqf1         = 1
  nqf2         = 1
  nqf3         = 1

  nk1         = 6
  nk2         = 6
  nk3         = 6
  nq1         = 3
  nq2         = 3
  nq3         = 3
/

Note: The k and q grids need to be commensurate, with the k grid at least of the size of the q grid. Since we chose a $6\times6 \times6$ k-point grid, the $3\times 3 \times3$ q-point grid used in the phonon calculation is appropriate; however a $6\times6 \times6$ q-point grid would be needed in order to interpolate the dynamical matrix and the electron-phonon matrix elements more accurately. For computational efficiency —see Phys. Rev. Research 3, 043022 (2021)— we recommend using the same k-point and q grids.

Note: In dvscf_dir = ’../phonon/save’ we specify the directory where the .dyn, .dvscf and patterns files are stored.

$RUN_PW < nscf.in > nscf.out
$RUN_EPW < epw1.in > epw1.out

In the first run, a nscf``calculation will be performed in the full coarse **k**-grid. Then, ``EPW will perform the following calculations:

Obtain the Maximally Localized Wannier Functions using Wannier90 as a library. You can find the Wannier function centers and spreads obtained in the epw1.out output file:
```
Wannier Function centers (cartesian, alat) and spreads (ang):

(   0.07779   0.07779   0.07779) :   2.22274
(   0.07779  -0.07779  -0.07779) :   2.22274
(  -0.07779   0.07779  -0.07779) :   2.22274
(  -0.07779  -0.07779   0.07779) :   2.22274
```
The full output from the Wannier90 run is in the pb.wout file. The input keywords for the Wannierization are in the block following wannierize = .true.; nbndsub corresponds to the number of Wannier functions (4, starting from Pb sp$^3$ orbitals as the initial guess), while bands_skipped = ... specifies the list of bands not wannierized (generally a set of bands lying at lower energies, such as semicore states in this example). For the other input keywords you can refer to the documentation page at https://docs.epw-code.org/doc/Inputs.html. It is also possible to add extra keywords that are read by Wannier90 using the input keyword wdata(index) with increasing index number. In this example, we are asking Wannier90 to plot the electronic band structure and specify the high-symmetry points that define the path. As a result, the code will produce the pb_band.dat, pb_band.gnu, and pb_band.kpt files.
With the wannier plot= .true. input keyword, we generate directly in EPW the cube files containing real-space Wannier functions, which can be plotted using visualization software such as VMD and VESTA; the plot of pb_00001.cube Wannier function in this example should look like this:

Compute the electron-phonon matrix elements on the (k, k + q) points for each irreducible q-point in the Brillouin zone, and unfold to the full Brillouin zone using symmetries. This is the most expensive part of the run. In the output you should see:

===================================================================
irreducible q point #    1
===================================================================

Symmetries of small group of q: 48
      in addition sym. q -> -q+G:

Number of q in the star =    1
List of q in the star:
      1   0.000000000   0.000000000   0.000000000
Imposing acoustic sum rule on the dynamical matrix

   q(    1 ) = (   0.0000000   0.0000000   0.0000000 )

...

Transform all quantities from reciprocal (Bloch) space to real (Wannier) space, and store the resulting matrices and additional information needed for restarting a calculation on disk (pb.epmatwp, crystal.fmt, vmedata.fmt, wigner.fmt, and epwdata.fmt files):
```
Writing Hamiltonian, Dynamical matrix and EP vertex in Wann rep to file
```

The Wannier-Fourier interpolation technique relies on the real-space decay of the Wannier functions and of the phonon perturbation. To examine how these quantities decay within the supercell defined by the coarse k- and q- grids, plot the files decay.H (Hamiltonian), decay.dynmat (dynamical matrix), and decay.epmate / decay.epmatp (electron–phonon matrix elements). For example, the following plots can be generated using gnuplot and the script below (note the logarithmic scale on the $y$ axis).

plot_decay.gnu

set encoding iso_8859_1
set logscale y
set format y "10^{%L}"
set pointsize 2

set xlabel "|R_e| (\305)"
set ylabel "max H_{nm} (Ry)"
plot "decay.H" u 1:2 w p pt 7
pause -1

set xlabel "|R_p| (\305)"
set ylabel "max D_{nm} (Ry)"
plot "decay.dynmat" u 1:2 w p pt 7
pause -1

set xlabel "|R_e| (\305)"
set ylabel "max g_{nm} (Ry)"
plot "decay.epmate" u 1:2 w p pt 7
pause -1

set xlabel "|R_p| (\305)"
set ylabel "max g_{nm} (Ry)"
plot "decay.epmatp" u 1:2 w p pt 7
pause -1

Run the script with gnuplot:

gnuplot plot_decay.gnu

These plots show that a larger supercell (i.e., a denser q-grid) is required to interpolate the phonon properties accurately, whereas the electronic contribution decays well. Recall that we are using a $6\times6\times6$ k-grid and only a $3\times3\times3$ q-grid.

Step 6: Restart from the real-space quantities obtained in the previous calculations and verify that the interpolated electronic and phonon band structures reproduce the direct DFT/DFPT results computed earlier.

Note: You should also check the quality of your electron-phonon matrix element interpolation by comparing with direct DFPT calculations. We will perform such check in Exercise 2.

Edit the file pb_band.kpt produced by Wannier90 so that the $\mathbf{k}$-points are specified in crystal units.

pb_band.kpt

crystal
000000    0.000000    0.000000   1.0
000000    0.050000    0.050000   1.0
000000    0.100000    0.100000   1.0
000000    0.150000    0.150000   1.0
000000    0.200000    0.200000   1.0
000000    0.250000    0.250000   1.0
000000    0.300000    0.300000   1.0
000000    0.350000    0.350000   1.0
...

Then move to the second folder, copy the necessary files for restart, and prepare the input files:

cd ../epw2
ln -s ../epw1/pb.epmatwp .
ln -s ../epw1/pb.ukk .
ln -s ../epw1/*.fmt .
ln -s ../epw1/pb_band.kpt .

epw2.in

&INPUTEPW
  prefix      = 'pb',
  outdir      = './'
  dvscf_dir   = '../phonon/save'

  etf_mem     = 0

  elph        = .true.
  epwread     = .true.

  wannierize  = .false.
  nbndsub     =  4
  bands_skipped = 'exclude_bands = 1-5'

  band_plot   = .true.
  fsthick     = 100    ! eV - should be large here for bandstructure plot

  filkf       = 'pb_band.kpt'
  filqf       = 'pb_band.kpt'

  nk1         = 6
  nk2         = 6
  nk3         = 6
  nq1         = 3
  nq2         = 3
  nq3         = 3
/

The flag band_plot = .true. in epw2.in instructs EPW to save the interpolated electronic bands (band.eig) and phonon frequencies (phband.freq) to file. In this example, the interpolation is performed along the same Brillouin-zone path used in the previous pw.x and matdyn.x calculations. The interpolation path is read from the file pb_band.kpt produced by Wannier90 in the previous step and modified above. This file is specified through the keywords filkf and filqf.

Run the EPW calculation:

$RUN_EPW < epw2.in > epw2.out

To extract the data for plotting, run the plotband.x utility from QuantumESPRESSO. Provide the input file (band.eig or phband.freq), the energy range (for example, -1 20), and the output file containing the data to plot (band.dat or freq.dat). The remaining prompts are not relevant; simply press ENTER when asked:

$QE/plotband.x
     Input file > band.eig
Reading    4 bands at     43 k-points
Range:   -0.4085   19.0142eV  Emin, Emax, [firstk, lastk] > -1 20
high-symmetry point:  0.0000 0.0000 0.0000   x coordinate   0.0000
high-symmetry point:  0.0000 1.0000 0.0000   x coordinate   1.0000
high-symmetry point: -0.5000 1.0000 0.0000   x coordinate   1.5000
high-symmetry point: -0.5000 0.5000 0.5000   x coordinate   2.2071
high-symmetry point:  0.0000 0.0000 0.0000   x coordinate   3.0731
high-symmetry point: -0.7500 0.7500 0.0000   x coordinate   4.1338
output file (gnuplot/xmgr) > band.dat
bands in gnuplot/xmgr format written to file band.dat
output file (ps) >
stopping ...

$QE/plotband.x
     Input file > phband.freq
Reading    3 bands at     43 k-points
Range:    0.0000   13.9958eV  Emin, Emax, [firstk, lastk] > 0 14
high-symmetry point:  0.0000 0.0000 0.0000   x coordinate   0.0000
high-symmetry point:  0.0000 1.0000 0.0000   x coordinate   1.0000
high-symmetry point: -0.5000 1.0000 0.0000   x coordinate   1.5000
high-symmetry point: -0.5000 0.5000 0.5000   x coordinate   2.2071
high-symmetry point:  0.0000 0.0000 0.0000   x coordinate   3.0731
high-symmetry point: -0.7500 0.7500 0.0000   x coordinate   4.1338
output file (gnuplot/xmgr) > freq.dat
bands in gnuplot/xmgr format written to file freq.dat
output file (ps) >
stopping ...

The electronic band structure obtained with pw.x can now be compared with the Wannier-interpolated bands produced by EPW. This can be done using gnuplot and the script below.

plot_bands.gnu

set terminal x11 enhanced
set encoding utf8
set ylabel "Energy (Ry)"
set xtics ("{/Symbol G}" 0, "X" 1, "W" 1.5, "L" 2.2071, "{/Symbol G}" 3.0731, "K" 4.1338)
set arrow from 1, graph 0 to 1, graph 1 nohead
set arrow from 1.5, graph 0 to 1.5, graph 1 nohead
set arrow from 2.2071, graph 0 to 2.2071, graph 1 nohead
set arrow from 3.0731, graph 0 to 3.0731, graph 1 nohead
plot "../band/bands.out.gnu" u 1:2 w l lw 3 title "pw.x", \
     "band.dat" u 1:2 w l lw 3 dt 2 title "EPW"
pause -1

Run the script with gnuplot:

gnuplot plot_bands.gnu

The phonon band structures can be compared in the same way using gnuplot and the script below.

plot_phbands.gnu

set terminal x11 enhanced
set encoding utf8
set ylabel "{/Symbol w} (meV)"
set xtics ("{/Symbol G}" 0, "X" 1, "W" 1.5, "L" 2.2071, "{/Symbol G}" 3.0731, "K" 4.1338)
set arrow from 1, graph 0 to 1, graph 1 nohead
set arrow from 1.5, graph 0 to 1.5, graph 1 nohead
set arrow from 2.2071, graph 0 to 2.2071, graph 1 nohead
set arrow from 3.0731, graph 0 to 3.0731, graph 1 nohead
plot "../phonon/pb.freq.gp" u 1:($2/8.06554) w l lw 3 lt 1 title "matdyn.x", \
     "" u 1:($3/8.06554) w l lw 3 lt 1 notitle, \
     "" u 1:($4/8.06554) w l lw 3 lt 1 notitle, \
     "freq.dat" u 1:2 w l lw 3 dt 2 title "EPW"
pause -1

Run the script with gnuplot:

gnuplot plot_phbands.gnu

Step 7: Compute the phonon linewidths $\gamma_{\mathbf{q}\nu}$ and the electron-phonon coupling strength $\lambda_{\mathbf{q}\nu}$ along the high-symmetry path.

This is done within the double-delta approximation by setting phonselfen = .true. and delta_approx = .true. in the input file.

epw3.in

&INPUTEPW
  prefix      = 'pb',
  outdir      = './'
  dvscf_dir   = '../phonon/save'

  etf_mem     = 0

  elph        = .true.
  epwread     = .true.

  wannierize  = .false.
  nbndsub     =  4
  bands_skipped = 'exclude_bands = 1-5'

  phonselfen  = .true.
  delta_approx= .true.

  fsthick     = 1     ! eV
  temps       = 0.075 ! K
  degaussw    = 0.2   ! eV

  filqf       = 'pb_band.kpt'
  nkf1         = 20
  nkf2         = 20
  nkf3         = 20

  nk1         = 6
  nk2         = 6
  nk3         = 6
  nq1         = 3
  nq2         = 3
  nq3         = 3
/

Note: The parameter fsthick determines the energy window around the Fermi level for which the electron-phonon matrix elements are interpolated. This can reduce significantly the cost of calculations; for example, only electronic states within a few phonon energies from the Fermi level will contribute to the phonon linewidth, therefore we can use fsthick = 1 (eV) in this example.

Note: Since calculating $\gamma_{\mathbf{q} \nu}$ and $\lambda_{\mathbf{q} \nu}$ requires an integration over k (see expressions below), we interpolate the electron-phonon matrix elements onto a denser $20\times20\times20$ k-point grid.

The phonon linewidth (i.e. imaginary part of the phonon self-energy) is defined as:

\[\gamma_{\mathbf{q} \nu} = 2\pi \omega_{\mathbf{q} \nu} \sum_{nm} \int_{\rm BZ} \frac{d\mathbf{k}}{\Omega_{\rm BZ}} |g_{nm\nu}(\mathbf{k,q})|^2 \delta(\varepsilon_{n\mathbf{k}}-\varepsilon_{\rm F})\delta(\varepsilon_{n\mathbf{k+q}}-\varepsilon_{\rm F}),\]

where $\varepsilon_{\rm F}$ is the Fermi level. The mode-resolved electron-phonon coupling strength is defined as:

\[\lambda_{\mathbf{q} \nu} = \frac{1}{N(\varepsilon_{\rm F})\omega_{{\mathbf{q}}\nu}} \sum_{nm} \int_{\rm BZ} \frac{d\mathbf{k}}{\Omega_{\rm BZ}} |g_{nm\nu}(\mathbf{k,q})|^2 \delta(\varepsilon_{n\mathbf{k}}-\varepsilon_{\rm F})\delta(\varepsilon_{n\mathbf{k+q}}-\varepsilon_{\rm F}) = \frac{\gamma_{\mathbf{q}\nu}}{2 \pi N(\varepsilon_{\rm F})\omega_{\mathbf{q} \nu}^2},\]

where $N(\varepsilon_{\rm F})$ is the density of states per spin at the Fermi level.

Run the EPW calculation:

$RUN_EPW < epw3.in > epw3.out

You can see that now the phonon linewidths and coupling strengths are printed in output for each phonon wavevector $\mathbf{q}$ and mode $\nu$ (lambda___ and gamma___). The sum of $\lambda_{\mathbf{q}\nu}$ over all phonon modes, $\lambda_\mathbf{q}$, is also written (lambda___( tot )). $\gamma_{\mathbf{q}\nu}$ and $\lambda_{\mathbf{q}\nu}$ are stored in the files linewidth.phself.XXXK and lambda.phself.XXXK, respectively. Inspect those files to familiarize yourself with the format and learn how to plot these quantities. For example, to plot the $\mathbf{q}$-dependent linewidth of the third phonon, you can type in gnuplot:

gnuplot
gnuplot> plot "linewidth.phself.0.075K" u 1:4 every 3::2 w l lw 2

and you should be able to produce a plot like this:

Exercise 2: Semiconducting 3C-SiC¶

In this exercise we examine the electron-phonon interaction in semiconducting 3C-SiC, a polar material in which electrons couple strongly to longitudinal optical phonons. In polar materials the electron-phonon matrix elements exhibit a long-range Fröhlich divergence proportional to $1/|\mathbf{q}|$ in the limit $|\mathbf{q}| \rightarrow 0$ (PRL 115, 176401 (2015)), together with an additional quadrupolar contribution in the long-wavelength limit (PRR 3, 043022 (2021)).

We begin by checking the quality of the Wannier interpolation, as in Exercise 1. We then show how to correctly interpolate the electron-phonon matrix elements in the case of polar materials, and compute the electron linewidth.

Step 1: Run a self-consistent calculation on a homogeneous $8\times8\times8$ k-point grid and a phonon calculation on a homogeneous $3\times3\times3$ q-point grid

Move to the exercise2 directory, and prepare the input files:

cd ../../exercise2
cd phonon

scf.in

&CONTROL
  calculation     = 'scf'
  prefix          = 'sic'
  pseudo_dir      = '../'
  outdir          = './'
/
&SYSTEM
  ibrav           = 2
  celldm(1)       = 8.237
  nat             = 2
  ntyp            = 2
  ecutwfc         = 30.0
/
&ELECTRONS
  diagonalization = 'david'
  mixing_beta     = 0.7
  conv_thr        = 1.0d-10
/
ATOMIC_SPECIES
  Si  28.0855     Si.pz-vbc.UPF
  C   12.01078    C.UPF
ATOMIC_POSITIONS alat
  Si  0.00  0.00  0.00
  C   0.25  0.25  0.25
K_POINTS automatic
  8 8 8 0 0 0

ph.in

&INPUTPH
  prefix   = 'sic'
  fildvscf = 'dvscf'
  ldisp    = .true
  fildyn   = 'sic.dyn.xml'
  nq1=3,
  nq2=3,
  nq3=3,
  tr2_ph   =  1.0d-16
/

Run the scf and ph calculations:

$RUN_PW < scf.in > scf.out
$RUN_PH < ph.in > ph.out

Note that the output of ph.out now contains also the Born effective charges $Z^\ast$ and the high-frequency dielectric constant $\epsilon_\infty$:

Electric Fields Calculation

...

End of electric fields calculation

     Dielectric constant in cartesian axis

     (       7.214365306       0.000000000       0.000000000 )
...

     Effective charges (d Force / dE) in cartesian axis

      atom      1   Si
 Ex  (        2.67023        0.00000       -0.00000 )
 ...

These quantities are automatically calculated for an insulating system and determine the frequency splitting between LO and TO phonon modes.

Step 2: Compute and plot the phonon dispersion using q2r.x and matdyn.x.

Prepare the input files:

q2r.in

&INPUT
  zasr='simple',
  fildyn='sic.dyn.xml',
  flfrc='sic333.fc'
/

matdyn.in

&INPUT
  asr              = 'simple',
  flfrc            = 'sic333.fc.xml'
  flfrq            = 'sic.freq'
  q_in_band_form   = .true.
  q_in_cryst_coord = .true.
/
  6
  0.000 0.000 0.000  30
  0.500 0.000 0.500  30
  0.500 0.250 0.750  30
  0.500 0.500 0.500  30
  0.000 0.000 0.000  30
  0.375 0.375 0.750  30

Run the calculations:

$RUN_Q2R < q2r.in > q2r.out
$RUN_MATDYN < matdyn.in > matdyn.out

And plot the phonon dispersion contained in the file sic.freq.gp gnuplot.

gnuplot
gnuplot> plot for [col=2:7] "sic.freq.gp" using 1:col w l lt 1 notitle

The result should look like this:

Note: Note that now there are three acoustic and three optical branches, separated by an energy gap, and that the so-called LO-TO splitting is present around $\Gamma$.

Gather the .dyn, .dvscf and patterns files into a save/ directory as in Exercise 1 using the pp.py script:

python3 $QE/../EPW/bin/pp.py

The script will ask you to enter the prefix of your calculation (in this case, sic):

Enter the prefix used for PH calculations (e.g. diam)
sic

Step 3: Compute the electron-phonon matrix elements along the selected high-symmetry lines using direct DFPT calculations.

Move to the ephline directory and copy the scf input file:

cd ../ephline
cp ../phonon/scf.in .

Prepare the ph input file:

ephline.in

&INPUTPH
  prefix   = 'sic'
  fildvscf = 'dvscf'
  ldisp    = .true.
  fildyn   = 'sic.dyn.xml'
  tr2_ph   =  1.0d-16
  qplot    = .true.
  q_in_band_form = .true.
  electron_phonon = 'prt'
  kx = 0.0
  ky = 0.0
  kz = 0.0
/
3
  1.0000 0.0000 0.0000 20 # X
  0.0000 0.0000 0.0000 20 # Gamma
 -0.5000 0.5000 0.5000 1  # L

Run the calculations:

$RUN_PW < scf.in > scf.out
$RUN_PH < ephline.in > ephline.out

The electron-phonon matrix elements are evaluated for $\mathbf{k}=\Gamma$ (the default) and for each $\mathbf{q}$ along the $X-\Gamma-L$ path. The results are written to the ephline.out file, with different lines corresponding to the different electronic bands and phonon modes. Because electron-phonon matrix elements are gauge dependent, the norm of the matrix elements is averaged over degenerate subspaces of bands and phonon modes. The resulting quantity is reported in the |g_sym| [meV] column.

Note: If you want to change the wavevector $\mathbf{k}$ of the initial state, you can do it by specifying the input flags of kx, ky, kz (in Cartesian $2\pi/a$ units).

To extract the electron-phonon matrix elements for the valence-band top ($n=m=4$) and selected phonon branches (for example $\nu=1,3,4,6$), filter the corresponding lines in ephline.out using grep (note the eight spaces between the indices):

grep "4        4        1" ephline.out > data1
grep "4        4        3" ephline.out > data3
grep "4        4        4" ephline.out > data4
grep "4        4        6" ephline.out > data6

Modes 2 and 5 are degenerate with modes 1 and 4, respectively, and are therefore not listed separately. These data will later be compared with the corresponding results obtained from EPW.

Step 4: Generate the wavefunctions on a uniform $\mathbf{k}$ grid and interpolate the electron-phonon matrix elements using EPW.

Move to the epw directory, and copy the .save/ directory from the previous scf run:

cd ../epw
cp -rp ../phonon/sic.save/ .

Prepare the input files:

nscf.in

&CONTROL
  calculation     = 'nscf'
  prefix          = 'sic'
  pseudo_dir      = '../'
  outdir          = './'
/
&SYSTEM
  ibrav           = 2
  celldm(1)       = 8.237
  nat             = 2
  ntyp            = 2
  ecutwfc         = 30.0
  nbnd            = 4
/
&ELECTRONS
  diagonalization = 'david'
  mixing_beta     = 0.7
  conv_thr        = 1.0d-10
/
ATOMIC_SPECIES
  Si  28.0855     Si.pz-vbc.UPF
  C   12.01078    C.UPF
ATOMIC_POSITIONS alat
  Si  0.00  0.00  0.00
  C   0.25  0.25  0.25
K_POINTS crystal
  216
  0.00000000  0.00000000  0.00000000  4.629630e-03
  0.00000000  0.00000000  0.16666667  4.629630e-03
  0.00000000  0.00000000  0.33333333  4.629630e-03
...

epw1.in

&INPUTEPW
  prefix      = 'sic'
  outdir      = './'
  dvscf_dir   = '../phonon/save'

  elph        = .true.
  epwwrite    = .true.
  lpolar      = .true.

  wannierize  = .true.
  nbndsub     =  4
  num_iter    = 300
  proj(1)     = 'Si:sp3'

  prtgkk      = .true.

  filqf       = 'path1.dat'
  nkf1        = 1
  nkf2        = 1
  nkf3        = 1

  nk1         = 6
  nk2         = 6
  nk3         = 6
  nq1         = 3
  nq2         = 3
  nq3         = 3
/

Note: To print the interpolated electron-phonon matrix elements to the output file, we set prtgkk = .true. in the EPW input file. We chose to print the electron-phonon matrix elements for the initial electronic states only at $\mathbf{k}$ = $\Gamma$ (nkf1=nkf2=nkf3=1).

And run the calculations:

$RUN_PW < nscf.in > nscf.out
$RUN_EPW < epw1.in > epw1.out

In this calculation we consider a path along the $X-\Gamma-L$ high-symmetry points in the BZ. This path is provided in the path1.dat file, which follows the same sequence of points as in ../ephline but with a denser sampling:

path1.dat

cartesian
0000000000  0.000000000  0.000000000  1.0
9800000000  0.000000000  0.000000000  1.0
9600000000  0.000000000  0.000000000  1.0
9400000000  0.000000000  0.000000000  1.0
9200000000  0.000000000  0.000000000  1.0
9000000000  0.000000000  0.000000000  1.0
8800000000  0.000000000  0.000000000  1.0
8600000000  0.000000000  0.000000000  1.0
....

Note: This file can be generated directly using Wannier90 called in library mode from EPW or with the following python script:

generate_XGL_path.py

import numpy as np
for ii in np.linspace(1,0,51):
    print( '{:16.10f}  0.000000000  0.000000000  1.0'.format(ii))
for ii in np.linspace(0+0.5/50,0.5,50):
    print( '{:16.10f}  {:16.10f}  {:16.10f}  1.0'.format(-ii,ii,ii))

Since 3C-SiC is not a metal, we set lpolar = .true. in order to correctly treat the long-range interaction. The strategy consists in subtracting the long-range component $g^\mathcal{L}$ from the full matrix element $g$ before interpolation and adding it back after interpolation (see e.g. npj Comput Mater 9, 156 (2023) for more details).

Note: An analogous strategy is implemented in q2r and matdyn to correctly interpolate the dynamical matrix including the long-range interactions which result in the LO-TO splitting.

In EPW, two long-range contributions can be considered: dipoles and quadrupoles. By default, only the dipoles will be considered when lpolar = .true.. To include quadrupoles, you need to provide a quadrupole.fmt file. The file must have that name and be located in the folder in which you run the calculation. For this example, we have taken the quadrupole values from Table II of Ref. PRR 3, 043022 (2021), so the quadrupole.fmt file (3 lines per atom) is as follows:

quadrupole.fmt

atom   dir       Qxx         Qyy         Qzz         Qyz         Qxz         Qxy
   1      0.000000    0.000000    0.000000    6.870000    0.000000    0.000000
   2      0.000000    0.000000    0.000000    0.000000    6.870000    0.000000
   3      0.000000    0.000000    0.000000    0.000000    0.000000    6.870000
   1      0.000000    0.000000    0.000000   -2.440000    0.000000    0.000000
   2      0.000000    0.000000    0.000000    0.000000   -2.440000    0.000000
   3      0.000000    0.000000    0.000000    0.000000    0.000000   -2.440000

Note: The quadrupole tensor can be obtained by fitting either the perturbed potential or the direct electron-phonon matrix elements (see for example here and here). Alternatively, it can be computed using perturbation theory with the Abinit code (see https://docs.abinit.org/topics/longwave). The latter is usually the simplest approach.

You can verify that quadrupoles are correctly included in the calculation by looking at the epw1.out file, where you should find:

------------------------------------
Quadrupole tensor is correctly read:
------------------------------------
atom   dir        Qxx       Qyy      Qzz        Qyz       Qxz       Qxy
  1        x       0.00000   0.00000   0.00000   6.87000   0.00000   0.00000
  1        y       0.00000   0.00000   0.00000   0.00000   6.87000   0.00000
  1        z       0.00000   0.00000   0.00000   0.00000   0.00000   6.87000
  2        x       0.00000   0.00000   0.00000  -2.44000   0.00000   0.00000
  2        y       0.00000   0.00000   0.00000   0.00000  -2.44000   0.00000
  2        z       0.00000   0.00000   0.00000   0.00000   0.00000  -2.44000

The electron-phonon matrix elements are evaluated for $\mathbf{k}=\Gamma$ (the default) and for each $\mathbf{q}$ along the selected Brillouin-zone path. The results are written to the epw1.out file, with different lines corresponding to the various electronic bands and phonon modes. Because these matrix elements are gauge dependent, the norm is averaged over degenerate subspaces of bands and phonon modes. The resulting (gauge-invariant) quantity is reported in the |g_sym| [meV] column. To extract the matrix elements for the valence-band top ($n=m=4$) and the different phonon branches, as in Step 4, filter the corresponding lines in epw1.out using grep (note the eight spaces between the indices):

grep "4        4        1" epw1.out > epwdata1
grep "4        4        3" epw1.out > epwdata3
grep "4        4        4" epw1.out > epwdata4
grep "4        4        6" epw1.out > epwdata6

You can now compare the direct DFPT results with the interpolated results using gnuplot. Create the following script:

plot_epw_vs_dfpt.gnu

set terminal x11 enhanced
set encoding utf8

set ylabel "|g|_{avg} (meV)" font ",20"
set xtics ("X" 0, "{/Symbol G}" 50, "L" 100) font ",20"
set ytics (0,300,600,900) font ",20"
set yrange [0:900]

set arrow from 50, graph 0 to 50, graph 1 nohead

plot "epwdata1" u 7 w l lw 2 lc rgb "red" title "EPW", \
     "epwdata3" u 7 w l lw 2 lc rgb "red" notitle, \
     "epwdata4" u 7 w l lw 2 lc rgb "red" notitle, \
     "epwdata6" u 7 w l lw 2 lc rgb "red" notitle, \
     "../ephline/data1" u ($0*2.5):7 w p pt 5 ps 1.5 lc rgb "black" title "ph.x", \
     "../ephline/data3" u ($0*2.5):7 w p pt 5 ps 1.5 lc rgb "black" notitle, \
     "../ephline/data4" u ($0*2.5):7 w p pt 5 ps 1.5 lc rgb "black" notitle, \
     "../ephline/data6" u ($0*2.5):7 w p pt 5 ps 1.5 lc rgb "black" notitle

pause -1

Note: We multiply the DFPT x-axis by 2.5 simply because this is the ratio of number of points computed with EPW.

And run the script with gnuplot:

gnuplot plot_epw_vs_dfpt.gnu

The resulting plot should look similar to the following figure:

Note: The figure also shows the results obtained without the quadrupole contribution (blue line). You can reproduce this curve by removing the quadrupole.fmt file and re-running the epw1.in calculation.

The interpolated results are not yet fully converged. A denser coarse $\mathbf{q}$ grid is required to accurately interpolate the electron-phonon matrix elements.

You can also inspect the interpolation when the long-range polar contribution is neglected by running a calculation with lpolar = .false.. This test should be performed in a separate directory, or followed by a fresh EPW calculation with lpolar = .true. to restore the default settings.

Step 5: Compute the linewidths of the valence states along the high-symmetry path.

The linewidths correspond to twice the imaginary part of the electron self-energy $\Sigma_{n\mathbf{k}}^{''}$:

\[\begin{split}\begin{gathered} \Sigma_{n\mathbf{k}}^{''}(\omega,T) = \pi \sum_{m\nu} \int_{\rm BZ} \frac{d\mathbf{q}}{\Omega_{\rm BZ}} | g_{mn,\nu}(\mathbf{k},\mathbf{q})|^2 \Big\{ [n_{\mathbf{q} \nu}(T) + f_{m\mathbf{k+q}}(T)]\delta(\omega - (\varepsilon_{m\mathbf{k+q}-\varepsilon_{\rm F}}) + \omega_{\mathbf{q} \nu}) \\ + [n_{\mathbf{q} \nu}(T) + 1 - f_{m\mathbf{k+q}}(T)]\delta(\omega - (\varepsilon_{m\mathbf{k+q}-\varepsilon_{\rm F}}) - \omega_{\mathbf{q} \nu}).\end{gathered}\end{split}\]

To perform this calculation, we set elecselfen = .true. in the epw input file. The calculation requires the $\mathbf{k}$-point path and a homogeneous $\mathbf{q}$ grid for the Brillouin-zone integration.

Prepare the input file:

epw2.in

&INPUTEPW
  prefix      = 'sic'
  amass(1)    = 28.0855
  amass(2)    = 12.0107
  outdir      = './'
  dvscf_dir   = '../phonon/save'

  etf_mem     = 0

  elph        = .true.
  epwwrite    = .false.
  epwread     = .true.
  lpolar      = .true.
  vme         = 'dipole'

  wannierize  = .false.
  nbndsub     =  4
  num_iter    = 300
  proj(1)     = 'Si:sp3'

  elecselfen  = .true.
  efermi_read = .true.
  fermi_energy= 9.6

  fsthick     = 7.0
  temps       = 20
  degaussw    = 0.05

  filkf       = 'path2.dat'
  nqf1        = 20
  nqf2        = 20
  nqf3        = 20

  nk1         = 6
  nk2         = 6
  nk3         = 6
  nq1         = 3
  nq2         = 3
  nq3         = 3
/

Note: Because the calculation uses a $\mathbf{k}$-point path rather than a homogeneous grid, the Fermi energy obtained from the interpolated bands is not reliable. To avoid this issue, we specify the Fermi energy explicitly in the input file by setting efermi_read = .true. and fermi_energy = 9.6 (slightly above the valence-band maximum).

And run the calculation:

$RUN_EPW < epw2.in > epw2.out

In the output, the progress of the $\mathbf{q}$-point integration is reported before the electron self-energy is printed for each $\mathbf{k}$ point:

     Progression iq (fine) =         50/      8000
     Progression iq (fine) =        100/      8000
     ...
     ...
     Average over degenerate eigenstates is performed
     Temperature:   20.000K
     WARNING: only the eigenstates within the Fermi window are meaningful
     ik =       1 coord.:    0.0000000   0.0000000   0.0000000
     -------------------------------------------------------------------
     E(   2 )=  -0.2406 eV   Re[Sigma]=      96.661664 meV Im[Sigma]=       0.240285 meV     Z=       0.860491 lam=       0.162127
     E(   3 )=  -0.2406 eV   Re[Sigma]=      96.661664 meV Im[Sigma]=       0.240285 meV     Z=       0.860491 lam=       0.162127
     E(   4 )=  -0.2406 eV   Re[Sigma]=      96.661664 meV Im[Sigma]=       0.240285 meV     Z=       0.860491 lam=       0.162127
     -------------------------------------------------------------------
...

The electron energies are now reported relative to the Fermi level $E_{\rm F}$. Since the parameter fsthick is set to 7 eV, the self-energy of the lowest valence band is not evaluated.

The linewidths can be plotted using the linewidth.elself file. For example, the imaginary part of the electron self-energy $\mathrm{Im}\,\Sigma$ for the highest valence band can be visualized with the following gnuplot script:

plot_linewidth.gnu

set terminal x11 enhanced
set encoding utf8

set ylabel "Im {/Symbol S} (meV)" font ",20"
set xtics ("{/Symbol G}" 1, "X" 31, "W" 61, "L" 91, "K" 121, "{/Symbol G}" 151, "K" 181)

set arrow from 31, graph 0 to 31, graph 1 nohead
set arrow from 61, graph 0 to 61, graph 1 nohead
set arrow from 91, graph 0 to 91, graph 1 nohead
set arrow from 121, graph 0 to 121, graph 1 nohead
set arrow from 151, graph 0 to 151, graph 1 nohead

set yrange [0:140]

plot "linewidth.elself.20.000K" u 1:4 every 3::2 w lp lw 2

pause -1

Run the script with:

gnuplot plot_linewidth.gnu

The resulting plot should resemble the following figure:

Note: The electron lifetime is related to the imaginary part of the self-energy through $\tau_{n\mathbf{k}} = \frac{\hbar}{2\,\mathrm{Im}\,\Sigma_{n\mathbf{k}}}$.

Note: You can inspect the contribution of individual phonon modes by setting iverbosity = 3 in the input file. In this case, linewidth.elself contains the mode-resolved linewidths.

Note: Try increasing the coarse $\mathbf{q}$ grid or using a random set of $\mathbf{q}$ points. You will observe that the linewidths are not yet fully converged. In polar materials, Brillouin-zone integrals converge more slowly because the Fröhlich interaction introduces a long-range $1/|\mathbf{q}|$ divergence in the electron-phonon matrix elements.