Electron-phonon coupling

Electron-phonon coupling is implemented for the LCAO mode.

Introduction

The electron-phonon interaction can be defined as

\[H_{el-ph} = \sum_{l,ij} g_{ij}^l c^{*}_i c_j ( a_l^{*} + a_l ) .\]

The phonon modes \(l\) are coupled to the electronic states \(i\), \(j\) via the electron-phonon coupling matrix

\[g_{ij}^l = \sqrt{ \frac{\hbar}{2 M \omega_l}} \langle i \vert \nabla_u V_{eff} \cdot \mathbf e_l \vert j \rangle .\]

\(\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)

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)

Zero matrix element inside/beyond the specified cutoffs.

This method is not tested. This method does not respect minimum image convention.

Parameters:

• cutmax (float) – Zero matrix elements for basis functions with a distance to the atomic gradient that is larger than the cutoff.

• cutmin (float) – Zero matrix elements where both basis functions have distances to the atomic gradient that is smaller than the cutoff.

bloch_matrix(kpts, qpts, c_kn, u_ql, omega_ql=None, kpts_from=None, spin=0, name='supercell')

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()

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)

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 or umklapp) 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)

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 or umklapp).

lcao_matrix(u_l, omega_l)

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')

Load supercell matrix from cache.

Parameters:

name (str) – User specified name of the cache.

set_basis_info(*args)

Store lcao basis info for atoms in reference cell in attribute.

Parameters:

args (tuple) – If the LCAO calculator is not available (e.g. if the supercell is loaded from file), the load_supercell_matrix member function provides the required info as arguments.

set_lcao_calculator(calc)

Set LCAO calculator for the calculation of the supercell matrix.