How to Calculate First-Order Raman Spectra

by Stefan Kontur & Pasquale Pavone for exciting nitrogen

Purpose: In this tutorial, you will learn how to calculate Raman shifts 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

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.

This tutorial require the use of the previous python scripts. Click above on the script name for downloading it. Assuming, e.g., that the script PLOT-raman-epsilon.py has been downloaded in the directory $HOME/Downloads, move it to the directory$EXCITINGTOOLS and make it executable:

$mv$HOME/Downloads/PLOT-raman-epsilon.py $EXCITINGTOOLS/PLOT-raman-epsilon.py$ chmod u+rwx $EXCITINGTOOLS/PLOT-raman-epsilon.py  Proceed in a similar way for the other two python scripts. Note: The symbol$ at beginnings of lines in code segments below indicates the shell prompt.

### 1. Introductory remarks

The Raman scattering efficiency $S$ can be expressed in terms of fluctuations of the dielectric function $\delta\varepsilon$ as (see Refs.[1], [2], and [3])

(1)
\begin{align} S_{i\rightarrow f} = \frac{V_s\; \omega_L \;\omega_s^3}{(4\pi)^2\,c^4}|{\bf e}_L \cdot \delta\varepsilon_{i\rightarrow f} \cdot {\bf e}_s|^2\;, \end{align}

where the fluctuations are caused by a transition between the two states $i$ and $f$. In Eq. (1) $\omega_L$ and $\omega_s$ are the angular frequencies, and ${\bf e}_L$ and ${\bf e}_s$ the polarizations of the laser and the scattered light, respectively. $V_s$ is the scattering volume, which contains $N$ unit cells of volume $V_c$, and $c$ the speed of light. Fluctuations in $\varepsilon$ due to phonons can be written in terms of quantum mechanical transition matrix elements between two vibrational states $|\mu\rangle$ and $|\nu\rangle$, times the derivative $\partial\varepsilon/\partial Q$ with respect to the normal coordinate $Q$ of the mode,

(2)
\begin{align} \left[ \delta\varepsilon_{i\rightarrow f}^{\alpha\beta}\right]_{\omega_s} \equiv \left[ \delta\varepsilon_{\mu\rightarrow \nu}^{\alpha\beta}\right]_{\omega_s} = \frac{\partial \varepsilon^{\alpha\beta}}{\partial Q} \langle\nu|Q|\mu\rangle\;, \end{align}

where the adiabatic approximation is inherent and terms beyond first order are omitted. The indeces $\alpha$ and $\beta$ refer to the optical components of the dielectric tensor. The complete expression including the temperature-dependent population of vibrational levels is then

(3)
\begin{align} S^{\alpha\beta}_{\omega_s} = \frac{N\,V_c \,\omega_L\, \omega_s^3 \;\sum_\mu \exp(-E_\mu/kT) \;\sum_\nu \left|\dfrac{\partial \varepsilon^{\alpha\beta}}{\partial Q} \langle\nu|Q|\mu\rangle\right|^2} {(4\pi)^2c^4 \sum_\mu \exp(-E_\mu/kT)}\;. \end{align}

Due to conservation of momentum only phonons at Γ contribute to the spectrum.

As can be seen in Eq. (2), two ingredients are necessary to compute Raman spectra: Nuclear 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:

### 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. A detailed description of how to achieve convergence with respect to these parameters can be found in Simple convergence tests. ##### 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.0312" 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) Input for computation 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" displ="0.01" degree="2" elaser="1.16" 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 and continue the work flow with it. 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 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 laser energy for the Raman spectrum is set by elaser = "1.16", 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 ##### 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"
displ="0.01"
degree="2"
elaser="1.16" elaserunit="eV"
temp="298.15"

<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.0312"
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 first


When the run has been compleded, move to the directory first.

$cd first  ### 3. Results ##### 3.1) Output files of the Raman calculation In addition to the basic output files already mentioned in Electronic-structure calculations, we requested a phonon calculation at Γ and obtain all the files from it (see Phonon properties of diamond-structure crystals 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 Electronic-structure calculations, 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 first folder 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.1) Analyzing the quality of the results Note: The numerical values and plots presented in this section might vary somewhat in your own results. This fact is due to the loose parameters used in the tutorial run and the (equivalent) results of the phonon eigenvectors for different diagonalization libraries and compilers. The variance of the results, up to equivalences, disappears when the calculations are converged. As can be seen from the introduction, Eq. (3), the intensity/efficiency of a raman line depends critically on the quality of the first derivatives of the components of the dielectric tensor for the given laser energy. So it is important to check the calculation with respect to it. The derivative of the dielectric tensor for a laser energy of 1.16 eV and the fourth phonon mode, written to file RAMAN.OUT [4], looks like the following [ 1.013 + 0.031 i 4.833 + 0.149 i 0.000 + 0.000 i ] [ 4.833 + 0.149 i 1.013 + 0.031 i 0.000 + 0.000 i ] [ 0.000 + 0.000 i 0.000 + 0.000 i -0.242 - 0.008 i ]  but the symmetry of the tensor should be [5] [ 0 a 0 ] [ a 0 0 ] [ 0 0 0 ]  Thus, the tensor has to be symmetric and have zero diagonal values (which is not the case for the numerical results we got). Next, we can plot$\varepsilon^{\alpha\beta}$as a function of the atomic displacement and look how the derivatives were obtained. Invoke the following script from the command line $ PLOT-raman-epsilon.py RAMAN_EPSILON_OC12_MOD004.OUT -0.05 0.05


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.05 and 0.05). The result 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.05 0.05  In this case, the result 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 in the exercise section related to the convergence of the calculation. ##### 3.2) First-order Raman spectrum The actual first-order Raman spectrum can be visualized with help of the script $ PLOT-raman-spectrum.py RAMAN_SPEC_OC12_MOD004.OUT 450 600


We display a simulated spectra in the range between 450 and 600 cm-1 for the fourth 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.

For silicon, a quantity that is often reported is the so-called Raman polarizability $a$. It can be expressed in terms of the derivatives of the dielectric function with respect to the total displacement $u$, $a = 1/4\,\pi\, V_c \left| \partial \varepsilon/\partial u \right|$. The derivatives are reported in the output file RAMAN.OUT. At a laser energy of 1.16 eV, for the derivative of the xy-component of 4.833 + 0.149 $i$ Bohr-1 and the cell volume of 265.3 Bohr-3, we get 102.1 Bohr2 = 28.7 Å2, which is already not too far from the experimental value of 23 $\pm$ 4 Å2 found in Ref.[6].

##### 3.3) 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 1 4  In this case we considered the contribution to the intensity of the Raman spectrum from the xy-component (it might be a different one in your case). 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 (4) \begin{align} \left|\frac{\partial\chi'}{\partial Q}\right|^2 = \left|\frac{1}{4\pi} V_s^{1/2} \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 just a 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 T2g (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.4) Converged results ##### Exercise: How to improve calculations • 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 first, set do = "skip" inside groundstate and getphonon = "fromfile" inside raman. Then, invoke excitingser 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 (as described in Section 2.5), 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-Tc 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 theu\$ that gives the minimum of the potential, written under the sub-header Phonon potential shifted into minimum.
5. For one of the three T2g phonons in the point group Oh, the other two being
[   0   0   a  ]
[   0   0   0  ]
[   a   0   0  ]

and
[   0   0   0  ]
[   0   0   a  ]
[   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).