Part 4: Conductivity of iron ¶

Content¶

 11 Conductivity using the constant relaxation-time approximation
 12 Conductivity using the self-energy relaxation-time approximation
 13 The effect of the lattice parameter, density functional, and pseudopotential

11 Conductivity using the constant relaxation-time approximation

$\uparrow$

¶

By the end of this section, you will be able to:

• Compute the conductivity of a simple metal in the constant relaxation-time approximation (CRTA) and compare to experiment.
• Test the convergence of conductivity with respect to $\textbf{k}$-point mesh.
• Perform a spin-polarised calculation accounting for the magnetic nature of iron and compare to the non-magnetic result.

11.1 Task¶

Compute the conductivity of body-centered-cubic (bcc) iron in the constant relaxation-time approximation (CRTA) with and without spin-polarization, converging with respect to the k-point mesh used for the transport coefficients

In many systems, the electronic and vibrational degrees of freedom (phonons) are treated separately, as electrons are much faster than the motion of nuclei. For conductivity, the electron-phonon coupling is prevalent. In metals, electrical transport arises from the motion of charge carriers (i.e., electrons and holes) under an applied electric field. The conductivity $\sigma$ quantifies how easily current flows through a material and is the inverse of the resistivity $\rho$ by $\sigma = 1/\rho$ (not to be confused with the density). From the linearized Boltzmann equation, the transport distribution function can be expressed as:

$$\tag{11.1} \mathcal{T}(\epsilon) = \frac{e^2}{N_\mathbf{k}\Omega} \sum_{n\mathbf{k}} \tau_{n\mathbf{k}} \, \mathbf{v}_{n\mathbf{k}} \otimes \mathbf{v}_{n\mathbf{k}} \, \delta(\epsilon_{n\mathbf{k}}-\epsilon), $$

where $\alpha,\beta$ are Cartesian coordinates, $e$ is the elementary charge, $N_\textbf{k}$ is number of $\textbf{k}$-points, $\Omega$ is the volume, $\tau_{n\mathbf{k}}$ is the relaxation time for state $(n,\mathbf{k})$ - band $n$ and $\textbf{k}$-point, $\mathbf{v}_{n\mathbf{k}}$ is the group velocity component, and $\epsilon_{n\mathbf{k}}$ is the electronic energy.

The Onsager transport coefficients $\mathcal{L}_{ij}$ are given as a $(2 \times 2)$ matrix (i.e., $i, j \in \{1, 2\} $), defined as [J.M. Ziman, Electrons and Phonons (2001)]: $$\tag{11.2} \mathcal{L}_{ij} = \int d\epsilon \, \mathcal{T}(\epsilon) \, (\epsilon-\mu)^{i+j-2} \left( -\frac{\partial f^0}{\partial \epsilon} \right). $$

Using $\mathcal{L}_{ij}$, we can compute the following transport coefficients:

$\begin{alignedat}{2} &\text{Electrical conductivity:} &&\quad \sigma = \mathcal{L}_{11} \quad &&\text{(11.3)} \\[6pt] &\text{Seebeck coefficient:} &&\quad S = \dfrac{\mathcal{L}_{11}}{\mathcal{L}_{12}} \quad &&\text{(11.4)} \\[6pt] &\text{Electronic thermal conductivity:} &&\quad \kappa_e = \tfrac{1}{T}(\mathcal{L}_{22} - \mathcal{L}_{21}\mathcal{L}_{11}^{-1}\mathcal{L}_{12}) \quad &&\text{(11.5)} \end{alignedat}$

There is also a lattice contribution to the thermal conductivity that we omit here [Phys. Rev. B 85, 024102 (2012)].

In the constant relaxation-time approximation (CRTA), we assume that all electrons at the Fermi surface share the same average scattering lifetime τ (also referred to as the relaxation time). This τ acts as an overall scaling factor of the average of the Fermi-surface electron velocities, making CRTA a quick and elegant way to see how the geometry of the bands alone controls conductivity. For metals, the electrons at the Fermi surface are already free to move as there is no bandgap and no additional considerations, so CRTA is a reasonable model. CRTA models the purely electronic component, before we move into the details of scattering later.

In this exercise, you will apply CRTA to iron, calculate the conductivity, and see the effect of $\textbf{k}$-point mesh and spin-polarization on the conductivity.

11.2 Constant relaxation-time approximation¶

11.2.1 Input¶

The input files to run this example are prepared in $TUTORIALS/electron-phonon/e11_conductivity_iron_crta/crta.

POSCAR

bcc iron
   1.0000000000000000
    -1.3704630709330003    1.3704630709330003    1.3704630709330003
     1.3704630709330003   -1.3704630709330003    1.3704630709330003
     1.3704630709330003    1.3704630709330003   -1.3704630709330003
  Fe
     1
Direct
     0.0000000000000000    0.0000000000000000    0.0000000000000000

INCAR

PREC = Accurate
EDIFF = 1e-8
ISMEAR = -15
SIGMA = 0.1

#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=0.0 #no additional carriers
ELPH_SCATTERING_APPROX=CRTA #constant-relaxation-time approximation does not require phonon calculations
TRANSPORT_RELAXATION_TIME=1e-14

Several of these INCAR tags, you will have encountered previously. We will skip over the standard settings and go to the electron-phonon specific ones. elph_run sets the electron-phonon calculation to run, ELPH_MODE = TRANSPORT means that a transport calculation is set up. Read our Wiki article to find out exactly what settings are used by default. The carrier density (e.g., holes and electrons) is defined by ELPH_SELFEN_CARRIER_DEN and the relaxation-time approximation by ELPH_SCATTERING_APPROX, in this case, the constant relaxation-time approximation (CRTA). The transport relaxation time in then set with TRANSPORT_RELAXATION_TIME.

KPOINTS

k-mesh
       0
Gamma
 12 12 12

KPOINTS_ELPH

k-mesh
0
Gamma
24 24 24

POTCAR

PAW_PBE Fe 06Sep2000

11.2.2 Calculation¶

Open a terminal, navigate to this example's directory and run the calculation by entering the following:

cd $TUTORIALS/electron-phonon/e11_*/crta
mpirun -np 4 vasp_std

Before we proceed, let us outline the different steps:

1. An electronic minimization to obtain the ground-state solution.
2. A non-self-consistent calculation on a different $\textbf{k}$-mesh KPOINTS_ELPH (or ELPH_KSPACING) - usually denser than the one used to obtain the ground-state.
3. Compute the electron group velocities using the commutator of the position operator with the Hamiltonian.
4. Evaluate the transport function on a Gauss-Legendre mesh and compute the Onsager coefficients.
5. Compute the transport coefficients.

All these steps are performed inside a single execution of vasp_std.

If you open the OUTCAR file, you can find details of the calculation. The first part, up until electron-phonon calculation is the electronic minimization, i.e., step 1. Various parameters of the calculation then follow, e.g.:

Transport:
  elph_transport = T
  elph_transport_nedos_plot =     -1
  elph_transport_emin_plot = min(eigenvalues)
  elph_transport_emax_plot = max(eigenvalues)
  elph_transport_dfermi_tol = .100E-05
  elph_transport_driver =   2
  elph_transport_nedos =    501
  transport_relaxation_time = .100E-13
  elph_velocity_mode =   2

Most of these tags are already defined in the input section; the remainder are the defaults or set by the ELPH_MODE=TRANSPORT tag, which offers a convenient way to set up a transport calculation. After comes the selection of electron bands and $\textbf{k}$-points that are used in the calculation (step 2):

 Transport energy range:   [   6.421:   7.461]  wich corresponds to    1.040 eV
 Number of selected k-points at which to compute the self-energy: [   412 /    413]
 Selected bands and k-points at which to compute the self-energy:
         3   4   5   6
     1   F   F   T   T
     2   F   F   T   T
     3   F   F   T   T
    ...

Finally, details of the calculation are shown, including the number of carriers, relaxation time, number of integration points, etc (steps 3-4):

Transport calculator N =  1
----------------------------------
 transport driver: 2 ! Gauss-Legendre integration
 Scattering approximation: constant relaxation-time approximation (CRTA)
 Static self-energy: F
 Transport number of points:   501
 temperature:    0.000 K
   Transport energy range:   [   6.940:   6.940]  wich corresponds to    0.000 eV
   Average relaxation time:    1.0000E-14 s
   Number of electrons:    3.8104E-01
   Number of holes:        3.8104E-01

In the following output, you can read the value of the conductivity by looking for sigma S/m (step 5):

Transport for self-energy accumulator N=     1
                 T K               mu eV           sigma S/m      mob cm^2/(V.s)       seebeck μV/K         peltier μV     kappa_e W/(m.K)
          0.00000000          6.93979457    9107828.46682731          0.00000000          0.00000000          0.00000000          0.00000000 Gauss-Legendre grids
        100.00000000          6.93996744    9129582.81044275          0.62456167          4.73450805        473.45080477         22.42750033 Gauss-Legendre grids
        200.00000000          6.94042241    9176143.35291047          2.74966584          9.26489540       1852.97908042         45.70863229 Gauss-Legendre grids
        300.00000000          6.94083401    9250240.03726177          4.26653127         12.89607103       3868.82130885         70.37075010 Gauss-Legendre grids
        400.00000000          6.94104058    9348336.51285969          5.24063646         15.54209909       6216.83963735         96.80282375 Gauss-Legendre grids
        500.00000000          6.94075176    9466692.72189808          5.93101115         17.08579625       8542.89812502        124.76327117 Gauss-Legendre grids

Before the end of the electron-phonon subroutines is reached:

------------------------ end of transport driver reached -----------------------------------------------

The key quantities reported are:

• sigma S/m: Electrical conductivity in Siemens per meter
• seebeck μV/K: Seebeck coefficient in microvolts per Kelvin
• kappa_e W/(m.K): Electronic thermal conductivity in Watts per meter-Kelvin

In this particular case for 300 K we obtained a conductivity of 9.2x10⁶ S m^-1. The value of the resistivity of iron, the how much it opposes the flow of current, is 97.1 nΩ m [CRC Handbook of Chemistry and Physics], which corresponds to a conductivity of 1.03x10⁷ S m^-1. Not bad for a first calculation! Confirm if you obtained the same result.

What about the electronic contribution to the thermal conductivity? How does our value compare with a reference?

The agreement with experiment is fortuitous. CRTA coincidentally leads to good agreement with experiment at 300 K. How can we know which is the correct value for the material we are considering? We will explore this question in the rest of the exercise.

First, let us investigate the convergence of the conductivity with $\textbf{k}$-point sampling!

11.3 Convergence with k-point sampling¶

In this section, we will see how the calculated conductivity depends on the density of the $\textbf{k}$-point mesh. As we increase the number of $\textbf{k}$-points, our results become more accurate—think of it as sharpening the resolution on a fuzzy picture of the Fermi surface! Denser meshes capture the fine details needed for reliable transport properties, but they also require more computational effort. Let's find out how much the conductivity changes as we increase the $\textbf{k}$-point density.

11.3.1 Input¶

The input files to run this example are prepared at $TUTORIALS/electron-phonon/e11_conductivity_iron_crta/kpoints/NxNxN.

INCAR

PREC = Accurate
EDIFF = 1e-8
ISMEAR = -15
SIGMA = 0.1

#transport
ELPH_MODE=TRANSPORT
ELPH_ISMEAR=-24
ELPH_SELFEN_CARRIER_DEN=0.0 #no additional carriers
ELPH_SCATTERING_APPROX=CRTA #constant-relaxation-time approximation
TRANSPORT_RELAXATION_TIME=1e-14

KPOINTS

k-mesh
       0
Gamma
 12 12 12

KPOINTS_ELPH

k-mesh
0
Gamma
10 10 10

k-mesh
0
Gamma
20 20 20

k-mesh
0
Gamma
30 30 30

k-mesh
0
Gamma
40 40 40

k-mesh
0
Gamma
50 50 50

k-mesh
0
Gamma
60 60 60

POTCAR

PAW_PBE Fe 06Sep2000

The INCAR file contains the same transport-related tags as the basic calculation, except for ELPH_ISMEAR, which selects the smearing method to determine the Fermi level and chemical potential before an electron-phonon calculation.

The key difference is that we use different $\textbf{k}$-point meshes to study convergence:

• KPOINTS: Used for the self-consistent ground-state DFT calculation
• KPOINTS_ELPH: for the non-self-consistent DFT calcualtion and computation of the transport properties.

Instead of KPOINTS_ELPH, one can use ELPH_KSPACING to specify the $\textbf{k}$-point mesh in the same spirit as KSPACING can be used instead of KPOINTS. In this part of the tutorial we will focus on the convergence of the mesh that will be used for the transport calculation.

The transport properties, especially electrical conductivity, are very sensitive to the $\textbf{k}$-point sampling because:

1. They depend on the Fermi-surface geometry
2. The Dirac delta function $\delta(\varepsilon_{n\mathbf{k}} - \varepsilon_F)$ picks out states exactly at the Fermi level

Therefore, convergence studies with respect to $\textbf{k}$-point density are crucial for accurate conductivity calculations within the approximations considered.

11.3.2 Calculation¶

Open a terminal and navigate to the $\textbf{k}$-point convergence directory:

cd $TUTORIALS/electron-phonon/e11_conductivity_iron_crta/kpoints
bash kpoints.sh

for k in 10 20 30 40 50 60; do
    cp ../crta/WAVECAR "${k}x${k}x${k}"
    cd "${k}x${k}x${k}"
    echo "Running calculation with ${k}x${k}x${k} k-mesh..."
    mpirun -np 4 vasp_std
    cd ../
done

This will run calculations with the same 12x12x12 mesh for the ground state calculation (KPOINTS), so we copy the WAVECAR from 11.2.2, then use multiple $\textbf{k}$-point meshes for the transport evaluation (KPOINTS_ELPH).

The transport coefficients are reported in the standard output and OUTCAR. Look for the section that begins with:

Transport for self-energy accumulator N=     1
                 T K               mu eV           sigma S/m      seebeck μV/K        kappa_e W/(m.K)

For each $\textbf{k}$-point mesh, you can extract the conductivity values at 300 K using:

cd $TUTORIALS/electron-phonon/e11_conductivity_iron_crta/
grep 300.00000 kpoints/*/OUTCAR | grep Gauss | awk -F'/' '{split($2,a,"x"); k=a[1]} {print k, $0}' | awk '{print $1, $5}'

It is easy to plot the conductivity convergence with respect to $\textbf{k}$-point grid, reading the transport properties using py4vasp:

The conductivity calculations require very dense $\textbf{k}$-point meshes because:

1. Fermi-surface sampling: The conductivity depends on electronic states exactly at the Fermi level. Accurate sampling of the Fermi surface is crucial.

2. Integration convergence: The transport integrals converge slowly with $\textbf{k}$-point density, particularly for metals with complex Fermi surfaces.

As the $\textbf{k}$-point mesh becomes denser, the conductivity calculation converges. Notice how the convergence is not completely smooth, with the 60x60x60 mesh showing a slight increase. Convergence can be slow, especially for:

• Complex Fermi surfaces: Metals with intricate Fermi surface topologies require very dense sampling
• Anisotropic materials: Different crystallographic directions may converge at different rates (not the case of bcc Iron)
• Light effective masses: Bands with high curvature (light effective masses) need denser $\textbf{k}$-meshes

For production calculations, you should monitor how the conductivity changes with $\textbf{k}$-point mesh sampling. For convergence, it is important to take a value close to the converged calculation, not the experiment. In the context of this tutorial and in the interest of time, we will consider 5% error from the converged value (cf. the larger 50x50x50 and 60x60x60 meshes) sufficient, for bcc iron means using a 30×30×30 mesh or denser in KPOINTS_ELPH.

11.4 Magnetism¶

Up to now, we have ignored the fact that iron is actually a magnetic material, with an on-site magnetic moment of around 2.2 $\mu_B$. That means our earlier calculations were “right for the wrong reasons”—we left out a key ingredient! In this section, we will rerun the conductivity calculation, but this time we will do it properly by performing a collinear spin calculation (ferromagnetic configuration).

Consider how you think including magnetism will affect the conductivity?

11.4.1 Input¶

Check out the input files in $TUTORIALS/electron-phonon/e11_conductivity_iron_crta/kpoints_mag/NxNxN!

INCAR

PREC = Accurate
EDIFF = 1e-8
ISMEAR = -15
SIGMA = 0.1
ISPIN = 2
MAGMOM = 6

#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=0.0 #no additional carriers
ELPH_SCATTERING_APPROX=CRTA #constant-relaxation-time approximation does not require phonon calculations
TRANSPORT_RELAXATION_TIME=1E-14

KPOINTS

k-mesh
       0
Gamma
12 12 12

KPOINTS_ELPH

k-mesh
0
Gamma
10 10 10

k-mesh
       0
Gamma
20 20 20

k-mesh
       0
Gamma
30 30 30

k-mesh
       0
Gamma
40 40 40

POTCAR

PAW_PBE Fe 06Sep2000

We have introduced several new tags in the INCAR:

ISPIN and MAGMOM, which turn on spin-polarized calculations and specify initial magnetic moment for each atom, respectively, KPAR, which determines the number of $\textbf{k}$-points that are treated in parallel, IBRION=6, which computes the finite atomic displacements, ELPH_POT_GENERATE calculates the electron-phonon potential from finite atomic displacements and LPHON_DISPERSION requests the calculation of the phonon dispersion along the $\textbf{q}$-point path supplied in file QPOINTS.

11.4.2 Calculation¶

Now we will repeat the previous calculation but this time on a spin-polarized calculation, thus including the effect of the usual ferromagnetic orientation in bcc iron:

cd "${TUTORIALS}/electron-phonon/e11_conductivity_iron_crta/kpoints_mag"
bash kpoints.sh

You will run calculations using 10x10x10, 20x20x20, and 30x30x30 $\textbf{k}$-point meshes, then compare to a pre-calculated 40x40x40.

for k in 10 20 30; do
    cd "${k}x${k}x${k}"
    echo "Running calculation with ${k}x${k}x${k} k-mesh..."
    mpirun -np 4 vasp_std
    cd ../
done

Now we will plot the results together with the non-magnetic case:

You can see that the non-magnetic result seems to converge more closely to the experimental value. However, the inclusion of magnetism is crucial for a physically accurate description of iron. Metals with intricate Fermi surface topologies require very dense sampling, hence the large $\textbf{k}$-point mesh required for convergence in each case.

11.5 Questions¶

1. What setting for the ELPH_MODE tag should you use to run a transport calculation?
2. Why does the electrical conductivity of iron converge slowly with respect to $\textbf{k}$-point mesh density?
3. If the non-magnetic solution is closer to experiment, can we ignore the magnetism of iron?

12 Conductivity using the self-energy relaxation-time approximation

$\uparrow$

¶

By the end of this section, you will be able to:

• Compute the conductivity of a simple metal in the self-energy relaxation-time approximation (SERTA).
• Test the convergence of conductivity with respect to $\textbf{k}$-point mesh and compare to CRTA.
• Visualize the scattering lifetimes as a function of the Kohn-Sham energies and learn how they influence the conductivity.

12.1 Task¶

Compute the conductivity of iron (2x2x2 cell) in the self-energy relaxation-time approximation (SERTA), converging with respect to the k-point mesh used for the transport coefficients

Now that we have an idea of the $\textbf{k}$-point mesh density required for accurate conductivity calculations, we can proceed to compute the relaxation times from first principles. This involves calculating the derivatives of the Kohn-Sham potential in a supercell w.r.t ionic displacements, so that as many phonons as possible can be considered. The transport calculation is then rerun with ELPH_SCATTERING_APPROX = SERTA [Phys. Rev. B 97 (2018) 121201(R)]:

$$\tag{12.1} \frac{1}{\tau^{\mathrm{SERTA}}_{n\mathbf{k}}} = \sum_{n'\nu\mathbf{k}'} \frac{2\pi}{\hbar} w_{n\mathbf{k},n'\mathbf{k}'} |g^{\nu}_{n\mathbf{k},n'\mathbf{k}'}|^2 \left[ (n_{\nu\mathbf{q}} + 1 - f_{n'\mathbf{k}'}) \, \delta(\varepsilon_{n\mathbf{k}} - \varepsilon_{n'\mathbf{k}'} - \hbar\omega_{\nu\mathbf{q}}) + (n_{\nu\mathbf{q}} + f_{n'\mathbf{k}'}) \, \delta(\varepsilon_{n\mathbf{k}} - \varepsilon_{n'\mathbf{k}'} + \hbar\omega_{\nu\mathbf{q}}) \right]$$

where ${\tau^{\mathrm{SERTA}}_{n\mathbf{k}}}$ is the relaxation time (or scattering time, or lifetime) for state $(n,\mathbf{k})$ using SERTA, $w_{n\mathbf{k},n'\mathbf{k}'}$ is the weight (set to 1 for SERTA), $g^{\nu}_{n\mathbf{k},n'\mathbf{k}'}$ is the electron-phonon coupling matrix element, $f_{n\mathbf{k}}$ is the population of the electronic state (Fermi-Dirac distribution), $n_{\nu\mathbf{q}}$ is the population of the phononic state (Bose-Einstein distribution), $\varepsilon_{n\mathbf{k}}$ is the energy of an electron band, $\omega_{\nu\mathbf{q}}$ is the phonon frequency, and $\delta$ is the Dirac delta function.

The scattering lifetime ${\tau^{\mathrm{SERTA}}_{n\mathbf{k}}}$ relates to the transport distribution function (cf. Eq. 11.1), which, in turn, allows the calculation of Onsager coefficients (cf. Eq. 11.2) and, thereby, the conductivity (cf. Eq. 11.3). Therefore, in SERTA, the energy of an electron band $\varepsilon_{n\mathbf{k}}$ influences the conductivity.

Let us first compute the derivative of the Kohn-Sham potential in the supercell, remembering to do a spin-polarized calculation.

12.2 Electron-phonon potential on a supercell¶

The first step is to do the supercell calculation. The input files are prepared at $TUTORIALS/electron-phonon/e12_conductivity_iron_serta/supercell.

12.2.1 Input¶

POSCAR

Fe 2x2x2 supercell of the conventional cell
   1.0
     5.4818522837320005    0.0000000000000000    0.0000000000000000
     0.0000000000000000    5.4818522837320005    0.0000000000000000
     0.0000000000000000    0.0000000000000000    5.4818522837320005
Fe
  16
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.5000000000000000  0.0000000000000000  0.0000000000000000
  0.0000000000000000  0.5000000000000000  0.0000000000000000
  0.5000000000000000  0.5000000000000000  0.0000000000000000
  0.0000000000000000  0.0000000000000000  0.5000000000000000
  0.5000000000000000  0.0000000000000000  0.5000000000000000
  0.0000000000000000  0.5000000000000000  0.5000000000000000
  0.5000000000000000  0.5000000000000000  0.5000000000000000
  0.2500000000000000  0.2500000000000000  0.2500000000000000
  0.7500000000000000  0.2500000000000000  0.2500000000000000
  0.2500000000000000  0.7500000000000000  0.2500000000000000
  0.7500000000000000  0.7500000000000000  0.2500000000000000
  0.2500000000000000  0.2500000000000000  0.7500000000000000
  0.7500000000000000  0.2500000000000000  0.7500000000000000
  0.2500000000000000  0.7500000000000000  0.7500000000000000
  0.7500000000000000  0.7500000000000000  0.7500000000000000

INCAR

#
PREC = Accurate
EDIFF = 1e-8
ISMEAR = 1
SIGMA = 0.1
ISPIN = 2
MAGMOM = 6.0*16
NBANDS = 96
#paralelization
KPAR = 4

#phonon dispersion
LPHON_DISPERSION = .TRUE.

#supercell
IBRION=6
POTIM=0.015
ELPH_POT_GENERATE=TRUE

KPOINTS

k-mesh
       0
Gamma
 2 2 2

POTCAR

PAW_PBE Fe 06Sep2000

All of the INCAR tags you have come across previously, with the exception of ELPH_POT_GENERATE, which calculates the electron-phonon potential from finite atomic displacements, and LPHON_DISPERSION which requests the calculation of the phonon dispersion along the $\textbf{q}$-point path supplied in QPOINTS file.

12.2.2 Calculation¶

To run the supercell calculation for electron-phonon coupling, open a terminal and navigate to the supercell directory:

cd $TUTORIALS/electron-phonon/e12_conductivity_iron_serta/supercell
python3 kpath.py -o QPOINTS

then run the VASP calculation

mpirun -np 4 vasp_std

This will generate the phelel_params.hdf5 file, necessary for the subsequent electron-phonon and transport calculations.

Take a look at the phonon dispersion:

12.3 Phonon limited electronic transport¶

Now we can finally compute the electronic transport, including scattering with phonons. The calculation looks exactly like the ones we have performed in the CRTA, with the difference that we need to pass on the phelel_params.hdf5 file and set ELPH_SCATTERING_APPROX=SERTA

12.3.1 Input¶

The input files are prepared at $TUTORIALS/electron_phonon/e12_conductivity_iron_serta/kpoints/NxNxN.

POSCAR

Fe primitive
   1.0000000000000000
     1.3704630709330001    1.3704630709330003   -1.3704630709329997
    -1.3704630709330001    1.3704630709329997    1.3704630709329997
     1.3704630709330001   -1.3704630709329997    1.3704630709330003
  Fe
     1
Direct
     0.0000000000000000    0.0000000000000000    0.0000000000000000

INCAR

PREC = Accurate
EDIFF = 1e-8
ISMEAR = -15
SIGMA = 0.1
ISPIN = 2
MAGMOM = 2.0

#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=0.0 #no additional carriers
ELPH_SCATTERING_APPROX=SERTA

KPOINTS

k-mesh
       0
Gamma
12 12 12

KPOINTS_ELPH

k-mesh
0
Gamma
10 10 10

k-mesh
0
Gamma
15 15 15

k-mesh
0
Gamma
20 20 20

k-mesh
0
Gamma
30 30 30

k-mesh
0
Gamma
40 40 40

POTCAR

PAW_PBE Fe 06Sep2000

ELPH_SCATTERING_APPROX=SERTA turns on the Self-Energy Relaxation-Time Approximation (SERTA). ELPH_SELFEN_CARRIER_DEN defines the number of additional carriers to be added in cm^-3.

There are four directories in e12_conductivity_iron_serta/kpoints/ with four different KPOINTS_ELPH meshes. You will only use the smallest two $\textbf{k}$-meshes. The others are too time-consuming for this tutorial, so we include their vaspout.h5 files for plotting later. If you have more time, you can run the 20x20x20, 30x30x30, and 40x40x40 calculations; it will take ~15 minutes to run the first two, and an hour for the third.

12.3.2 Calculation¶

Navigate to the kpoints directory and copy the supercell phelel_params.hdf5 file across to each directory. You can do this and then run the calculation for different $\textbf{k}$-point samplings (10x10x10 and 15x15x15) in the same step with the following:

cd $TUTORIALS/electron-phonon/e12_conductivity_iron_serta/kpoints/
bash kpoints.sh

for k in 10 15; do
    cp ../supercell/phelel_params.hdf5 "${k}x${k}x${k}"
    cd "${k}x${k}x${k}"
    mpirun -np 4 vasp_std
    cd ../
done

Once the calculation is finalized you can find the conductivity as a function of $\textbf{k}$-point sampling:

cd $TUTORIALS/electron-phonon/e12_conductivity_iron_serta/
grep 300.00000 kpoints/*/OUTCAR | grep Gauss | awk -F'/' '{split($2,a,"x"); k=a[1]} {print k, $0}' | awk '{print $1, $5}' 

To save time, we have pre-computed calcuations using 20x20x20, 30x30x30, and 40x40x40 $\textbf{k}$-points meshes. You can plot these with your results using the python script below:

Our calculated conductivity is now 40% more than the experimental value! The inclusion of phonon scattering has worsened the agreement with experiment. Sometimes, "more physics" just means "more problems"!

Recalling equations 11.1, 11.2, and 11.3, you can see how the lifetimes calculated using Eq. 12.1 connect to the transport distribution function (Eq. 11.1):

$$\tag{11.1} \mathcal{T}(\epsilon) = \frac{e^2}{N_\mathbf{k}\Omega} \sum_{n\mathbf{k}} \tau_{n\mathbf{k}} \, \mathbf{v}_{n\mathbf{k}} \otimes \mathbf{v}_{n\mathbf{k}} \, \delta(\epsilon_{n\mathbf{k}}-\epsilon), $$

which allows calculation of the Onsager coefficients:

$$\tag{11.2} \mathcal{L}_{ij} = \int d\epsilon \, \mathcal{T}(\epsilon) \, (\epsilon-\mu)^{i+j-2} \left( -\frac{\partial f^0}{\partial \epsilon} \right), $$

and therefore the conductivity:

$$\tag{11.3} \sigma = \mathcal{L}_{11}, $$

A change in scattering lifetimes should therefore change the conductivity. Let us try to understand what is going on by plotting the scattering lifetimes as a function of the Kohn-Sham energies (i.e., the energy of the electronic bands) at 300 K.

The purple line is if you assume that all carriers have the same lifetime $\tau$, i.e., 10^-14 s; this constant relaxation time is the approximation made in CRTA. In this figure, you can see that the carriers have a wide range of lifetimes; there is a lot of deviation in $\tau$, particularly at higher carrier energies. Varying the lifetimes in this way is the approximation made in SERTA, where $\tau$ contains many values and is not a single number (in contrast to CRTA, which fixes the lifetime independent of carrier energy).

Mostly, the lifetimes are less than 2x10^-14 s, but there are many carriers in the upper-right portion of the plot that go up to ~9x10^-14 s. The longer the lifetime, the higher the conductivity. By including these carriers, the conductivity increases a lot, hence the greater disagreement with experiment for SERTA. CRTA fails to consider these higher energies, so it matches experiment by coincidence, rather than a good description of the physics. SERTA predicts a higher conductivity by including the higher energy and higher mobility carriers, which CRTA neglects.

Try changing from the 40x40x40 $\textbf{k}$-mesh to another one to see how $\textbf{k}$-mesh convergence changes the plot.

In the next section, you will consider how the choice of pseudopotential, density functional, and lattice parameter can influence the conductivity.

12.4 Questions¶

1. Does CRTA or SERTA use averaged lifetimes?
2. How does including a range of lifetimes affect the conductivity?
3. How do you define the k- and q-point meshes?

13 The effect of the lattice parameter, density functional, and pseudopotential

$\uparrow$

¶

By the end of this section, you will have

• Computed the conductivity of iron at the experimental lattice parameter.
• Understand how different density functionals affect the electronic structure and transport properties.
• Compare conductivity calculations using different iron pseudopotentials (standard, GW).

13.1 Task¶

Generate the electron-phonon potential, then calculate the conductivity using SERTA (10x10x10, 20x20x20, 30x30x30 q-point meshes) for bcc iron at the experimental lattice parameter, while varying both the pseudopotential (standard and GW) the and density functional (CA and PBE)

The electrical conductivity of metals is intrinsically connected to their electronic band structure, which in turn depends on the crystal lattice. If we change the lattice parameter (and thus the volume) from the calculated PBE value to the experimental value [CRC Handbook of Chemistry and Physics], several effects occur:

1. Band width changes: Changing the lattice parameter alters the overlap between atomic orbitals, affecting the band width of the conduction bands.

2. Fermi surface modification: The size and shape of the Fermi surface change with volume, directly affecting the number of charge carriers and their velocities.

3. Density of states variations: The electronic density of states at the Fermi level $N(\varepsilon_F)$ changes with volume, influencing the conductivity.

The choice of pseudopotential can significantly affect the calculated electronic properties and transport coefficients. Different pseudopotentials vary in:

1. Core-valence partitioning: How many electrons are treated explicitly vs. frozen in the core
2. Cutoff energy requirements: Harder pseudopotentials require higher plane-wave cutoffs
3. Accuracy of the band structure: Particularly important near the Fermi level for transport

We will calculate the conductivity with each potential to assess:

• The difference that a choice of a more or less accurate pseudopotential can make in the final result.
• What relative error can be expected due to the choice of a different exchange-correlation functional.

13.2 Input¶

The first step is as before: to run the supercell calculation. However, previously, a PBE lattice parameter was used; this time, the experimental lattice parameter will be used. The input files to run this example are prepared in $TUTORIALS/electron-phonon/e13_conductivity_iron_exp/, which is a parent directory with the structure:

tree -d

├── exp_gw_lda
│   ├── serta_exp_gw_lda_10x10x10
│   ├── serta_exp_gw_lda_20x20x20
│   ├── serta_exp_gw_lda_30x30x30
│   └── supercell_exp_gw_lda
├── exp_gw_pbe
│   ├── serta_exp_gw_pbe_10x10x10
│   ├── serta_exp_gw_pbe_20x20x20
│   ├── serta_exp_gw_pbe_30x30x30
│   └── supercell_exp_gw_pbe
├── exp_std_lda
│   ├── serta_exp_std_lda_10x10x10
│   ├── serta_exp_std_lda_20x20x20
│   ├── serta_exp_std_lda_30x30x30
│   └── supercell_exp_std_lda
└── exp_std_pbe
    ├── serta_exp_std_pbe_10x10x10
    ├── serta_exp_std_pbe_20x20x20
    ├── serta_exp_std_pbe_30x30x30
    └── supercell_exp_std_pbe

where exp_PP_DF means that the calculations use the experimental lattice parameter, standard or GW POTCARs, and the LDA or PBE density functional. serta_exp_PP_DF_NxNxN then defines that a SERTA calculation is performed with an NxNxN $\textbf{k}$-point mesh (cf. KPOINTS_ELPH file) is used, and supercell_exp_PP_DF. You will only run the smallest $\textbf{k}$-point mesh 10x10x10, as the other calculations take a long time. They have been run previously and their vaspout.h5 supplied for comparing against later. To conserve space, we do not replicate duplicate files and instead label directories according to the above numbering.

POSCAR

Fe 2x2x2
   1.0
     5.7199999999999998    0.0000000000000000    0.0000000000000000
     0.0000000000000000    5.7199999999999998    0.0000000000000000
     0.0000000000000000    0.0000000000000000    5.7199999999999998
Fe
  16
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.5000000000000000  0.0000000000000000  0.0000000000000000
  0.0000000000000000  0.5000000000000000  0.0000000000000000
  0.5000000000000000  0.5000000000000000  0.0000000000000000
  0.0000000000000000  0.0000000000000000  0.5000000000000000
  0.5000000000000000  0.0000000000000000  0.5000000000000000
  0.0000000000000000  0.5000000000000000  0.5000000000000000
  0.5000000000000000  0.5000000000000000  0.5000000000000000
  0.2500000000000000  0.2500000000000000  0.2500000000000000
  0.7500000000000000  0.2500000000000000  0.2500000000000000
  0.2500000000000000  0.7500000000000000  0.2500000000000000
  0.7500000000000000  0.7500000000000000  0.2500000000000000
  0.2500000000000000  0.2500000000000000  0.7500000000000000
  0.7500000000000000  0.2500000000000000  0.7500000000000000
  0.2500000000000000  0.7500000000000000  0.7500000000000000
  0.7500000000000000  0.7500000000000000  0.7500000000000000

Fe primitive
   1.0000000000000000
     1.4299999999999999    1.4300000000000002   -1.4299999999999997
    -1.4299999999999999    1.4299999999999997    1.4299999999999997
     1.4299999999999999   -1.4299999999999997    1.4300000000000002
  Fe
     1
Direct
     0.0000000000000000    0.0000000000000000    0.0000000000000000

INCAR

#
PREC = Accurate
EDIFF = 1e-8
ISMEAR = 1
SIGMA = 0.1
ISPIN = 2
MAGMOM = 6.0*16
NBANDS = 96

#supercell
IBRION=6
ELPH_POT_GENERATE=TRUE

#
PREC = Accurate
EDIFF = 1e-8
ISMEAR = 1
SIGMA = 0.1
ISPIN = 2
MAGMOM = 6.0*16

#supercell
IBRION=6
ELPH_POT_GENERATE=TRUE

PREC = Accurate
EDIFF = 1e-8
ISMEAR = -15
SIGMA = 0.1
ISPIN = 2
MAGMOM = 2.0

#transport
ELPH_MODE=TRANSPORT
ELPH_SELFEN_CARRIER_DEN=0.0 #no additional carriers
ELPH_SCATTERING_APPROX=SERTA

KPOINTS

Automatic
0
Gamma
2 2 2

Automatic
0
Gamma
12 12 12

KPOINTS_ELPH

Automatic
0
Gamma
10 10 10

Automatic
0
Gamma
20 20 20

Automatic
0
Gamma
30 30 30

POTCAR

PAW Fe_sv_GW 05Dec2013

PAW_PBE Fe_sv_GW 05Dec2013

PAW Fe 03Mar1998

PAW_PBE Fe 06Sep2000

All INCAR tags you have previously encountered. Note that the density functional defaults to the one defined by the POTCAR file, so it does not need to be defined in the INCAR file in this case.

The phelel_params.hdf5 is already provided, so there is no need to repeat the supercell calculation.

13.3 Conductivity using the experimental lattice parameter¶

To begin with, navigate to the exercise directory:

cd $TUTORIALS/electron-phonon/e13_conductivity_iron_exp/

In this exercise, you will compare the conductivity using SERTA for 2 different DFT functionals (CA (LDA) and PBE (GGA)), using 2 different POTCARs (Standard and GW), and with 3 different $\textbf{k}$-point meshes (10x10x10, 20x20x20, and 30x30x30), totalling 16 calculations (12 SERTA, 4 supercell). For the purposes of this tutorial, we have run the supercell calculations previously, so the phelel_params.hdf5 files are already prepared. You will run the 10x10x10 $\textbf{k}$-point meshes for (LDA and PBE) exchange-correlation functionals and (Standard and GW) pseudopotentials, totalling 4 calculations. You can run these conductivity calculations efficiently using the following command:

bash conductivity.sh

for calc in exp_std_pbe exp_std_lda exp_gw_pbe exp_gw_lda;
do
    # run supercell
    cd $calc/supercell_$calc
    #mpirun -np 4 vasp_std | tee vasp.out
    cd -

    # run transport
    for mesh in 10 #20 30;
    do
        cp $calc/supercell_$calc/phelel_params.hdf5 $calc/serta_${calc}_${mesh}x${mesh}x${mesh}
        cd $calc/serta_${calc}_${mesh}x${mesh}x${mesh}
        mpirun -np 4 vasp_std | tee vasp.out
        cd -
    done
done

Once the calculations are finished, you can plot the results using the following script:

The transport calculations in this section are performed using 10×10×10, 20×20×20, and 30×30×30 $\textbf{k}$-point meshes. Your CRTA calculation from earlier is also included. The decay looks quite linear and (cf. the R² values that have been printed out, showing a reasonable fit) so we can extrapolate the conductivities to the "infinite" $\textbf{k}$-point limit.

Comparing the two pseudopotentials (PP) used, standard std (blue and purple) and GW (red and cyan), you can see that the std pseudopotential underestimates the conductivity relative to experiment with increasing $\textbf{k}$-point density, while using GW pseudopotentials overestimates. At the "infinite" $\textbf{k}$-point limit, this does not change, though GW lie closer to experiment. The GW pseudopotentials are more accurate and were prepared to consider unoccupied states, reproducing the all-electron scattering properties at high energies.

In fact, the PP is the dominant effect in this case. Moving from an LDA to a GGA: CA (blue and red) and PBE (purple and cyan), respectively (see this article for more details on the functionals used here), makes a smaller difference, with both underestimating relative to experiment, worsening as convergence is approached. Little difference is seen at the infinite limit. LDA often performs well for metals, so the matching for LDA and PBE is expected. However, if you were to use LDA for relaxing the structure, you would get a poor lattice parameter, so it is important to carefully choose your density functional in any case.

You can also see how big the difference that the lattice parameter can make, compare your SERTA calculations (PBE + std PP) at the calculated lattice parameter in this exercise (green) to the experimental lattice parameter from Example 12 (purple). The difference between the calculated and experimental lattice parameters is around 50%, which is enough to change the prediction from an underestimation to an overestimation.

The trends become even clearer when looking at the "infinite" $\textbf{k}$-point limit (the third plot) (Note the "SERTA magnetic - a (calc.)" value is the value for a 40x40x40 $\textbf{k}$-mesh, as this is already converged, so extrapolation can only distrort the value with the unconverged poinst). Here, using the std PP underestimates the conductivity by half (blue and purple), though still the correct order of magnitude. It is important to note that a factor of 2 in this case is relatively small on the order of 10⁷. Ideally, using the experimental lattice parameter would give the value closest to experiment; however, this is not always the case, as seen here. There is no established practice on whether it is better to use the experimental lattice parameter, or one specific to a density functional (it also changes between density functionals).

This tells us that the CRTA line (gray) that you calculated in Example 11 lies close to experiment for all the wrong reasons. It assumes a constant relaxation-time so it does not include as physical a description as SERTA (green) does, the $\textbf{k}$-mesh is not converged and agreement worsens as it improves, the pseudopotential (std) is inappropriate for this problem, and the lattice parameter does not match the experiment. All of these factors can be important for achieving agreement with experiment. It is easy to match experiment for the wrong reasons, making it essential to achieve numerical convergence by carefully testing in order to determine how large the error is for each factor. If you have time, try calculating the CRTA using all you have learned in this exercise and see how much it differs from experiment using converged settings.

There are still other factors that could be considered for comparison with experiment, including the supercell size, different magnetic phases, the impact of defects, the $\textbf{k}$-point mesh of the self-consistency calculation, and electron-electron scattering (e.g., post-DFT methods).

13.4 Questions¶

1. How can you obtain the "infinite" k-point limit for the conductivity?
2. What states do GW POTCARs consider that standard POTCARs fail to?

Well done! You have finished Part 4: Conductivity of metals!¶

In this tutorial, you have learned:

1. Basic conductivity calculations: How to compute the electrical conductivity of metals using the constant relaxation-time approximation (CRTA)

2. Convergence studies: The importance of dense k-point meshes for accurate transport property calculations

3. Lattice parameter effects: How volume changes affect the electronic structure and transport properties

4. Pseudopotential dependence: The role of pseudopotential choice in determining the reliability of conductivity calculations

Key takeaways:

• Transport properties require much denser k-point sampling than ground-state calculations
• The CRTA provides a good starting point for conductivity calculations
• Lattice parameter variations have significant effects on conductivity through band structure changes

Next steps: In Part 5, you will learn how to compute the mobilities of semiconductors limited by phonon scattering.

Go to Top $\uparrow$