# Electron-phonon interaction {#sec:dfpteph}

The electron-phonon (eph) interaction is an important mechanism that
affects a broad range of materials properties. For example, it is one
the main scattering mechanisms in charge and heat transport, and it
determines the lifetime of electrons and phonons. DFPT provides a way to
obtain scattering rates due to electron-phonon interactions given by
electron-phonon interaction (EPI) matrix,
$g_{\mathbf{k+q}m,\mathbf{k}n}^{\mathbf{q}\nu}$, which is the central
quantity that characterizes eph interactions. Once EPI matrix is known,
related physical properties can be calculated through various
formalisms. In this section, we demonstrate simulating eph interactions
in RESCU by looking at aluminum example and calculating electrical
consuctivity.

EPI matrix for phonons with frequency $\omega_{\mathbf{q}\nu}$ and
polarization vector $\eta_{\mathbf{q}\nu}$ and electronic bloch states
of $| \mathbf{k}+\mathbf{q}m \rangle$ and $| \mathbf{k}n \rangle$ is
given by [\[SS96\]](#SS96){.citation}:

[$$g_{\mathbf{k}+\mathbf{q}m,\mathbf{k}n}^{\mathbf{q}\nu} = 
\sum_{R,\mu}
\frac{\eta_{\mathbf{q}\nu}(R\mu)}{(M_{R} \omega_{\mathbf{q}\nu})^{1/2}} 
\langle \mathbf{k}+\mathbf{q}m|
\frac{\delta V_{eff}}{\delta R_{\mu}}
| \mathbf{k}n \rangle$$]{label="elphgmat"}

where $M_{R}$ is the nuclei mass. Therefore, following quantities are
required for computing EPI matrix:

-   Ground state wavefunctions
-   Phonon modes
-   Perturbed self-consistent potentials for atomic displacements
    perturbations

In RESCU these quantities are obtained using the following workflow:

-   Performing a ground state calculation to obtain the ground state
    density and unperturbed wavefunctions
-   Performing a phonon (and dielectric in case of non-metallic systems)
    DFPT calculation to obtain dynamical matrix (dielectric constant and
    Born effective charges) and Perturbed self-consistent potentials
-   Performing a DFPT electron phonon calculation to obtain EPI matrix
    and calculate other eph related quantities

:::: important
::: title
Important
:::

Electron-phonon calculations often require much denser k- and q-sampling
compared to ground state and phonon (dielectric) DFPT calculations.
RESCU uses a nonself-consistent calculation to compute quantities on
dense k-mesh, and Fourier interpolation to compute dynamical matrix and
perturbed self-consistent potentials on dense q-mesh, to speed up the
calculations at eph level [\[B20\]](#B20){.citation}.
::::

:::: warning
::: title
Warning
:::

k- and q-mesh used in this example are only for demonstration and does
note generate accurate physical results
::::

We begin by calculating the ground state of Al using the following input
file:

``` matlab
info.calculationType = 'self-consistent'
info.savepath = './results/al_real_scf'
element(1).species = 'Al'
element(1).path = './Al_AtomicData.mat' 
atom.element = [1]
atom.xyz = [0 0 0]
domain.latvec = 7.665*[0.5,0.5,0.0;0.0,0.5,0.5;0.5,0.0,0.5]
domain.lowres = 0.2
functional.libxc = true
functional.list = {'XC_LDA_X','XC_LDA_C_PW'}
kpoint.gridn = [3 3 3] 
kpoint.sampling = 'gauss'   % required for dfpt for metals
kpoint.sigma = 0.025        % required for dfpt for metals
mixing.tol = 1e-8*[1 1]     % advised for dfpt
option.saveWavefunction = 1 % required for dfpt phonon (and dielectric)
diffop.method = 'fft'       % required for dfpt
domain.fourierInit = false  % required for dfpt
```

RESCU can be executed as

``` shell
rescu -i al_real_scf
```

Next, we perform a DFPT phonon calculation as follows. Note that this
time we also ask RESCU to save perturbed self-consistent potentials
which are to be used later in electron-phonon step:

``` matlab
info.calculationType = 'dfpt-phonon'
info.savepath = './results/al_real_phonon'
rho.in = './results/al_real_scf.mat'
psi.in = './results/al_real_scf.h5'
mixing.tol = 1e-5*[1 1]       % advised for dfpt
option.saveWavefunction = 0   % No need to save perturbed wavefunctions
option.savePotential = 1      % MUST be saved to be used in eph calc
symmetry.timereversal   = 1 
symmetry.spacesymmetry  = 0   % required for eph calc
symmetry.pointsymmetry  = 0   % required for eph calc
```

:::: important
::: title
Important
:::

Perturbed self-consistent potentials for atomic displacements
perturbation are required in electron-phonon calculations, hence they
must be stored using:

``` 
option.savePotential = 1
```
::::

:::: note
::: title
Note
:::

RESCU currently supports timereversal symmetry for calculating perturbed
self-consistent potentials. It must be ensured other symmetry settings
are deactivated if phonon results are to be used later for an
electron-phonon calculation:

``` matlab
symmetry.timereversal   = 1 
symmetry.spacesymmetry  = 0   
symmetry.pointsymmetry  = 0  
```
::::

RESCU is executed by

``` shell
rescu -i al_real_phonon
```

The important outputs of this step are already explained in sections
`dmandbe`{.interpreted-text role="ref"} and
`DFPTMetal`{.interpreted-text role="ref"}. It should be noted that for
system with an electronic band-gap we must perform a dielectric
calculation before taking the next step. This provides additional
quantities which are required to calculate LO-TO effect corresponding to
phonons in long wave-length limit, e.g see
`sec:dfptbs`{.interpreted-text role="ref"}, [\[B20\]](#B20){.citation}
and [\[GL97\]](#GL97){.citation}.

Now we can proceed with the electron-phonon calculation using the
following input file

``` matlab
info.calculationType = 'dfpt-elph'
info.savepath = './results/al_real_elph'
rho.in = './results/al_real_scf.mat' 
dfpt.phononData = './results/al_real_phonon.mat' % path to DFPT phonon data
kpoint.gridn = [5,5,5]   
dfpt.qpointGridn = [5,5,5]                 
dfpt.elphTemperature = linspace(5,300,6)    % required for elph
dfpt.elphZcut = 0.01                        % required for elph  
units.energy = 'eV'

dfpt.saveEPImat = 0             % Optional, default value 0    
dfpt.saveElectronVelocity = 0   % Optional, default value 0             
dfpt.saveElectronSE = 0         % Optional, default value 0       
dfpt.saveElectronTau   = 0      % Optional, default value 0    
```

Description of important keywords follows:

-   `info.calculationType` determines the type of the calculation;;
-   `info.savepath` points to the file which stores the output data of
    the elph calculation;
-   `rho.in` this DFPT calculation loads the ground state electronic
    density at the beginning of the computation, and this keyword points
    to the file which contains the output of the ground state
    calculation, so that RESCU can load corresponding electronic density
    data;
-   `dfpt.phononData` points to the path of the file containing the
    dynamical matrix (and Born effective charges);
-   `kpoint.grid` determines the k-sampling of the reciprocal space used
    for nonself-consistent calculation of ground state wavefunctions and
    electronic structure;
-   `dfpt.qpointGridn` determines the q-sampling of the reciprocal space
    used in interplation of dynamical matrix for phonons and perturbed
    self-consistent potentials;

:::: note
::: title
Note
:::

Current version of RESCU requires an identical, $\Gamma$-centered, and
odd `kpoint.gridn` and `dfpt.qpointGridn` for DFPT calculations;
::::

-   `dfpt.elphZcut` small imaginary parameter used in evaluating
    self-energy integrals. Here, it is combined with `units.energy`, so
    that it can be set in eV units;
-   `dfpt.elphTemperature` a row vector containing temperatures at which
    eph related variables are calculated;

::: {#para:dfptsave}
Some of the large eph related outputs can be saved using following
keywords:
:::

-   `dfpt.saveEPImat` optional keyword to save EPI matrix to
    `./results/al_real_elph.h` file. Default value is set to `0`.
    Results can be found in `dfpt/EPImat` field;
-   `dfpt.saveElectronVelocity` optional keyword to save electron
    velocities to `./results/al_real_elph.h` file. Default value is set
    to `0`. Results can be found in `dfpt/electronVelocitiy` field;
-   `dfpt.saveElectronSE` optional keyword to save electron self-energy
    to `./results/al_real_elph.h` file. Default value is set to `0`.
    Results can be found in `dfpt/electronSE` field;
-   `dfpt.saveElectronTau` optional keyword to save electron lifetime to
    `./results/al_real_elph.h` file. Default value is set to `0`.
    Results can be found in `dfpt/electronTau` field.

:::: warning
::: title
Warning
:::

Saving elph variables in `h5` file may require a large disk space.
::::

RESCU is executed as follows

``` shell
rescu -i al_real_elph
```

Important output variables can be found in `dfpt` structure of
`./results/bn_real_phonon.mat` file:

``` matlab
>> load results/al_real_elph.mat
>> dfpt

      dfpt = 

        struct with fields:

                          BornEC: [3x3 double]
                             IFC: [3x3x27 double]
                             asr: [1 1 1]
                            desc: [1x1 struct]
                dielectricTensor: [3x3 double]
                          dynMat: [3x3x125 double]
          electricalConductivity: [3x3x6 double]
                 elphTemperature: [5 64 123 182 241 300]
                  elphTransCoeff: [3x3x6 double]
                      phononFreq: [3x125 double]
                      phononMode: {125x1 cell}
                    qpointDirect: [125x3 double]
                     qpointGridn: [5 5 5]
                            type: 'elph'   
```

RESCU computes generalized transport coefficients using relaxation time
approximation (RTA) of linearized Boltzmann transport equations (LBTE),
where electrons lifetimes are evaluated from electrons self-energy due
to eph interactions [\[B20\]](#B20){.citation}:

[$$\begin{aligned}
\begin{split}
   \Sigma_{n\mathbf{k}}^{FM}(\omega, \varepsilon_{F}, T) = & \sum_{m, \nu} \int_{BZ} \frac{d\mathbf{q}}{\Omega_{BZ}} |g_{mn\nu}(\mathbf{k}, \mathbf{q})|^{2} \\
   & \times \Big[ \frac{n_{\mathbf{q}\nu}(T) + f_{m\mathbf{k+q}}(\varepsilon_{F},T)}{\omega-\varepsilon_{m\mathbf{k+q}}+\omega_{\mathbf{q}\nu} + i\eta} \\
   & + \frac{n_{\mathbf{q}\nu}(T) + 1 - f_{m\mathbf{k+q}}(\varepsilon_{F},T)}{\omega-\varepsilon_{m\mathbf{k+q}} - \omega_{\mathbf{q}\nu} + i\eta}\Big]
\end{split}
\end{aligned}$$]{label="eleSE"}

where $f_{m\mathbf{k+q}}(\varepsilon_{F},T)$ and $n_{\mathbf{q}\nu}(T)$
correspond to the Fermi-Dirac and Bose-Einstein distributions. In the
$\eta \rightarrow  0^{+}$ limit, the imaginary part evaluated at the
Kohn-Sham energy can be obtained by

[$$\begin{aligned}
\begin{split}
\lim_{\eta \rightarrow 0^{+}} \Im\{ \Sigma_{n\mathbf{k}}^{FM}(\varepsilon_{F})\} = &
\pi \sum_{m, \nu} \int_{BZ}  \frac{d\mathbf{q}}{\Omega_{BZ}} |g_{mn\nu}(\mathbf{k}, \mathbf{q})|^{2}
\\
& \times \Big[
(
n_{\mathbf{q}\nu} + f_{m\mathbf{k+q}})
\delta(\varepsilon_{n\mathbf{k}} - \varepsilon_{m\mathbf{k+q}} + \omega_{\mathbf{q}\nu}
)
\\
& +
(
n_{\mathbf{q}\nu} + 1 - f_{m\mathbf{k+q}}
)
\delta
(\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k+q}}-\omega_{\mathbf{q}\nu}
)
\Big]
\end{split}
\end{aligned}$$]{label="eleSEiamg"}

from which electrons lifetimes are calculated:

[$$\begin{aligned}
\frac{1}{\tau_{n\mathbf{k}}} = 2 \lim_{\eta \rightarrow 0^{+}} \Im\{ \Sigma_{n\mathbf{k}}^{FM}(\varepsilon_{F})\}.\end{aligned}$$]{label="eleLT"}

These coefficients are stored in `elphTransCoeff` as a
$3\times 3\times N_{temperature}$ array in atomic units:

``` matlab
>> dfpt.elphTransCoeff(:,:,1)

      ans =

         1.0e+04 *

          1.4984   -0.0020   -0.0020
         -0.0020    1.4984   -0.0020
         -0.0020   -0.0020    1.4984   
```

Also, electrical conductivity is stored as a
$3\times 3\times N_{temperature}$ array in (Siemens cm^-1^) units in
`electricalConductivity`:

``` matlab
>> dfpt.electricalConductivity(:,:,1)

      ans =

         1.0e+07 *

          1.2244   -0.0016   -0.0016
         -0.0016    1.2244   -0.0016
         -0.0016   -0.0016    1.2244   
```

For small data sizes, RESCU also outputs the following quantities in
text file format:

-   Electrons lifetimes in `results/Lifetime_electron`

``` text
# Electrons lifetime (femto-seconds) with SERTA approximation
# nbands = 10
# T =    5.0 (K)      Efermi = -4.338795e-02 (Ha)
# spin =  1
# kpoint = [0.0000000   0.0000000   0.0000000]
+3.830255e-27 +6.056664e-28 +4.693947e-28 +7.562906e-28 +1.244860e-30 +6.569030e-31 +4.699217e-31 +1.647781e-28 +1.298661e-28 +2.911374e-31 
# kpoint = [0.2000000   0.0000000   0.0000000]
...
```

-   Electrical conductivity in `results/Conductivity_electron`

``` text
# Electrical conductivity (Siemens/cm^-1) with SERTA approximation
# T =    5.0 (K)      Efermi = -4.338795e-02 (Ha)
+1.224412e+07   -1.645853e+04   -1.646234e+04   
-1.645853e+04   +1.224413e+07   -1.647936e+04   
-1.646234e+04   -1.647936e+04   +1.224416e+07      
```

-   Electrons self-energy in `results/Selfenergy_electron`

``` text
# Electrons self-energy (Hartree) with SERTA approximation
# nbands = 10
# T =    5.0 (K)      Efermi = -4.338795e-02 (Ha)
# spin =  1
# kpoint = [0.0000000   0.0000000   0.0000000]
-2.845214e-35 -7.637874e-38 +2.644990e-36 -4.830218e-37 -3.199617e-35 -6.232497e-37 -1.452827e-35 -3.868223e-37 +2.390173e-34 -2.350065e-34 
-6.871185e-35 -4.453474e-34 -6.396189e-34 -6.225507e-34 -2.519247e-35 -1.775418e-36 -6.435553e-35 -2.252705e-36 +8.060548e-34 -1.004852e-33 
# kpoint = [0.2000000   0.0000000   0.0000000]   
```

As mentioned `above<para:dfptsave>`{.interpreted-text role="ref"}, these
quantities as well as EPI matrix can also be saved in `h5` format. For
example EPI matrix information can be accessed using:

``` matlab
% (nbnad,nband) EPI matrix for [i_mode,i_qpt,i_kpt_ispin] = ind2sub([nmode,nqpt,nkpt,nspin],10)
>> h5info('results/al_real_elph.h5','/dfpt/EPImat10')    
      ans = 

        struct with fields:

            Filename: '/home/saeedb/Projects/nano/rescumat/doc-examples/al-eph/results/al_real_elph.h5'
                Name: 'EPImat10'
            Datatype: [1x1 struct]
           Dataspace: [1x1 struct]
           ChunkSize: []
           FillValue: 0
             Filters: []
          Attributes: [1x1 struct]

% real and imaginary parts
>> data = h5read('results/al_real_elph.h5','/dfpt/EPImat10');
>> size(data)        

      ans =

          10    20
% (nband,nband) transition rates
>> g = data(:,1:10) + 1j*data(:,11:20)   
>> size(g)

   ans =

       10    10
```

::: {#citations}

[B20]{#B20 .citation-label}

:   Guillaume Brunin et al.
    `Phonon-limited electron mobility in Si, GaAs, and GaP with exact treatment of dynamical quadrupoles. <10.1103/PhysRevB.102.094308>`{.interpreted-text
    role="doi"} Physical Review B 102.9 (2020), p. 094308.

[GL97]{#GL97 .citation-label}

:   Xavier Gonze and Changyol Lee.
    `Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory. <10.1103/PhysRevB.55.10355>`{.interpreted-text
    role="doi"} Physical Review B 55.16 (1997), p. 10355.

[SS96]{#SS96 .citation-label}

:   S Yu Savrasov and D Yu Savrasov.
    `Electron-phonon interactions and related physical properties of metals from linear-response theory. <10.1103/PhysRevB.54.16487>`{.interpreted-text
    role="doi"} Physical Review B 54.23 (1996), p. 16487.
:::
