Coverage for gpaw/elph/raman_data.py: 18%
78 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1"""
2Calculates Raman matrices from Raman tensor
3"""
4import numpy as np
5from typing import Tuple
7from ase.units import invcm
8from ase.utils.filecache import MultiFileJSONCache
9from gpaw.typing import ArrayND
12def gaussian(w, sigma):
13 g = 1.0 / (np.sqrt(2.0 * np.pi) * sigma) * np.exp(-(w**2) / (2 * sigma**2))
14 return g
17class RamanData:
18 def __init__(self, raman_name="Rlab", gridspacing=1.0) -> None:
19 """Raman Spectroscopy data.
21 Parameters
22 ----------
23 raman_name: str
24 Name of Rlab file cache. Default 'Rlab'
25 gridspacing: float
26 grid spacing in cm^-1, default 1 rcm
27 """
28 # JSON cache
29 self.raman_cache = MultiFileJSONCache(raman_name)
30 self.wph_w = self.raman_cache["phonon_frequencies"] # eV
31 self.wrl_w = self.raman_cache["frequency_grid"] # eV
33 gridmax = self.wph_w[-1] + 6.3e-3 # highest freq + 50rcm
34 self.gridev_w = np.linspace(
35 0.0,
36 gridmax,
37 num=int(gridmax / (gridspacing * invcm) + 1),
38 ) # eV
40 def calculate_raman_intensity(self, d_i, d_o, T=300, sigma=5.0):
41 """
42 Calculates Raman intensities from Raman tensor.
44 Returns bare `|R|^2` and and Bose occupation weighted values
46 Parameters
47 ----------
48 d_i: int
49 Incoming polarization
50 d_o: int
51 Outgoing polarization
52 T: float
53 Temperature in Kelvin
54 sigma: float
55 Gaussian line width in rcm, default 5rcm
56 """
58 KtoeV = 8.617278e-5
59 sigmaev = sigma * invcm # 1eV = 8.065544E3 rcm
61 # grid for spectrum with extra space and sigma/5 spacing
62 if self.gridev_w is None:
63 self.gridev_w = np.linspace(
64 0.0,
65 self.wph_w[-1] * 1.1,
66 num=int(self.wph_w[-1] / (sigmaev / 5) + 100),
67 )
69 int_bare = np.zeros_like(self.gridev_w, dtype=float)
70 int_occ = np.zeros_like(self.gridev_w, dtype=float)
72 # Note: Resonant Raman does not have a frequency grid
73 R_lw = self.raman_cache[f"{'xyz'[d_i]}{'xyz'[d_o]}"]
75 for l in range(len(self.wph_w)):
76 occ = 1.0 / (np.exp(self.wph_w[l] / (KtoeV * T)) - 1.0) + 1.0
77 delta_w = gaussian((self.gridev_w - self.wph_w[l]), sigmaev)
78 R2_w = np.abs(R_lw[l])**2
79 R2d_w = R2_w * delta_w
80 int_bare += R2d_w
81 int_occ += occ / self.wph_w[l] * R2d_w
83 return int_bare, int_occ
85 def calculate_raman_spectrum(self,
86 entries, T=300, sigma=5.0
87 ) -> Tuple[ArrayND, ArrayND]:
88 """
89 Calculates Raman intensities from Raman tensor.
91 Returns Raman shift in eV, bare `|R|^2`
92 and Bose occupation weighted `|R|^2` values
94 Parameters
95 ----------
96 entries: str, list
97 Sting or list of strings with desired polarisaitons
98 For example: ["xx", "yy", "xy", "yx"]
99 T: float
100 Temperature in Kelvin
101 sigma: float
102 Gaussian line width in rcm, default 5rcm
103 """
104 if isinstance(entries, str):
105 entries = [
106 entries,
107 ]
109 spectrum_w = np.zeros_like(self.gridev_w)
111 polnum = {"x": 0, "y": 1, "z": 2}
112 for entry in entries:
113 d_i = polnum[entry[0]]
114 d_o = polnum[entry[1]]
115 (int_bare, _) = self.calculate_raman_intensity(d_i, d_o, T, sigma)
116 spectrum_w += int_bare
118 return self.gridev_w, spectrum_w
120 @classmethod
121 def plot_raman(
122 cls,
123 figname,
124 grid_w,
125 spectra_nw,
126 labels_n=None,
127 relative=False,
128 ):
129 """Plots a given Raman spectrum.
131 Parameters
132 ----------
133 figname: str
134 Filename for figure.
135 grid_w:
136 Frequency grid of spectrum in eV
137 spectra_nw: np.ndarray
138 Raman spectra to plot
139 labels_n: list
140 Labels for the legend
141 relative: bool
142 If true, normalize each spectrum to 1
143 Default is False
144 """
146 from scipy import signal
147 import matplotlib.pyplot as plt
149 if not isinstance(spectra_nw, np.ndarray):
150 spectra_nw = np.asarray(spectra_nw)
152 ylabel = "Intensity (arb. units)"
153 if relative:
154 ylabel = "I/I_max"
155 spectra_nw /= np.max(spectra_nw, axis=1)[:, np.newaxis]
156 else:
157 spectra_nw /= np.max(spectra_nw)
159 nspec = spectra_nw.shape[0]
160 if labels_n is None:
161 labels_n = nspec * [None]
163 grid_w = cls.eVtorcm(grid_w)
165 for i in range(nspec):
166 peaks = signal.find_peaks(spectra_nw[i])[0]
167 locations = np.take(grid_w, peaks)
168 intensities = np.take(spectra_nw[i], peaks)
170 plt.plot(grid_w, spectra_nw[i], label=labels_n[i])
172 for j, loc in enumerate(locations):
173 if intensities[j] / np.max(intensities) > 0.05:
174 plt.axvline(x=loc, color="grey", linestyle="--")
176 plt.minorticks_on()
177 if labels_n[0] is not None:
178 plt.legend()
179 plt.title("Raman intensity")
180 plt.xlabel("Raman shift (cm$^{-1}$)")
181 plt.ylabel(ylabel)
182 if relative:
183 plt.yticks([])
185 plt.savefig(figname, dpi=300)
186 plt.clf()
188 @classmethod
189 def eVtorcm(cls, energy):
190 return energy / invcm