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