Part 1: Bandgap renormalization from perturbation theory ¶
By the end of this tutorial, you will be able to:
- Compute the electron-phonon potential from supercells
- Calculate the phonon-induced bandgap renormalization
- Compare the phonon-induced bandgap at the zero-point and a range of temperatures.
1.1 Task¶
Calculate the zero-point and temperature-dependent phonon-induced bandgap renormalization for a 16-atom diamond supercell.
The computation of electron-phonon interactions is at a bare minimum a two-step procedure. In the first calculation, a supercell is used to compute the electron-phonon potential, while in the second calculation, a primitive unit cell is used to compute the bandgap renormalization due to electron-phonon interactions. This workflow allows for an efficient and accurate evaluation of electron-phonon effects by separating the computationally intensive supercell step from the final property calculations in the unit cell.
In the first step in this exercise, you will compute the electron-phonon potential of diamond using supercells and finite differences. Diamond, a covalent insulator with a wide indirect band gap, serves as an excellent model system to study electron-phonon physics. Its strong covalent bonds lead to high-frequency phonons and significant electron-phonon coupling. Understanding its bandgap renormalization is crucial for accurately predicting its optical and electronic properties at different temperatures.
The relevant quantities from this calculation are saved in an HDF5 file named phelel_params.hdf5. During the second step of this exercise, you will then use this information to calculate $g_{mn\nu}(\mathbf{k}, \mathbf{q})$ in a unit cell to compute the electron-phonon interactions and the resulting bandgap renormalization. In the end, you will compute the renormalization of the direct and indirect (fundamental) band gap of diamond at both zero and finite temperature.
1.2 Computing the electron-phonon potential¶
1.2.1 Input¶
diamond 2x2x2
1.0000000000
3.5618615596 -3.5618615596 0.0000000000
3.5618615596 0.0000000000 3.5618615596
0.0000000000 -3.5618615596 3.5618615596
C
16
Direct
0.0000000000 0.0000000000 0.0000000000
0.5000000000 0.0000000000 1.0000000000
0.0000000000 0.5000000000 1.0000000000
0.5000000000 0.5000000000 0.0000000000
0.0000000000 1.0000000000 0.5000000000
0.5000000000 1.0000000000 0.5000000000
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.5000000000 0.5000000000
0.3750000000 0.3750000000 0.3750000000
0.8750000000 0.3750000000 0.3750000000
0.3750000000 0.8750000000 0.3750000000
0.8750000000 0.8750000000 0.3750000000
0.3750000000 0.3750000000 0.8750000000
0.8750000000 0.3750000000 0.8750000000
0.3750000000 0.8750000000 0.8750000000
0.8750000000 0.8750000000 0.8750000000
system = diamond
ismear = 0; sigma = 0.1; efermi = midgap
prec = accurate
ediff = 1e-7
encut = 400
ibrion = 6
potim = 0.015
elph_pot_generate = True
k-mesh
0
Gamma
1 1 1
PAW_PBE C_s 06Sep2000
The POSCAR file defines a supercell structure of diamond. In this case, we have a 2x2x2 replication of the primitive unit cell. A supercell is required in order to use finite atomic displacements.
Note that for the sake of this tutorial, we only use a relatively small 2x2x2 supercell since larger cells would be too computationally demanding. In practice, however, the supercell should be carefully chosen by conducting convergence tests by testing your property of interest, e.g. bandgap renormalization for incrementally increasing supercell sizes, i.e., 2x2x2, 3x3x3, 4x4x4.
The INCAR file contains two important tags that are required to compute the electron-phonon potential:
- IBRION, which activates the finite-difference driver in VASP
- ELPH_POT_GENERATE, which computes the electron-phonon potential.
1.2.2 Calculation¶
Open a terminal, navigate to this example's directory (supercell) and run VASP by entering the following:
cd $TUTORIALS/electron-phonon/e01_*/supercell
mpirun -np 4 vasp_gam
VASP first performs an electronic minimization to obtain the ground-state solution in the supercell. Afterwards, VASP displaces atoms one at a time in each of the three Cartesian directions and computes the corresponding change of the Kohn-Sham potential. The results are written to the file phelel_params.hdf5. Notice that we can run this calculation using the gamma-only version of VASP since we have chosen only a single $\mathbf{k}$-point in the KPOINTS file.
Setting ELPH_POT_GENERATE also generates the CONTCAR_ELPH output file. This file is similar to the regular CONTCAR file but instead contains the primitive-cell structure. Open the POSCAR and CONTCAR_ELPH files side by side and compare them:
vi -O POSCAR CONTCAR_ELPH
Can you spot the difference in the lattice vectors and how they relate to the primitive cell and supercell?
You should use the CONTCAR_ELPH as the input POSCAR for the following electron-phonon calculation. This is going to be the next step in this exercise. If you do not like the automatic choice of primitive cell, you can supply your own choice via the ELPH_POT_LATTICE tag.
VASP writes the electron-phonon potential from the supercell calculation to the phelel_params.hdf5 file on a real-space grid for the unit-cell calculation. The electron-phonon potential is calculated on the fine FFT grid (supercell) but is output as an FFT grid (unit cell), as defined by ELPH_POT_FFT_MESH. By default, this ELPH_POT_FFT_MESH grid has to have the exact same dimensions as the FFT grid (unit cell) that is used during the electron-phonon calculation.
How can you influence the size of the FFT grid that VASP uses for the unit-cell calculation?
The good news is that if you do not change ENCUT and PREC between the supercell and the unit-cell calculations (recommended!), then the two grids will match automatically. No further input is required.
However, in some cases, the default FFT grid (unit cell) for the electron-phonon potential from the fine FFT grid (supercell) does not give electron-phonon properties (e.g., bandgap renormalization) that are sufficiently converged. Changing ENCUT to define a new FFT grid (unit cell) is overly complicated and we save this discussion for the end of 2.3.2. Instead, a new FFT grid (unit cell) can then be defined using ELPH_POT_FFT_MESH during a new supercell calculation, i.e., you need to repeat the supercell calculation.
1.3 Bandgap renormalization¶
1.3.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron_phonon/e01_diamond/el-ph.
Make sure that you copy the required files from the supercell directory of the previous step into the el-ph directory as explained below!
system = diamond
ismear = 0; sigma = 0.1; efermi = midgap
prec = accurate
ediff = 1e-7
encut = 400
#run electron-phonon calculation
elph_mode = renorm
elph_selfen_temps = 0 100 200 300 400 500 600
Automatically generated mesh
0
Gamma
4 4 4
Automatically generated mesh
0
Gamma
10 10 10
PAW_PBE C_s 06Sep2000
The INCAR file is largely the same as for the supercell calculation.
Notice that the ENCUT and PREC parameters are exactly the same.
However, this time, we request an electron-phonon calculation to compute the band-structure renormalization by setting ELPH_MODE equal to renorm.
This tag sets reasonable defaults for other tags related to electron-phonon interactions, in this case, renormalization.
In particular,
- ELPH_RUN to start an electron-phonon calculation
- ELPH_SELFEN_GAPS to compute the renormalization automatically for the band gap.
- ELPH_SELFEN_FAN and ELPH_SELFEN_DW to include both the FM and the DW contribution in the self-energy.
In addition, we specify a list of temperatures in K via ELPH_SELFEN_TEMPS. VASP will calculate the bandgap renormalization at all temperatures in a single run and report the results.
The KPOINTS file contains a few $\mathbf{k}$-points in order to perform a typical self-consistent electronic minimization. A different (denser) $\mathbf{k}$-point mesh is used for the computation of electron-phonon physics, which is specified in the KPOINTS_ELPH file. Quantities such as the electron-phonon matrix elements, electronic energies, and phonon frequencies are computed on this denser electron-phonon $\mathbf{k}$-point mesh during the calculation of the self-energy.
1.3.2 Calculation¶
First, open a terminal and navigate to this example's directory (unit cell):
cd $TUTORIALS/electron-phonon/e01_*/el-ph
We need to copy the phelel_params.hdf5 file from the previous step (supercell directory) into this directory so that VASP can read it.
Similarly, we also need to copy the CONTCAR_ELPH file and rename it to POSCAR.
Then we have all the input files required to run VASP.
To copy the files and run VASP, you can type the following:
cp ../supercell/phelel_params.hdf5 .
cp ../supercell/CONTCAR_ELPH POSCAR
mpirun -np 4 vasp_std
VASP reports the output of the electron-phonon calculation in three places:
- The standard output. Here you will see a bare minimum summary of the calculation, including the temperature-dependent bandgap renormalization for the direct and indirect (fundamental) band gaps.
- Towards the end of the OUTCAR file. This output is much more verbose than the standard output, containing a lot of additional information.
- The phelel_params.hdf5 file. It contains machine-readable binary data that can be read by post-processing tools. Here, we use py4vasp!
Inspect the numbers for the bandgap renormalization in the OUTCAR using:
grep -A 8 "Direct gap" $TUTORIALS/electron-phonon/e01_*/el-ph/OUTCAR
grep -A 8 "Fundamental gap" $TUTORIALS/electron-phonon/e01_*/el-ph/OUTCAR
%%bash
grep -A 8 "Direct gap" e01_*/el-ph/OUTCAR
%%bash
grep -A 8 "Fundamental gap" e01_*/el-ph/OUTCAR
Alternatively, you can print the bandgap out using py4vasp:
from py4vasp import Calculation
calc = Calculation.from_path("./e01_diamond/el-ph")
bandgap = calc.electron_phonon.bandgap[0]
bandgap
calc.electron_phonon.bandgap[0] is an array because in a single execution of VASP, you can compute the bandgap using different convergence settings, e.g., the number of bands in Example 2.3 and the smearing parameter in Example 2.4.
Lattice vibrations exist even at zero temperature due to quantum fluctuations in the ground state. This is known as zero-point motion and leads to the zero-point renormalization (ZPR). This means that even for calculations in the ground-state, electron-phonon physics are important to get the correct size of the bandgap.
The bandgap renormalization has a negative sign, which means that the quasi-particle gap (renormalized gap) is smaller than the Kohn-Sham (KS) gap. In other words, the ZPR shrinks the band gap, which is often the case.
Now, plot the renormalized direct gap of diamond as a function of temperature using py4vasp:
from py4vasp import Calculation
calc = Calculation.from_path("./e01_diamond/el-ph")
bandgap = calc.electron_phonon.bandgap[0]
bandgap.plot("direct_renorm")
A higher temperature means that there are more phonon states available for scattering processes. In this case, the electronic states are renormalized by these additional scattering processes such that the band gap shrinks.
1.4 Questions¶
- Why is there a renormalization even at zero temperature?
- What does CONTCAR_ELPH contain?
By the end of this tutorial, you will be able to do the following:
- perform convergence studies with respect to the number of $\mathbf{q}$-points, $N_q$, the number of intermediate states, $N_b$, and the smearing parameter, $\delta$
- understand each type of convergence parameter and how they influence the self-energy
2.1 Task¶
Perform convergence studies of the electron-phonon self-energy with respect to the number of q-points, $N_q$, the number of intermediate states, $N_b$, and the smearing parameter, $\delta$, for cd C
For this exercise, we will look only at the zero-point renormalization (ZPR). Since the temperature dependence enters the self-energy $\Sigma$ only via the fermionic and bosonic distribution functions, a convergence study of the ZPR is usually sufficient to obtain good values also at finite temperature. There are three important convergence parameters that appear in the electron self-energy. In this exercise, you will learn how to control these parameters to obtain good numerical accuracy for high-quality calculations.
Let us look again at the Fan-Migdal (FM) self-energy. At $T=0$, it takes a simple form: $$\tag{2.1} \Sigma^{\text{FM}}_{n\mathbf{k}}(T=0) = \frac{1}{N_q} \sum_{q}^{N_q} \sum_{m}^{N_b} \sum_{\nu} \frac{|g_{mn\nu}(\mathbf{k}, \mathbf{q})|^2}{\varepsilon_{n\mathbf{k}} - \varepsilon_{m\mathbf{k}+\mathbf{q}} \pm \omega_{\nu\mathbf{q}} + \mathrm{i}\delta}. $$ The first convergence parameter is the number of $\mathbf{q}$-points, $N_q$, or more precisely, the $\mathbf{q}$-point mesh that is used to integrate this self-energy. You will perform electron-phonon calculations on 4x4x4, 6x6x6, 8x8x8 and 10x10x10 $\mathbf{q}$-point meshes to see how the ZPR changes with a growing number of $\mathbf{q}$-points.
The second convergence parameter is the number of intermediate states, $N_b$, included in the sum over states. This sum converges slowly as scattering processes with high-energy states become less and less physically relevant. You will compute the ZPR as a function of the number of intermediate states that are included in the self-energy.
The third convergence parameter is the magnitude of the imaginary shift in the denominator, $\delta$. This parameter is a remnant of many-body perturbation theory where it ensures the correct pole structure of the propagators. However, in our numerical context, we can think of it as a smearing parameter that helps us deal with (near) degenerate states during the integration. In general, it should be chosen as small as possible because large values introduce unphysical contributions to the self-energy. Note, however, that too small values might introduce numerical noise due to the finite precision of floating-point operations on common hardware. The default value of $0.01\ \text{eV}$ is usually a good compromise and starting point. You will compute the ZPR as a function of the $\delta$ parameter to see how different values can affect the result.
2.2 Convergence with respect to q-point mesh¶
2.2.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron_phonon/e02_diamond_convergence/q-points.
Each subdirectory specifies a different $\mathbf{q}$-point mesh via the KPOINTS_ELPH file.
Check them out!
diamond
1.0000000000000000
1.7809307798000000 -1.7809307798000000 0.0000000000000000
1.7809307798000000 0.0000000000000000 1.7809307798000000
0.0000000000000000 -1.7809307798000000 1.7809307798000000
C
2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.7500000000000000 0.7500000000000000 0.7500000000000000
system = diamond
ismear = 0; sigma = 0.1; efermi = midgap
prec = accurate
ediff = 1e-7
encut = 400
# run electron-phonon calculation
elph_mode = renorm
elph_selfen_temps = 0
Automatically generated mesh
0
Gamma
4 4 4
PAW_PBE C_s 06Sep2000
Click to see q-mesh-4x4x4/KPOINTS_ELPH!
Automatically generated mesh
0
Gamma
4 4 4
Click to see q-mesh-6x6x6/KPOINTS_ELPH!
K points
0
Gamma
6 6 6
Click to see q-mesh-8x8x8/KPOINTS_ELPH!
K points
0
Gamma
8 8 8
Click to see q-mesh-10x10x10/KPOINTS_ELPH!
K points
0
Gamma
10 10 10
The POSCAR and phelel_params.hdf5 are already provided, so there is no need to run another supercell calculation. Notice how we explicitly select only the temperature 0 K in the INCAR file via ELPH_SELFEN_TEMPS.
2.2.2 Calculation¶
In order to check whether the calculation is converged with respect to the number of $\mathbf{q}$-points, perform electron-phonon calculations at different $\mathbf{q}$-point densities. Remember that the number of k-points used for electron-phonon calculations is specified in the KPOINTS_ELPH. This also determines the number of $\mathbf{q}$-points used in the integration of the self-energy.
Open a terminal and navigate to the $\mathbf{q}$-points directory of this exercise:
cd $TUTORIALS/electron-phonon/e02_*/q-points
Next, run an electron-phonon calculation for each $\mathbf{q}$-point mesh by looping over the directories like so:
for q in 4 6 8 10; do
cd "${q}x${q}x${q}"
mpirun -np 4 vasp_std
cd ../
done
Congratulations on completing your first convergence study!
Let's take a look at the convergence with respect to the $\mathbf{q}$-point mesh:
import plotly.graph_objects as go
from py4vasp import Calculation
zpr_values = []
for path in ["4x4x4", "6x6x6", "8x8x8", "10x10x10"]:
calc = Calculation.from_path(f"./e02_diamond_convergence/q-points/{path}/")
gap_dict = calc.electron_phonon.bandgap[0].to_dict()
gap = gap_dict["fundamental"][0]
renorm_gap = gap_dict["fundamental_renorm"][0][0] # first temperature T = 0K
zpr_values.append((renorm_gap - gap) * 1000) # convert to meV
fig = go.Figure(go.Scatter(
x=[4, 6, 8, 10],
y=zpr_values,
mode='lines+markers'
))
fig.update_xaxes(
tickvals=[4, 6, 8, 10],
ticktext=["4x4x4", "6x6x6", "8x8x8", "10x10x10"],
)
fig.update_layout(
xaxis_title="q-points",
yaxis_title="ZPR (meV)",
width=640,
height=480,
)
fig.show()
In this particular case, we are seemingly able to converge the ZPR with only a 10x10x10 $\mathbf{q}$-point mesh. However, this is not always the case. In fact, it is often necessary to go to much higher $\mathbf{q}$-point densities. You can try to plot the direct gap instead of the fundamental one, and then you will see that the 10x10x10 mesh is not enough to converge it. In practice, it is always good to double check these convergence studies since every material and every band gap is different.
2.3 Convergence with respect to number of states¶
2.3.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron_phonon/e02_diamond_convergence/bands.
Check them out below!
diamond
1.0000000000000000
1.7809307798000000 -1.7809307798000000 0.0000000000000000
1.7809307798000000 0.0000000000000000 1.7809307798000000
0.0000000000000000 -1.7809307798000000 1.7809307798000000
C
2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.7500000000000000 0.7500000000000000 0.7500000000000000
system = diamond
ismear = 0; sigma = 0.1; efermi = midgap
prec = accurate
ediff = 1e-7
encut = 400
#run electron-phonon calculation
elph_mode = renorm
elph_selfen_temps = 0
elph_nbands_sum = 20 40 60 80 100 120 140 160 180 200
Automatically generated mesh
0
Gamma
4 4 4
Automatically generated mesh
0
Gamma
10 10 10
PAW_PBE C_s 06Sep2000
Just like in the previous step, we provide a phelel_params.hdf5 file, to skip the supercell calculation. However, this time, it is not necessary to run separate VASP calculations. Notice the INCAR tag ELPH_NBANDS_SUM, which specifies the maximum number of bands, $N_b$, as a list. For each such number, the self-energy will be evaluated, all in a single VASP calculation.
2.3.2 Calculation¶
Go to your terminal and navigate to the bands folder for this example.
Once there, execute VASP to perform an electron-phonon calculation:
cd $TUTORIALS/electron-phonon/e02_*/bands
mpirun -np 4 vasp_std
from py4vasp import Calculation
import plotly.graph_objects as go
calc = Calculation.from_path("./e02_diamond_convergence/bands/")
band_counts = []
zpr_values = []
for bandgap in calc.electron_phonon.bandgap:
gap_dict = bandgap.to_dict()
# Extracting the number of bands and ZPR values
nbands_sum = gap_dict['metadata']['nbands_sum'][()]
gap = gap_dict['direct'][0]
renormalized_gap = gap_dict['direct_renorm'][0, 0]
# Append the values to the lists
band_counts.append(nbands_sum)
zpr_values.append((renormalized_gap - gap) * 1000) # Convert to meV
fig = go.Figure(go.Scatter(
x=band_counts,
y=zpr_values,
mode='lines+markers'
))
# Update x-axis with band counts
fig.update_xaxes(tickvals=band_counts, ticktext=[str(int(x)) for x in band_counts])
fig.update_layout(
xaxis_title="N<sub>b</sub>",
yaxis_title="ZPR (meV)",
width=640,
height=480,
)
fig.show()
Notice how even at 200 bands, this system is barely converged. One issue with this type of convergence is the non-monotonic behavior that makes it difficult to judge when the system is truly converged. Depending on the system and the used pseudopotential (POTCAR), convergence could require a few hundreds to well over a thousand intermediate states. It is therefore very important to increase the number of bands to check if convergence is reached.
When the meta tag ELPH_MODE is set to renorm, VASP already uses all available bands during the electron-phonon calculation by setting ELPH_NBANDS to -2.
This sets the number of bands equal to the number of plane-wave coefficients in the calculation.
What could you do to increase the number of available bands further?
What would that mean for your FFT meshes?
Would you have to re-run your supercell calculations?
Click to see the answer!
Yes, you would have to re-run the supercell calculation.
Since the number of bands in the unit cell calculation is limited by the number of plane waves, the only way to increase it is by raising the plane-wave cutoff via ENCUT. However, this means that the FFT grid dimensions (unit cell) will likely change, and they will no longer match the grid in the phelel_params.hdf5 file from the supercell calculation. In this case, you must re-run your supercell calculation with the higher ENCUT. This will produce a finer real-space grid in the phelel_params.hdf5 for the unit-cell calculation to match the FFT grid (unit cell). You should keep the same ENCUT and PREC between the supercell and unit-cell calculations. You will then be able to increase the number of bands until convergence is reached.
Alternatively, you can also control the FFT grid (unit cell) via NGX, NGY, and NGZ.
In general, it is therefore recommended to start with a higher cutoff during the initial supercell calculation to be able to properly converge the unit cell with respect to the number of bands later on. This avoids having to re-run the expensive supercell calculation multiple times.
2.4 Convergence with respect to the smearing parameter¶
2.4.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron_phonon/e02_diamond_convergence/delta.
Check them out below!
diamond
1.0000000000000000
1.7809307798000000 -1.7809307798000000 0.0000000000000000
1.7809307798000000 0.0000000000000000 1.7809307798000000
0.0000000000000000 -1.7809307798000000 1.7809307798000000
C
2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.7500000000000000 0.7500000000000000 0.7500000000000000
system = diamond
ismear = 0; sigma = 0.1; efermi = midgap
prec = accurate
ediff = 1e-7
encut = 400
#run electron-phonon calculation
elph_mode = renorm
elph_selfen_temps = 0
elph_selfen_delta = 0.0001 0.001 0.01 0.1 1
Automatically generated mesh
0
Gamma
6 6 6
Automatically generated mesh
0
Gamma
10 10 10
PAW_PBE C_s 06Sep2000
Just like for the case of band-convergence, VASP provides a convenient method for specifying the smearing parameter as a list via ELPH_SELFEN_DELTA.
2.4.2 Calculation¶
The workflow is exactly the same as in the previous step. Navigate to the example folder and execute a single VASP electron-phonon calculation:
cd $TUTORIALS/electron-phonon/e02_*/delta
mpirun -np 4 vasp_std
from py4vasp import Calculation
import plotly.graph_objects as go
calc = Calculation.from_path("./e02_diamond_convergence/delta/")
delta_values = []
zpr_values = []
for bandgap in calc.electron_phonon.bandgap:
gap_dict = bandgap.to_dict()
# Extracting the number of bands and ZPR values
delta = gap_dict['metadata']['selfen_delta']
gap = gap_dict['direct'][0]
renormalized_gap = gap_dict[f'direct_renorm'][0, 0]
# Append the values to the lists
delta_values.append(delta)
zpr_values.append((renormalized_gap - gap) * 1000) # Convert to meV
fig = go.Figure(go.Scatter(
x=delta_values,
y=zpr_values,
mode='lines+markers'
))
# Update x-axis with band counts
fig.update_xaxes(tickvals=delta_values, ticktext=delta_values)
fig.update_layout(
xaxis_title="δ",
yaxis_title="ZPR (meV)",
width=640,
height=480,
xaxis_type="log"
)
fig.show()
Judging from the result, can you tell if the default value of $0.01\ \text{eV}$ for the smearing parameter is a good choice for this system? Take another look at the formula for the zero-point renormalization. What other quantities or parameters could influence the best choice for the smearing parameter?
Lowering the smearing parameter from $0.01\ \text{eV}$ down to $0.1\ \text{meV}$ did not lead to a significant change in the renormalization. Therefore, you can conclude that this default value is sufficient for diamond.
However, for more complicated systems, the delta parameter must be chosen more carefully.
In general, the smearing parameter should be chosen in accordance with the $\mathbf{q}$-point mesh used for integrating the self-energy.
A denser $\mathbf{q}$-mesh introduces more states that are close in energy, which leads to small energies in the denominator of the self-energy.
We recommend conducting a double-convergence study where you change both the $\mathbf{q}$-mesh and smearing parameter at the same time. For each $\mathbf{q}$-point mesh (as in Example 2.2), you would run the ELPH_SELFEN_DELTA with an array of values (as in Example 2.4). In practice, this means adding e.g., ELPH_SELFEN_DELTA = 0.0001 0.001 0.01 0.1 1 to the INCAR file in $TUTORIALS/electron_phonon/e02_diamond_convergence/q-points.
This is particularly important for polar materials, where the electron-phonon matrix elements corresponding to the longitudinal optical modes diverge in the long-wavelength limit due to long-range electrostatic interactions [J. Chem. Phys. 143 (2015) 102813]. In this case, the required $\mathbf{q}$-meshes can become very dense and the smearing parameter should be chosen accordingly.
2.5 Questions¶
- What are $\mathbf{q}$-points?
- What should you do if you don't have enough bands in your electron-phonon calculation to achieve convergence?
- Why does the smearing parameter affect the zero-point renormalization?
By the end of this tutorial, you will be able to do the following for a polar compound:
- Compute the Born effective charges and dielectric tensor
- Compute the electron-phonon potential
- Compute the phonon-induced bandgap renormalization
- Determine the zero-point renormalization via extrapolation
3.1 Task¶
Calculate the Born effective charges and dielectric tensor for a primitive MgO unit cell, followed by the electron-phonon potential for a 2x2x2 MgO supercell, and the zero-point bandgap renormalization for a 2x2x2 MgO supercell, extrapolating to converge the q-point mesh.
MgO is a polar material, which makes the computation of the zero-point renormalization (ZPR) more challenging. In polar materials, long-range electrostatic interactions cause the electron-phonon matrix elements to diverge for certain phonon modes. This is similar to the phenomenon of LO-TO (longitudinal-optical - transverse-optical) splitting for phonons in polar materials. This means you need much denser $\mathbf{q}$-point meshes and careful convergence tests to accurately capture the ZPR. Just like for phonons, a special correction scheme exists that incorporates these long-range interactions into the electron-phonon matrix elements computed by VASP.
The long-range (dipole only) part of the electron-phonon matrix element for a polar material is given by: $$\tag{3.1} g_{mn\nu}^{\text{LR}}(\mathbf{k}, \mathbf{q}) = i \frac{4 \pi}{V_{\text{pc}}} \sum_{\kappa} \sqrt{\frac{\hbar}{2 m_\kappa \omega_{\nu\mathbf{q}}}} \sum_{\mathbf{G}} \frac{ (\mathbf{q} + \mathbf{G}) \cdot \mathbf{Z}^\star_\kappa \cdot \mathbf{e}_{\kappa,\nu \mathbf{q}} }{ (\mathbf{q} + \mathbf{G}) \cdot \boldsymbol{\epsilon}_\infty \cdot (\mathbf{q} + \mathbf{G}) } \langle \psi_{m\mathbf{k}+\mathbf{q}} | e^{i(\mathbf{q} + \mathbf{G}) \cdot (\mathbf{r} + \mathbf{\tau}_\kappa)} | \psi_{n\mathbf{k}} \rangle. $$ where $\mathbf{Z}^*_\kappa$ is the Born effective charge tensor, $\boldsymbol{\epsilon}_\infty$ is the high-frequency dielectric tensor, $m_\kappa$ is the mass of atom $\kappa$, $\omega_{\nu\mathbf{q}}$ is the phonon frequency, and $\mathbf{e}_{\kappa\alpha,\nu}(\mathbf{q})$ is the phonon eigenvector.
The phonon frequencies and eigenvectors, as well as the electronic Bloch orbitals can all be computed by VASP during the electron-phonon calculation. However, the Born effective charges and the dielectric tensor need to be calculated before and given as input.
In this exercise, you will
- Compute the Born effective charges and the dielectric tensor from a linear-response calculation in the unit cell.
- Compute the electron-phonon potential from a supercell while supplying the dielectric properties.
- Compute the bandgap ZPR using dipole corrections and study the q-point convergence.
3.2 Computing dielectric properties¶
3.2.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron_phonon/e03_MgO/nac.
Check them out!
Note: The term "NAC" means non-analytic-term correction and is used in the treatment of the non-analytic behavior of the dynamical matrix at the Brillouin-zone center [Phys. Rev. B 55 (1997) 10355]. Since the Born effective charges and dielectric tensor are used for NAC and a similar concept exists for electron-phonon interactions, we make use of this term here.
MgO
1.0000000000000000
0.0000000000000000 2.1229207162079331 2.1229207162079331
2.1229207162079331 0.0000000000000000 2.1229207162079331
2.1229207162079331 2.1229207162079331 0.0000000000000000
Mg O
1 1
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.4999999999999999 0.4999999999999999 0.4999999999999999
# dielectric INCAR
PREC = Accurate
EDIFF = 1.E-08
ISMEAR = 0
SIGMA = 0.05
EFERMI = MIDGAP
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
ENCUT = 400
# Compute dielectric tensor and Born effective charges
LEPSILON = .TRUE.
Automatically generated mesh
0
Gamma
6 6 6
0 0 0
PAW_PBE Mg 13Apr2007
PAW_PBE O_s 07Sep2000
Setting the INCAR tag LEPSILON instructs VASP to compute the response due to an electric field using density-functional perturbation theory. This includes the Born effective charge tensor and the static ion-clamped dielectric tensor. You can learn more about this type of calculation on the VASP wiki.
3.2.2 Calculation¶
To run this example, open a terminal, navigate to the nac directory and execute VASP:
cd $TUTORIALS/electron-phonon/e03_*/nac
mpirun -np 4 vasp_std
After the calculation has finished, inspect the OUTCAR file and search for the information towards the end. What you are looking for begins with these lines:
MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in DFT)
and
BORN EFFECTIVE CHARGES (including local field effects) (in |e|, cummulative output)
You can find them using:
grep -A 6 "MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in DFT)" OUTCAR
and
grep -A 10 "BORN EFFECTIVE CHARGES (including local field effects) (in |e|, cummulative output)" OUTCAR
%%bash
grep -A 6 "MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in DFT)" e03_*/nac/OUTCAR
grep -A 10 "BORN EFFECTIVE CHARGES (including local field effects) (in |e|, cummulative output)" e03_*/nac/OUTCAR
Alternatively, the output can be obtained from the vaspout.h5 file using py4vasp:
from py4vasp import Calculation
calc = Calculation.from_path("./e03_MgO/nac")
calc.dielectric_tensor.print()
calc.born_effective_charge.print()
The .print() function is a convenient and reliable way to visualize the data. When scripts and workflows are involved, a much more convenient way to handle the data is to use the .read() function, which you will use in 3.3.1.
In the next step of this exercise, you will need to specify these tensors as input in the INCAR file.
3.3 Computing the electron-phonon potential¶
3.3.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron_phonon/e03_MgO/supercell.
The supercell calculation is similar to the simpler case of diamond, which was covered in previous exercises.
For the polar material MgO, however, you also need to specify the Born effective charges and the dielectric tensor of MgO in the INCAR file.
To help you with the setup, the provided INCAR file already has the correct tags and the format.
Check it out and then fill in the missing information X X X from the previous step in the PHON_DIELECTRIC and the PHON_BORN_CHARGES tags.
You can do that with the script below:
import py4vasp
import sys
from py4vasp import Calculation
def format_val(v):
return f"{v:14.8f}"
calc = Calculation.from_path("./e03_MgO/nac")
born_charges = calc.born_effective_charge.read()['charge_tensors']
dielectric_tensor = calc.dielectric_tensor.read()['clamped_ion']
dielec_str = "PHON_DIELECTRIC = \\\n" + " ".join(map(format_val, dielectric_tensor.flatten()))
born_str = "PHON_BORN_CHARGES = \\\n" + " \\\n".join([" ".join(map(format_val, b.flatten())) for b in born_charges])
print(dielec_str)
print(born_str)
with open("e03_MgO/supercell/INCAR", "a") as incar:
incar.write(f"\n{dielec_str}\n{born_str}\n")
Mg8 O8
1.0
0.0000000000000000 4.2458414324158662 4.2458414324158662
4.2458414324158662 0.0000000000000000 4.2458414324158662
4.2458414324158662 4.2458414324158662 0.0000000000000000
Mg O
8 8
direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.0000000000000000 0.5000000000000000 Mg
0.0000000000000000 0.5000000000000000 0.0000000000000000 Mg
0.0000000000000000 0.5000000000000000 0.5000000000000000 Mg
0.5000000000000000 0.0000000000000000 0.0000000000000000 Mg
0.5000000000000000 0.0000000000000000 0.5000000000000000 Mg
0.5000000000000000 0.5000000000000000 0.0000000000000000 Mg
0.5000000000000000 0.5000000000000000 0.5000000000000000 Mg
0.2500000000000000 0.2500000000000000 0.2500000000000000 O
0.2500000000000000 0.2500000000000000 0.7500000000000000 O
0.2500000000000000 0.7500000000000000 0.2500000000000000 O
0.2500000000000000 0.7500000000000000 0.7500000000000000 O
0.7500000000000000 0.2500000000000000 0.2500000000000000 O
0.7500000000000000 0.2499999999999999 0.7500000000000000 O
0.7500000000000000 0.7500000000000000 0.2499999999999999 O
0.7500000000000001 0.7500000000000000 0.7500000000000000 O
# INCAR file for MgO supercell calculations
system = MgO
ismear = 0; sigma = 0.01; efermi = midgap
prec = accurate
ediff = 1e-7
encut = 400
ibrion = 6
potim = 0.015
elph_pot_generate = True
# fill with dielectric tensor
# phon_dielectric = \
# X X X \
# X X X \
# X X X
# fill with Born effective charges
# phon_born_charges = \
# X X X \
# X X X \
# X X X \
# X X X \
# X X X \
# X X X
k-mesh
0
Gamma
2 2 2
PAW_PBE Mg 13Apr2007
PAW_PBE O_s 07Sep2000
This time, we use a 2x2x2 k-point mesh to help with convergence in such a small supercell. Thus, make sure that you use the standard version of VASP for this calculation! In practice, you would want to increase the supercell size until you reach convergence.
3.3.2 Calculation¶
Open a terminal, navigate to this example's directory (supercell) and run VASP by entering the following:
cd $TUTORIALS/electron-phonon/e03_*/supercell
mpirun -np 4 vasp_std
Once again, all relevant quantities are written to the binary phelel_params.hdf5 file, which is then used in the electron-phonon calculation. In particular, phelel_params.hdf5 now also contains the Born effective charges and the dielectric tensor.
3.4 Bandgap renormalization¶
3.4.1 Input¶
The input files to run this example are prepared at $TUTORIALS/electron_phonon/e03_MgO/el-ph.
We will compute the renormalization at different $\mathbf{q}$-point densities, so there is a subdirectory for each $\mathbf{q}$-mesh: 4x4x4, 6x6x6, 8x8x8, 10x10x10.
Since the Born-effective charges and the dielectric tensor of MgO were set already in the supercell calculation, there are no additional inputs required in the INCAR file compared to the diamond example.
First, copy the phelel_params.hdf5 from the previous supercell step to each subdirectory in the el-ph directory.
Similarly, copy the CONTCAR_ELPH and use it as the POSCAR for the electron-phonon calculations.
3.4.2 Calculation¶
Open a terminal, navigate to the el-ph directory and run an electron-phonon calculation with VASP in each subdirectory.
Do not forget to copy the required input files from the supercell calculation (phelel_params.hdf5 and CONTCAR_ELPH)!
You can do everything in one simple script by looping over the directories:
cd $TUTORIALS/electron-phonon/e03_*/el-ph
for q in 4 6 8 10; do
cd "${q}x${q}x${q}"
cp ../../supercell/phelel_params.hdf5 .
cp ../../supercell/CONTCAR_ELPH POSCAR
mpirun -np 4 vasp_std
cd ../
done
Let us once again inspect the OUTCAR files from the four different calculations. MgO only has a direct bandgap so there is no difference in the output between the "direct" and the "fundamental" gap. Therefore, simply search for the direct gap in the output:
grep -A 2 "Direct gap" */OUT*
The quantity of interest is the ZPR, which is written as KS-QP gap (meV) but is the quasiparticle (renormalized) gap minus the Kohn-Sham gap.
You may already be able to tell that this system does not seem to converge well with the number of $\mathbf{q}$-points.
To see it more clearly, plot the ZPR as a function of the $\mathbf{q}$-point density using py4vasp:
import plotly.graph_objects as go
from py4vasp import Calculation
zpr_values = []
for path in ["4x4x4","6x6x6","8x8x8","10x10x10"]:
calc = Calculation.from_path(f"./e03_MgO/el-ph/{path}/")
gap_dict = calc.electron_phonon.bandgap[0].to_dict()
gap = gap_dict["fundamental"][0]
renorm_gap = gap_dict["fundamental_renorm"][0][0] # first temperature T = 0K
zpr_values.append((renorm_gap - gap) * 1000) # convert to meV
fig = go.Figure(go.Scatter(
x=[4, 6, 8, 10],
y=zpr_values,
mode='lines+markers'
))
fig.update_xaxes(
tickvals=[4, 6, 8, 10],
ticktext=["4x4x4", "6x6x6", "8x8x8", "10x10x10"],
)
fig.update_layout(
xaxis_title="q-points",
yaxis_title="ZPR (meV)",
width=640,
height=480,
)
fig.show()
From this plot, it is clear that the ZPR is not yet converged. However, even if you increase the q-point density (maybe give it a try), the ZPR will not converge. This is because the electron-phonon matrix elements associated with the longitudinal optical phonon modes diverge at long wavelengths (as $\mathbf{q} \to \mathbf{0}$). Fortunately, this divergence is integrable, so the (non-adiabatic) self-energy is still finite and the ZPR is well defined.
Given all this, how can you still obtain the ZPR from these finite $\mathbf{q}$-point grids? Hint: The answer lies in the long-range behavior of the electron-phonon matrix, i.e., when $\mathbf{q} \to \mathbf{0}$.
Click to see the answer!
Let us once again inspect the equation of the long-range electron-phonon matrix element. For $\mathbf{q} \to \mathbf{0}$, the denominator diverges as $$ g_{mn\nu}^{\text{LR}}(\mathbf{k}, \mathbf{q}) \propto \frac{|\mathbf{q}|}{|\mathbf{q}|^2} = \frac{1}{|\mathbf{q}|} $$ Therefore, when we integrate the self-energy for $\mathbf{q} \to \mathbf{0}$ it behaves as $$ \int d^3 q |g_{mn\nu}^{\text{LR}}(\mathbf{k}, \mathbf{q})|^2 \propto \int d^3 q \frac{1}{|\mathbf{q}|^2} $$ This function is integrable around the origin and well-behaved for small $\mathbf{q}$, so the integral in the self-energy can be estimated by an iterative refinement of the q-mesh. The self-energy can then be extrapolated towards the $N_q \to \infty$ limit.
In order to extrapolate the ZPR to the limit of an infinitely dense $\textbf{q}$-point mesh ($N_q \to \infty$), plot the ZPR as a function of the inverse grid-spacing. The mathematically inclined reader might want to prove that this is indeed the right thing to do!
import plotly.graph_objects as go
from py4vasp import Calculation
import numpy as np
q_values = np.array([4, 6, 8, 10])
zpr_values = []
for path in ["4x4x4", "6x6x6", "8x8x8", "10x10x10"]:
calc = Calculation.from_path(f"./e03_MgO/el-ph/{path}/")
gap_dict = calc.electron_phonon.bandgap[0].to_dict()
gap = gap_dict["fundamental"][0]
renorm_gap = gap_dict["fundamental_renorm"][0][0] # first temperature T = 0K
zpr_values.append((renorm_gap - gap) * 1000) # convert to meV
inv_q = 1 / q_values
# linearly extrapolate the last two data points of x and y
coefficients = np.polyfit(inv_q[-2:], zpr_values[-2:], 1)
# get the extrapolated value at x = 0
zpr_at_0 = np.polyval(coefficients, 0)
linex = [inv_q[-1], 0]
liney = [zpr_values[-1], zpr_at_0]
fig = go.Figure(go.Scatter(
x=inv_q,
y=zpr_values,
mode='lines+markers',
name='ZPR samples',
))
fig.add_scatter(
x=linex,
y=liney,
mode='lines',
line=dict(dash='dash', color='gray'),
name='Extrapolation'
)
fig.update_xaxes(
tickvals=[*inv_q, 0],
ticktext=["4x4x4", "6x6x6", "8x8x8", "10x10x10", "Continuum limit"],
)
fig.update_layout(
xaxis_title="q-points",
yaxis_title="ZPR (meV)",
width=640,
height=480,
)
fig.show()
Important: In this tutorial, we demonstrate the above idea of extrapolation in the linear regime of the ZPR. However, in practice you will need much denser $\mathbf{q}$-point grids than the ones we used here to enter the linear regime. Keep in mind that it will likely be necessary to go to 32x32x32 or even higher $\mathbf{q}$-point densities to obtain a reliable extrapolation towards infinity.
3.5 Questions¶
- Why do you need to calculate the Born effective charges and the dielectric tensor?
- Why is a dense $\mathbf{q}$-point mesh not sufficient to converge the ZPR for a polar material?