Raman and Infrared spectra with the Time-Dependent Self Consistent Harmonic Approximation

An introduction

The infrared spectrum is related to the dipole-dipole response function:

where the average \(\left<M(t)M(0)\right>\) is the quantum average at finite temperature.

Exploiting the TD-SCHA formalism introduced in the previous lecture, this response function can be written as:

where \(\boldsymbol{G}(\omega)\) is the TD-SCHA green function, while the \(\boldsymbol{r}\) and \(\boldsymbol{q}\) are vectors that quantify the perturbation and response, respectively.

In particular, if we neglect two-phonon effects (nonlinear coupling with light), we get that

where \({\mathcal Z}_{\alpha a}\) is the Born effective charge of atom \(a\), with polarization \(\alpha\), and \(G_{ab}(\omega)\) is the one-phonon green function, (its imaginary part is precisely the spectral function).

Indeed, we need to compute the effective charges. This can be done directly by quantum espresso using linear response theory (ph.x).

Prepare the infrared response

We need to attach the Raman Tensor and effective charges computed inside ir_raman_complete.pho to the final dynamical matrix, we will use this to initialize the response function calculation, as in Eq.1.

To attach the content of an espresso ph calculation (only Dielectric tensor, Raman Tensor, and Born effective charges) to a specific dynamical matrix, use

dyn.ReadInfoFromESPRESSO("ir_raman_complete.pho")

If you save the dynamical matrix in quantum espresso format, before the frequencies and the diagonalization, there will be the Dielectric tensor

Dielectric Tensor:

     1.890128098000           0.000000000000           0.000000000000
     0.000000000000           1.912811137000           0.000000000000
     0.000000000000           0.000000000000           1.916728724000

Followed by the effective charges and the Raman tensor.

Submitting the IR calculation

With the following script, we submit a TD-SCHA calculation for the IR.

import numpy as np
import cellconstructor as CC, cellconstructor.Phonons
import sscha, sscha.Ensemble
import tdscha, tdscha.DynamicalLanczos as DL

# Load the starting dynamical matrix
dyn_start = CC.Phonons.Phonons("start_dyn_ice")

# Load the ensemble
temperature = 0 # K
population = 2
n_configs = 10000

ensemble = sscha.Ensemble.Ensemble(dyn_start, temperature)
ensemble.load("data", population, n_configs)

# Load the final dynamical matrix
final_dyn = CC.Phonons.Phonons("final_dyn_ice")
final_dyn.ReadInfoFromESPRESSO("ir_raman_complete.pho")

# Update the ensemble weights for the final dynamical matrix
ensemble.update_weights(final_dyn, temperature)

# Setup the TD-SCHA calculation with the Lanczos algorithm
lanczos = DL.Lanczos(ensemble)
lanczos.ignore_v3 = True
lanczos.ignore_v4 = True

# If you have julia-enabled tdscha installed uncomment
# lanczos.mode = DL.MODE_FAST_JULIA
# for a x10-x15 speedup.

lanczos.init()


# Setup the IR response
polarization = np.array([1,0,0])  # Polarization of light
lanczos.prepare_ir(pol_vec = polarization)


# Run the algorithm
n_iterations = 1000
lanczos.run_FT(n_iterations)
lanczos.save_status("ir_xpol")

Congratulations! You ran your first TD-SCHA calculation. You can plot the results by using:

tdscha-plot.py ir_xpol.npz

The script tdscha-plot.py is automatically installed with the tdscha package.

Additionally, tdscha-plot.py takes three more parameters: the range of the frequencies to be displayed and the smearing.

Deep dive into the calculation

Let us dive a bit into the calculation. The beginning of the script should be almost self-explanatory, as we are just loading dynamical matrices, dielectric tensors, and effective charges.

The line

ensemble.update_weights(final_dyn, temperature)

deserves special attention. Here, we are changing the weights of the configurations inside the ensemble to simulate the specified dynamical matrix and temperature, even if they differ from those used to generate the ensemble. This is useful to compute the spectrum at several temperatures without extracting and calculating a new ensemble each time.

# Setup the TD-SCHA calculation with the Lanczos algorithm
lanczos = DL.Lanczos(ensemble)
lanczos.ignore_v3 = True
lanczos.ignore_v4 = True
lanczos.init()

Then we initialize the Lanczos algorithm for the tdscha, passing the ensemble.

The ignore_v3 and ignore_v4 are flags that, if set to True, the 3-phonon and 4-phonon scattering will be ignored during the calculation.

As you can see from the output, our IR signal had very sharp peaks because we ignored any phonon-phonon scattering process that may give rise to a finite lifetime.

By setting only ignore_v4 to True, we reproduce the behavior of the bubble approximation. Notably, while the four-phonon scattering is exceptionally computationally and memory demanding in free energy hessian calculations, within the Lanczos algorithm, accounting for the four-phonon scattering is only a factor two more expensive than using just the third order, without requiring any additional memory.

# Setup the IR response
polarization = np.array([1,0,0])  # Polarization of light
lanczos.prepare_ir(pol_vec = polarization)

Here we tell the Lanczos which kind of calculation we want to do. In other words, we set the \(\boldsymbol{r}\) and \(\boldsymbol{q}\) vectors in Eq.1 for the Lanczos calculation. - prepare_ir - prepare_raman - prepare_mode - prepare_perturbation

The names are intuitive; besides the Raman and IR, prepare_mode allows you to study the response function of a specific phonon mode, and prepare_perturbation enables defining a custom perturbation function.

# Run the algorithm
n_iterations = 1000
lanczos.run_FT(n_iterations)
lanczos.save_status("ir_xpol")

Here we start the calculation of the response function. The number of iterations indicates how many Lanczos steps are required. Each step adds a new pole to the green function. Therefore, many steps are necessary to converge broad spectrum features, while much less if the peaks are sharp. We save the status in such a way that we can get it back later.

Last, the commented line

lanczos.mode = DL.MODE_FAST_JULIA

This line only works if Julia and PyCall are correctly set up in the PC; in that case, run the script with python-jl instead of python. It will exploit a massive speedup of a factor between 10x and 15x. The calculation can also be run in parallel using mpirun before calling the Python executable (or python-jl). In this case, to work correctly, you should have mpi4py installed and working.

Analyze the output

In the last part, we employed the script tdscha-plot.py to display the simulation result. This is a quick way to show the results of a calculation.

Here, we will dive deeper into the calculation output file to extract the response function and get the results.

The Lanczos algorithm provides a set of coefficients \(a_i\), \(b_i\), and \(c_i\) through which the green function is evaluated thanks to a continued fraction:

Each iteration of the algorithm adds a new set of coefficients written in the standard output. Thanks to this expression, we only need the series of coefficients to compute the dynamical Green function at any frequency and with any smearing. The Green function can be computed with:

green_function = lanczos.get_green_function_continued_fraction(frequencies, smearing=smearing)

Here frequencies is an array in Rydberg. The response function is the opposite of the imaginary part of the green function; thus, to reproduce the plot, we have:

import tdscha, tdscha.DynamicalLanczos
import cellconstructor as CC, cellconstructor.Units
import numpy as np
import matplotlib.pyplot as plt

# Load the result of the previous calculation
lanczos = tdscha.DynamicalLanczos.Lanczos()
lanczos.load_status("ir_xpol_v4")

# Get the green function
W_START = 0
W_END = 3700
N_W = 10000
SMEARING = 10

frequencies = np.linspace(W_START, W_END, N_W)

# Convert in RY
frequencies_ry = frequencies / CC.Units.RY_TO_CM
smearing_ry = SMEARING / CC.Units.RY_TO_CM

# Compute the green function
green_function = lanczos.get_green_function_continued_fraction(frequencies_ry,
        smearing=smearing_ry)

# Get the response function
ir_response_function = - np.imag(green_function)

# Plot the data
plt.plot(frequencies, ir_response_function)
plt.show()

The previous script plots the data, precisely like plot-tdscha.py; however, now you have full access to the response function, both its imaginary and real parts.

Unpolarize Raman and IR

In the previous section, we saw how to compute Raman and IR with specific polarization of the incoming and outcoming radiation, and on oriented crystals (single crystals). However, the most common situation is a powder sample probed with unpolarized light.

In this case, we need to look at the Raman and IR response for unpolarized samples. While this is just the average of the IR signal’s x, y, and z, the Raman is more complex. In particular, unpolarized Raman signal can be computed from the so-called invariants, where the perturbations in the polarizations are the following:

The total Intensity of unpolarized Raman is:

The tdscha code implements a way to compute each perturbation separately. For example, the Raman response related to \(I_A\) is calculated with

lanczos.prepare_raman(unpolarized=0)

While the \(I_{B_i}\) is computed using index \(i\). For example, to compute \(I_{B_5}\):

# To compute I_B5 we do
lanczos.prepare_raman(unpolarized=5)

To get the total spectrum, you need to add the scattering factor \(n(\omega) + 1\) and sum all these perturbation with the correct prefactor (45 for \(I_A\) and 7 for the sum of all \(I_B\)).

To reset a calculation and start a new one, you can use

lanczos.reset()

which may be called before preparing the perturbation.

You should employ a supercell size sufficiently big to converge the simulation properly. In this case, the 1x1x1 supercell is too tiny to converge the calculation and get meaningful results.