How to Calculate First-Order Raman Spectra

by Stefan Kontur & Pasquale Pavone for exciting neon

Purpose: In this tutorial, you will learn how to calculate Raman polarizability tensors and frequency-dependent first-order scattering intensities using exciting. You will be guided through the example of silicon.

---

0. Prerequisites

In this tutorial it is assumed that you have already set the relevant environment variables. Otherwise, please have a look at How to Set Environment Variables for Tutorial Scripts. Before starting this tutorial, we recommend you to work through the tutorials

• Lattice dynamics of diamond and zincblende-structure crystals, using method = "sc".
• Excited states from TDDFT
• Excited states from BSE

Here, is a list of the scripts which are relevant for this tutorial with a short description.

• EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
• PLOT-raman-epsilon.py: Python visualization tool for the quality of dielectric function fits.
• PLOT-raman-spectrum.py: Python visualization tool for Raman spectra.
• PLOT-raman-resonance.py: Python visualization tool for the resonance behavior of Raman intensities.

Note: The symbol $ at beginnings of lines in code segments below indicates the shell prompt.

---

1. Theoretical background

---

☛ Theory of Raman scattering

1.1) Theory of first-order Raman scattering

The Raman scattering efficiency $S_{\textrm{i}\rightarrow \textrm{f}}^\phantom{M}$ can be expressed in terms of fluctuations of the dielectric function $\delta\pmb{\varepsilon}$ (see Refs. [1], [2], and [3]) as

(1)

\begin{align} S_{\textrm{i}\rightarrow \textrm{f}}^\phantom{M} = \frac{V_\textrm{s}\, \omega_\textrm{laser}^\phantom{I} \,\omega_\textrm{s}^3}{(4\pi)^2\,c^4_\phantom{I}}|\pmb{e}_\textrm{laser}^\phantom{I} \cdot \delta\pmb{\varepsilon}_{\textrm{i}\rightarrow \textrm{f}}^\phantom{I} \cdot \pmb{e}_\textrm{s}|^2\;, \end{align}

where the fluctuations are caused by a transition between the two states (i) and (f). In Eq. (1), $\omega_\textrm{laser}^\phantom{I}$ and $\omega_\textrm{s}$ are the angular frequencies, and $\pmb{e}_\textrm{laser}^\phantom{I}$ and $\pmb{e}_\textrm{s}$ the polarizations of the laser and the scattered light, respectively. $V_\textrm{s}$ is the scattering volume, which contains $N_\textrm{s}$ unit cells of volume $\Omega$, and $c$ the speed of light.

Due to conservation of momentum only phonons at Γ contribute to the first-order Raman spectrum. We simplify Eq. (1) even further by considering only a phonon mode with label $\nu$. The fluctuations in $\pmb{\varepsilon}$ due to this phonon mode can be written in terms of quantum-mechanical transition matrix elements between two vibrational states $|\sigma\rangle$ and $|\gamma\rangle$, times the derivative $\partial\pmb{\varepsilon}/\partial Q_\nu$ w.r.t. the normal coordinate $Q_\nu$ of the phonon mode $\nu$ ,

(2)

\begin{align} \left[ \delta\pmb{\varepsilon}_{\textrm{i}\rightarrow \textrm{f}}\right]_{\omega_\textrm{s}}^\nu \equiv \left[ \delta\pmb{\varepsilon}_{\sigma\rightarrow \gamma}\right]_{\omega_\textrm{s}}^\nu = \frac{\partial\pmb{\varepsilon}}{\partial Q_\nu} \langle\gamma|\hat{Q}_\nu|\sigma\rangle\;, \end{align}

where the adiabatic approximation is inherent. The complete expression including the temperature-dependent population of vibrational levels is then

(3)

\begin{align} S_{ij}^{\nu,\omega_\textrm{s}} = \frac{N_\textrm{s}\,\Omega\,\omega_\textrm{laser}^\phantom{I}\,\omega_\textrm{s}^3\,\sum_\sigma \exp(-E_\sigma/kT)\,\sum_\gamma \left|\dfrac{\partial\varepsilon_{ij}}{\partial Q_\nu} \langle\gamma|\hat{Q}_\nu|\sigma\rangle \right|^2} {(4\pi)^2c^4_\phantom{I}\,\sum_\sigma \exp(-E_\sigma/kT)}\;, \end{align}

where $i$ and $j$ labels the Cartesian (optical) components.

The macroscopic dielectric tensor $\pmb{\varepsilon}$ is related to the macroscopic dielectric susceptibility $\pmb{\chi}$ by

(4)

\begin{align} \pmb{\varepsilon} = \pmb{1}+4\pi\,\pmb{\chi}\;. \end{align}

---

☛ Work-flow of a Raman calculation

1.2) Work-flow of a Raman calculation

As can be seen above in Eq. (2), two ingredients are necessary to compute Raman spectra: Vibrational wavefunctions and derivatives of the dielectric function for every phonon mode. They are obtained from frozen phonon calculations of the potential and optical properties, so you have to input parameters for

• distorting the equilibrium geometry according to the displacement pattern of the mode,
• ground-state computations and forces,
• the calculation of optical spectra, and
• for the final calculation of the Raman spectrum.

Schematically, the computation involves the following steps:

---

☛ Phonon eigenvectors, actual displacements, and normal coordinates

1.3) Phonon eigenvectors, actual displacements, and normal coordinates

The actual structure distortion due to a given phonon mode $\nu$ at the Γ point is expressed in terms of the phonon amplitude $u_\nu$ and the normalized phonon eigenvectors $w^\nu_{\kappa, k}$ as

(5)

\begin{align} U^\nu_{\kappa, k} = u_\nu\, \sqrt{\frac{\mu_\nu}{M_\kappa}} w^\nu_{\kappa, k}\;, \end{align}

where $U^\nu_{\kappa, k}$ is the actual displacement in the Cartesian direction $k$ of the atom $\kappa$ in the unit cell, $M_\kappa$ is the mass of atom $\kappa$, and $\mu_\nu$ is the effective mass of the phonon mode $\nu$ defined by

(6)

\begin{align} \frac{1}{\mu_\nu}= \sum_{\kappa, k} \frac{|w^\nu_{\kappa, k}|^2}{M_\kappa}\;. \end{align}

Notice that the phonon eigenvectors $w^\nu_{\kappa, k}$ are normalized as

(7)

\begin{align} \sum_{\kappa, k} |w^\nu_{\kappa, k}|^2 = 1\;. \end{align}

Finally, the normal coordinate of the phonon mode $\nu$ is defined in terms of the phonon amplitude $u_\nu$ as

(8)

\begin{align} Q_\nu=\sqrt{\mu_\nu}\,u_\nu \;. \end{align}

---

☛ Raman mode polarizability and Raman tensors

1.4) Raman mode polarizability and Raman tensors

The frequency-dependent Raman polarizability tensor associated with the phonon mode $\nu$ is defined as

(9)

\begin{align} a_{ij}^\nu(\omega) = \Omega \sum_{\kappa, k} \frac{\partial\chi_{ij}(\omega)}{\partial U_{\kappa, k}} \sqrt{\frac{\mu}{M_\kappa}} \,w^\nu_{\kappa, k}\;, \end{align}

where $\mu$ is the reduced mass defined as

(10)

\begin{align} \frac{1}{\mu}= \sum_{\kappa} \frac{1}{M_\kappa}\;. \end{align}

Using the definition of the Raman tensor:

(11)

\begin{align} P_{ijk}^\kappa \equiv \Omega \frac{\partial\chi_{ij}(\omega)}{\partial U_{\kappa, k}}\;, \end{align}

the Raman mode polarizability tensor can be given (omitting the frequency dependence) as

(12)

\begin{align} a_{ij}^\nu = \sum_{\kappa, k} P_{ijk} \sqrt{\frac{\mu}{M_\kappa}}\:w^\nu_{\kappa, k}\;, \end{align}

An alternative expression of the Raman mode polarizability tensor in terms of the derivative of the dielectric tensor w.r.t. the phonon amplitude is

(13)

\begin{align} a_{ij}^\nu = \frac{\Omega}{4\pi} \sqrt{\frac{\mu}{\mu_\nu}} \frac{\partial\varepsilon_{ij}}{\partial u_\nu}\;. \end{align}

In this tutorial, the Raman mode polarizability will be calculated using this last expression.

---

2. Set up an input file for a Raman calculation

Setting up an input file for Raman computations involves the following steps.

---

2.1) Define structural parameters

First, the parameter defining your structure should be defined. In this example, we consider silicon. The structure part of the input file looks then like the following.

...
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="10.20">
         <basevect>0.5 0.5 0.0</basevect>
         <basevect>0.5 0.0 0.5</basevect>
         <basevect>0.0 0.5 0.5</basevect>
      </crystal>
 
      <species speciesfile="Si.xml">
         <atom coord="0.00 0.00 0.00" />
         <atom coord="0.25 0.25 0.25" />
      </species>
 
   </structure>
...

The chosen lattice parameter of 10.20 refers to the converged equilibirum parameter for silicon, as obtained with the LDA exchange-correlation functional.

---

2.2) Parameters for the ground-state calculation

Then, parameters for performing the SCF have to be specified. In this example, we use the LDA functional of Perdew and Wang (PW), other choices can be found in the description of the attribute xctype of the element groundstate.

...
   <groundstate 
      do="fromscratch" 
      ngridk="4 4 4" 
      rgkmax="5.0"
      xctype="LDA_PW"/>
...

Note that this example of input contains values for the parameters ngridk and rgkmax which allow for fast calculations but should not be considered realistic. The convergence w.r.t. the computational parameters (e.g., the value of the attributes ngridk and rgkmax) of all target physical quantities (in this tutorial, e.g., phonon frequencies, Raman tensors and intensities) should be investigated in order to achieve the desired precision.

---

2.3) Set up the calculation of optical spectra

...
   <xs 
      xstype="TDDFT"
      ngridk="8 8 8"
      ngridq="8 8 8"
      vkloff="0.093 0.279 0.461"
      rgkmax="5.0"
      nempty="30"
      scissor="0.0316"
      broad="0.00367"
      dfoffdiag="true"
      tevout="true">
 
      <energywindow 
         intv="0.0 1.0" 
         points="1000"/>
 
      <tddft fxctype="RPA"/>
 
      <qpointset>
         <qpoint> 0.0 0.0 0.0 </qpoint>
      </qpointset>
   </xs>
...

For an detailed explanation of this input we refer to the tutorial Excited states from TDDFT and the Input Reference for the xs element. Like in the ground-state part, we use here parameters which allow for a quick calculation, rather than yielding a well converged result. In particular, we use RPA to obtain the dielectric function for silicon, which is not the optimal choice (see also the execises section on how to improve the calculation).

---

2.4) Set up the calculation of Raman intensities

The calculation of Raman intensities is triggered by the presence of the subelement raman inside the element properties. The main attributes of the raman element are shown in the following example.

...
   <properties>
      <raman 
         getphonon="fromscratch"
         mode="4"
         nstep="5"
         ninter="600"
         displ="0.01"
         degree="2"
         elaser="1.1" elaserunit="eV"
         temp="298.15"
         broad="3.0">
 
         <energywindow 
            intv="0.0 0.005" 
            points="4000" />
 
      </raman>
   </properties>
...

First of all, we need to calculate the relevant phonon modes, as indicated by the presence of the attribute getphonon = "fromscratch". This calculation does not take very long, as we are looking at the Γ-point modes only. The element phonons can be added optionally to finetune the calculation (i.e., for tuning the default parameters). However, the phonon calculation will always be done for a "supercell" that consists of just one primitive unit cell and all the phonons have vanishing wavevector q=0. Out of the 6 phonon modes we choose the 4th one (mode = "4"), i.e., one of the degenerate optical modes (in our example, it should be the one polarized along the z-direction).

Next, we have to create distorted geometries according to the considered phonon mode and to compute total energies, forces, and, more involving, the dielectric functions for them. The generation of these distorted geometries is tuned by the attributes nstep = "5" and displ = "0.01". Their meaning is the following: 5 geometries are constructed along the phonon mode, with a stepsize in the phonon amplitude of 0.01 Bohr. This stepsize refers to the phonon amplitude of the wavelike displacement of all atoms in the unit cell. Thus, for unit cells with a large number of atoms, large values should be used. As in the case of phonons, the stepsize should be checked carefully. Too low or too large values result in numerical instabilities. One geometry is very close to the equilibrium one, the other ones are distorted by one time, two times, etc., the stepsize in both positive and negative directions.

The potential along the phonon mode and the value of the components of the dielectric function are fitted to the values from the frozen phonon calculations. The fit function is always a polynomial, whose degree we set to the value of 2 (i.e., a harmonic oscillator for the phonon problem) by degree = "2", its default value. The value of niter is the number of intervals in the numerical (FE) solution of the oscillator problem.

The laser energy for the Raman spectrum is set by elaser = "1.1", the default units of this energy (Hartree) can be changed to standard electron volts by setting elaserunit = "eV". In any case, a complete resonance spectrum over the whole range of the available optical spectra (see the element energywindow inside the xs block) is computed. The energy range of the Raman spectrum is set by the element energywindow inside the raman block and should not be confused with other occurrences of this element, e.g., inside xs. The energy range in this example is set between 0 and 0.005 Hartree, i.e., roughly 0 to 1100 cm^-1. Room temperature is chosen for this run by setting (temp = "298.15"), which is also the default value.

---

2.5) Other options to get phonon modes

---

☛ Show options

The considerations in this section are useful either if results of a previous phonon calculations are available or if the phonon displacement patterns are easy to construct directly from symmetry considerations. In this case, there are the following alternative options:

• Setting getphonon = "fromfile" allows the code to read the dynamical matrices from a collection of DYN_fileext.OUT files (see Lattice dynamics of diamond and zincblende-structure crystals, with method = "sc", for details) and proceeds as above. In this case exciting searches DYN_fileext.OUT files in the same directory as where input.xml is located.
• The choice getphonon = "symvec" tries to construct symmetry vectors and uses them instead of the true eigenvectors. In cases of low symmetry and/or many atoms in the unit cell the result becomes ambiguous and the procedure fails. getphonon = "symveccheck" stops right after the construction of symmetry vectors and can be used for a dry run.
• getphonon = "readinput" reads the eigenvector of a single mode from the input file. The eigenvector does not need to be normalized in the input. The order of the (complex) components of the eigenvector follows the order of atoms inside the element structure and lists three items (Cartesian components) for every atom. The way to give this input is shown in the following example for T_2g mode of silicon which is polarized along the z direction:

...
   <properties >
      <raman getphonon="readinput" ... >
 
         <eigvec comp=" 0.0 0.0" />
         <eigvec comp=" 0.0 0.0" />
         <eigvec comp=" 1.0 0.0" />
         <eigvec comp=" 0.0 0.0" />
         <eigvec comp=" 0.0 0.0" />
         <eigvec comp="-1.0 0.0" />
 
      </raman>
   </properties>
...

---

2.6) The basic input file for a Raman calculation

We are now ready to run the actual Raman calculation. Create a subdirectory si-raman and move inside it to run the calculations.

$ mkdir si-raman
$ cd si-raman

The complete input file looks like the following:

<input>
   <title>Silicon Raman RPA</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="10.20">
         <basevect>0.50 0.50 0.00</basevect>
         <basevect>0.50 0.00 0.50</basevect>
         <basevect>0.00 0.50 0.50</basevect>
      </crystal>
      <species speciesfile="Si.xml">
         <atom coord="0.00 0.00 0.00"></atom>
         <atom coord="0.25 0.25 0.25"></atom>
      </species>
   </structure>
 
   <groundstate 
      do="fromscratch" 
      ngridk="4 4 4" 
      rgkmax="5.0"
      xctype="LDA_PW"/>
 
   <properties>
      <raman 
         getphonon="fromscratch"
         mode="4"
         nstep="5"
         ninter="600"
         displ="0.01"
         degree="2"
         elaser="1.1" elaserunit="eV"
         temp="298.15"
         broad="3.0">
 
         <energywindow 
            intv="0.0 0.005" 
            points="4000" />
 
      </raman>
   </properties>
 
   <xs 
      xstype="TDDFT"
      ngridk="8 8 8"
      ngridq="8 8 8"
      vkloff="0.093 0.279 0.461"
      rgkmax="5.0"
      nempty="30"
      scissor="0.0316"
      broad="0.00367"
      dfoffdiag="true"
      tevout="true">
 
      <energywindow 
         intv="0.0 1.0" 
         points="1000" />
 
      <tddft fxctype="RPA"/>
 
      <qpointset>
         <qpoint> 0.0 0.0 0.0 </qpoint>
      </qpointset>
   </xs>
 
</input>

Please, remember that after creating the input file input.xml, the path for the species file must be set using the command

$ SETUP-excitingroot.sh

To perform actual calculations, you may run the script EXECUTE-single.sh and choose a name for the directory in which you run with the calculations.

$ EXECUTE-single.sh mode-4

This run will take about 5 to 10 minutes, depending on the velocity of the processors in your computer.

---

3. Results

When the run has been completed, move to the directory mode-4 in order to analyze the results.

$ cd mode-4

---

3.1) Output files of the Raman calculation

In addition to the basic output files already mentioned in How to start an exciting calculation, we requested a phonon calculation at Γ and obtain all the files from it (see Lattice dynamics of diamond and zincblende-structure crystals, with method = "sc", for details). The ground-state and optical runs can be found in subdirectories, that are named MOD004_DISP01, MOD004_DISP02, etc., and refer to the different distorted geometries for all modes requested (in our case, five displacements for the 4th mode). Inside these directories you can find all output files as described in How to start an exciting calculation, Excited states from TDDFT, and Excited states from BSE.

Output files which are specific to the Raman calculation are described below and can be found in the folder mode-4 of our example.

filename description INFO_RAMAN.OUT the file containing info about what steps where taken during the run and what time they took RAMAN.OUT the main output file containing data for the potential of the phonon mode, changes of the dielectric function, fitting results, and scattering efficiencies from all contributions RAMAN_SYM.OUT file containing symmetry information which is useful in the context of Raman scattering RAMAN_SPEC_OCxy_MODmod.OUT first-order Raman spectra for optical components xy and modes mod for the laser energy fiven in the input file RAMAN_RESONANCE_OCxy_MODmod.OUT resonance behavior of optical components xy and modes mod across the whole range of energies in the optical spectra RAMAN_EPSILON_OCxy_MODmod.OUT file containing results of the fits to data for the dielectric function (real and imaginary parts) for the laser energy given in the input file RAMAN_POTENTIAL_MODmod.OUT file containing results of the fits to data for the potential along mode mod

---

3.2) Checking the eigenvector of the phonon mode

First, one has to check if the phonon mode which has been calculated (in this case the mode 4) is the correct one. This step is very important in case the phonon mode we have chosen is a degenerate mode (as in our example). The polarization of the phonon mode affects the symmetry of the resulting Raman tensor.

By running the following command:

$ grep 'Mode   4' -A9 INFO_RAMAN.OUT | head -n10

we obtain on the screen

Mode   4
 Eigenvalue:       526.09420 cm^-1

      Atom   Polarization    Eigenvector
         1              1  -0.0000000   0.0000000
         1              2   0.0000000   0.0000000
         1              3   0.7071068   0.0000000
         2              1  -0.0000000   0.0000000
         2              2   0.0000000   0.0000000
         2              3  -0.7071068   0.0000000

which confirms that we have chosen the optical phonon polarized along the z-direction (OPT-001).

---

Note: The order of the phonon modes can depend on the choice of the computational parameters. For instance, the optical mode polarized along the z-direction calculated with rgkmax = "5.0" is the 4th one. Nevertheless, it could be the 5th or the 6th one, if a different value (e.g., "7.0") is chosen. This is a very frequent case for degenerate phonon modes. Furthermore, this often happens by using different diagonalization libraries and compilers. For these reasons, it is strongly suggested to check the actual polarization of the phonon mode which has been calculated.

---

3.3) Derivatives of the dielectric tensor

The derivatives of the dielectric tensor w.r.t. the phonon amplitude are a fundamental ingredient for the calculation of the Raman tensors (see Eq. (13)). The values of this derivative for a laser energy of 1.1 eV and the 4th phonon mode are written inside the file RAMAN.OUT [4], and can be obtained by typing

$ grep 'Derivatives' -A5 RAMAN.OUT | tail -n6

The output on the screen should be (notice that the derivatives are given in units of Bohr^-1)

 Derivatives of the dielectric function deps/du: 
 Re                                         Im

  (       -0.249      -4.645       0.000 )       (       -0.008      -0.141       0.000 ) 
  (       -4.645      -0.249       0.000 )       (       -0.141      -0.008       0.000 ) 
  (        0.000       0.000       0.028 )       (        0.000       0.000       0.001 )

In this output, we notice the two blocks (3$\times$3 matrices) corresponding to the real and imaginary part of the derivative, respectively. By symmetry considerations [5], the resulting tensor should be a complex matrix of the following type

[   0   a   0  ]
[   a   0   0  ]       for an optical phonon polarized along the z-direction
[   0   0   0  ]

with a very small imaginary part (as it was indeed found in our calculations).

However, the tensor has to be symmetric and have zero diagonal values, but this is not the case for the numerical results we got. This is due to our choice of not fully optimized computational parameters.

In order to see how these derivatives have been obtained, we can plot the components ε_ij (in the figures below indicated as ε_αβ) as a function of the phonon amplitude u. Invoke the following script from the command line

$ PLOT-raman-epsilon.py RAMAN_EPSILON_OC12_MOD004.OUT -0.04 0.04

In this way, you will obtain a plot of the optical component "12" (or, equivalently "xy") of the 4th mode as a function of the phonon amplitude (in the range between -0.04 and 0.04 Bohr). The result (PLOT.png) can be seen in the next figure.

The fit looks quite fine, nevertheless, notice that the value of this dielectric-tensor component should be zero at zero phonon amplitude due to the cubic symmetry of the crystal. We can look now at the component "11" (or, equivalently "xx").

$ PLOT-raman-epsilon.py RAMAN_EPSILON_OC11_MOD004.OUT -0.04 0.04

In this case, the result stored in PLOT.png is

Clearly, the minimum of this function is not at the equilibrium geometry and its first derivative is not zero. This problem will be addressed below in the discussion about the convergence of the calculations.

---

3.4) First-order Raman spectrum

The actual first-order Raman spectrum can be visualized using the command

$ PLOT-raman-spectrum.py RAMAN_SPEC_OC12_MOD004.OUT 500 550

We display a simulated spectra in the range between 500 and 550 cm^-1 for the 4th mode and the "12" optical component, broadened by 3.0 cm^-1, as requested by the attribute broad = "3.0".

The full first-order Raman spectrum is then a superposition of the spectra of all modes for the optical component of interest. In our case of silicon, the three optical modes are degenerate and the spectrum looks just as above.

---

3.5) Raman polarizability

The Raman mode polarizability tensor is calculated using Eq. (13). This equation can be further simplified in the case of Silicon. First, if $M$ is the atomic mass of Silicon, then

(14)

\begin{align} \mu=\frac{M}{2}\;. \end{align}

Furthermore, for the OPT-001 phonon mode the effective mass is

(15)

\begin{align} \mu_{{\tt OPT\textrm{-}001}}^\phantom{M}=M\;. \end{align}

Therefore, the only non vanisching component of the frequency-dependent Raman mode polarizability tensor for the OPT-001 phonon mode can be written as

(16)

\begin{align} a_{12}^\nu(\omega)=a_{{\tt OPT\textrm{-}001}}^{\tt RPA}\!(\omega) = \frac{\Omega}{4\sqrt{2}\pi}\,\frac{\partial \varepsilon_{12}^\phantom{M}(\omega)}{\partial u}\,, \end{align}

where we explicitly include the label RPA denoting the use of the random-phase approximation in the calculations.

For $\omega$=$\omega_\textrm{laser}$ the components defined in Eq. (16) have a small, non vanishing, imaginary parts. In the literature, the following real quantity

(17)

\begin{align} a_\textrm{R}^{\tt RPA}(\omega_\textrm{laser}^\phantom{M}) \equiv \left| a_{{\tt OPT\textrm{-}001}}^{\tt RPA}\!(\omega_\textrm{laser}^\phantom{M}) \right| = \frac{\Omega}{4\sqrt{2}\pi}\,\left|\frac{\partial \varepsilon_{12}^\phantom{M}(\omega_\textrm{laser}^\phantom{M})}{\partial u}\right|\,, \end{align}

is referred as Raman polarizability of the given phonon mode.

For $\omega_\textrm{laser}$ = 1.1 eV, using the derivative $\partial\varepsilon/\partial u$ reported in the output file RAMAN.OUT (expressed in units of Å^-1) and a f.c.c. lattice constant of 10.20 Bohr = 5.3976 Å, we obtain

• Raman polarizability = $a_\textrm{R}^{\tt RPA}$(1.1 eV) = 19.43 Ų

---

3.6) Resonance behavior

The main result of our calculation is the resonance behavior of the Raman signal. We can produce an appropriate plot invoking the tool script

$ PLOT-raman-resonance.py RAMAN_RESONANCE_OC12_MOD004.OUT 0 5

In this case we considered the contribution to the intensity of the Raman spectrum from the "xy"-component. The obtained plot looks like the following.

What we see is a logarithmic plot of the squared Raman susceptibility, which is defined in this context as

(18)

\begin{align} \left|\frac{\partial\chi}{\partial Q}\right|^2 = \left|\frac{1}{4\pi}\frac{\partial\varepsilon}{\partial Q}\right|^2\;. \end{align}

In absence of any resonance (i.e., if the laser energy is not close to an optical transition of the material), the squared Raman susceptibility is almost constant. In our case of silicon, resonances of the Raman signal appear at about 3.4 eV and beyond. Clearly, the quality of this result depends also on the method chosen to compute the dielectric function, RPA in our case.

An experiment on the Raman resonance of the T_2g (or Γ_25') phonon in silicon can be found in Ref. [7]. Furthermore, the theoretical results from this tutorial can be compared with the work of Ref. [8].

---

3.7) Improving the results by a better choice of the computational parameters

---

☛ Show more precise results

More precise results can be obtained by an improved choice of the computational parameters. Notice that in this case the computational times could be drastically increased!

The results obtained in this section correspond to the following settings:

• (rgkmax = "7.0" and ngridk = "24 24 24") for both ground-state and optical properties;
• nstep = "7", degree = "4", and ninter = "1800" inside raman.

Notice that an increased value of ngridk is fundamental for getting the correct symmetry of the dielectric tensor and its derivatives.

In our specific case, moreover, the optical phonon mode of silicon with eigendisplacements along the z-direction (OPT-001) is fully determined by symmetry:

Atom   Polarization    Eigenvector      Non-normalized eigenvector
  1         1            0       0               0      0
  1         2            0       0               0      0
  1         3         sqrt(2)/2  0               1      0
  2         1            0       0               0      0
  2         2            0       0               0      0
  2         3        -sqrt(2)/2  0              -1      0

Therefore, we can avoid the quite expensive calculation of the dynamical matrix at the Γ point, by further changing the input file as follows.

• Set the attribute getphonon of the raman element to "readinput":

      <raman getphonon="readinput" ... >

• Inside raman, just after the energywindow block, add the block defining the non-normalized eigenvector as:

         <eigvec comp=" 0.0 0.0" />
         <eigvec comp=" 0.0 0.0" />
         <eigvec comp=" 1.0 0.0" />
         <eigvec comp=" 0.0 0.0" />
         <eigvec comp=" 0.0 0.0" />
         <eigvec comp="-1.0 0.0" />

Notice that with this choice of parameters the output files corresponding to the optical phonon mode of silicon with eigendisplacements along the z-direction are labeled with the extension MODE001 (only one mode is taken into account).

The calculation can be performed with these fine settings using with the command

$ time exciting_smp &

After the calculation is completed (be aware that it lasts much longer…), the values obtained for the derivative of the dielectric tensor for a laser energy of 1.1 eV and the chosen phonon mode can be obtained by typing inside the outcoming directory the command

$ grep 'Derivatives' -A5 RAMAN.OUT | tail -n6

The values should look like the following (notice that the derivatives are given in units of Bohr^-1):

 Derivatives of the dielectric function deps/du: 
 Re                                         Im

  (        0.001      -3.821       0.000 )       (        0.000      -0.119       0.000 ) 
  (       -3.821       0.001       0.000 )       (       -0.119       0.000       0.000 ) 
  (        0.000       0.000       0.001 )       (        0.000       0.000       0.000 )

In the following we show the other converged results.

---

• Raman polarizability, for the component "xy"

Phonon mode Laser frequency [eV] Raman polarizability [Å2] Description Reference OPT-001 1.1 15.98 TDDFT-RPA This work OPT-001 1.1 15.87 RPA Ref.[8] OPT-001 1.1 19.75 BSE Ref.[8] OPT-001 1.1 23 ± 4 Experiment Ref.[6]

---

• Dielectric function, component "xy":

$ PLOT-raman-epsilon.py  RAMAN_EPSILON_OC12_MOD001.OUT -0.04 0.04

---

• Dielectric function, component "xx":

$ PLOT-raman-epsilon.py  RAMAN_EPSILON_OC11_MOD001.OUT -0.04 0.04

---

• First-order Raman spectrum:

$ PLOT-raman-spectrum.py RAMAN_SPEC_OC12_MOD001.OUT 490 540

---

• Resonance behavior:

$ PLOT-raman-resonance.py RAMAN_RESONANCE_OC12_MOD001.OUT 1 4

---

Exercise: Further improvements

• Rerun the calculation for different laser energies in the attribute elaser and check the behavior of the dielectric function again. Do you see some dependence of the quality of fits and derivatives on the laser energy? You can restart from the data you have just computed: Move inside the directory mode-4, set do = "skip" inside groundstate and getphonon = "fromfile" inside raman. Then, invoke exciting_smp from the command line.
• Converge the parameters relevant for ground-state and optics calculations. Run the calculation for better fit functions, by setting nstep and degree to higher values. These runs have to be done every time in a new directory, otherwise old data will be read in.
• Suggestion for further work on the topic: RPA is in general not the best choice for a material such as silicon, where excitonic effects are significant. Try using BSE instead and investigate how this changes the computed resonance behavior.

---

Bibliography

1. Knoll and Ambrosch-Draxl, in Proceedings of the International Workshop on Anharmonic Properties of High-T_c Cuprates, edited by Mihailovich, Ruani, Kaldis, and Müller (World Scientific, Singapore, 1995), p. 220 (PDF).

2. Ambrosch-Draxl, Auer, Kouba, and Sherman, Phys. Rev. B 65, 064501 (2002) (link).

3. Kontur and Draxl, in preparation.

4. This information can be found in the RAMAN.OUT file after the fitted polynomial coefficients for the potential (force) and the components of the dielectric function, under the header START CALCULATION. The derivatives are computed for the $u$ that gives the minimum of the potential, written under the sub-header Phonon potential shifted into minimum.

5. For one of the three T_2g phonons in the point group O_h, the other two being

[   0   0   a  ]
[   0   0   0  ]       for an optical phonon polarized along the y-direction
[   a   0   0  ]

and

[   0   0   0  ]
[   0   0   a  ]       for an optical phonon polarized along the x-direction
[   0   a   0  ]

respectively. Depending on the appearance of the phonon eigenvector also linear combinations can be obtained. See Bilbao crystallographic server for other irreducible representations and point groups.

6. Wagner and Cardona, Solid State Commun. 48, 301 (1983) (link).

7. Compaan and Trodahl, Phys. Rev. B 29, 793 (1984) (link).

8. Gillet, Giantomassi, and Gonze, Phys. Rev. B 88, 094305 (2013) (link).