by Stefan Kontur & Pasquale Pavone for exciting nitrogen
Purpose: In this tutorial, you will learn how to calculate Raman shifts and frequencydependent firstorder scattering intensities using exciting. You will be guided through the example of silicon.
Table of Contents

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.
 EXECUTEsingle.sh: (Bash) shell script for running a single exciting calculation.
 PLOTramanepsilon.py: Python visualization tool for the quality of dielectric function fits.
 PLOTramanspectrum.py: Python visualization tool for Raman spectra.
 PLOTramanresonance.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 PLOTramanepsilon.py has been downloaded in the directory $HOME/Downloads, move it to the directory $EXCITINGTOOLS and make it executable:
$ mv $HOME/Downloads/PLOTramanepsilon.py $EXCITINGTOOLS/PLOTramanepsilon.py
$ chmod u+rwx $EXCITINGTOOLS/PLOTramanepsilon.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)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)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 temperaturedependent population of vibrational levels is then
(3)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,
 groundstate computations and forces,
 the calculation of optical spectra, and
 for the final calculation of the Raman spectrum.
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 exchangecorrelation functional.
2.2) Parameters for the groundstate 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 groundstate 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 siraman and move inside it to run the calculations.
$ mkdir siraman
$ cd siraman
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" 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.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> </input>
Please, remember that after creating the input file input.xml, the path for the species file must be set using the command
$ SETUPexcitingroot.sh
To perform actual calculations, you may run the script EXECUTEsingle.sh and choose a name for the directory in which you run with the calculations.
$ EXECUTEsingle.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 Electronicstructure calculations, we requested a phonon calculation at Γ and obtain all the files from it (see Phonon properties of diamondstructure crystals for details). The groundstate 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 Electronicstructure 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.

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
$ PLOTramanepsilon.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 dielectrictensor 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).
$ PLOTramanepsilon.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) Firstorder Raman spectrum
The actual firstorder Raman spectrum can be visualized with help of the script
$ PLOTramanspectrum.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 firstorder 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 socalled 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 xycomponent of 4.833 + 0.149 $i$ Bohr^{1} and the cell volume of 265.3 Bohr^{3}, we get 102.1 Bohr^{2} = 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
$ PLOTramanresonance.py RAMAN_RESONANCE_OC12_MOD004.OUT 1 4
In this case we considered the contribution to the intensity of the Raman spectrum from the xycomponent (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)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 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.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 groundstate 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.
[ 0 0 a ]
[ 0 0 0 ]
[ a 0 0 ]
[ 0 0 0 ]
[ 0 0 a ]
[ 0 a 0 ]