Quantum ESPRESSO: Calculating IR and Raman Spectra
**************************************************

Quantum ESPRESSO Spectroscopy Crystal IR PLAMS Raman

Infrared (IR) and Raman spectroscopy are powerful techniques used to
get information about the vibrational modes of molecules and solids
and also provide valuable insights into the structure, bonding, and
dynamics of materials. IR spectroscopy measures the absorption of
infrared light, which excites vibrational modes that have a change in
dipole moment. On the other hand, Raman spectroscopy measures the
scattering of light, which excites vibrational modes that have a
change in polarizability. By analyzing the frequencies and intensities
of the absorbed or scattered light, we can identify the vibrational
modes present in the material and gain information about its
properties.

Beryllium oxide (BeO), also known as bromellite, is a wide band gap
(7.46 eV) semiconductor. It has a wurtzite crystal structure (space
group P6_3mc, No. 186) and exhibits interesting optical properties,
making it relevant for applications in catalysis, high-power
electronics, and neutron moderators in nuclear reactors. We obtained
the system's geometry from The Materials Project (ID: mp-2542). From
the computational point of view, this system is relatively simple
(only 4 atoms per unit cell) and has a small number of electrons (20
electrons per unit cell), which leads to faster calculations, making
it suitable for illustrative purposes. BeO has a hexagonal structure
(C_{6v} (6mm) symmetry). Some vibrational modes are active in both IR
and Raman spectroscopy, while others are only Raman active. This
selectivity arises from the relationship between the symmetry of the
vibrational mode and the selection rules governing IR and Raman
activity. Theoretical results for comparison can be found in The
Computational Raman Database (mpid: mp-2542), and the experimental
Raman spectrum is available in the RRUFF Project (RRUFF ID: X050194).

Tip:

  The following paper is a highly recommended resource for anyone
  performing first-principles Raman spectroscopy calculations. It
  provides a solid foundation and practical guidance for accurate and
  efficient computations, complemented by a valuable database of
  theoretical results and corresponding experimental data.G. Petretto,
  S. Dwaraknath, H. P.C. Miranda, *et al.* High-throughput density-
  functional perturbation theory phonons for inorganic materials. *Sci
  Data* **5**, 180065 (2018).
  https://doi.org/10.1038/s41597-023-01988-5


Computational Methodology
=========================

The calculation of IR and Raman spectra using Quantum Espresso (QE)
involves a series of steps that are handled internally by AMS, making
the process transparent to the user. These steps use the QE tools
**pw.x**, **ph.x**, and **dynmat.x**. For detailed information about
these tools, refer to the official Quantum Espresso documentation
online.

Starting from an optimized geometry, a self-consistent field (SCF)
calculation is performed using **pw.x** to determine the electronic
ground state of the system. This step solves the Kohn-Sham equations
iteratively to obtain the converged electron density and ground state
energy. Next, **ph.x** is employed to calculate the phonon frequencies
and eigenvectors at a specific q-point (typically \Gamma point) in the
Brillouin zone using density functional perturbation theory (DFPT) to
compute the dynamical matrix, which captures the interatomic force
constants. Finally, **dynmat.x** is used to post-process the phonon
data, obtaining the IR and Raman intensities, frequencies, and
vibrational modes. This step involves computing the Born effective
charges and the Raman tensors, which describe the coupling of the
vibrational modes to the electric field and the change in
polarizability, respectively. These quantities are essential for
determining the IR and Raman activity of the vibrational modes.

IR and Raman spectra can be calculated using Quantum Espresso through
AMS via two main interfaces: the graphical user interface (GUI) and
the Python library, PLAMS. The GUI presents a user-friendly, visual
approach that is ideal for users performing routine calculations on
simple systems. It simplifies the process of setting up input files
and visualizing results. On the other hand, PLAMS offers a powerful
and flexible Python environment suitable for complex workflows,
parameter sweeps, and integration with other scientific libraries.
This makes it particularly advantageous for advanced users, high-
throughput calculations, or situations that require automation and
scripting, providing greater control and customization of the
computational process. In this tutorial, we will explore both options.
They are completely independent, so you can jump directly to the
section that interests you the most.


Using the GUI
=============

We will use the following structure as basis for the subsequent
calculations:

   Be       0.0000000000       1.5555413800       0.0007009000
   Be       1.3471383600       0.7777706900       2.1887103500
    O       0.0000000000       1.5555413800       1.6512462200
    O       1.3471383600       0.7777706900       3.8392556700
   VEC1    2.6942767100     0.0000000000     0.0000000000
   VEC2   -1.3471383600     2.3333120800     0.0000000000
   VEC3    0.0000000000     0.0000000000     4.3760188900

Hint:

  To build a custom crystal, check out the Tutorial on Building
  Crystals and Slabs

Setting Up the Calculation:

   **1.** Open AMSinput
   **2.** Copy & paste the coordinates above into AMSinput.
   **3.** Switch to [image: QuantumESPRESSOPanel][image]
   **4.** Select **Task → GeometryOptimization**
   **5.** Change the **Wavefunction energy cutoff** to **25 Ryd**
   **6.** Change the **K-Points N** to **4 4 3**

You should get something like this:

   [image]

In the steps above, we specified a *GeometryOptimization* task to
ensure our BeO structure is at a local energy minimum, preventing any
spurious imaginary frequencies in the phonon spectrum. Next, for the
electronic structure part, we set the plane-wave energy cutoff
(ecutwfc) to 25 Ry. While this value is lower than typically
recommended for production calculations (40 Ry), we use it here for
demonstration purposes to reduce computational cost. Then, the
Brillouin zone sampling is done using a Monkhorst-Pack k-point grid of
size *4 × 4 × 3*, ensuring a balanced sampling in all directions.

To specify the pseudopotentials, go to the Pseudopotentials
configuration panel by selecting **Model → Pseudopotentials**. In this
panel, use the file selection areas next to each element under
"Pseudopotentials used:". Clicking on the folder icon will open a file
browser window. Navigate to the "QE" directory and select the
appropriate pseudopotential file for each element.

   [image]

Notice that we have specified the pseudopotentials by selecting the
files directly instead of using a pseudopotential family. This is due
to Quantum Espresso's limitation of only supporting Norm-Conserving
(NC) and LDA-type pseudopotentials for Raman spectra calculations,
which are not included in our standard pseudopotential families. For
this example, we selected the pseudopotential files: *Be.pz-
mt_fhi.UPF* for beryllium and *O.pz-mt_fhi.UPF* for oxygen. The
filenames of these pseudopotentials indicate that they use the Perdew-
Zunger LDA functional (*pz*), the Martins-Troullier generation method
(*mt*), and that they originate from the FHI-aims library (*fhi*).

Finally, we explicitly request the calculation of NormalModes (which
yields the IR spectrum) and the Raman spectrum. To request the IR
spectrum, go to **Properties → IR (Frequencies), VCD** and check the
box next to **Frequencies:** by selecting *Yes*. For the Raman
spectrum, navigate to **Properties → Raman, VROA** and select *Yes* in
the **Raman:** checkbox.

   [image]


Computation & Results
---------------------

We are now ready to start the calculation

   **1.** Select **File → Save As** and give it an appropriate name (e.g. "BeO")
   **2.** Select **File → Run**
   **3.** Wait for the calculation to finish (this should take a few minutes). Click Yes when asked to load the optimized geometry in AMSinput.

When you run a calculation, AMSjobs should come to the forefront, and
your job will appear at the top of the list. On the right side, you'll
see the job running, indicated by a gear icon. While the simulation is
in progress, the AMSjobs window will display the progress based on the
logfile. Once the job is complete, the job status icon will change to
a full circle to indicate that it is ready.

   [image]

The AMSspectra module allows you to analyze both the IR and Raman
spectra:

   In AMSjobs (with your job selected) or AMSinput
   Start the AMSspectra module by clicking (**SCM → Spectra**)

You should get something like this (note that the figure below uses a
width of 10.0\text{cm}^{-1}, which may differ from your default
value):

   [image]

On the left side, there is a 3D representation of the BeO molecule,
displaying the movement of atoms according to the selected normal
mode. The right panel presents the calculated Raman spectrum. At the
bottom, you will find a table that provides detailed information about
each vibrational mode, including its symmetry, frequency, and Raman
intensity.

Note:

  You can click on a vibrational mode in the bottom table or the Raman
  spectrum to see the atoms moving according to that mode's
  vibrational frequency.

Note:

  To view the overlay arrows that indicate the displacement vectors
  for a specific vibrational mode, click on **Play → Displacement
  Vectors**. You can adjust the size of these arrows by using Ctrl + L
  to make them larger or Ctrl + K to make them smaller. Note that the
  figure above uses larger arrows than the default size.

Note:

  You can change several parameters related to the spectrum and its
  broadening using the text boxes at the bottom of the spectrum. You
  can modify the scaling factor and offset (useful for DFT functional
  correction factors), as well as the shape (such as Gaussian or
  Lorentzian) and the width.

To see the IR spectrum go to **Spectra → Vibration**. This will
refresh the window to display the IR spectrum results in a layout
similar to that of the Raman spectrum.

   [image]

As we said before, our system BeO has a hexagonal structure with
symmetry C_{6v} (6mm) symmetry. So, we can verify if the transitions
found in our calculations are in agreement with the table of
characters of the point group C_{6v}


Character table for point group C_{6v}
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| C_{6v}      | E           | 2C_6(z)     | 2C_3(z)     | C_2(z)      | 3σ_v        | 3σ_d        | lin, rot    | quadrat     |
|=============|=============|=============|=============|=============|=============|=============|=============|=============|
| A_1         | +1          | +1          | +1          | +1          | +1          | +1          | z           | x^2+y^2,    |
|             |             |             |             |             |             |             |             | z^2         |
+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| A_2         | +1          | +1          | +1          | +1          | -1          | -1          | R_z         | --          |
+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| B_1         | +1          | -1          | +1          | -1          | +1          | -1          | --          | --          |
+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| B_2         | +1          | -1          | +1          | -1          | -1          | +1          | --          | --          |
+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| E_1         | +2          | +1          | -1          | -2          | 0           | 0           | (x, y)      | (xz, yz)    |
|             |             |             |             |             |             |             | (R_x, R_y)  |             |
+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| E_2         | +2          | -1          | -1          | +2          | 0           | 0           | --          | (x^2-y^2,   |
|             |             |             |             |             |             |             |             | xy)         |
+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+

Note:

  **Infrared selection rules:**If a vibration results in a change in
  the dipole moment, it is IR-active. These vibrational modes have the
  same symmetry as the x, y, and z axes.**Raman selection rules:**If a
  vibration results in a change in the polarizability, it is Raman-
  active. These vibrational modes have the same symmetry as quadratic
  functions like xy, xz, and x^2-y^2.

From the characters table above, we can anticipate the IR and Raman
activity of the BeO vibrational modes. This table indicates that modes
with A_1 and E_1 symmetry should be IR-active, while modes with A_1,
E_1, and E_2 symmetry should be Raman-active. Notice this analysis is
consistent with the results obtained from our calculations, as shown
in the IR and Raman spectra plots above (see the columns symmetry and
intensity at the bottom table).


Using Python (PLAMS)
====================

The full PLAMS workflow is now maintained as a scripting example:

* Quantum ESPRESSO IR and Raman spectra for BeO

It contains both downloadable notebook/script files and the rendered
step-by-step Python workflow.

See also:

* Quantum ESPRESSO phonon band structure and DOS for BeO for the
  related phonon workflow.


Further Challenge
=================

To further challenge your understanding, let's calculate the IR and
Raman spectra of a larger, and more complex system Bi_2Y_2O_6 (10
atoms and 52 electrons). You can find theoretical results for
comparison in The Computational Raman Database (mpid: mp-768221).

   Y       0.000000029801245     -0.000000017205758     11.584484628486793
   Y      -0.000000000000000     -0.000000000000000      4.118754113258331
   Bi      0.000000029801245     -0.000000017205758     14.817148949424018
   Bi     -0.000000029801246      0.000000017205756      7.351418334652482
   O       1.327040764604285     -1.891211026034423      8.229234902990678
   O       1.653083765594469     -0.170635326573514      5.740657998219140
   O      -0.974316410415964     -1.346294872078748      5.740657998219141
   O      -2.301357234622741     -0.203645500808650      8.229234902990676
   O       0.974316380614717      2.094856578460344      8.229234902990676
   O      -0.678767355178506      1.516930198652259      5.740657998219140
   VEC1     2.980     1.721     4.977
   VEC2    -2.980     1.721     4.977
   VEC3     0.000    -3.441     4.977

   [image]

   [image]
