Raman spectroscopy for extended systems¶

This tutorial shows how to calculate the Raman spectrum of diamond using the electron-phonon coupling based Raman methods as described in Raman spectroscopy.

Phonons¶

The phonon frequencies enter the Raman calculation by defining the line positions in the final spectrum. Their associated mode vectors are necessary to project the electron-phonon matrix from Cartesian coordinates in the modes.

It is crucial to use tightly relaxed structures and well-converged parameters to obtain accurate phonon frequencies. In this case we choose the separate out the phonon from the potential calculation for illustrative purposes: (phonon.py)

import numpy as np
from ase.build import bulk
from ase.phonons import Phonons
from ase.units import invcm
from gpaw import GPAW
from gpaw.mpi import world

atoms = bulk('C', 'diamond', a=3.567)
calc = GPAW(mode={'name': 'pw', 'ecut': 500},
            kpts=(5, 5, 5), xc='PBE',
            symmetry={'point_group': False},
            convergence={'density': 0.5e-5},
            txt='phonons.txt')

# Phonon calculation
ph = Phonons(atoms, calc, supercell=(1, 1, 1), delta=0.01)
ph.run()

# To display results (optional)
ph.read(method='frederiksen', acoustic=True)
frequencies = ph.band_structure([[0, 0, 0], ])[0]  # Get frequencies at Gamma
if world.rank == 0:
    # save phonon frequencies for later use
    np.save("vib_frequencies.npy", frequencies)
    print('  i    cm^-1')
    print('------------')
    for i, fr in enumerate(frequencies):
        print('{:3d}  {:4.2f}'.format(i, fr / invcm))

As exercise, check the dependence of the phonon frequencies with the calculation mode, supercell size and convergence parameters.

This given set of parameters yield this result:

  i    cm^-1
------------
  -0.03
  -0.00
  -0.00
1304.46
1304.61
1304.73

Effective Potential¶

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, of 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 (2x2x2) 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. (elph.py)

from ase.build import bulk
from gpaw import GPAW
from gpaw.elph.electronphonon import ElectronPhononCoupling

atoms = bulk('C', 'diamond', a=3.567)
calc = GPAW(mode='lcao', basis='dzp',
            kpts=(3, 3, 3), xc='PBE',
            symmetry={'point_group': False},
            txt='elph.txt')

# Displacements for potential
elph = ElectronPhononCoupling(atoms, calc, supercell=(2, 2, 2),
                              calculate_forces=False)
elph.run()

This calculation merely dumped the effective potential at various displacements onto the harddrive. We now need to calculate the actual derivative and the matrix elements with the Kohn-Sham wave-functions.

For this we first need to complete a ground-state calculation for the supercell with all bands converged. This calculation needs to be done in LCAO mode with parallelisation over domains disabled. (supercell_matrix.py)

from ase.build import bulk
from gpaw import GPAW
from gpaw.raman.elph import EPC

supercell = (2, 2, 2)
atoms = bulk('C', 'diamond', a=3.567)
atoms_sc = atoms * supercell
calc = GPAW(mode='lcao', basis='dzp',
            kpts=(5, 5, 5), xc='PBE',
            symmetry={'point_group': False},
            convergence={'bands': 'nao', 'density': 1e-5},
            parallel={'domain': 1},
            txt='gs_super.txt')
atoms_sc.calc = calc
atoms_sc.get_potential_energy()

# Supercell matrix
elph = EPC(atoms, supercell=supercell)
elph.calculate_supercell_matrix(calc, include_pseudo=True)

The calculate_supercell_matrix() method will then compute the gradients and calculate the matrix elements. The results are saved in pckl files in a basis of LCAO orbitals and supercell indices.

Momentum matrix¶

The Raman tensor does not only depend on the electron-phonon matrix, but also on the transition probability between electronic states, which in this case is expressed in the momentum matrix elements. Those can be extracted directly from a converged LCAO calculation of the primitive cell. (momentum_matrix.py)

from ase.build import bulk
from gpaw import GPAW
from gpaw.raman.dipoletransition import get_momentum_transitions

atoms = bulk('C', 'diamond', a=3.567)
calc = GPAW(mode='lcao', basis='dzp',
            kpts=(5, 5, 5), xc='PBE',
            symmetry={'point_group': False},
            convergence={'bands': 'nao', 'density': 1e-5},
            parallel={'band': 1},
            txt='mom.txt')

atoms.calc = calc
atoms.get_potential_energy()
atoms.calc.write("gs.gpw", mode="all")
get_momentum_transitions(calc.wfs, savetofile=True)

It it convienient to save the calculation, as we will need it for the next step.

Phonon mode projected electron-phonon matrix¶

With all the above calculations finished we can extract the electron-phonon matrix in the Bloch basis of the primitive cell projected onto the phonon modes: (elph_matrix.py)

from ase.phonons import Phonons
from gpaw import GPAW
from gpaw.raman.elph import EPC

# Load pre-computed ground state calculation (primitive cell)
calc = GPAW('gs.gpw', parallel={'band': 1})
atoms = calc.atoms

# Load results from phonon and electron-phonon coupling calculations
phonon = Phonons(atoms, supercell=(1, 1, 1))
elph = EPC(atoms, supercell=(2, 2, 2))

# Construct electron-phonon matrix of Bloch functions
elph.get_elph_matrix(calc, phonon)

This will save the electron-phonon matrix as a numpy file to the disk. The optional load_gx_as_needed tag prevents from all supercell pickle files being read at once. Instead they are loaded and processed one-by-one. This can save lots of memory for larger systems with hundreds of atoms, where the supercell matrix can be over 100GiB large.

Note: This part has not been tested properly for parallel runs and should be done in serial mode only.

Raman spectrum¶

With all ingredients provided we can now commence with the computation of the Raman tensor. This can take some time, if calculation of the non-resonant terms is included. (raman_intensities.py)

# web-page: Raman_all.png
import numpy as np
from gpaw import GPAW
from gpaw.raman.raman import (calculate_raman, calculate_raman_intensity,
                              plot_raman)

# Load pre-computed calculation
calc = GPAW('gs.gpw', parallel={'domain': 1, 'band': 1})
atoms = calc.atoms

# laser frequency 633 nm approx 1.958676 eV
w_l = 1.958676
suffix = '632nm'

# use previously saved phonon frequencies
w_ph = np.load('vib_frequencies.npy')

# Scan through all polarisations
pollist = []
for d_i in (0, 1, 2):
    for d_o in (0, 1, 2):
        # Calculate mode resolved Raman tensor for given direction
        calculate_raman(calc, w_ph, w_l, d_i, d_o, resonant_only=True,
                        suffix=suffix)
        if calc.wfs.kd.comm.rank == 0:
            # Calculate Raman intensity
            calculate_raman_intensity(w_ph, d_i, d_o, suffix=suffix)
            pollist.append('{}{}_{}'.format('xyz'[d_i], 'xyz'[d_o], suffix))

# And plot
if calc.wfs.kd.comm.rank == 0:
    plot_raman(figname='Raman_all.png', RIsuffix=pollist)

This part consists of three steps. First we calculate the Raman tensor entry for a given set of polarisations, which are saved in files of the form Rlab_zy.npy. From those the intensities, summed over all phonon modes, are calculated on a grid and saved in files of the form RI_zy.npy. Lastly, the spectrum is plotted.

Note: This part has not been tested properly for parallel runs and should be done in serial mode only.