Part 3: Electron-phonon matrix elements and using VASP with phelel ¶

Content¶

 8 ELPH_DRIVER tag to obtain the electron-phonon matrix elements + visualization
 9 Bandstructure of C diamond
 10 Dielectric properties for polar materials: MgO

8 Obtaining the electron-phonon matrix elements

$\uparrow$

¶

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

• Calculate the electron-phonon matrix elements.
• Plot the matrix elemenets along a $\textbf{q}$-path.
• Relate the resulting plot to the electronic and phonon band structures.

8.1 Task¶

Use VASP to plot averaged electron phonon matrix elements along a Γ-X-U-K-Γ-L-W-X q-path for C diamond (primitive) using a 4x4x4 k-mesh

The interaction between electrons and lattice vibrations (phonons) leads to a modification of the electronic energies, described by the electron self-energy. Using second-order perturbation theory, the electron self-energy can be separated into two contributions: the Fan-Migdal (FM) and Debye-Waller (DW) terms. In practice, these terms are evaluated in terms of the electron-phonon matrix elements $g_{mn\mathbf{k}, \nu \mathbf{q}}$, which can be defined in terms of the Kohn-Sham orbitals $\psi_{n\mathbf{k}}$ and the electron-phonon potential $\partial_{\nu \mathbf{q}} V_{\text{KS}}$:

$$\tag{8.1} g_{mn\mathbf{k}, \nu \mathbf{q}} = \left\langle \psi_{m\mathbf{k}+\mathbf{q}} \left| \partial_{\nu \mathbf{q}} V_{\text{KS}} \right| \psi_{n\mathbf{k}} \right\rangle$$

where $(\nu, \textbf{q})$ are the index $\nu$ and wavevector $\textbf{q}$ of a phonon mode, and $m,n$ are electronic bands with wavevectors $\textbf{k}, \textbf{k} + \textbf{q}$.

These matrix elements are usually computed on-the-fly and discarded, improving the efficiency. However, in some cases, you may want to obtain the electron-phonon matrix elements directly, for example, to analyze or plot them individually.

In this tutorial, you will learn how to obtain and plot the electron-phonon matrix elements.

8.2 Input¶

The input files to run this example are prepared at $TUTORIALS/electron_phonon/e08_matrix_elements.

POSCAR

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

INCAR

system = diamond

ismear = 0; sigma = 0.1; efermi = midgap
prec = accurate
ediff = 1e-7
encut = 400

elph_run = True
elph_driver = mels
elph_selfen_ikpt = 1
elph_selfen_band_start = 1
elph_selfen_band_stop = 4
elph_nbands = 4

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_driver = mels means that the electron-phonon matrix elements are printed to a Vaspelph.h5 file. elph_selfen_ikpt determines which $\textbf{k}$-points the electron self-energy is calculated for, in this case, only the $\Gamma$. elph_selfen_band_start and elph_selfen_band_stop define the initial and final electron bands that are used in the electron-phonon calculation, in this case, 1 2 3 4. Finally, elph_nbands determines the number of bands to compute on the dense $\textbf{k}$-point grid for the electron-phonon driver.

KPOINTS

Automatically generated mesh
       0
Gamma
 4 4 4

KPOINTS_ELPH

kpath for band structure
21
line
reciprocal
    0.00000000     0.00000000     0.00000000  Gamma
    0.50000000     0.00000000     0.50000000  X

    0.50000000     0.00000000     0.50000000  X
    0.62500000     0.25000000     0.62500000  U

    0.37500000     0.37500000     0.75000000  K
    0.00000000     0.00000000     0.00000000  Gamma

    0.00000000     0.00000000     0.00000000  Gamma
    0.50000000     0.50000000     0.50000000  L

    0.50000000     0.50000000     0.50000000  L
    0.50000000     0.25000000     0.75000000  W

    0.50000000     0.25000000     0.75000000  W
    0.50000000     0.00000000     0.50000000  X

POTCAR

PAW_PBE C_s 06Sep2000

We provide a phelel_params.hdf5 containing the electron-phonon potential to skip the supercell calculation, and also a script elph_mels.py to plot the matrix elements. There are also a few new tags: elph_selfen_ikpt, , which computes the electron self-energy for a list of $\textbf{k}$-points specified by their index, and elph_selfen_band_start and elph_selfen_band_stop, which compute the electron self-energy due to electron-phonon coupling for bands between the defined start and stop. ELPH_DRIVER=mels writes the electron-phonon matrix elements to disk.

8.3 The electron-phonon matrix elements¶

The calculation is quick to run, navigate to the exercise directory and then run the calculations with the following:

cd $TUTORIALS/electron-phonon/e08_matrix_elements
mpirun -np 4 vasp_std

When the phonon-induced electron self-energy is computed, the electron-phonon matrix elements are calculated on the fly and discarded afterwards. This allows the implementation of the self-energy to be extremely efficient, computing only the matrix elements that are really required during the calculation. However, in certain cases, you may want to obtain the electron-phonon matrix elements directly, for example, to analyze or plot them individually.

Once the calculation has finished, you can extract the electron-phonon matrix elements using elph_mels.py and plot them in comparison to the electronic and phonon band structure with the following script:

When ELPH_DRIVER tag is equal to mels, the electron-phonon matrix elements are written to disk. In this case, VASP does not compute the self-energy expression but instead accumulates all matrix elements that meet the specified selection rules and writes them to the vaspelph.h5 HDF5 file. The layout of this structured file is detailed on the corresponding VASP wiki page.

Here, we show how it is possible to compute the matrix elements along a path in the first Brillouin zone, extract that information from the file and plot it in a similar fashion as an electronic band structure (1st plot) and phonon band structure (2nd plot). The electronic state $n,\textbf{k}$ at the $\Gamma$-point, and the electronic state at $m,\textbf{k+q}$ along the $\textbf{q}$-path interact via the electron-phonon potential by emitting or adsorbing a phonon with a branch index $\nu$ and a wavevector $\textbf{q}$.

We then plot (3rd plot) the electronic state $\varepsilon_{m,\textbf{k+q}}$ with $\textbf{k} = \Gamma$, the phonon state $\omega_{\nu,\textbf{q}}$, and the matrix element $g_{mn\mathbf{k}, \nu \mathbf{q}}$ along the $\textbf{q}$-path above. The average of $g_{mn\mathbf{k}, \nu \mathbf{q}}$ is plotted at $\mathbf{k} = \Gamma$-point over all degenerate states at band $m=n=4$ (which includes $n=2,3,4$ due to degeneracy), cf. the electronic band structure plot. For more details of the averaging, look in elph_mels.py or see Eq. 97 of Phys. Rev. B 100, 174304 (2019)).

You can follow the path of a phonon index $\nu=6$ from the phonon band structure in the $g_{mn\mathbf{k}, \nu \mathbf{q}}$ plot (4th plot). At a crossing between indices in the phonon band structure, there is a sudden change in the derivative, hence the discontinuities in $g_{mn\mathbf{k}, \nu \mathbf{q}}$, if you keep the same $\nu$. Alternatively, you could follow the physical state throughout the crossing, then you would change the phonon index $\nu$, but there would be no discontinuity.

8.4 Questions¶

1. How do you follow a phonon state through the phonon band structure and electron-phonon matrix elements?
2. Why are there discontinuities if you follow the phonon band index? Are they physical?

9 Bandstructure of diamond

$\uparrow$

¶

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

• Use velph to generate a velph.toml file and prepare directories for band structure and DOS calculations.
• Plot the band structure and DOS for C diamond.

9.1 Task¶

Using velph to set up a VASP calculation for primitive C diamond with a k-spacing of 0.6 and calculate the bands structure and DOS for a 2x2x2 supercell.

Unlike the manual approach in Part 1, where we had preprepared input files in separate directories, the velph workflow is managed through a single configuration file (velph.toml) that contains all the information needed to create the different calculations and directories automatically. The velph.toml file will be generated with the velph init command, which allows you to choose supercell size, $\textbf{k}$-point sampling and electronic cutoff energies that are kept consistent across calculations. In Part 1 of this tutorial, you learned how to perform electron-phonon calculations using VASP via the ELPH_POT_GENERATE tag. While this approach is simple and straightforward, it also has its shortcomings.

The phelel package provides a high-level interface to VASP that automates many of the tedious steps in electron-phonon calculations. The velph command-line tool, part of the phelel ecosystem, can automatically:

• Set up calculations of the Born effective charges and dielectric tensor used for the non-analytic term correction (NAC)
• Generate supercells and compute electron-phonon potentials
• Prepare primitive cell electron-phonon calculations
• Handle file management and maintain consistent settings between calculations

This integrated workflow reduces the chance of errors in setting up calculations and transferring data between steps. In addition, phelel offers greater flexibility in choosing the atomic displacements. It also employs more symmetry operations when constructing the electron-phonon potential, which can lead to reduced run times.

In this tutorial, we will demonstrate the complete VASP+phelel workflow for C, starting from just the primitive cell structure and a POTCAR file.

9.2 Input¶

The input files to run this example are prepared at $TUTORIALS/electron_phonon/e09_phelel_band.

For the velph workflow, we start with just two essential files: a primitive cell structure (POSCAR) and the pseudopotentials (POTCAR). Check them out!

POSCAR

diamond
   1.0000000000000000
     1.7809307797928542   -1.7809307797928542    0.0000000000000000
     1.7809307797928542    0.0000000000000000    1.7809307797928546
     0.0000000000000000   -1.7809307797928542    1.7809307797928535
  C
     2
Direct
     0.0000000000000000    0.0000000000000000    0.0000000000000000
     0.7500000000000002    0.7499999999999998    0.7499999999999999

POTCAR

 PAW_PBE C_s 06Sep2000

template.toml

[vasp.incar]
encut = 400
[vasp.phelel.incar]
elph_prepare = true
lwap = false

The POSCAR file defines the primitive unit cell structure of C diamond. This is the same structure that was used in Part 1, Exercise 2 for the dielectric calculation.

The POTCAR file contains the pseudopotentials for C atoms. This is the same POTCAR used in Part 1.

The template.toml file is used to provide input settings that are otherwise set by default in velph, e.g., ENCUT.

9.3 Generating the workflow configuration¶

Now we use velph to setup the calculational parameters. Take some time to familiarize yourself with the command-line options of velph via the following help commands:

cd $TUTORIALS/electron-phonon/e09_phelel_band
velph --help
velph hints

The first step is to initialize the project and generate the velph.toml configuration file. This is done via velph init, so once again, familiarize yourself with the available command-line options of this subcommand:

velph init --help

In addition to the POSCAR file, velph init also requires you to specify a project directory and the supercell size. For the project directory, we can simply choose the example directory and for the supercell size, we use the --max-num-atoms command-line option. Now, Navigate to the example directory and execute the following:

velph init POSCAR . --max-num-atoms=100 --kspacing=0.6 --force --template-toml template.toml

We set the --kspacing in this example to call KSPACING, which determines the spacing between $\textbf{k}$-points in reciprocal space, generating a KPOINTS file. A setting of --kspacing=0.6 gives the same 2x2x2 $\textbf{k}$-point sampling in the supercell as in the previous MgO example. Otherwise, a default $\textbf{k}$-point mesh will be chosen for all calculations.

velph will now take the input geometry and options to create a velph.toml configuration file. This file contains all input settings for the following steps and calculations. Take some time to inspect the file.

It is possible to adjust all the settings within the velph.toml file. Most of them map one-to-one to INCAR tags. Any subsequent command, for example velph phelel init, will use these settings to generate VASP input files.

[vasp.el_bands.incar]
ibrion = -1
nsw = 0
lorbit = 11
nedos = 5001
ismear = 0
sigma = 0.05
ediff = 1e-08
encut = 400
prec = "accurate"
lreal = false
lwave = false
lcharg = false
[vasp.el_bands.kpoints]
mesh = [5, 5, 5]
[vasp.el_bands.kpoints_opt]
line = 51
[vasp.el_bands.kpoints_dense]
mesh = [61, 61, 61]

Take a look inside with vi velph.toml. You can see that there are many lines with [vasp.*.incar] or [vasp.*.kpoints] that provide the settings for the INCAR or KPOINTS files. At the end of the file, you can also find details of the symmetry [symmetry], unit cell [unitcell.*], and lattice [lattice]. For our current calculations, you can find the details under [vasp.el_bands.incar]. You can generate the files for calculating the band structure with the following command:

cd $TUTORIALS/electron-phonon/e09_*
velph el_bands generate

Now that the calculation is prepared, you can run it with the following:

cd $TUTORIALS/electron-phonon/e09_*/el_bands/bands
mpirun -np 4 vasp_std

Congratulations, you have done your first calculation set up by velph. Take a look at the band structure:

The velph el_bands generate command that you ran earlier produced two sets of files; bands you have already calculated, now try dos to calculate the density of states (DOS):

cd $TUTORIALS/electron-phonon/e09_*/el_bands/dos
mpirun -np 4 vasp_std

Once the calculation is finished, plot the DOS:

and then use velph to plot them side-by-side:

In this exercise, you have learned the basic introduction to phelel and the velph.toml file. You have used the following workflow:

1. Initialize velph.toml configuration file.
2. Generate directories el_bands/bands and el_bands/dos.
3. Run band structure and DOS calculations.

which has allowed you to run two calculations beginning from only POSCAR and POTCAR files.

In the next exercise, you will go through the full VASP+phelel workflow for C diamond to calculate the bandgap renormalization.

9.4 Questions¶

1. Where can you find the INCAR settings for your calculation when using velph?
2. What are the strengths of this structure?

10 Dielectric properties for polar materials: MgO

$\uparrow$

¶

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

• Compute dielectric properties for polar materials using velph.
• Calculate bandgap renormalization using the integrated workflow.
• Compare results between manual and automated approaches.

10.1 Task¶

Use velph to generate a velph.toml file (kspacing=0.6), calculate the Born effective charges and the dielectric tensor, automatically generate the supercell (2x2x2), compute the electron-phonon potential, then perform the final electron-phonon calculation in the primitive cell and compute the bandgap renormalization for MgO over a range of temperatures from 0-700 K.

You have learned how to perform a basic band structure calculation in VASP using phelel. This approach is a simple introduction to how phelel sets up a calculation. Unlike the manual approach in Example 3, where we had ready-prepared input files in separate directories, the velph workflow is managed through a single configuration file (velph.toml) that contains all the information needed to create the different calculations and directories automatically. The velph.toml file will be generated with the velph init command, which allows you to choose supercell size, $\textbf{k}$-point sampling and electronic cutoff energies that are kept consistent across calculations.

In this tutorial, we will demonstrate how to generate the Born effective charges and dielectric tensor for MgO, starting from just the primitive cell structure and a POTCAR file. Using these, we will calculate the zero-point bandgap renormalization. We will compare the results with those obtained in Part 1 to show that both approaches yield identical results.

10.2 Input¶

The input files to run this example are prepared at $TUTORIALS/electron_phonon/e10_phelel_MgO.

For the velph workflow, we start with just two essential files: a primitive cell structure (POSCAR) and the pseudopotentials (POTCAR). Check them out!

POSCAR

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

POTCAR

PAW_PBE Mg 13Apr2007
PAW_PBE O_s 07Sep2000

The POSCAR file defines the primitive unit cell structure of MgO. This is the same structure that was used in Part 1, Exercise 3 for the dielectric calculation.

The POTCAR file contains the pseudopotentials for Mg and O atoms. This is the same POTCAR used in Part 1.

10.3 Generating the workflow configuration¶

We use velph to setup the calculational parameters. The first step is to initialize the project, as before, and generate the velph.toml configuration file, done via velph init. Now, Navigate to the example directory and execute the following:

cd $TUTORIALS/electron-phonon/e10_phelel_MgO
velph init POSCAR . --max-num-atoms=100 --kspacing=0.6 --force --template-toml template.toml

10.4 Computing dielectric properties using velph¶

10.4.1 Understanding the configuration¶

The velph.toml configuration file contains all the predefined parameters to run the VASP calculation.

Inspect the vasp.nac section in the velph.toml configuration file:

cd $TUTORIALS/electron-phonon/e10_phelel_MgO
grep -A 14 "\[vasp.nac\]" velph.toml

In the vasp.nac.incar subsection, you will find some basic convergence tags as well as LEPSILON = True. Additionally, there is vasp.nac.kpoints subsection which defines the $\textbf{k}$-point mesh used in the calculation. Since we chose a very coarse $\textbf{k}$-point mesh for tutorial purposes during the initialization (velph init --kspacing=0.6), the $\textbf{k}$-point mesh is only 4x4x4. In order to get results consistent with the MgO calculation in part 1 of the tutorials, increase the $\textbf{k}$-point mesh to 6x6x6 by editing the corresponding entry in the velph.toml file.

With velph, we eliminate the need to manually create input files for each step and ensure consistent parameters across all calculations. Now, it is time to generate the VASP input files. To do so, run velph nac generate, which will create a new directory called nac containing all required VASP input files:

cd $TUTORIALS/electron-phonon/e10_*
velph nac generate

10.4.2 Calculating the dielectric properties with VASP¶

Navigate to the nac directory and verify that the INCAR file contains the expected tags and that the KPOINTS file has indeed 6x6x6 $\textbf{k}$-points. Afterwards, run VASP to compute the Born effective charges and the dielectric tensor.

cd $TUTORIALS/electron-phonon/e10_*/nac
mpirun -np 4 vasp_std

With this, the NAC part of the electron-phonon workflow is completed.

This produces the usual output as we have already seen in the dielectric calculation in Example 3. Try printing out the dielectric properties using py4vasp:

Let's compare to Example 3:

You can see that they give identical values; the calculations are the same, and the only difference is the workflow used.

velph now automatically knows where to look for the Born effective charges and the dielectric tensor and will include them in the next steps.

10.5 Supercell electron-phonon potential using velph¶

10.5.1 Generating the supercell calculation¶

Let's use velph to automatically set up the supercell calculation. Once again, you can use the velph hints command to learn about this step and how to carry it out. The subcommand corresponding to the supercell calculation is called velph phelel (since it uses phelel to perform the finite differences in the supercell), and you can learn more about it by typing velph phelel --help.

After you have familiarized yourself with the available documentation, it is time to initialize the supercell calculation and generate the VASP input files. Make sure that you are at the example root directory where the velph.toml file is located. Then, start by running velph phelel init. This command finds the required atomic displacements and generates a new directory, phelel, for the supercell calculation. It contains a single configuration file phelel_disp.yml that specifies exactly which atoms are to be displaced and what symmetries are applied.

First initialize the calculation:

cd $TUTORIALS/electron-phonon/e10_*
velph phelel init

and then generate the input files based on this information, run velph phelel generate:

velph phelel generate

Within the phelel directory, velph automatically created a set of subdirectories, each corresponding to a specific atomic displacement in the supercell. Each subdirectory contains all the necessary input files for that particular displacement. In addition, there is also a subdirectory for the unperturbed geometry. This organization allows you to run independent calculations for each displaced structure.

Depending on the symmetry settings that were chosen during the initialization process, the number of subdirectories can be different. Navigate to each of the new subdirectories in the phelel directory and confirm that the input files were correctly generated. The only difference between the input files from different directories should be the POSCAR, containing the atomic displacements. Notice that there is a new INCAR tag in each of the directories: ELPH_PREPARE. This tag will write all the necessary information (such as the local potential) to the vaspout.h5.

In the next step, you will run VASP for each of these displacements and then use phelel to process the results.

10.5.2 Running the supercell calculation with VASP and phelel¶

First, navigate to the phelel directory and then run VASP inside each of the subdirectories:

cd $TUTORIALS/electron-phonon/e10_*/phelel
for d in disp-*
do
   cd $d
   mpirun -np 4 vasp_std
   cd ..
done

VASP+phelel presently offers more flexibility over the choice of displacements than the VASP-only approach.

After all supercell calculations are finished, it is now time to execute phelel to gather the results, perform finite differences and generate the phelel_params.hdf5. However, instead of running phelel manually, you can make use of velph to run it for you. To do this, make once again sure that you are back at the example root directory where the velph.toml file is located. Once there, simply run velph phelel differentiate:

cd $TUTORIALS/electron-phonon/e10_*
velph phelel differentiate

This command will supply phelel with the correct settings and tell it where the VASP calculations are located. It will then perform the finite differences and produce the phelel_params.hdf5 file. After everything is done, it is now time for the final step, the electron-phonon calculation.

10.6 Primitive cell electron-phonon calculation using velph¶

10.6.1 Preparing the electron-phonon calculation¶

Let's use velph to prepare the final electron-phonon calculation. We will compute the zero-point renormalization (ZPR) and temperature dependence. Before you run any commands, take some time to study the velph.toml file once again. This time, focus on the section [vasp.selfenergy.incar]. This section contains the INCAR tags used during the electron-phonon calculation. You can change the settings here before generating the input files.

For this tutorial, we want to compute the bandgap renormalization as a function of temperature at a fixed $\textbf{k}$-point sampling. To do this, change some of the input parameters in the velph.toml file to suit our needs.

First, we want to set the dense $\textbf{k}$-point mesh to 8x8x8. Like almost everything else in this tutorial, this is will not yield converged results, but it will be sufficient to produce reasonable values. Find the [vasp.selfenergy.kpoints_dense] section in the velph.toml file and change the mesh to 8x8x8.

Next, some of the INCAR tags given here are not optimal for our purposes. Got to the [vasp.selfenergy.incar] section and change the list of temperatures (ELPH_SELFEN_TEMPS) to span from 0 to 700 K:

elph_selfen_temps = [0, 100, 200, 300, 400, 500, 600, 700]

Since we are not interested in studying the convergence with respect to the number of intermediate states in this exercise, remove the ELPH_NBANDS_SUM tag. Instead, we want to use all available bands, so set ELPH_NBANDS = -2.

With all settings adjusted, it is time to generate the VASP input files. Make once again sure that you are in the exercise root directory, where the velph.toml is located. Once there, run velph selfenergy generate:

cd $TUTORIALS/electron-phonon/e10_*
velph selfenergy generate

10.6.2 Running the calculation¶

velph has created a new directory called selfenergy in which the input files for the final step were generated. Navigate to the selfenergy directory and confirm for yourself that the INCAR file contains the desired tags. Check also the KPOINTS_ELPH to see if it has the chosen 8x8x8 k-point mesh. Once everything is in place, execute the final calculation to compute the bandgap renormalization:

cd $TUTORIALS/electron-phonon/e10_*/selfenergy
mpirun -np 4 vasp_std

Let's examine the computed bandgap renormalization using py4vasp:

In Example 03, you calculate only the first point on the plot, focusing instead on converging $\textbf{k}$-meshes. Here, you can see that the point at 0 K is identical to your new values. What you can now see is how the ZPR changes with temperature. It decreases with increasing temperature. We hope that this convinces you that the velph procedure sets up the same calculations with an alternative workflow.

The VASP-only workflow will likely suit someone who is familiar with running calculations in VASP. The VASP+phelel workflow will likely suit someone who is familiar with phonopy already and has started using VASP. At present, VASP+phelel offers more flexibility. Whether VASP-only or VASP+phelel suits you better depends largely on your personal preference.

10.7 Questions¶

1. How does the VASP+phelel workflow assist you in calculations?
2. When might you want to use the VASP-only or the VASP+phelel workflow?

Congratulations! You have completed Part 3 of the electron-phonon tutorial!¶

You have successfully:

• Plotted electron-phonon matrix elements and followed a phonon band.
• Set up and executed a complete VASP+phelel workflow using velph.
• Computed Born effective charges and dielectric tensor automatically.
• Generated supercell electron-phonon potentials with minimal manual intervention.
• Calculated bandgap renormalization using the integrated workflow.

Go to Top $\uparrow$