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