Electron-phonon coupling¶
Electron-phonon coupling is implemented for the LCAO mode.
Introduction¶
The electron-phonon interaction can be defined as
The phonon modes \(l\) are coupled to the electronic states \(i\), \(j\) via the electron-phonon coupling matrix
\(\omega_l\) and \(\mathbf e_l\) are the frequency and mass-scaled polarization vector of the \(l\) th phonon. \(M\) is an effective mass and nabla_u denotes the gradient wrt atomic displacements.
The implementation supports calculations of the el-ph coupling in both finite and periodic systems, i.e. expressed in a basis of molecular orbitals or Bloch states.
The implementation is based on finite-difference calculations of the the atomic gradients of the effective potential expressed on a real-space grid. The el-ph couplings are obtained from LCAO representations of the atomic gradients of the effective potential and the electronic states.
The current implementation supports spin-paired and spin-polarized computations.
A short example is given below. Another worked out example can be found in the tutorial for Raman calculations here.
Example¶
In a typical application one would compute the phonon modes separately as those need very different convergence settings. (phonon.py)
from ase.build import bulk
from ase.phonons import Phonons
from gpaw import GPAW, FermiDirac
a = 3.567
atoms = bulk('C', 'diamond', a=a)
calc = GPAW(mode='lcao',
basis='dzp',
kpts=(5, 5, 5),
xc='PBE',
occupations=FermiDirac(0.01),
symmetry={'point_group': False},
convergence={'energy': 2e-5, 'density': 1e-5},
txt='phonons.txt')
atoms.calc = calc
# Phonon calculator
ph = Phonons(atoms, atoms.calc, supercell=(1, 1, 1), delta=0.01)
ph.run()
The corresponding calculation of the effective potential changes can be done simultaneously. (elph.py)
from ase.build import bulk
from gpaw import GPAW, FermiDirac
from gpaw.elph.electronphonon import ElectronPhononCoupling
a = 3.567
atoms = bulk('C', 'diamond', a=a)
calc = GPAW(mode='lcao',
basis='dzp',
kpts=(5, 5, 5),
xc='PBE',
occupations=FermiDirac(0.01),
symmetry={'point_group': False},
convergence={'bands': 'nao'},
txt='elph.txt')
atoms.calc = calc
atoms.get_potential_energy()
calc.write("scf.gpw", 'all')
elph = ElectronPhononCoupling(atoms, atoms.calc, calculate_forces=False)
elph.set_lcao_calculator(atoms.calc)
elph.run()
elph.calculate_supercell_matrix()
The last line in the above script constructs the electron-phonon matrix in terms of LCAO orbitals (and cell repetitions) and saves it as ASE JSON cache to the \(supercell\) directory.
After both calculations are finished the final electron-phonon matrix can be constructed with a ‘simple’ script. (construct_matrix.py)
import numpy as np
from ase.phonons import Phonons
from gpaw import GPAW
from gpaw.elph.electronphonon import ElectronPhononCoupling
calc = GPAW("scf.gpw")
kpts = calc.get_ibz_k_points()
qpts = [[0, 0, 0], ]
# Phonon calculation, We'll read the forces from the elph.run function
# This only looks at gamma point phonons
ph = Phonons(atoms=calc.atoms, supercell=(1, 1, 1))
ph.read()
frequencies, modes = ph.band_structure(qpts, modes=True)
# Find el-ph matrix in the LCAO basis
elph = ElectronPhononCoupling(calc.atoms, calc=None, calculate_forces=False)
elph.set_lcao_calculator(calc)
elph.load_supercell_matrix()
# Find the bloch expansion coefficients
g_sqklnn = []
for s in range(calc.wfs.nspins):
c_kn = []
for k in range(calc.wfs.kd.nibzkpts):
C_nM = calc.wfs.collect_array('C_nM', k, s)
c_kn.append(C_nM)
c_kn = np.array(c_kn)
# And we finally find the electron-phonon coupling matrix elements!
# elph.g_xNNMM = elph.g_xsNNMM[:, s]
g_qklnn = elph.bloch_matrix(kpts, qpts, c_kn, u_ql=modes, spin=s)
g_sqklnn.append(g_qklnn)
# np.save("gsqklnn.npy", np.array(g_sqklnn))
Code¶
- class gpaw.elph.electronphonon.ElectronPhononCoupling(atoms, calc=None, supercell=(1, 1, 1), name='elph', delta=0.01, calculate_forces=False)[source]¶
Class for calculating the electron-phonon coupling in an LCAO basis.
The derivative of the effective potential wrt atomic displacements is obtained from a finite difference approximation to the derivative by doing a self-consistent calculation for atomic displacements in the +/- directions. These calculations are carried out in the
run
member function.The subsequent calculation of the coupling matrix in the basis of atomic orbitals (or Bloch-sums hereof for periodic systems) is handled by the
calculate_matrix
member function.Initialize with base class args and kwargs.
- Parameters
atoms (Atoms) – The atoms to work on.
calc (GPAW) – Calculator for the supercell finite displacement calculation.
supercell (tuple, list) – Size of supercell given by the number of repetitions (l, m, n) of the small unit cell in each direction.
name (str) – Name to use for files (default: ‘elph’).
delta (float) – Magnitude of displacements.
calculate_forces (bool) – If true, also calculate and store the dynamical matrix.
- apply_cutoff(cutmax=None, cutmin=None)[source]¶
Zero matrix element inside/beyond the specified cutoffs.
This method is not tested. This method does not respect minimum image convention.
- bloch_matrix(kpts, qpts, c_kn, u_ql, omega_ql=None, kpts_from=None, spin=0, name='supercell')[source]¶
Calculate el-ph coupling in the Bloch basis for the electrons.
This function calculates the electron-phonon coupling between the specified Bloch states, i.e.:
______ mnl / hbar ^ g = /------- < m k + q | e . grad V | n k > kq \/ 2 M w ql q ql
In case the
omega_ql
keyword argument is not given, the bare matrix element (in units of eV / Ang) without the sqrt prefactor is returned.Phonon frequencies and mode vectors must be given in ase units.
- Parameters
kpts (np.ndarray or tuple) – k-vectors of the Bloch states. When a tuple of integers is given, a Monkhorst-Pack grid with the specified number of k-points along the directions of the reciprocal lattice vectors is generated.
qpts (np.ndarray or tuple) – q-vectors of the phonons.
c_kn (np.ndarray) – Expansion coefficients for the Bloch states. The ordering must be the same as in the
kpts
argument.u_ql (np.ndarray) – Mass-scaled polarization vectors (in units of 1 / sqrt(amu)) of the phonons. Again, the ordering must be the same as in the corresponding
qpts
argument.omega_ql (np.ndarray) – Vibrational frequencies in eV.
kpts_from (list[int] or int) – Calculate only the matrix element for the k-vectors specified by their index in the
kpts
argument (default: all).spin (int) – In case of spin-polarised system, define which spin to use (0 or 1).
- calculate_gradient()[source]¶
Calculate gradient of effective potential and projector coefs.
This function loads the generated pickle files and calculates finite-difference derivatives.
- calculate_supercell_matrix(name='supercell', filter=None, include_pseudo=True)[source]¶
Calculate matrix elements of the el-ph coupling in the LCAO basis.
This function calculates the matrix elements between LCAOs and local atomic gradients of the effective potential. The matrix elements are calculated for the supercell used to obtain finite-difference approximations to the derivatives of the effective potential wrt to atomic displacements.
- Parameters
name (str) – User specified name of the generated JSON cache. Default is ‘supercell’.
filter (str) – Fourier filter atomic gradients of the effective potential. The specified components (
normal
orumklapp
) are removed (default: None).include_pseudo (bool) – Include the contribution from the psedupotential in the atomic gradients. If
False
, only the gradient of the effective potential is included (default: True).
- fourier_filter(V1t_xG, components='normal', criteria=0)[source]¶
Fourier filter atomic gradients of the effective potential.
This method is not tested.
- Parameters
V1t_xG (np.ndarray) – Array representation of atomic gradients of the effective potential in the supercell grid.
components (str) – Fourier components to filter out (
normal
orumklapp
).
- lcao_matrix(u_l, omega_l)[source]¶
Calculate the el-ph coupling in the electronic LCAO basis.
For now, only works for Gamma-point phonons.
This method is not tested.
- Parameters
u_l (np.ndarray) – Mass-scaled polarization vectors (in units of 1 / sqrt(amu)) of the phonons.
omega_l (np.ndarray) – Vibrational frequencies in eV.
- load_supercell_matrix(name='supercell')[source]¶
Load supercell matrix from cache.
- Parameters
name (str) – User specified name of the cache.