Coverage for gpaw/raman/raman.py: 4%
208 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1"""
2Calculates Raman matrices and intensities.
4i -> j -> m -> n
5i, n are valence; j, m are conduction, also i=n in the end
6see https://doi.org/10.1038/s41467-020-16529-6
7"""
9import numpy as np
10from ase.units import invcm, Hartree
13def lorentzian(w, gamma):
14 l = 0.5 * gamma / (np.pi * (w**2 + 0.25 * gamma**2))
15 return l
18def gaussian(w, sigma):
19 g = 1. / (np.sqrt(2. * np.pi) * sigma) * np.exp(-w**2 / (2 * sigma**2))
20 return g
23def calculate_raman(calc, w_ph, w_in, d_i, d_o, resonant_only=False,
24 gamma_l=0.1, momname='mom_skvnm.npy',
25 elphname='gsqklnn.npy', gridspacing=1.0, suffix=None):
26 """
27 Calculates the first order Raman tensor contribution for a given
28 polarization.
30 Parameters
31 ----------
32 calc: GPAW
33 Converged ground state calculation
34 w_ph: np.ndarray, str
35 Array of phonon frequencies in eV, or name of file with them
36 w_in: float
37 Laser energy in eV
38 d_i: int
39 Incoming polarization
40 d_o: int
41 Outgoing polarization
42 gamma_l: float
43 Line broadening in eV, (0.1eV=806rcm)
44 momname: str
45 Name of momentum file
46 elphname: str
47 Name of electron-phonon file
48 gridspacing: float
49 grid spacing in cm^-1
50 suffix: str
51 suffix for Rlab file name. Defaults to XXXnm if not given.
52 """
53 # those won't make sense here
54 assert calc.wfs.gd.comm.size == 1
55 assert calc.wfs.bd.comm.size == 1
56 kd = calc.wfs.kd
58 print(f"Calculating Raman spectrum: Laser frequency = {w_in} eV")
59 if suffix is None:
60 suffix = f"{int(1239.841 / w_in)}nm"
62 # Phonon frequencies
63 if isinstance(w_ph, str):
64 w_ph = np.load(w_ph)
65 assert max(w_ph) < 1. # else not eV units
66 nmodes = len(w_ph)
68 # Set grid
69 w_max = np.round(np.max(w_ph) / invcm + 50, -1) # max of grid in rcm
70 ngrid = int(w_max / gridspacing) + 1
71 w = np.linspace(0., w_max, num=ngrid) * invcm # in eV
73 # Load files
74 mom_skvnm = np.load(momname, mmap_mode='c')
75 g_sqklnn = np.load(elphname, mmap_mode='c') # [s,q=0,k,l,n,m]
77 # Define a few more variables
78 opt = 'optimal' # mode for np.einsum. not sure what is fastest
79 nspins = g_sqklnn.shape[0]
80 nk = g_sqklnn.shape[2]
81 assert nmodes == g_sqklnn.shape[3]
82 assert mom_skvnm.shape[0] == nspins
83 assert mom_skvnm.shape[1] == nk
84 assert mom_skvnm.shape[-1] == g_sqklnn.shape[-1]
86 # valence is ket in ij and bra in mn
87 # ab is in and out polarization
88 # l is the phonon mode and w is the raman shift
90 # XXX: The below can probably be made better by lambdas a lot
91 def _term1(f_vc, E_vc, mom_dnn, elph_lnn, nc, nv):
92 term1_l = np.zeros((elph_lnn.shape[0]), dtype=complex)
93 t1_ij = f_vc * mom_dnn[0, :nv, nc:] / (w_in - E_vc)
94 # t1_ij = f_vc * mom_dnn[0, nc:, :nv].T / (w_in - E_vc)
95 for l in range(nmodes):
96 t1_xx = elph_lnn[l]
97 # t1_mn = (f_vc * mom_dnn[1, :nv, nc:] / (w_in - w_ph[l] - E_vc)).T
98 t1_mn = f_vc.T * mom_dnn[1, nc:, :nv] / (w_in - w_ph[l] - E_vc.T)
99 term1_l[l] += np.einsum('sj,jm,ms', t1_ij, t1_xx[nc:, nc:], t1_mn,
100 optimize=opt)
101 term1_l[l] -= np.einsum('is,ni,sn', t1_ij, t1_xx[:nv, :nv], t1_mn,
102 optimize=opt)
103 return term1_l
105 def _term2(f_vc, E_vc, mom_dnn, elph_lnn, nc, nv):
106 term2_lw = np.zeros((nmodes, ngrid), dtype=complex)
107 t2_ij = f_vc * mom_dnn[0, :nv, nc:] / (w_in - E_vc)
108 t2_xx = mom_dnn[1]
109 for l in range(nmodes):
110 for wi in range(ngrid):
111 t2_mn = f_vc.T * elph_lnn[l][nc:, :nv] / (w[wi] - E_vc.T)
112 term2_lw[l, wi] += np.einsum('sj,jm,ms', t2_ij,
113 t2_xx[nc:, nc:], t2_mn,
114 optimize=opt)
115 term2_lw[l, wi] -= np.einsum('is,ni,sn', t2_ij,
116 t2_xx[:nv, :nv], t2_mn,
117 optimize=opt)
118 return term2_lw
120 def _term3(f_vc, E_vc, mom_dnn, elph_lnn, nc, nv):
121 term3_lw = np.zeros((nmodes, ngrid), dtype=complex)
122 for l in range(nmodes):
123 t3_xx = elph_lnn[l]
124 for wi in range(ngrid):
125 t3_ij = f_vc * mom_dnn[1, :nv, nc:] / (-w_in + w[wi] - E_vc)
126 t3_mn = (f_vc.T * mom_dnn[0, :nv, nc:] / (-w_in - w_ph[l] +
127 w[wi] - E_vc.T))
128 term3_lw[l, wi] += np.einsum('sj,jm,ms', t3_ij,
129 t3_xx[nc:, nc:], t3_mn,
130 optimize=opt)
131 term3_lw[l, wi] -= np.einsum('is,ni,sn', t3_ij,
132 t3_xx[:nv, :nv], t3_mn,
133 optimize=opt)
134 return term3_lw
136 def _term4(f_vc, E_vc, mom_dnn, elph_lnn, nc, nv):
137 term4_lw = np.zeros((nmodes, ngrid), dtype=complex)
138 t4_xx = mom_dnn[0]
139 for l in range(nmodes):
140 for wi in range(ngrid):
141 t4_ij = f_vc * mom_dnn[1, :nv, nc:] / (-w_in + w[wi] - E_vc)
142 t4_mn = (f_vc.T * elph_lnn[l, nc:, :nv] / (w[wi] - E_vc.T))
143 term4_lw[l, wi] += np.einsum('sj,jm,ms', t4_ij,
144 t4_xx[nc:, nc:], t4_mn,
145 optimize=opt)
146 term4_lw[l, wi] -= np.einsum('is,ni,sn', t4_ij,
147 t4_xx[:nv, :nv], t4_mn,
148 optimize=opt)
149 return term4_lw
151 def _term5(f_vc, E_vc, mom_dnn, elph_lnn, nc, nv):
152 term5_l = np.zeros((nmodes), dtype=complex)
153 t5_xx = mom_dnn[0]
154 for l in range(nmodes):
155 t5_ij = f_vc * elph_lnn[l, :nv, nc:] / (-w_ph[l] - E_vc)
156 t5_mn = (f_vc.T * mom_dnn[1, nc:, :nv] / (w_in - w_ph[l] - E_vc.T))
157 term5_l[l] += np.einsum('sj,jm,ms', t5_ij, t5_xx[nc:, nc:], t5_mn,
158 optimize=opt)
159 term5_l[l] -= np.einsum('is,ni,sn', t5_ij, t5_xx[:nv, :nv], t5_mn,
160 optimize=opt)
161 return term5_l
163 def _term6(f_vc, E_vc, mom_dnn, elph_lnn, nc, nv):
164 term6_lw = np.zeros((nmodes, ngrid), dtype=complex)
165 t6_xx = mom_dnn[1]
166 for l in range(nmodes):
167 t6_ij = f_vc * elph_lnn[l, :nv, nc:] / (-w_ph[l] - E_vc)
168 for wi in range(ngrid):
169 t6_mn = (f_vc.T * mom_dnn[0, nc:, :nv] / (-w_in - w_ph[l]
170 + w[wi] - E_vc.T))
171 term6_lw[l, wi] += np.einsum('sj,jm,ms', t6_ij,
172 t6_xx[nc:, nc:], t6_mn,
173 optimize=opt)
174 term6_lw[l, wi] -= np.einsum('is,ni,sn', t6_ij,
175 t6_xx[:nv, :nv], t6_mn,
176 optimize=opt)
177 return term6_lw
179 # -------------------------------------------------------------------------
181 print("Evaluating Raman sum")
182 raman_lw = np.zeros((nmodes, ngrid), dtype=complex)
184 # Loop over kpoints - this is parallelised
185 for kpt in calc.wfs.kpt_u:
186 print(f"Rank {kd.comm.rank}: s={kpt.s}, k={kpt.k}")
188 # Check if we need to add timer-add time reversed kpoint
189 if (calc.symmetry.time_reversal and not
190 np.allclose(kd.ibzk_kc[kpt.k], [0., 0., 0.])):
191 add_time_reversed = True
192 # Currently broken.
193 raise NotImplementedError
194 else:
195 add_time_reversed = False
197 # Limit sums to relevant bands, partially occupied bands are a pain.
198 # So, in principal, partially occupied bands can be initial and
199 # final states, but we need to restrict to a positive deltaE if we
200 # allow this.
201 f_n = kpt.f_n / kpt.weight
202 assert np.isclose(max(f_n), 1.0, atol=0.1)
204 if 1:
205 vs = np.arange(0, len(f_n))
206 cs = np.arange(0, len(f_n))
207 nc = 0
208 nv = len(f_n)
209 else:
210 vs = np.where(f_n >= 0.1)[0]
211 cs = np.where(f_n < 0.9)[0]
212 nv = max(vs) + 1 # VBM+1 index
213 nc = min(cs) # CBM index
215 # Precalculate f * (1-f) term
216 f_vc = np.outer(f_n[vs], 1. - f_n[cs])
218 # Precalculate E-E term
219 E_vc = np.empty((len(vs), len(cs)), dtype=complex)
220 for n in range(len(vs)):
221 E_vc[n] = (kpt.eps_n[cs] - kpt.eps_n[n]) * Hartree + 1j * gamma_l
223 # set weights for negative energy transitions zero
224 if 0:
225 neg = np.where(E_vc[n].real <= 0.)[0]
226 f_vc[n, neg] = 0.
228 # Obtain appropriate part of mom and g arrays
229 mom_dnn = np.ascontiguousarray(mom_skvnm[kpt.s, kpt.k, [d_i, d_o]])
230 assert mom_dnn.shape[0] == 2
231 g_lnn = np.ascontiguousarray(g_sqklnn[kpt.s, 0, kpt.k])
233 # Raman contribution of this k-point
234 this_lw = np.zeros((nmodes, ngrid), dtype=complex)
236 # Resonant term
237 term1_l = _term1(f_vc, E_vc, mom_dnn, g_lnn, nc, nv)
238 # print("Term1: ", np.max(np.abs(term1_l)))
239 this_lw += term1_l[:, None]
241 if not resonant_only:
242 term2_lw = _term2(f_vc, E_vc, mom_dnn, g_lnn, nc, nv)
243 # print("Term2: ", np.max(np.abs(term2_lw)))
244 this_lw += term2_lw
246 term3_lw = _term3(f_vc, E_vc, mom_dnn, g_lnn, nc, nv)
247 # print("Term3: ", np.max(np.abs(term3_lw)))
248 this_lw += term3_lw
250 term4_lw = _term4(f_vc, E_vc, mom_dnn, g_lnn, nc, nv)
251 # print("Term4: ", np.max(np.abs(term4_lw)))
252 this_lw += term4_lw
254 term5_l = _term5(f_vc, E_vc, mom_dnn, g_lnn, nc, nv)
255 # print("Term5: ", np.max(np.abs(term5_l)))
256 this_lw += term5_l[:, None]
258 term6_lw = _term6(f_vc, E_vc, mom_dnn, g_lnn, nc, nv)
259 # print("Term6: ", np.max(np.abs(term6_lw)))
260 this_lw += term6_lw
262 # At the moment we only allow no symmetry, so all k-points same
263 # weight, or time_reversal only. In the later case r -> 2*Re(r)
264 # because gdd-> (gdd)^* for k-> -k
265 if add_time_reversed:
266 raman_lw += 2. * this_lw.real * kpt.weight
267 else:
268 raman_lw += this_lw * kpt.weight
270 # Collect parallel contributions
271 kd.comm.sum(raman_lw)
273 if kd.comm.rank == 0:
274 print('Raman intensities per mode')
275 print('--------------------------')
276 ff = ' Phonon {} with energy = {:4.2f} rcm: {:.4f}'
277 for l in range(nmodes):
278 print(ff.format(l, w_ph[l] / invcm, abs(np.sum(raman_lw[l]))))
280 # Save Raman intensities to disk
281 raman = np.vstack((w / invcm, raman_lw))
282 np.save("Rlab_{}{}_{}.npy".format('xyz'[d_i], 'xyz'[d_o], suffix),
283 raman)
286def calculate_raman_tensor(calc, w_ph, w_in, resonant_only=False,
287 gamma_l=0.1, momname='mom_skvnm.npy',
288 elphname='gsqklnn.npy'):
289 for i in range(3):
290 for j in range(3):
291 calculate_raman(calc, w_ph, w_in, d_i=i, d_o=j,
292 resonant_only=resonant_only, gamma_l=gamma_l,
293 momname=momname, elphname=elphname)
296def calculate_raman_intensity(w_ph, d_i, d_o, suffix, T=300):
297 """
298 Calculates Raman intensities from Raman tensor.
300 Parameters
301 ----------
302 w_ph: np.ndarray, str
303 Array of phonon frequencies in eV, or name of file with them
304 d_i: int
305 Incoming polarization
306 d_o: int
307 Outgoing polarization
308 suffix: str
309 Suffix for Rlab and RI files. Usually XXXnm
310 """
311 KtoeV = 8.617278e-5
312 # Phonon frequencies
313 if isinstance(w_ph, str):
314 w_ph = np.load(w_ph)
315 assert max(w_ph) < 1. # else not eV units
317 # Load raman matrix elements R_lab
318 tmp = np.load("Rlab_{}{}_{}.npy".format('xyz'[d_i], 'xyz'[d_o], suffix))
319 w = tmp[0].real # already in rcm
320 raman_lw = tmp[1:]
321 assert raman_lw.shape[0] == len(w_ph)
323 int_bare = np.zeros_like(w, dtype=float)
324 int_occ = np.zeros_like(w, dtype=float)
325 for l in range(len(w_ph)):
326 occ = 1. / (np.exp(w_ph[l] / (KtoeV * T)) - 1.) + 1.
327 delta = gaussian(w=(w - w_ph[l] / invcm), sigma=5.)
328 R = raman_lw[l] * raman_lw[l].conj()
329 int_bare += R.real * delta
330 int_occ += occ / w_ph[l] * R.real * delta
332 raman = np.vstack((w, int_bare, int_occ))
333 np.save("RI_{}{}_{}.npy".format('xyz'[d_i], 'xyz'[d_o], suffix), raman)
336def plot_raman(figname, RIsuffix, relative=False, w_min=None, w_max=None):
337 """Plots a given Raman spectrum.
339 Parameters
340 ----------
341 figname: str
342 Filename for figure.
343 RIsuffix: str, list
344 Suffix of Raman intensity files to use for plotting. For example
345 "0_1_455nm" for RI_0_1_455nm.npy
346 """
347 from scipy import signal
348 import matplotlib.pyplot as plt
349 import matplotlib.colors as colors
350 import matplotlib.cm as cmx
352 if isinstance(RIsuffix, str):
353 legend = False
354 RI_name = [f"RI_{RIsuffix}.npy"]
355 else: # assume list
356 legend = True
357 RI_name = [f"RI_{name}.npy" for name in RIsuffix]
358 cm = plt.get_cmap('inferno')
359 cNorm = colors.Normalize(vmin=0, vmax=len(RI_name))
360 scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)
362 ylabel = "Intensity (arb. units)"
363 if w_min is None:
364 w_min = 0.
366 peaks = None
367 for i, name in enumerate(RI_name):
368 RI = np.load(name) # load raman intensity
369 if w_max is None:
370 w_max = np.max(RI[0])
371 r = RI[1][np.logical_and(RI[0] >= w_min, RI[0] <= w_max)]
372 w = RI[0][np.logical_and(RI[0] >= w_min, RI[0] <= w_max)]
373 if relative:
374 ylabel = "I/I_max"
375 r /= np.max(r)
377 if peaks is None:
378 peaks = signal.find_peaks(r[np.logical_and(w >= w_min, w <= w_max)
379 ])[0]
380 locations = np.take(w[np.logical_and(w >= w_min, w <= w_max)],
381 peaks)
382 intensities = np.take(r[np.logical_and(w >= w_min, w <= w_max)],
383 peaks)
385 if legend:
386 cval = scalarMap.to_rgba(i)
387 plt.plot(w, r, color=cval, label=RIsuffix[i])
388 else:
389 plt.plot(w, r)
391 for i, loc in enumerate(locations):
392 if intensities[i] / np.max(intensities) > 0.05:
393 plt.axvline(x=loc, color="grey", linestyle="--")
395 plt.minorticks_on()
396 if legend:
397 plt.legend()
398 plt.title("Raman intensity")
399 plt.xlabel("Raman shift (cm$^{-1}$)")
400 plt.ylabel(ylabel)
401 if relative:
402 plt.yticks([])
403 plt.savefig(figname, dpi=300)
404 plt.clf()