by Stefan Kontur & Pasquale Pavone for exciting neon
Purpose: In this tutorial, you will learn how to calculate Raman polarizability tensors 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
 Lattice dynamics of diamond and zincblendestructure 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.
 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.
Note: The symbol $ at beginnings of lines in code segments below indicates the shell prompt.
1. Theoretical background
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. 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 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) 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 zdirection).
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
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" 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
$ 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 mode4
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 mode4 in order to analyze the results.
$ cd mode4
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 zincblendestructure crystals, with method = "sc", 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 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 mode4 of our example.

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 zdirection (OPT001).
Note: The order of the phonon modes can depend on the choice of the computational parameters. For instance, the optical mode polarized along the zdirection 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 zdirection
[ 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
$ PLOTramanepsilon.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 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.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) Firstorder Raman spectrum
The actual firstorder Raman spectrum can be visualized using the command
$ PLOTramanspectrum.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 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.
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)Furthermore, for the OPT001 phonon mode the effective mass is
(15)Therefore, the only non vanisching component of the frequencydependent Raman mode polarizability tensor for the OPT001 phonon mode can be written as
(16)where we explicitly include the label RPA denoting the use of the randomphase 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)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 Å^{2}
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
$ PLOTramanresonance.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)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
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 mode4, set do = "skip" inside groundstate and getphonon = "fromfile" inside raman. Then, invoke exciting_smp 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, 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 ] for an optical phonon polarized along the ydirection
[ a 0 0 ]
[ 0 0 0 ]
[ 0 0 a ] for an optical phonon polarized along the xdirection
[ 0 a 0 ]