Electron-phonon coupling

At the heart of the electron-phonon coupling is the calculation of the gradient of the effective potential, which is done using finite displacements just like the phonons. Those two calculations can run simultaneous, if the required set of parameters coincide. The electron-phonon matrix is quite sensitive to self-interaction of a displaced atom with its periodic images, so a sufficiently large supercell needs to be used. The (3x3x3) supercell used in this example is a bit too small. When using a supercell for the calculation you have to consider, that the atoms object needs to contain the primitive cell, not the supercell, while the parameters for the calculator object need to be good for the supercell, not the primitive cell. Convergence parameters and k-points need to be tuned to ensure convergence. (effective_potential.py)

from ase.build import bulk
from gpaw import GPAW, FermiDirac
from gpaw.elph import DisplacementRunner

atoms = bulk('Si', 'diamond', a=5.431)

calc = GPAW(mode='lcao', h=0.18, basis='dzp',
            kpts=(4, 4, 4),
            xc='PBE',
            occupations=FermiDirac(0.01),
            symmetry={'point_group': False},
            convergence={'energy': 2e-5, 'density': 1e-5},
            txt='displacement.txt')

elph = DisplacementRunner(atoms=atoms, calc=calc,
                          supercell=(3, 3, 3), name='elph',
                          calculate_forces=True)
elph.run()

This script executes \(2*3*N\) displacements and saves the change in total energy and effective potential into a file cache in the directory \(elph\). The phonon/effective potential calculation can take quite some time, but can be distributed over several images. The DisplacementRunner() class is based on ASEs ase.phonon.Displacement class, which allows to select atoms to be displaced using the set_atoms function.

Note

If you use the calculate_forces=True option, you can load the phonon part of the cache using Phonons(..., name='elph', center_refcell=True), with center_refcell=True being crucial. Else calculated phonon energies will be wrong, as ase.phonons.Phonons defaults to center_refcell=False.

In the second step we map the gradient of the effective potential to the LCAO orbitals of the supercell (supercell.py).

from ase.build import bulk
from gpaw.elph import Supercell

atoms = bulk('Si', 'diamond', a=5.431)

sc = Supercell(atoms, supercell_name="supercell", supercell=(3, 3, 3))
calcdict = {'kpts': (4, 4, 4), 'basis': 'dzp', 'txt': 'supercell.txt'}
sc.calculate_supercell_matrix(calcdict, fd_name='elph')

This does not require a full GPAW calculation and should not take too long. You need to specify two options as a dictonary, which will be used to create a GPAW object internally: The size of the LCAO basis and the size of the k-point grid. The last line invokes the calculation of the supercell matrix, which is stored in the \(supercell\) file cache. It is also possible to pass a full GPAW calculator instead of a dictonary.

Note

If you give a full GPAW calculator instead of a dictonary to calculate_supercell_matrix you need to consider a few things. calculate_supercell_matrix currently only works for k-point parallelization, so you need to include parallel={'domain': 1, 'band': 1} in your script. Also, the real space grids used in the finite displacement calculations needs to be the same as the one used in the supercell matrix calculation. If you use planewave mode for the finite displacement calculation you should set the required grid manually, for example by adding gpts=(nx, ny, nz) where nx, ny, nz need be substituted with the required number of grid points in each direction. You can use python3 -m gpaw.elph.gpts to get help with this.

After both calculations are finished the final electron-phonon matrix can be constructed. (gmatrix.py)

from gpaw import GPAW
from gpaw.elph import ElectronPhononMatrix

calc = GPAW("scf.gpw")
atoms = calc.atoms

elph = ElectronPhononMatrix(atoms, 'supercell', 'elph')
q = [[0., 0., 0.], [1. / 11., 1. / 11., 1. / 11.]]
g_sqklnn = elph.bloch_matrix(calc, k_qc=q,
                             savetofile=True, prefactor=False)

For this we need to perform another ground state calculation (scf.py) to obtain the wave function used to project the electron-phonon coupling matrix into. The electron-phonon matrix is computed for a list of q values, which need to be commensurate with the k-point mesh chosen. The ElectronPhononMatrix class doesn’t like parallelization so much and should be done separately from the rest.

Exercise

The electron-phonon matrix can be used to calculate acoustic and optical deformation potentials without much effort. The optical deformation potential is identical to the bare electron-phonon coupling matrix Ref. [1]:

\[D_{ODP} = \langle m \vert \nabla_u V_{eff} \cdot \mathbf e_l \vert n \rangle .\]

For the silicon VBM at Gamma, for the LO phonons at Gamma, we expect a deformation potential of about \(\vert M \vert \approx 3.6\) eV/Angstrom. As the optical phonons are degenerate, as well as the three highest valence bands, we need to sum over all 27 contributions here.

Converge the deformation potential of Si with respect to the supercell size, k-points, grid spacing and SCF convergence parameters.

References

[1]

26. Li, P. Graziosi, and N. Neophytou, “Deformation potential extraction and computationally efficient mobility calculations in silicon from first principles”, Physical Review B 104, 195201 (2021)