Coverage for /builds/debichem-team/python-ase/ase/vibrations/albrecht.py : 89.86%
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
1import sys
2import numpy as np
3from itertools import combinations_with_replacement
5import ase.units as u
6from ase.parallel import parprint, paropen
7from ase.vibrations.resonant_raman import ResonantRaman
8from ase.vibrations.franck_condon import FranckCondonOverlap
9from ase.vibrations.franck_condon import FranckCondonRecursive
12class Albrecht(ResonantRaman):
13 def __init__(self, *args, **kwargs):
14 """
15 Parameters
16 ----------
17 all from ResonantRaman.__init__
18 combinations: int
19 Combinations to consider for multiple excitations.
20 Default is 1, possible 2
21 skip: int
22 Number of first transitions to exclude. Default 0,
23 recommended: 5 for linear molecules, 6 for other molecules
24 nm: int
25 Number of intermediate m levels to consider, default 20
26 """
27 self.combinations = kwargs.pop('combinations', 1)
28 self.skip = kwargs.pop('skip', 0)
29 self.nm = kwargs.pop('nm', 20)
30 approximation = kwargs.pop('approximation', 'Albrecht')
32 ResonantRaman.__init__(self, *args, **kwargs)
34 self.set_approximation(approximation)
36 def set_approximation(self, value):
37 approx = value.lower()
38 if approx in ['albrecht', 'albrecht b', 'albrecht c', 'albrecht bc']:
39 if not self.overlap:
40 raise ValueError('Overlaps are needed')
41 elif not approx == 'albrecht a':
42 raise ValueError('Please use "Albrecht" or "Albrecht A/B/C/BC"')
43 self._approx = value
45 def calculate_energies_and_modes(self):
46 if hasattr(self, 'im_r'):
47 return
49 ResonantRaman.calculate_energies_and_modes(self)
51 # single transitions and their occupation
52 om_Q = self.om_Q[self.skip:]
53 om_v = om_Q
54 ndof = len(om_Q)
55 n_vQ = np.eye(ndof, dtype=int)
57 l_Q = range(ndof)
58 ind_v = list(combinations_with_replacement(l_Q, 1))
60 if self.combinations > 1:
61 if not self.combinations == 2:
62 raise NotImplementedError
64 for c in range(2, self.combinations + 1):
65 ind_v += list(combinations_with_replacement(l_Q, c))
67 nv = len(ind_v)
68 n_vQ = np.zeros((nv, ndof), dtype=int)
69 om_v = np.zeros((nv), dtype=float)
70 for j, wt in enumerate(ind_v):
71 for i in wt:
72 n_vQ[j, i] += 1
73 om_v = n_vQ.dot(om_Q)
75 self.ind_v = ind_v
76 self.om_v = om_v
77 self.n_vQ = n_vQ # how many of each
78 self.d_vQ = np.where(n_vQ > 0, 1, 0) # do we have them ?
80 def get_energies(self):
81 self.calculate_energies_and_modes()
82 return self.om_v
84 def _collect_r(self, arr_ro, oshape, dtype):
85 """Collect an array that is distributed."""
86 if len(self.myr) == self.ndof: # serial
87 return arr_ro
88 data_ro = np.zeros([self.ndof] + oshape, dtype)
89 if len(arr_ro):
90 data_ro[self.slize] = arr_ro
91 self.comm.sum(data_ro)
92 return data_ro
94 def Huang_Rhys_factors(self, forces_r):
95 """Evaluate Huang-Rhys factors derived from forces."""
96 return 0.5 * self.unitless_displacements(forces_r)**2
98 def unitless_displacements(self, forces_r, mineigv=1e-12):
99 """Evaluate unitless displacements from forces
101 Parameters
102 ----------
103 forces_r: array
104 Forces in cartesian coordinates
105 mineigv: float
106 Minimal Eigenvalue to consider in matrix inversion to handle
107 numerical noise. Default 1e-12
109 Returns
110 -------
111 Unitless displacements in Eigenmode coordinates
112 """
113 assert(len(forces_r.flat) == self.ndof)
115 if not hasattr(self, 'Dm1_q'):
116 self.eigv_q, self.eigw_rq = np.linalg.eigh(
117 self.im_r[:, None] * self.H * self.im_r)
118 # there might be zero or nearly zero eigenvalues
119 self.Dm1_q = np.divide(1, self.eigv_q,
120 out=np.zeros_like(self.eigv_q),
121 where=np.abs(self.eigv_q) > mineigv)
122 X_r = self.eigw_rq @ np.diag(self.Dm1_q) @ self.eigw_rq.T @ (
123 forces_r.flat * self.im_r)
125 d_Q = np.dot(self.modes_Qq, X_r)
126 s = 1.e-20 / u.kg / u.C / u._hbar**2
127 d_Q *= np.sqrt(s * self.om_Q)
129 return d_Q
131 def omegaLS(self, omega, gamma):
132 omL = omega + 1j * gamma
133 omS_Q = omL - self.om_Q
134 return omL, omS_Q
136 def init_parallel_excitations(self):
137 """Init for paralellization over excitations."""
138 n_p = len(self.ex0E_p)
140 # collect excited state forces
141 exF_pr = self._collect_r(self.exF_rp, [n_p], self.ex0E_p.dtype).T
143 # select your work load
144 myn = -(-n_p // self.comm.size) # ceil divide
145 rank = self.comm.rank
146 s = slice(myn * rank, myn * (rank + 1))
147 return n_p, range(n_p)[s], exF_pr
149 def meA(self, omega, gamma=0.1):
150 """Evaluate Albrecht A term.
152 Returns
153 -------
154 Full Albrecht A matrix element. Unit: e^2 Angstrom^2 / eV
155 """
156 self.read()
158 if not hasattr(self, 'fcr'):
159 self.fcr = FranckCondonRecursive()
161 omL = omega + 1j * gamma
162 omS_Q = omL - self.om_Q
164 n_p, myp, exF_pr = self.init_parallel_excitations()
165 exF_pr = np.where(np.abs(exF_pr) > 1e-2, exF_pr, 0)
167 m_Qcc = np.zeros((self.ndof, 3, 3), dtype=complex)
168 for p in myp:
169 energy = self.ex0E_p[p]
170 d_Q = self.unitless_displacements(exF_pr[p])
171 energy_Q = energy - self.om_Q * d_Q**2 / 2.
172 me_cc = np.outer(self.ex0m_pc[p], self.ex0m_pc[p].conj())
174 wm_Q = np.zeros((self.ndof), dtype=complex)
175 wp_Q = np.zeros((self.ndof), dtype=complex)
176 for m in range(self.nm):
177 fco_Q = self.fcr.direct0mm1(m, d_Q)
178 e_Q = energy_Q + m * self.om_Q
179 wm_Q += fco_Q / (e_Q - omL)
180 wp_Q += fco_Q / (e_Q + omS_Q)
181 m_Qcc += np.einsum('a,bc->abc', wm_Q, me_cc)
182 m_Qcc += np.einsum('a,bc->abc', wp_Q, me_cc.conj())
183 self.comm.sum(m_Qcc)
185 return m_Qcc # e^2 Angstrom^2 / eV
187 def meAmult(self, omega, gamma=0.1):
188 """Evaluate Albrecht A term.
190 Returns
191 -------
192 Full Albrecht A matrix element. Unit: e^2 Angstrom^2 / eV
193 """
194 self.read()
196 if not hasattr(self, 'fcr'):
197 self.fcr = FranckCondonRecursive()
199 omL = omega + 1j * gamma
200 omS_v = omL - self.om_v
201 nv = len(self.om_v)
202 om_Q = self.om_Q[self.skip:]
203 nQ = len(om_Q)
205 # n_v:
206 # how many FC factors are involved
207 # nvib_ov:
208 # delta functions to switch contributions depending on order o
209 # ind_ov:
210 # Q indicees
211 # n_ov:
212 # # of vibrational excitations
213 n_v = self.d_vQ.sum(axis=1) # multiplicity
215 nvib_ov = np.empty((self.combinations, nv), dtype=int)
216 om_ov = np.zeros((self.combinations, nv), dtype=float)
217 n_ov = np.zeros((self.combinations, nv), dtype=int)
218 d_ovQ = np.zeros((self.combinations, nv, nQ), dtype=int)
219 for o in range(self.combinations):
220 nvib_ov[o] = np.array(n_v == (o + 1))
221 for v in range(nv):
222 try:
223 om_ov[o, v] = om_Q[self.ind_v[v][o]]
224 d_ovQ[o, v, self.ind_v[v][o]] = 1
225 except IndexError:
226 pass
227 # XXXX change ????
228 n_ov[0] = self.n_vQ.max(axis=1)
229 n_ov[1] = nvib_ov[1]
231 n_p, myp, exF_pr = self.init_parallel_excitations()
233 m_vcc = np.zeros((nv, 3, 3), dtype=complex)
234 for p in myp:
235 energy = self.ex0E_p[p]
236 d_Q = self.unitless_displacements(exF_pr[p])[self.skip:]
237 S_Q = d_Q**2 / 2.
238 energy_v = energy - self.d_vQ.dot(om_Q * S_Q)
239 me_cc = np.outer(self.ex0m_pc[p], self.ex0m_pc[p].conj())
241 fco1_mQ = np.empty((self.nm, nQ), dtype=float)
242 fco2_mQ = np.empty((self.nm, nQ), dtype=float)
243 for m in range(self.nm):
244 fco1_mQ[m] = self.fcr.direct0mm1(m, d_Q)
245 fco2_mQ[m] = self.fcr.direct0mm2(m, d_Q)
247 wm_v = np.zeros((nv), dtype=complex)
248 wp_v = np.zeros((nv), dtype=complex)
249 for m in range(self.nm):
250 fco1_v = np.where(n_ov[0] == 2,
251 d_ovQ[0].dot(fco2_mQ[m]),
252 d_ovQ[0].dot(fco1_mQ[m]))
254 em_v = energy_v + m * om_ov[0]
255 # multiples of same kind
256 fco_v = nvib_ov[0] * fco1_v
257 wm_v += fco_v / (em_v - omL)
258 wp_v += fco_v / (em_v + omS_v)
259 if nvib_ov[1].any():
260 # multiples of mixed type
261 for n in range(self.nm):
262 fco2_v = d_ovQ[1].dot(fco1_mQ[n])
263 e_v = em_v + n * om_ov[1]
264 fco_v = nvib_ov[1] * fco1_v * fco2_v
265 wm_v += fco_v / (e_v - omL)
266 wp_v += fco_v / (e_v + omS_v)
268 m_vcc += np.einsum('a,bc->abc', wm_v, me_cc)
269 m_vcc += np.einsum('a,bc->abc', wp_v, me_cc.conj())
270 self.comm.sum(m_vcc)
272 return m_vcc # e^2 Angstrom^2 / eV
274 def meBC(self, omega, gamma=0.1,
275 term='BC'):
276 """Evaluate Albrecht BC term.
278 Returns
279 -------
280 Full Albrecht BC matrix element.
281 Unit: e^2 Angstrom / eV / sqrt(amu)
282 """
283 self.read()
285 if not hasattr(self, 'fco'):
286 self.fco = FranckCondonOverlap()
288 omL = omega + 1j * gamma
289 omS_Q = omL - self.om_Q
291 # excited state forces
292 n_p, myp, exF_pr = self.init_parallel_excitations()
293 # derivatives after normal coordinates
294 exdmdr_rpc = self._collect_r(
295 self.exdmdr_rpc, [n_p, 3], self.ex0m_pc.dtype)
296 dmdq_qpc = (exdmdr_rpc.T * self.im_r).T # unit e / sqrt(amu)
297 dmdQ_Qpc = np.dot(dmdq_qpc.T, self.modes_Qq.T).T # unit e / sqrt(amu)
299 me_Qcc = np.zeros((self.ndof, 3, 3), dtype=complex)
300 for p in myp:
301 energy = self.ex0E_p[p]
302 S_Q = self.Huang_Rhys_factors(exF_pr[p])
303 # relaxed excited state energy
304 # # n_vQ = np.where(self.n_vQ > 0, 1, 0)
305 # # energy_v = energy - n_vQ.dot(self.om_Q * S_Q)
306 energy_Q = energy - self.om_Q * S_Q
308 # # me_cc = np.outer(self.ex0m_pc[p], self.ex0m_pc[p].conj())
309 m_c = self.ex0m_pc[p] # e Angstrom
310 dmdQ_Qc = dmdQ_Qpc[:, p] # e / sqrt(amu)
312 wBLS_Q = np.zeros((self.ndof), dtype=complex)
313 wBSL_Q = np.zeros((self.ndof), dtype=complex)
314 wCLS_Q = np.zeros((self.ndof), dtype=complex)
315 wCSL_Q = np.zeros((self.ndof), dtype=complex)
316 for m in range(self.nm):
317 f0mmQ1_Q = (self.fco.directT0(m, S_Q) +
318 np.sqrt(2) * self.fco.direct0mm2(m, S_Q))
319 f0Qmm1_Q = self.fco.direct(1, m, S_Q)
321 em_Q = energy_Q + m * self.om_Q
322 wBLS_Q += f0mmQ1_Q / (em_Q - omL)
323 wBSL_Q += f0Qmm1_Q / (em_Q - omL)
324 wCLS_Q += f0mmQ1_Q / (em_Q + omS_Q)
325 wCSL_Q += f0Qmm1_Q / (em_Q + omS_Q)
327 # unit e^2 Angstrom / sqrt(amu)
328 mdmdQ_Qcc = np.einsum('a,bc->bac', m_c, dmdQ_Qc.conj())
329 dmdQm_Qcc = np.einsum('ab,c->abc', dmdQ_Qc, m_c.conj())
330 if 'B' in term:
331 me_Qcc += np.multiply(wBLS_Q, mdmdQ_Qcc.T).T
332 me_Qcc += np.multiply(wBSL_Q, dmdQm_Qcc.T).T
333 if 'C' in term:
334 me_Qcc += np.multiply(wCLS_Q, mdmdQ_Qcc.T).T
335 me_Qcc += np.multiply(wCSL_Q, dmdQm_Qcc.T).T
337 self.comm.sum(me_Qcc)
338 return me_Qcc # unit e^2 Angstrom / eV / sqrt(amu)
340 def electronic_me_Qcc(self, omega, gamma):
341 self.calculate_energies_and_modes()
343 approx = self.approximation.lower()
344 assert(self.combinations == 1)
345 Vel_Qcc = np.zeros((len(self.om_Q), 3, 3), dtype=complex)
346 if approx == 'albrecht a' or approx == 'albrecht':
347 Vel_Qcc += self.meA(omega, gamma) # e^2 Angstrom^2 / eV
348 # divide through pre-factor
349 with np.errstate(divide='ignore'):
350 Vel_Qcc *= np.where(self.vib01_Q > 0,
351 1. / self.vib01_Q, 0)[:, None, None]
352 # -> e^2 Angstrom / eV / sqrt(amu)
353 if approx == 'albrecht bc' or approx == 'albrecht':
354 Vel_Qcc += self.meBC(omega, gamma) # e^2 Angstrom / eV / sqrt(amu)
355 if approx == 'albrecht b':
356 Vel_Qcc += self.meBC(omega, gamma, term='B')
357 if approx == 'albrecht c':
358 Vel_Qcc = self.meBC(omega, gamma, term='C')
360 Vel_Qcc *= u.Hartree * u.Bohr # e^2 Angstrom^2 / eV -> Angstrom^3
362 return Vel_Qcc # Angstrom^2 / sqrt(amu)
364 def me_Qcc(self, omega, gamma):
365 """Full matrix element"""
366 self.read()
367 approx = self.approximation.lower()
368 nv = len(self.om_v)
369 V_vcc = np.zeros((nv, 3, 3), dtype=complex)
370 if approx == 'albrecht a' or approx == 'albrecht':
371 if self.combinations == 1:
372 # e^2 Angstrom^2 / eV
373 V_vcc += self.meA(omega, gamma)[self.skip:]
374 else:
375 V_vcc += self.meAmult(omega, gamma)
376 if approx == 'albrecht bc' or approx == 'albrecht':
377 if self.combinations == 1:
378 vel_vcc = self.meBC(omega, gamma)
379 V_vcc += vel_vcc * self.vib01_Q[:, None, None]
380 else:
381 vel_vcc = self.meBCmult(omega, gamma)
382 V_vcc = 0
383 elif approx == 'albrecht b':
384 assert(self.combinations == 1)
385 vel_vcc = self.meBC(omega, gamma, term='B')
386 V_vcc = vel_vcc * self.vib01_Q[:, None, None]
387 if approx == 'albrecht c':
388 assert(self.combinations == 1)
389 vel_vcc = self.meBC(omega, gamma, term='C')
390 V_vcc = vel_vcc * self.vib01_Q[:, None, None]
392 return V_vcc # e^2 Angstrom^2 / eV
394 def summary(self, omega=0, gamma=0,
395 method='standard', direction='central',
396 log=sys.stdout):
397 """Print summary for given omega [eV]"""
398 if self.combinations > 1:
399 return self.extended_summary()
401 om_v = self.get_energies()
402 intensities = self.get_absolute_intensities(omega, gamma)[self.skip:]
404 if isinstance(log, str):
405 log = paropen(log, 'a')
407 parprint('-------------------------------------', file=log)
408 parprint(' excitation at ' + str(omega) + ' eV', file=log)
409 parprint(' gamma ' + str(gamma) + ' eV', file=log)
410 parprint(' approximation:', self.approximation, file=log)
411 parprint(' Mode Frequency Intensity', file=log)
412 parprint(' # meV cm^-1 [A^4/amu]', file=log)
413 parprint('-------------------------------------', file=log)
414 for n, e in enumerate(om_v):
415 if e.imag != 0:
416 c = 'i'
417 e = e.imag
418 else:
419 c = ' '
420 e = e.real
421 parprint('%3d %6.1f %7.1f%s %9.1f' %
422 (n, 1000 * e, e / u.invcm, c, intensities[n]),
423 file=log)
424 parprint('-------------------------------------', file=log)
425 parprint('Zero-point energy: %.3f eV' %
426 self.vibrations.get_zero_point_energy(),
427 file=log)
429 def extended_summary(self, omega=0, gamma=0,
430 method='standard', direction='central',
431 log=sys.stdout):
432 """Print summary for given omega [eV]"""
433 self.read(method, direction)
434 om_v = self.get_energies()
435 intens_v = self.intensity(omega, gamma)
437 if isinstance(log, str):
438 log = paropen(log, 'a')
440 parprint('-------------------------------------', file=log)
441 parprint(' excitation at ' + str(omega) + ' eV', file=log)
442 parprint(' gamma ' + str(gamma) + ' eV', file=log)
443 parprint(' approximation:', self.approximation, file=log)
444 parprint(' observation:', self.observation, file=log)
445 parprint(' Mode Frequency Intensity', file=log)
446 parprint(' # meV cm^-1 [e^4A^4/eV^2]', file=log)
447 parprint('-------------------------------------', file=log)
448 for v, e in enumerate(om_v):
449 parprint(self.ind_v[v], '{0:6.1f} {1:7.1f} {2:9.1f}'.format(
450 1000 * e, e / u.invcm, 1e9 * intens_v[v]),
451 file=log)
452 parprint('-------------------------------------', file=log)
453 parprint('Zero-point energy: %.3f eV' %
454 self.vibrations.get_zero_point_energy(),
455 file=log)