Coverage for /builds/ase/ase/ase/vibrations/resonant_raman.py: 65.27%
311 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""Resonant Raman intensities"""
5import sys
6from pathlib import Path
8import numpy as np
10import ase.units as u
11from ase.parallel import paropen, parprint, world
12from ase.vibrations import Vibrations
13from ase.vibrations.raman import Raman, RamanCalculatorBase
16class ResonantRamanCalculator(RamanCalculatorBase, Vibrations):
17 """Base class for resonant Raman calculators using finite differences.
18 """
20 def __init__(self, atoms, ExcitationsCalculator, *args,
21 exkwargs=None, exext='.ex.gz', overlap=False,
22 **kwargs):
23 """
24 Parameters
25 ----------
26 atoms: Atoms
27 The Atoms object
28 ExcitationsCalculator: object
29 Calculator for excited states
30 exkwargs: dict
31 Arguments given to the ExcitationsCalculator object
32 exext: string
33 Extension for filenames of Excitation lists (results of
34 the ExcitationsCalculator).
35 overlap : function or False
36 Function to calculate overlaps between excitation at
37 equilibrium and at a displaced position. Calculators are
38 given as first and second argument, respectively.
40 Example
41 -------
43 >>> from ase.calculators.h2morse import (H2Morse,
44 ... H2MorseExcitedStatesCalculator)
45 >>> from ase.vibrations.resonant_raman import ResonantRamanCalculator
46 >>>
47 >>> atoms = H2Morse()
48 >>> rmc = ResonantRamanCalculator(atoms, H2MorseExcitedStatesCalculator)
49 >>> rmc.run()
51 This produces all necessary data for further analysis.
52 """
53 self.exobj = ExcitationsCalculator
54 if exkwargs is None:
55 exkwargs = {}
56 self.exkwargs = exkwargs
57 self.overlap = overlap
59 super().__init__(atoms, *args, exext=exext, **kwargs)
61 def _new_exobj(self):
62 # XXXX I have to duplicate this because there are two objects
63 # which have exkwargs, why are they not unified?
64 return self.exobj(**self.exkwargs)
66 def calculate(self, atoms, disp):
67 """Call ground and excited state calculation"""
68 assert atoms == self.atoms # XXX action required
69 forces = self.atoms.get_forces()
71 if self.overlap:
72 """Overlap is determined as
74 ov_ij = int dr displaced*_i(r) eqilibrium_j(r)
75 """
76 ov_nn = self.overlap(self.atoms.calc,
77 self.eq_calculator)
78 if world.rank == 0:
79 disp.save_ov_nn(ov_nn)
81 disp.calculate_and_save_exlist(atoms)
82 return {'forces': forces}
84 def run(self):
85 if self.overlap:
86 # XXXX stupid way to make a copy
87 self.atoms.get_potential_energy()
88 self.eq_calculator = self.atoms.calc
89 Path(self.name).mkdir(parents=True, exist_ok=True)
90 fname = Path(self.name) / 'tmp.gpw'
91 self.eq_calculator.write(fname, 'all')
92 self.eq_calculator = self.eq_calculator.__class__(restart=fname)
93 try:
94 # XXX GPAW specific
95 self.eq_calculator.converge_wave_functions()
96 except AttributeError:
97 pass
98 Vibrations.run(self)
101class ResonantRaman(Raman):
102 """Base Class for resonant Raman intensities using finite differences.
103 """
105 def __init__(self, atoms, Excitations, *args,
106 observation=None,
107 form='v', # form of the dipole operator
108 exkwargs=None, # kwargs to be passed to Excitations
109 exext='.ex.gz', # extension for Excitation names
110 overlap=False,
111 minoverlap=0.02,
112 minrep=0.8,
113 comm=world,
114 **kwargs):
115 """
116 Parameters
117 ----------
118 atoms: ase Atoms object
119 Excitations: class
120 Type of the excitation list object. The class object is
121 initialized as::
123 Excitations(atoms.calc)
125 or by reading form a file as::
127 Excitations('filename', **exkwargs)
129 The file is written by calling the method
130 Excitations.write('filename').
132 Excitations should work like a list of ex obejects, where:
133 ex.get_dipole_me(form='v'):
134 gives the velocity form dipole matrix element in
135 units |e| * Angstrom
136 ex.energy:
137 is the transition energy in Hartrees
138 approximation: string
139 Level of approximation used.
140 observation: dict
141 Polarization settings
142 form: string
143 Form of the dipole operator, 'v' for velocity form (default)
144 and 'r' for length form.
145 overlap: bool or function
146 Use wavefunction overlaps.
147 minoverlap: float ord dict
148 Minimal absolute overlap to consider. Defaults to 0.02 to avoid
149 numerical garbage.
150 minrep: float
151 Minimal representation to consider derivative, defaults to 0.8
152 """
154 if observation is None:
155 observation = {'geometry': '-Z(XX)Z'}
157 kwargs['exext'] = exext
158 Raman.__init__(self, atoms, *args, **kwargs)
159 assert self.vibrations.nfree == 2
161 self.exobj = Excitations
162 if exkwargs is None:
163 exkwargs = {}
164 self.exkwargs = exkwargs
165 self.observation = observation
166 self.dipole_form = form
168 self.overlap = overlap
169 if not isinstance(minoverlap, dict):
170 # assume it's a number
171 self.minoverlap = {'orbitals': minoverlap,
172 'excitations': minoverlap}
173 else:
174 self.minoverlap = minoverlap
175 self.minrep = minrep
177 def read_exobj(self, filename):
178 return self.exobj.read(filename, **self.exkwargs)
180 def get_absolute_intensities(self, omega, gamma=0.1, delta=0, **kwargs):
181 """Absolute Raman intensity or Raman scattering factor
183 Parameter
184 ---------
185 omega: float
186 incoming laser energy, unit eV
187 gamma: float
188 width (imaginary energy), unit eV
189 delta: float
190 pre-factor for asymmetric anisotropy, default 0
192 References
193 ----------
194 Porezag and Pederson, PRB 54 (1996) 7830-7836 (delta=0)
195 Baiardi and Barone, JCTC 11 (2015) 3267-3280 (delta=5)
197 Returns
198 -------
199 raman intensity, unit Ang**4/amu
200 """
201 alpha2_r, gamma2_r, delta2_r = self._invariants(
202 self.electronic_me_Qcc(omega, gamma, **kwargs))
203 return 45 * alpha2_r + delta * delta2_r + 7 * gamma2_r
205 @property
206 def approximation(self):
207 return self._approx
209 @approximation.setter
210 def approximation(self, value):
211 self.set_approximation(value)
213 def read_excitations(self):
214 """Read all finite difference excitations and select matching."""
215 if self.overlap:
216 return self.read_excitations_overlap()
218 disp = self._eq_disp()
219 ex0_object = disp.read_exobj()
220 eu = ex0_object.energy_to_eV_scale
221 matching = frozenset(ex0_object)
223 def append(lst, disp, matching):
224 exo = disp.read_exobj()
225 lst.append(exo)
226 matching = matching.intersection(exo)
227 return matching
229 exm_object_list = []
230 exp_object_list = []
231 for a, i in zip(self.myindices, self.myxyz):
232 mdisp = self._disp(a, i, -1)
233 pdisp = self._disp(a, i, 1)
234 matching = append(exm_object_list,
235 mdisp, matching)
236 matching = append(exp_object_list,
237 pdisp, matching)
239 def select(exl, matching):
240 mlst = [ex for ex in exl if ex in matching]
241 assert len(mlst) == len(matching)
242 return mlst
244 ex0 = select(ex0_object, matching)
245 exm = []
246 exp = []
247 r = 0
248 for a, i in zip(self.myindices, self.myxyz):
249 exm.append(select(exm_object_list[r], matching))
250 exp.append(select(exp_object_list[r], matching))
251 r += 1
253 self.ex0E_p = np.array([ex.energy * eu for ex in ex0])
254 self.ex0m_pc = (np.array(
255 [ex.get_dipole_me(form=self.dipole_form) for ex in ex0]) *
256 u.Bohr)
257 exmE_rp = []
258 expE_rp = []
259 exF_rp = []
260 exmm_rpc = []
261 expm_rpc = []
262 r = 0
263 for a, i in zip(self.myindices, self.myxyz):
264 exmE_rp.append([em.energy for em in exm[r]])
265 expE_rp.append([ep.energy for ep in exp[r]])
266 exF_rp.append(
267 [(em.energy - ep.energy)
268 for ep, em in zip(exp[r], exm[r])])
269 exmm_rpc.append(
270 [ex.get_dipole_me(form=self.dipole_form)
271 for ex in exm[r]])
272 expm_rpc.append(
273 [ex.get_dipole_me(form=self.dipole_form)
274 for ex in exp[r]])
275 r += 1
276 # indicees: r=coordinate, p=excitation
277 # energies in eV
278 self.exmE_rp = np.array(exmE_rp) * eu
279 self.expE_rp = np.array(expE_rp) * eu
280 # forces in eV / Angstrom
281 self.exF_rp = np.array(exF_rp) * eu / 2 / self.delta
282 # matrix elements in e * Angstrom
283 self.exmm_rpc = np.array(exmm_rpc) * u.Bohr
284 self.expm_rpc = np.array(expm_rpc) * u.Bohr
286 def read_excitations_overlap(self):
287 """Read all finite difference excitations and wf overlaps.
289 We assume that the wave function overlaps are determined as
291 ov_ij = int dr displaced*_i(r) eqilibrium_j(r)
292 """
293 ex0 = self._eq_disp().read_exobj()
294 eu = ex0.energy_to_eV_scale
295 rep0_p = np.ones((len(ex0)), dtype=float)
297 def load(disp, rep0_p):
298 ex_p = disp.read_exobj()
299 ov_nn = disp.load_ov_nn()
300 # remove numerical garbage
301 ov_nn = np.where(np.abs(ov_nn) > self.minoverlap['orbitals'],
302 ov_nn, 0)
303 ov_pp = ex_p.overlap(ov_nn, ex0)
304 ov_pp = np.where(np.abs(ov_pp) > self.minoverlap['excitations'],
305 ov_pp, 0)
306 rep0_p *= (ov_pp.real**2 + ov_pp.imag**2).sum(axis=0)
307 return ex_p, ov_pp
309 def rotate(ex_p, ov_pp):
310 e_p = np.array([ex.energy for ex in ex_p])
311 m_pc = np.array(
312 [ex.get_dipole_me(form=self.dipole_form) for ex in ex_p])
313 r_pp = ov_pp.T
314 return ((r_pp.real**2 + r_pp.imag**2).dot(e_p),
315 r_pp.dot(m_pc))
317 exmE_rp = []
318 expE_rp = []
319 exF_rp = []
320 exmm_rpc = []
321 expm_rpc = []
322 exdmdr_rpc = []
323 for a, i in zip(self.myindices, self.myxyz):
324 mdisp = self._disp(a, i, -1)
325 pdisp = self._disp(a, i, 1)
326 ex, ov = load(mdisp, rep0_p)
327 exmE_p, exmm_pc = rotate(ex, ov)
328 ex, ov = load(pdisp, rep0_p)
329 expE_p, expm_pc = rotate(ex, ov)
330 exmE_rp.append(exmE_p)
331 expE_rp.append(expE_p)
332 exF_rp.append(exmE_p - expE_p)
333 exmm_rpc.append(exmm_pc)
334 expm_rpc.append(expm_pc)
335 exdmdr_rpc.append(expm_pc - exmm_pc)
337 # select only excitations that are sufficiently represented
338 self.comm.product(rep0_p)
339 select = np.where(rep0_p > self.minrep)[0]
341 self.ex0E_p = np.array([ex.energy * eu for ex in ex0])[select]
342 self.ex0m_pc = (np.array(
343 [ex.get_dipole_me(form=self.dipole_form)
344 for ex in ex0])[select] * u.Bohr)
346 if len(self.myr):
347 # indicees: r=coordinate, p=excitation
348 # energies in eV
349 self.exmE_rp = np.array(exmE_rp)[:, select] * eu
350 self.expE_rp = np.array(expE_rp)[:, select] * eu
351 # forces in eV / Angstrom
352 self.exF_rp = np.array(exF_rp)[:, select] * eu / 2 / self.delta
353 # matrix elements in e * Angstrom
354 self.exmm_rpc = np.array(exmm_rpc)[:, select, :] * u.Bohr
355 self.expm_rpc = np.array(expm_rpc)[:, select, :] * u.Bohr
356 # matrix element derivatives in e
357 self.exdmdr_rpc = (np.array(exdmdr_rpc)[:, select, :] *
358 u.Bohr / 2 / self.delta)
359 else:
360 # did not read
361 self.exmE_rp = self.expE_rp = self.exF_rp = np.empty(0)
362 self.exmm_rpc = self.expm_rpc = self.exdmdr_rpc = np.empty(0)
364 def read(self, *args, **kwargs):
365 """Read data from a pre-performed calculation."""
366 self.vibrations.read(*args, **kwargs)
367 self.init_parallel_read()
368 if not hasattr(self, 'ex0E_p'):
369 if self.overlap:
370 self.read_excitations_overlap()
371 else:
372 self.read_excitations()
374 self._already_read = True
376 def get_cross_sections(self, omega, gamma):
377 """Returns Raman cross sections for each vibration."""
378 I_v = self.intensity(omega, gamma)
379 pre = 1. / 16 / np.pi**2 / u._eps0**2 / u._c**4
380 # frequency of scattered light
381 omS_v = omega - self.om_Q
382 return pre * omega * omS_v**3 * I_v
384 def get_spectrum(self, omega, gamma=0.1,
385 start=None, end=None, npts=None, width=20,
386 type='Gaussian',
387 intensity_unit='????', normalize=False):
388 """Get resonant Raman spectrum.
390 The method returns wavenumbers in cm^-1 with corresponding
391 Raman cross section.
392 Start and end point, and width of the Gaussian/Lorentzian should
393 be given in cm^-1.
394 """
396 self.type = type.lower()
397 assert self.type in ['gaussian', 'lorentzian']
399 frequencies = self.get_energies().real / u.invcm
400 intensities = self.get_cross_sections(omega, gamma)
401 if width is None:
402 return [frequencies, intensities]
404 if start is None:
405 start = min(self.om_Q) / u.invcm - 3 * width
406 if end is None:
407 end = max(self.om_Q) / u.invcm + 3 * width
409 if not npts:
410 npts = int((end - start) / width * 10 + 1)
412 prefactor = 1
413 if self.type == 'lorentzian':
414 intensities = intensities * width * np.pi / 2.
415 if normalize:
416 prefactor = 2. / width / np.pi
417 else:
418 sigma = width / 2. / np.sqrt(2. * np.log(2.))
419 if normalize:
420 prefactor = 1. / sigma / np.sqrt(2 * np.pi)
421 # Make array with spectrum data
422 spectrum = np.empty(npts)
423 energies = np.linspace(start, end, npts)
424 for i, energy in enumerate(energies):
425 energies[i] = energy
426 if self.type == 'lorentzian':
427 spectrum[i] = (intensities * 0.5 * width / np.pi /
428 ((frequencies - energy)**2 +
429 0.25 * width**2)).sum()
430 else:
431 spectrum[i] = (intensities *
432 np.exp(-(frequencies - energy)**2 /
433 2. / sigma**2)).sum()
434 return [energies, prefactor * spectrum]
436 def write_spectrum(self, omega, gamma,
437 out='resonant-raman-spectra.dat',
438 start=200, end=4000,
439 npts=None, width=10,
440 type='Gaussian'):
441 """Write out spectrum to file.
443 Start and end
444 point, and width of the Gaussian/Lorentzian should be given
445 in cm^-1."""
446 energies, spectrum = self.get_spectrum(omega, gamma,
447 start, end, npts, width,
448 type)
450 # Write out spectrum in file. First column is absolute intensities.
451 outdata = np.empty([len(energies), 3])
452 outdata.T[0] = energies
453 outdata.T[1] = spectrum
455 with paropen(out, 'w') as fd:
456 fd.write('# Resonant Raman spectrum\n')
457 if hasattr(self, '_approx'):
458 fd.write(f'# approximation: {self._approx}\n')
459 for key in self.observation:
460 fd.write(f'# {key}: {self.observation[key]}\n')
461 fd.write('# omega={:g} eV, gamma={:g} eV\n'
462 .format(omega, gamma))
463 if width is not None:
464 fd.write('# %s folded, width=%g cm^-1\n'
465 % (type.title(), width))
466 fd.write('# [cm^-1] [a.u.]\n')
468 for row in outdata:
469 fd.write('%.3f %15.5g\n' %
470 (row[0], row[1]))
472 def summary(self, omega, gamma=0.1,
473 method='standard', direction='central',
474 log=sys.stdout):
475 """Print summary for given omega [eV]"""
476 self.read(method, direction)
477 hnu = self.get_energies()
478 intensities = self.get_absolute_intensities(omega, gamma)
479 te = int(np.log10(intensities.max())) - 2
480 scale = 10**(-te)
481 if not te:
482 ts = ''
483 elif te > -2 and te < 3:
484 ts = str(10**te)
485 else:
486 ts = f'10^{te}'
488 if isinstance(log, str):
489 log = paropen(log, 'a')
491 parprint('-------------------------------------', file=log)
492 parprint(' excitation at ' + str(omega) + ' eV', file=log)
493 parprint(' gamma ' + str(gamma) + ' eV', file=log)
494 parprint(' method:', self.vibrations.method, file=log)
495 parprint(' approximation:', self.approximation, file=log)
496 parprint(' Mode Frequency Intensity', file=log)
497 parprint(f' # meV cm^-1 [{ts}A^4/amu]', file=log)
498 parprint('-------------------------------------', file=log)
499 for n, e in enumerate(hnu):
500 if e.imag != 0:
501 c = 'i'
502 e = e.imag
503 else:
504 c = ' '
505 e = e.real
506 parprint('%3d %6.1f%s %7.1f%s %9.2f' %
507 (n, 1000 * e, c, e / u.invcm, c, intensities[n] * scale),
508 file=log)
509 parprint('-------------------------------------', file=log)
510 parprint('Zero-point energy: %.3f eV' %
511 self.vibrations.get_zero_point_energy(),
512 file=log)
515class LrResonantRaman(ResonantRaman):
516 """Resonant Raman for linear response
518 Quick and dirty approach to enable loading of LrTDDFT calculations
519 """
521 def read_excitations(self):
522 eq_disp = self._eq_disp()
523 ex0_object = eq_disp.read_exobj()
524 eu = ex0_object.energy_to_eV_scale
525 matching = frozenset(ex0_object.kss)
527 def append(lst, disp, matching):
528 exo = disp.read_exobj()
529 lst.append(exo)
530 matching = matching.intersection(exo.kss)
531 return matching
533 exm_object_list = []
534 exp_object_list = []
535 for a in self.indices:
536 for i in 'xyz':
537 disp1 = self._disp(a, i, -1)
538 disp2 = self._disp(a, i, 1)
540 matching = append(exm_object_list,
541 disp1,
542 matching)
543 matching = append(exp_object_list,
544 disp2,
545 matching)
547 def select(exl, matching):
548 exl.diagonalize(**self.exkwargs)
549 mlist = list(exl)
550# mlst = [ex for ex in exl if ex in matching]
551# assert(len(mlst) == len(matching))
552 return mlist
553 ex0 = select(ex0_object, matching)
554 exm = []
555 exp = []
556 r = 0
557 for a in self.indices:
558 for i in 'xyz':
559 exm.append(select(exm_object_list[r], matching))
560 exp.append(select(exp_object_list[r], matching))
561 r += 1
563 self.ex0E_p = np.array([ex.energy * eu for ex in ex0])
564# self.exmE_p = np.array([ex.energy * eu for ex in exm])
565# self.expE_p = np.array([ex.energy * eu for ex in exp])
566 self.ex0m_pc = (np.array(
567 [ex.get_dipole_me(form=self.dipole_form) for ex in ex0]) *
568 u.Bohr)
569 self.exF_rp = []
570 exmE_rp = []
571 expE_rp = []
572 exmm_rpc = []
573 expm_rpc = []
574 r = 0
575 for a in self.indices:
576 for i in 'xyz':
577 exmE_rp.append([em.energy for em in exm[r]])
578 expE_rp.append([ep.energy for ep in exp[r]])
579 self.exF_rp.append(
580 [(em.energy - ep.energy)
581 for ep, em in zip(exp[r], exm[r])])
582 exmm_rpc.append(
583 [ex.get_dipole_me(form=self.dipole_form) for ex in exm[r]])
584 expm_rpc.append(
585 [ex.get_dipole_me(form=self.dipole_form) for ex in exp[r]])
586 r += 1
587 self.exmE_rp = np.array(exmE_rp) * eu
588 self.expE_rp = np.array(expE_rp) * eu
589 self.exF_rp = np.array(self.exF_rp) * eu / 2 / self.delta
590 self.exmm_rpc = np.array(exmm_rpc) * u.Bohr
591 self.expm_rpc = np.array(expm_rpc) * u.Bohr