# The Raman spectrum and non-linear optical susceptibility {#sec:dfptraman}

Within DFPT framework, 3rd order non-linear response functions can be
calculated from the 1st order perturbation wavefunctions using \"2n+1\"
theorem \[[VGG05](#VGG05){.citation}, [BGCG01](#BGCG01){.citation}\]. In
this section, we present how to compute the non-linear optical
susceptibility and the Raman tensor based on \"2n+1\" theorem. The
former is related to 3rd order derivatives of total energy with respect
to three electric field directions, and the latter is calculated from
3rd order derivatives of total energy with respect to two atomic
displacements and an electric field direction. By combining this two
non-linear response functions we can compute the Raman susceptibilities
including the non-analytical correction due to LO-TO effect. At last,
these susceptibilities are used to generate the Raman intensities and
spectra.

:::: note
::: title
Note
:::

Raman related properties are only well-defined for non-metallic systems.
::::

RESCU currently requires three types of wavefunctions to perform a raman
calculations:

-   Unperturbed wavefunctions which are computed and stored at
    `dfptprep`{.interpreted-text role="ref"}
-   Perturbed wavefunctions with respect to electric field perturbation.
    Although the response of the system to electric field perturbation
    is calculated in `ioncd`{.interpreted-text role="ref"}, saving the
    perturbed wavefunctions was left as an option
-   Perturbed wavefunctions with respect to atomic displacements.
    Although the response of the system to atomic displacements is
    calculated in `dmandbe`{.interpreted-text role="ref"}, saving the
    perturbed wavefunctions was left as an option

In addition, currently RESCU only support time-reversal symmetry for
calculating the perturbed wavefunctions and self-consistent potentials.
Therefor, to prepare the input data for a Raman calculation, it is
necessary to run `ioncd`{.interpreted-text role="ref"} and
`dmandbe`{.interpreted-text role="ref"} examples with following
modifications to the input files:

``` matlab
option.saveWavefunction = 1
symmetry.timereversal = 1      
symmetry.spacesymmetry = 0    
symmetry.pointsymmetry = 0                                  
```

:::: warning
::: title
Warning
:::

Saving unperturbed as well as perturbed wavefunctions may require a
large disk space.
::::

The Raman spectroscopy depends on following properties of the used
light, which are determined as the input of calculations:

-   Direction of incoming photon which is set by either
    `dfpt.qpointKdirect` in direct coordinates or by
    `dfpt.qpointCartesian` in cartesian coordinates.
-   Polarization of incoming photon that is determined by
    `dfpt.ramanIncPhotonPol`
-   Polarization of outgoing photon that is determined by
    `dfpt.ramanOutPhotonPol`

RESCU assumes that measurement is performed with a 90 degree setup
regarding the direction of outgoing photon [\[PWLM\]](#PWLM){.citation}.
Therefore, these keywords work together to define the incoming and
outgoing photons in Raman spectroscopy measurements. They accept a
$N_{directions}\times 3$ matrix where each row corresponds to a
spectroscopy measurement. By default, RESCU assumes three Raman
measurements with ($\hat{z},\hat{x},\hat{y}$) for incoming photon and
($\hat{z},\hat{x},\hat{y}$) for outgoing photon polarization values.

:::: note
::: title
Note
:::

Use `dfpt.qpointCartesian` (or `dfpt.qpointKdirect`),
`dfpt.ramanIncPhotonPol` and `dfpt.ramanOutPhotonPol` to specify the
properties of photons used in the Raman spectroscopy measurement.
::::

:::: important
::: title
Important
:::

Currently only LDA exchange-correlation functional can be used with
\"2n+1\" theorem. Although previous steps of DFPT calculations can still
be performed with another XC-functional type, it should be directly set
to LDA for the Raman calculation. Final results should be checked
carefully for possible inconsistencies.
::::

Following input file can be used for DFPT Raman simulation:

``` matlab
info.calculationType = 'dfpt-raman'
info.savepath = './results/bn_real_raman'
rho.in = './results/bn_real_scf.mat'
psi.in = './results/bn_real_scf.h5'
dfpt.phononData = './results/bn_real_phonon.mat'
dfpt.dielectricData = './results/bn_real_dielectric.mat'
dfpt.asr = [1,1,1]
dfpt.qpointCartesian = [1 0 0;1 0 0]   % Two measurements both having incoming photon in x-direction
dfpt.ramanIncPhotonPol = [0,1,0;0,0,1] % Incoming polarizations for 1st and 2nd measurements in y and z directions
dfpt.ramanOutPhotonPol = [0,1,0;0,1,0] % Outgoing polarizations for 1st and 2nd measurements in y and y directions
units.length = 'Angstrom'              
dfpt.ramanIncPhotonWav = 2660          % Wavelength of incoming photon used for the Raman spectroscopy
dfpt.ramanTemperature = 300            % Temperature of the Raman spectroscopy in Kelvin
dfpt.opticDamping = 1e-4
```

The keywords are explained as

-   `info.calculationType` is a string determining type of the
    calculation;
-   `info.savepath` points to the file where the output data of
    calculation is stored;
-   `rho.in` this DFPT calculation loads the ground state electronic
    density at the beginning of the computation, and this keyword points
    to the file which contains the output of the ground state
    calculation, so that RESCU can load corresponding electronic density
    data;
-   `psi` this DFPT calculation also loads the ground state
    wavefunctions at the beginning of the computation, and this keyword
    points to the file which contains the stored ground state
    wavefunctions;
-   `dfpt.dielectricData` points to the path of the file containing the
    first-order perturbed wavefunctions with respect to external
    electric fields and the ion-clamped dielectric tensor;
-   `dfpt.phononData` points to the path of the file containing the
    first-order perturbed wavefunctions with respect to atomic
    displacements and the dynamical matrix and Born effective charges;
-   `dfpt.asr` it is a $1\times 3$ vector where the first and third
    entries determines if acoustic sum rule will be applied to the
    dynamical matrix and the Raman tensor, respectively. The second
    entry determines if charge neutrality will be enforced on Born
    effective charges. RESCU assumes that all sum rules are active by
    default;
-   `dfpt.qpointCartesian` each row determines a direction of incident
    electric field for which the Raman related properties are
    calculated. By default, RESCU uses the three Cartesian directions;
-   `dfpt.ramanIncPhotonPol` polarization of incoming photon in the
    Raman measurement in Cartesian coordinates. It is a
    $N_{direction}\times 3$ matrix where each row corresponds to a
    direction given by a row of `dfpt.qpointCartesian`. RESCU assumes
    ($\hat{z},\hat{x},\hat{y}$) as default values.
-   `dfpt.ramanOutPhotonPol` polarization of outgoing photon in the
    Raman measurement in Cartesian coordinates. It is a
    $N_{direction}\times 3$ matrix where each row corresponds to a
    direction given by a row of `dfpt.qpointCartesian`. RESCU assumes
    ($\hat{z},\hat{x},\hat{y}$) as default values.
-   `dfpt.ramanIncPhotonWav` wavelength if incoming photon used in the
    Raman measurement. It can be set in Angstrom if used with
    `units.length` keyword. RESCU assumes a wavelength of 4500 Angstrom
    (450 nm) as the default value.
-   `dfpt.ramanTemperature` temperature of the Raman measurement in
    Kelvin which is used for evaluating phonon distributions according
    to Bose-Einstein statistics. The default value is 300 Kelvin.
-   `dfpt.opticDamping` damping factor for Lorentzian broadening of the
    Raman spectrum in Hartree unit (1 Hartree $\simeq 6579.7$ THz).

RESCU is executed as follows

``` shell
rescu -i bn_real_raman
```

The output data will be stored in `./results/bn_real_raman.mat` where
`dfpt` and `raman` fields hold the main variables of the simulation.

``` matlab
>> load results/bn_real_raman.mat
>> dfpt

      dfpt = 

        struct with fields:

                     BornEC: [3x3x2 double]
                        asr: [1 1 1]
                    dTensor: [3x3x3 double]
                       desc: [1x1 struct]
           dielectricTensor: [3x3 double]
                     dynMat: [6x6x2 double]
             nonlinOpticSus: [3x3x3 double]
               opticDamping: 1.0000e-04
                 phononFreq: [6x2 double]
                 phononMode: {2x1 cell}
               qpointDirect: [2x3 double]
          ramanIncPhotonPol: [2x3 double]
          ramanIncPhotonWav: 2660
          ramanOutPhotonPol: [2x3 double]
           ramanTemperature: 300
                ramanTensor: [3x3x6 double]
                       type: 'raman'

>> raman

      raman = 

        struct with fields:

                   Suscep: {2x1 cell}
                   Tensor: {2x1 cell}
               depolRatio: {[6x1 double]  [6x1 double]}
                     desc: [1x1 struct]
                freqShift: [1x1000 double]
                intensity: {[6x1 double]  [6x1 double]}
          intensityAvgDpl: {[6x1 double]  [6x1 double]}
          intensityAvgPol: {[6x1 double]  [6x1 double]}
          intensityAvgTot: {[6x1 double]  [6x1 double]}
                 spectrum: {[1x1000 double]  [1x1000 double]}
           spectrumAvgDpl: {[1x1000 double]  [1x1000 double]}
           spectrumAvgPol: {[1x1000 double]  [1x1000 double]}
           spectrumAvgTot: {[1x1000 double]  [1x1000 double]}
```

Description of important outputs follows:

-   `Suscep` it is a $N_{direction}\times 1$ cell-array where each cell
    entry stores the ($3\times 3\times 3N_{atom}$) Raman susceptibility
    for corresponding direction calculated from the Raman tensor and
    phonons modes in atomic units.

``` matlab
>> raman.Suscep{1}(:,:,1)

      ans =

         1.0e-07 *

        -0.0000 + 0.0000i  -0.5275 - 0.0000i  -0.5266 - 0.0000i
        -0.5275 - 0.0000i  -0.0000 - 0.0000i  -0.4732 + 0.0000i
        -0.5266 - 0.0000i  -0.4732 + 0.0000i  -0.0000 - 0.0000i   
```

-   `Tensor` is a $N_{direction}\times 1$ cell-array where each cell
    entry stores the ($3\times 3\times 3N_{atom}$) Raman tensor for
    corresponding direction including non-analytical effect of LO
    phonons in atomic units. Sum rule is applied if set by `dfpt.asr`.

``` matlab
>> raman.Tensor{1}(:,:,1)

      ans =

        -0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i
         0.0000 - 0.0000i  -0.0000 - 0.0000i  -0.0878 + 0.0000i
         0.0000 - 0.0000i  -0.0878 + 0.0000i  -0.0000 - 0.0000i   
```

-   `intensity` polarization dependent Raman intensities for all
    available phonon modes. It is a $N_{direction}\times 1$ cell-array
    where each entry stores a $3N_{atom}\times 1$ column vector
    corresponds to intensities of available modes.

``` matlab
>> raman.intensity{1}

      ans =

         1.0e-22 *

               0
               0
               0
          0.1461
          0.0408
          0.0122   
```

-   `intensityAvgPol` polarizated averaged Raman intensities
    ($I_{\parallel}$) for phonon modes where it is assumed that the
    incoming and scattered photons have the same polarizations. It is a
    $N_{direction}\times 1$ cell-array where each entry stores a
    $3N_{atom}\times 1$ column vector corresponds to intensities of
    available modes.

``` matlab
>> raman.intensityAvgPol{1}

      ans =

         1.0e-14 *

               0
               0
               0
          0.4077
          0.4078
          0.3471   
```

-   `intensityAvgDpl` depolarized averaged Raman intensities
    ($I_{\perp}$) for phonon modes where the incoming and scattered
    photons have perpendicular polarizations. It is a
    $N_{direction}\times 1$ cell-array where each entry stores a
    $3N_{atom}\times 1$ column vector corresponds to intensities of
    available modes.

``` matlab
>> raman.intensityAvgDpl{1}

      ans =

         1.0e-14 *

               0
               0
               0
          0.3058
          0.3058
          0.2603   
```

-   `intensityAvgTot` total (sum of polarized and depolarized) averaged
    Raman intensities for phonon modes. It is a $N_{direction}\times 1$
    cell-array where each entry stores a $3N_{atom}\times 1$ column
    vector corresponds to intensities of available modes.

``` matlab
>> raman.intensityAvgTot{1}

      ans =

         1.0e-14 *

               0
               0
               0
          0.7135
          0.7136
          0.6074
```

-   `depolRatio` depolarization ratio of averaged intensities, i.e.
    $I_{\parallel}/I_{\perp}$. It is a $N_{direction}\times 1$
    cell-array where each entry stores a $3N_{atom}\times 1$ column
    vector corresponds to depolarization ratio of available modes.
-   `freqShift` frequencies of Raman measurement in cm^-1^ (x-axis) at
    which the Raman spectrum is calculated.
-   `spectrum` $N_{direction}\times 1$ cell-array where each entry
    stores polarization dependent Raman spectrum, and can be plotted
    against `freqShift`.
-   `spectrumAvgPol` $N_{direction}\times 1$ cell-array where each entry
    stores polarized averaged Raman spectrum. It can be plotted against
    `freqShift`.
-   `spectrumAvgDpl` $N_{direction}\times 1$ cell-array where each entry
    stores depolarized averaged Raman spectrum. It can be plotted
    against `freqShift`.
-   `spectrumAvgTot` $N_{direction}\times 1$ cell-array where each entry
    stores total averaged Raman spectrum. It can be plotted against
    `freqShift`.

Also `dfpt` field holds following important information:

-   `BornEC` Born effective charges.
-   `dynMat` updated dynamical matrix for each incoming direction with
    acoustic sum rule applied if set in the input.
-   `phononFreq` phonons frequencies in Hartree.
-   `phononMode` phonons eigenvectors stored as a
    $N_{directions}\times 1$ cell-array.
-   `ramanTensor` $3\times 3\times 3N_{atom}$ array which stores the raw
    Raman tensor before adding the non-analytical LO phonons
    contribution and sum rule (if applicable).
-   `nonlinOpticSus` non-linear optical susceptibility. It is a
    $3\times 3\times 3$ array referring to three Cartesian coordinates
    [\[VGG05\]](#VGG05){.citation}.
-   `dTensor` it is defined as half on the non-linear optical
    susceptibility and given in pmV^-1^ [\[VGG05\]](#VGG05){.citation}.

RESCU\'s plotting utility can be used to visualize the Raman spectrum:

``` shell
rescu -p ./results/bn_real_raman.mat
```

For 3 first given photon directions in `dfpt.qpointCartesian` (or
`dfpt.qpointKdirect`) polarization dependent and averaged Raman spectra
are plotted and saved to corresponding files.

`fig:bn_rampol`{.interpreted-text role="numref"} shows polarization
dependent Raman spectrum which can also be found in
`bn_real_raman_plot_polarized.fig`.

<figure>
<img src="figs/bn_rampol.png" id="fig:bn_rampol"
alt="figs/bn_rampol.png" />
<figcaption>Polarization dependent Raman spectrum of boron
nitrite.</figcaption>
</figure>

`fig:bn_ramavg`{.interpreted-text role="numref"} shows averaged Raman
spectra which can also be found in `bn_real_raman_plot_averaged.fig`.

<figure>
<img src="figs/bn_ramavg.png" id="fig:bn_ramavg"
alt="figs/bn_ramavg.png" />
<figcaption>Averaged Raman spectrum of boron nitrite.</figcaption>
</figure>

::: {#citations}

[BGCG01]{#BGCG01 .citation-label}

:   Stefano Baroni, Stefano De Gironcoli, Andrea Dal Corso, and Paolo
    Giannozzi.
    `Phonons and related crystal properties from density-functional perturbation theory. <10.1103/RevModPhys.73.515>`{.interpreted-text
    role="doi"} Reviews of Modern Physics 73.2 (2001), p. 515.

[PWLM]{#PWLM .citation-label}

:   S. A. Prosandeev, U. Waghmare, I. Levin, and J. Maslar
    `First-order Raman spectra of $\mathrm{A}{\mathrm{B}}_{1∕2}^{\ensuremath{'}}{\mathrm{B}}_{1∕2}^{\ensuremath{''}}{\mathrm{O}}_{3}$ double perovskites<10.1103/PhysRevB.71.214307>`{.interpreted-text
    role="doi"} Physical Review B 71, 214307 (2005)

[VGG05]{#VGG05 .citation-label}

:   M Veithen, Xavier Gonze, and Ph Ghosez.
    `Nonlinear optical susceptibilities, Raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory. <10.1103/PhysRevB.71.125107>`{.interpreted-text
    role="doi"} Physical Review B 71.12 (2005), p. 125107.
:::
