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