.. _tutorial_raman: Calculating Raman intensities ============================= PyNAO can be used with the ASE ``vibrations`` module to calculate both static and resonant intensities of the Raman vibration modes. It uses the Placzek approximation implemented in ASE by Prof. M. Walter. When using this module for your work, please cite *M. Walter and M. Moseler, Ab Initio Wavelength-Dependent Raman Spectra: Placzek Approximation and Beyond, J. Chem. Theory Comput. 2020, 16, 1, 576–586* In this example, we will calculate the Raman intensities of a CO2 molecule Calculating the static Raman intensities ---------------------------------------- .. code:: ipython3 import numpy as np from ase import Atoms from ase.calculators.siesta import Siesta from ase.calculators.siesta.siesta_lrtddft import RamanCalculatorInterface from ase.vibrations.raman import StaticRamanCalculator from ase.vibrations.placzek import PlaczekStatic from ase.units import Ry, eV, Ha CO2 = Atoms('CO2', positions=[[-0.009026, -0.020241, 0.026760], [1.167544, 0.012723, 0.071808], [-1.185592, -0.053316, -0.017945]], cell=[20, 20, 20]) # set-up the Siesta parameters CO2.calc = Siesta( mesh_cutoff=250 * Ry, basis_set='DZP', pseudo_qualifier='gga', xc="PBE", energy_shift=(25 * 10**-3) * eV, fdf_arguments={ 'SCFMustConverge': False, 'COOP.Write': True, 'WriteDenchar': True, 'PAO.BasisType': 'split', 'DM.Tolerance': 1e-4, 'DM.MixingWeight': 0.01, "MD.NumCGsteps": 0, "MD.MaxForceTol": (0.02, "eV/Ang"), 'MaxSCFIterations': 10000, 'DM.NumberPulay': 4, 'XML.Write': True, "WriteCoorXmol": True, "DM.UseSaveDM": True,}) name = 'co2' rm = StaticRamanCalculator(CO2, RamanCalculatorInterface, name=name, delta=0.011, exkwargs=dict(label="siesta", jcutoff=7, iter_broadening=0.15, xc_code='LDA,PZ', tol_loc=1e-6, tol_biloc=1e-7, krylov_options={"tol": 1.0e-5, "atol": 1.0e-5})) # save dipole moments from DFT calculation in order to get # infrared intensities as well rm.ir = True rm.run() pz = PlaczekStatic(CO2, name=name) e_vib = pz.get_energies() pz.summary() .. parsed-literal:: Writing co2.eq.pckl, dipole moment = (-0.000052 0.000020 -0.000097) Total number of iterations: 30 Writing co2.0x-.pckl, dipole moment = (-0.022999 -0.000500 -0.000805) Total number of iterations: 30 Writing co2.0x+.pckl, dipole moment = (0.022946 0.000541 0.000611) Total number of iterations: 30 Writing co2.0y-.pckl, dipole moment = (-0.000533 -0.004412 -0.000116) Total number of iterations: 30 Writing co2.0y+.pckl, dipole moment = (0.000488 0.004502 -0.000077) Total number of iterations: 30 Writing co2.0z-.pckl, dipole moment = (-0.000745 -0.000018 -0.004569) Total number of iterations: 30 Writing co2.0z+.pckl, dipole moment = (0.000678 0.000041 0.004399) Total number of iterations: 30 Writing co2.1x-.pckl, dipole moment = (0.011596 0.000281 0.000263) Total number of iterations: 30 Writing co2.1x+.pckl, dipole moment = (-0.011393 -0.000232 -0.000437) Total number of iterations: 30 Writing co2.1y-.pckl, dipole moment = (0.000231 0.002235 -0.000086) Total number of iterations: 30 Writing co2.1y+.pckl, dipole moment = (-0.000311 -0.002257 -0.000108) Total number of iterations: 30 Writing co2.1z-.pckl, dipole moment = (0.000287 0.000057 0.002160) Total number of iterations: 30 Writing co2.1z+.pckl, dipole moment = (-0.000439 0.000007 -0.002387) Total number of iterations: 30 Writing co2.2x-.pckl, dipole moment = (0.011357 0.000273 0.000251) Total number of iterations: 30 Writing co2.2x+.pckl, dipole moment = (-0.011644 -0.000239 -0.000450) Total number of iterations: 30 Writing co2.2y-.pckl, dipole moment = (0.000278 0.002249 -0.000085) Total number of iterations: 30 Writing co2.2y+.pckl, dipole moment = (-0.000258 -0.002249 -0.000105) Total number of iterations: 30 Writing co2.2z-.pckl, dipole moment = (0.000340 0.000061 0.002157) Total number of iterations: 30 Writing co2.2z+.pckl, dipole moment = (-0.000338 0.000010 -0.002365) Total number of iterations: 30 ------------------------------------- Mode Frequency Intensity # meV cm^-1 [0.1A^4/amu] ------------------------------------- 0 0.0 0.0 11.20 1 0.0 0.0 12.43 2 17.1 137.6 1.93 3 20.7 166.9 6.13 4 35.2 284.1 8.87 5 76.7 618.4 0.07 6 78.0 628.9 0.35 7 173.3 1397.4 212.10 8 309.9 2499.5 0.01 ------------------------------------- .. code:: ipython3 from ase.vibrations.infrared import InfraRed # finite displacement for vibrations ir = InfraRed(CO2, name=name) ir.run() ir.summary() .. parsed-literal:: ------------------------------------- Mode Frequency Intensity # meV cm^-1 (D/Å)^2 amu^-1 ------------------------------------- 0 34.2i 275.6i 0.0026 1 29.7i 239.4i 0.0014 2 17.1 137.6 0.0001 3 20.7 166.9 0.0004 4 35.2 284.1 0.0071 5 76.7 618.4 0.5261 6 78.0 628.9 0.5142 7 173.3 1397.4 0.0008 8 309.9 2499.5 13.9984 ------------------------------------- Zero-point energy: 0.355 eV Static dipole moment: 0.001 D Maximum force on atom in `equilibrium`: 1.0566 eV/Å Summary ------- Let’s summarized our results for Infra-red and Raman intensities in a clear table. +-----------------------+---------+---------+-------------+--------------+ | Quantity | Method | 1 | 2 | 3 asymmetric | | | | bending | stretching | streching | +=======================+=========+=========+=============+==============+ | :math:`\omega` | exp | 667.00 | 1330.00 | 2349.00 | | (cm\ :math:`^{-1}`) +---------+---------+-------------+--------------+ | | ASE | 623.65 | 1397.40 | 2499.50 | | +---------+---------+-------------+--------------+ | | QE | 608.45 | 1271.13 | 2223.67 | +-----------------------+---------+---------+-------------+--------------+ | IR :math:`(D/A)^{2}` | ASE | 0.52 | 0.00 | 14.00 | | amu\ :math:`^{-1}` +---------+---------+-------------+--------------+ | | QE | 0.45 | 0.02 | 12.33 | +-----------------------+---------+---------+-------------+--------------+ | Raman (:math:`A^{4}` | ASE | 0.02 | 21.21 | 0.00 | | amu\ :math:`^{-1}`) +---------+---------+-------------+--------------+ | | QE | 0.00 | 23.82 | 0.00 | +-----------------------+---------+---------+-------------+--------------+ The intensities we calculated are given by the ASE rows, experimental values are given for the infra-red intensities values as well as Quantum Expresso calculations. More details can be found in Chapter 7 of M. Barbry PhD thesis. Resonant Raman calculations --------------------------- As explain in M. Walter paper, the Placzek is good enough to get the resonant Raman intensities, for this one can pass the parameter ``omega`` to the ``RamanCalculatorInterface`` object. But first lets find the frequency we want to calculate for. In order to do this, we will calculate the polarizability of the CO2 molecule. .. code:: ipython3 from ase.calculators.siesta.siesta_lrtddft import SiestaLRTDDFT LRTDDFT = SiestaLRTDDFT(label="siesta", jcutoff=7, iter_broadening=0.15, xc_code='LDA,PZ', tol_loc=1e-6, tol_biloc=1e-7, krylov_options={"tol": 1.0e-5, "atol": 1.0e-5}) # run siesta LRTDDFT.get_ground_state(CO2, mesh_cutoff=250 * Ry, basis_set='DZP', pseudo_qualifier='gga', xc="PBE", energy_shift=(25 * 10**-3) * eV, fdf_arguments={ 'SCFMustConverge': False, 'COOP.Write': True, 'WriteDenchar': True, 'PAO.BasisType': 'split', 'DM.Tolerance': 1e-4, 'DM.MixingWeight': 0.01, "MD.NumCGsteps": 0, "MD.MaxForceTol": (0.02, "eV/Ang"), 'MaxSCFIterations': 10000, 'DM.NumberPulay': 4, 'XML.Write': True, "WriteCoorXmol": True, "DM.UseSaveDM": True,}) # Get polarizability freq = np.arange(0.0, 25.0, 0.5) pmat = LRTDDFT.get_polarizability(freq) .. parsed-literal:: Total number of iterations: 1756 .. code:: ipython3 import matplotlib.pyplot as plt h = 9 w = 16*h/9 fig = plt.figure(1, figsize=(w, h)) iax = 1 for i in range(3): for j in range(3): ax = fig.add_subplot(3, 3, iax) ax.plot(freq, pmat[i, j, :].imag, linewidth=3) iax += 1 fig.tight_layout() plt.show() .. image:: output_6_0.png Let’s focus on the excitation around 15 eV .. code:: ipython3 # Get polarizability freq = np.arange(14.0, 16.0, 0.01) pmat2 = LRTDDFT.get_polarizability(freq) .. parsed-literal:: Total number of iterations: 8596 .. code:: ipython3 import matplotlib.pyplot as plt h = 9 w = 16*h/9 fig = plt.figure(1, figsize=(w, h)) iax = 1 for i in range(3): for j in range(3): ax = fig.add_subplot(3, 3, iax) ax.plot(freq, pmat2[i, j, :].imag, linewidth=3) iax += 1 fig.tight_layout() plt.show() .. image:: output_9_0.png .. code:: ipython3 freq_max_int = freq[np.argmax(pmat2[2, 2, :])] print("excitation frequency along zz is ", freq_max_int) .. parsed-literal:: excitation frequency along zz is 14.639999999999986 .. code:: ipython3 name = 'co2-res' rm = StaticRamanCalculator(CO2, RamanCalculatorInterface, name=name, delta=0.011, exkwargs=dict(omega=freq_max_int, label="siesta", jcutoff=7, iter_broadening=0.15, xc_code='LDA,PZ', tol_loc=1e-6, tol_biloc=1e-7, krylov_options={"tol": 1.0e-5, "atol": 1.0e-5}, verbose=5)) # save dipole moments from DFT calculation in order to get # infrared intensities as well rm.ir = True rm.run() pz = PlaczekStatic(CO2, name=name) e_vib = pz.get_energies() pz.summary() .. parsed-literal:: Writing co2-res.eq.pckl, dipole moment = (-0.000026 0.000020 -0.000117) Total number of iterations: 45 Writing co2-res.0x-.pckl, dipole moment = (-0.023009 -0.000500 -0.000805) Total number of iterations: 46 Writing co2-res.0x+.pckl, dipole moment = (0.022946 0.000541 0.000611) Total number of iterations: 46 Writing co2-res.0y-.pckl, dipole moment = (-0.000533 -0.004412 -0.000116) Total number of iterations: 49 Writing co2-res.0y+.pckl, dipole moment = (0.000488 0.004502 -0.000077) Total number of iterations: 49 Writing co2-res.0z-.pckl, dipole moment = (-0.000745 -0.000018 -0.004569) Total number of iterations: 49 Writing co2-res.0z+.pckl, dipole moment = (0.000678 0.000041 0.004399) Total number of iterations: 49 Writing co2-res.1x-.pckl, dipole moment = (0.011596 0.000281 0.000263) Total number of iterations: 47 Writing co2-res.1x+.pckl, dipole moment = (-0.011393 -0.000232 -0.000437) Total number of iterations: 46 Writing co2-res.1y-.pckl, dipole moment = (0.000231 0.002235 -0.000086) Total number of iterations: 49 Writing co2-res.1y+.pckl, dipole moment = (-0.000311 -0.002257 -0.000108) Total number of iterations: 49 Writing co2-res.1z-.pckl, dipole moment = (0.000287 0.000057 0.002160) Total number of iterations: 49 Writing co2-res.1z+.pckl, dipole moment = (-0.000439 0.000007 -0.002387) Total number of iterations: 49 Writing co2-res.2x-.pckl, dipole moment = (0.011357 0.000273 0.000251) Total number of iterations: 46 Writing co2-res.2x+.pckl, dipole moment = (-0.011644 -0.000239 -0.000450) Total number of iterations: 47 Writing co2-res.2y-.pckl, dipole moment = (0.000278 0.002249 -0.000085) Total number of iterations: 49 Writing co2-res.2y+.pckl, dipole moment = (-0.000258 -0.002249 -0.000105) Total number of iterations: 49 Writing co2-res.2z-.pckl, dipole moment = (0.000340 0.000061 0.002157) Total number of iterations: 49 Writing co2-res.2z+.pckl, dipole moment = (-0.000338 0.000010 -0.002365) Total number of iterations: 49 ------------------------------------- Mode Frequency Intensity # meV cm^-1 [100A^4/amu] ------------------------------------- 0 0.0 0.0 0.32 1 0.0 0.0 0.40 2 17.1 137.6 0.35 3 20.7 166.9 0.20 4 35.2 284.1 0.25 5 76.7 618.4 0.02 6 78.0 628.9 0.02 7 173.3 1397.4 318.36 8 309.9 2499.3 0.02 ------------------------------------- .. code:: ipython3 res = 100*318.36 static = 0.1*212.10 print(static, res, res/static) .. parsed-literal:: 21.21 31836.0 1500.9900990099009 Raman intensities enhancement ----------------------------- Let’s analyse the enhancement of the resonant Raman intensity compared to the static one for the mode at :math:`\omega = 1397.4` cm\ :math:`^{-1}`. The intensities are (taking into account the scaling factor in the summary) - static: :math:`I_{sta} = 21.21` A\ :math:`^{4}`/amu - resonant: :math:`I_{res} = 31836.0` A\ :math:`^{4}`/amu Leading to an enhancement of around 1500