Coverage for /builds/kinetik161/ase/ase/vibrations/albrecht.py: 89.82%

1import sys

2from itertools import combinations_with_replacement

4import numpy as np

6import ase.units as u

7from ase.parallel import paropen, parprint

8from ase.vibrations.franck_condon import (FranckCondonOverlap,

9 FranckCondonRecursive)

10from ase.vibrations.resonant_raman import ResonantRaman

13class Albrecht(ResonantRaman):

14 def __init__(self, *args, **kwargs):

15 """

16 Parameters

17 ----------

18 all from ResonantRaman.__init__

19 combinations: int

20 Combinations to consider for multiple excitations.

21 Default is 1, possible 2

22 skip: int

23 Number of first transitions to exclude. Default 0,

24 recommended: 5 for linear molecules, 6 for other molecules

25 nm: int

26 Number of intermediate m levels to consider, default 20

27 """

28 self.combinations = kwargs.pop('combinations', 1)

29 self.skip = kwargs.pop('skip', 0)

30 self.nm = kwargs.pop('nm', 20)

31 approximation = kwargs.pop('approximation', 'Albrecht')

33 ResonantRaman.__init__(self, *args, **kwargs)

35 self.set_approximation(approximation)

37 def set_approximation(self, value):

38 approx = value.lower()

39 if approx in ['albrecht', 'albrecht b', 'albrecht c', 'albrecht bc']:

40 if not self.overlap:

41 raise ValueError('Overlaps are needed')

42 elif not approx == 'albrecht a':

43 raise ValueError('Please use "Albrecht" or "Albrecht A/B/C/BC"')

44 self._approx = value

46 def calculate_energies_and_modes(self):

47 if hasattr(self, 'im_r'):

48 return

50 ResonantRaman.calculate_energies_and_modes(self)

52 # single transitions and their occupation

53 om_Q = self.om_Q[self.skip:]

54 om_v = om_Q

55 ndof = len(om_Q)

56 n_vQ = np.eye(ndof, dtype=int)

58 l_Q = range(ndof)

59 ind_v = list(combinations_with_replacement(l_Q, 1))

61 if self.combinations > 1:

62 if not self.combinations == 2:

63 raise NotImplementedError

65 for c in range(2, self.combinations + 1):

66 ind_v += list(combinations_with_replacement(l_Q, c))

68 nv = len(ind_v)

69 n_vQ = np.zeros((nv, ndof), dtype=int)

70 om_v = np.zeros((nv), dtype=float)

71 for j, wt in enumerate(ind_v):

72 for i in wt:

73 n_vQ[j, i] += 1

74 om_v = n_vQ.dot(om_Q)

76 self.ind_v = ind_v

77 self.om_v = om_v

78 self.n_vQ = n_vQ # how many of each

79 self.d_vQ = np.where(n_vQ > 0, 1, 0) # do we have them ?

81 def get_energies(self):

82 self.calculate_energies_and_modes()

83 return self.om_v

85 def _collect_r(self, arr_ro, oshape, dtype):

86 """Collect an array that is distributed."""

87 if len(self.myr) == self.ndof: # serial

88 return arr_ro

89 data_ro = np.zeros([self.ndof] + oshape, dtype)

90 if len(arr_ro):

91 data_ro[self.slize] = arr_ro

92 self.comm.sum(data_ro)

93 return data_ro

95 def Huang_Rhys_factors(self, forces_r):

96 """Evaluate Huang-Rhys factors derived from forces."""

97 return 0.5 * self.unitless_displacements(forces_r)**2

99 def unitless_displacements(self, forces_r, mineigv=1e-12):

100 """Evaluate unitless displacements from forces

101

102 Parameters

103 ----------

104 forces_r: array

105 Forces in cartesian coordinates

106 mineigv: float

107 Minimal Eigenvalue to consider in matrix inversion to handle

108 numerical noise. Default 1e-12

109

110 Returns

111 -------

112 Unitless displacements in Eigenmode coordinates

113 """

114 assert len(forces_r.flat) == self.ndof

115

116 if not hasattr(self, 'Dm1_q'):

117 self.eigv_q, self.eigw_rq = np.linalg.eigh(

118 self.im_r[:, None] * self.H * self.im_r)

119 # there might be zero or nearly zero eigenvalues

120 self.Dm1_q = np.divide(1, self.eigv_q,

121 out=np.zeros_like(self.eigv_q),

122 where=np.abs(self.eigv_q) > mineigv)

123 X_r = self.eigw_rq @ np.diag(self.Dm1_q) @ self.eigw_rq.T @ (

124 forces_r.flat * self.im_r)

125

126 d_Q = np.dot(self.modes_Qq, X_r)

127 s = 1.e-20 / u.kg / u.C / u._hbar**2

128 d_Q *= np.sqrt(s * self.om_Q)

129

130 return d_Q

131

132 def omegaLS(self, omega, gamma):

133 omL = omega + 1j * gamma

134 omS_Q = omL - self.om_Q

135 return omL, omS_Q

136

137 def init_parallel_excitations(self):

138 """Init for paralellization over excitations."""

139 n_p = len(self.ex0E_p)

140

141 # collect excited state forces

142 exF_pr = self._collect_r(self.exF_rp, [n_p], self.ex0E_p.dtype).T

143

144 # select your work load

145 myn = -(-n_p // self.comm.size) # ceil divide

146 rank = self.comm.rank

147 s = slice(myn * rank, myn * (rank + 1))

148 return n_p, range(n_p)[s], exF_pr

149

150 def meA(self, omega, gamma=0.1):

151 """Evaluate Albrecht A term.

152

153 Returns

154 -------

155 Full Albrecht A matrix element. Unit: e^2 Angstrom^2 / eV

156 """

157 self.read()

158

159 if not hasattr(self, 'fcr'):

160 self.fcr = FranckCondonRecursive()

161

162 omL = omega + 1j * gamma

163 omS_Q = omL - self.om_Q

164

165 n_p, myp, exF_pr = self.init_parallel_excitations()

166 exF_pr = np.where(np.abs(exF_pr) > 1e-2, exF_pr, 0)

167

168 m_Qcc = np.zeros((self.ndof, 3, 3), dtype=complex)

169 for p in myp:

170 energy = self.ex0E_p[p]

171 d_Q = self.unitless_displacements(exF_pr[p])

172 energy_Q = energy - self.om_Q * d_Q**2 / 2.

173 me_cc = np.outer(self.ex0m_pc[p], self.ex0m_pc[p].conj())

174

175 wm_Q = np.zeros((self.ndof), dtype=complex)

176 wp_Q = np.zeros((self.ndof), dtype=complex)

177 for m in range(self.nm):

178 fco_Q = self.fcr.direct0mm1(m, d_Q)

179 e_Q = energy_Q + m * self.om_Q

180 wm_Q += fco_Q / (e_Q - omL)

181 wp_Q += fco_Q / (e_Q + omS_Q)

182 m_Qcc += np.einsum('a,bc->abc', wm_Q, me_cc)

183 m_Qcc += np.einsum('a,bc->abc', wp_Q, me_cc.conj())

184 self.comm.sum(m_Qcc)

185

186 return m_Qcc # e^2 Angstrom^2 / eV

187

188 def meAmult(self, omega, gamma=0.1):

189 """Evaluate Albrecht A term.

190

191 Returns

192 -------

193 Full Albrecht A matrix element. Unit: e^2 Angstrom^2 / eV

194 """

195 self.read()

196

197 if not hasattr(self, 'fcr'):

198 self.fcr = FranckCondonRecursive()

199

200 omL = omega + 1j * gamma

201 omS_v = omL - self.om_v

202 nv = len(self.om_v)

203 om_Q = self.om_Q[self.skip:]

204 nQ = len(om_Q)

205

206 # n_v:

207 # how many FC factors are involved

208 # nvib_ov:

209 # delta functions to switch contributions depending on order o

210 # ind_ov:

211 # Q indicees

212 # n_ov:

213 # # of vibrational excitations

214 n_v = self.d_vQ.sum(axis=1) # multiplicity

215

216 nvib_ov = np.empty((self.combinations, nv), dtype=int)

217 om_ov = np.zeros((self.combinations, nv), dtype=float)

218 n_ov = np.zeros((self.combinations, nv), dtype=int)

219 d_ovQ = np.zeros((self.combinations, nv, nQ), dtype=int)

220 for o in range(self.combinations):

221 nvib_ov[o] = np.array(n_v == (o + 1))

222 for v in range(nv):

223 try:

224 om_ov[o, v] = om_Q[self.ind_v[v][o]]

225 d_ovQ[o, v, self.ind_v[v][o]] = 1

226 except IndexError:

227 pass

228 # XXXX change ????

229 n_ov[0] = self.n_vQ.max(axis=1)

230 n_ov[1] = nvib_ov[1]

231

232 n_p, myp, exF_pr = self.init_parallel_excitations()

233

234 m_vcc = np.zeros((nv, 3, 3), dtype=complex)

235 for p in myp:

236 energy = self.ex0E_p[p]

237 d_Q = self.unitless_displacements(exF_pr[p])[self.skip:]

238 S_Q = d_Q**2 / 2.

239 energy_v = energy - self.d_vQ.dot(om_Q * S_Q)

240 me_cc = np.outer(self.ex0m_pc[p], self.ex0m_pc[p].conj())

241

242 fco1_mQ = np.empty((self.nm, nQ), dtype=float)

243 fco2_mQ = np.empty((self.nm, nQ), dtype=float)

244 for m in range(self.nm):

245 fco1_mQ[m] = self.fcr.direct0mm1(m, d_Q)

246 fco2_mQ[m] = self.fcr.direct0mm2(m, d_Q)

247

248 wm_v = np.zeros((nv), dtype=complex)

249 wp_v = np.zeros((nv), dtype=complex)

250 for m in range(self.nm):

251 fco1_v = np.where(n_ov[0] == 2,

252 d_ovQ[0].dot(fco2_mQ[m]),

253 d_ovQ[0].dot(fco1_mQ[m]))

254

255 em_v = energy_v + m * om_ov[0]

256 # multiples of same kind

257 fco_v = nvib_ov[0] * fco1_v

258 wm_v += fco_v / (em_v - omL)

259 wp_v += fco_v / (em_v + omS_v)

260 if nvib_ov[1].any():

261 # multiples of mixed type

262 for n in range(self.nm):

263 fco2_v = d_ovQ[1].dot(fco1_mQ[n])

264 e_v = em_v + n * om_ov[1]

265 fco_v = nvib_ov[1] * fco1_v * fco2_v

266 wm_v += fco_v / (e_v - omL)

267 wp_v += fco_v / (e_v + omS_v)

268

269 m_vcc += np.einsum('a,bc->abc', wm_v, me_cc)

270 m_vcc += np.einsum('a,bc->abc', wp_v, me_cc.conj())

271 self.comm.sum(m_vcc)

272

273 return m_vcc # e^2 Angstrom^2 / eV

274

275 def meBC(self, omega, gamma=0.1,

276 term='BC'):

277 """Evaluate Albrecht BC term.

278

279 Returns

280 -------

281 Full Albrecht BC matrix element.

282 Unit: e^2 Angstrom / eV / sqrt(amu)

283 """

284 self.read()

285

286 if not hasattr(self, 'fco'):

287 self.fco = FranckCondonOverlap()

288

289 omL = omega + 1j * gamma

290 omS_Q = omL - self.om_Q

291

292 # excited state forces

293 n_p, myp, exF_pr = self.init_parallel_excitations()

294 # derivatives after normal coordinates

295 exdmdr_rpc = self._collect_r(

296 self.exdmdr_rpc, [n_p, 3], self.ex0m_pc.dtype)

297 dmdq_qpc = (exdmdr_rpc.T * self.im_r).T # unit e / sqrt(amu)

298 dmdQ_Qpc = np.dot(dmdq_qpc.T, self.modes_Qq.T).T # unit e / sqrt(amu)

299

300 me_Qcc = np.zeros((self.ndof, 3, 3), dtype=complex)

301 for p in myp:

302 energy = self.ex0E_p[p]

303 S_Q = self.Huang_Rhys_factors(exF_pr[p])

304 # relaxed excited state energy

305 # # n_vQ = np.where(self.n_vQ > 0, 1, 0)

306 # # energy_v = energy - n_vQ.dot(self.om_Q * S_Q)

307 energy_Q = energy - self.om_Q * S_Q

308

309 # # me_cc = np.outer(self.ex0m_pc[p], self.ex0m_pc[p].conj())

310 m_c = self.ex0m_pc[p] # e Angstrom

311 dmdQ_Qc = dmdQ_Qpc[:, p] # e / sqrt(amu)

312

313 wBLS_Q = np.zeros((self.ndof), dtype=complex)

314 wBSL_Q = np.zeros((self.ndof), dtype=complex)

315 wCLS_Q = np.zeros((self.ndof), dtype=complex)

316 wCSL_Q = np.zeros((self.ndof), dtype=complex)

317 for m in range(self.nm):

318 f0mmQ1_Q = (self.fco.directT0(m, S_Q) +

319 np.sqrt(2) * self.fco.direct0mm2(m, S_Q))

320 f0Qmm1_Q = self.fco.direct(1, m, S_Q)

321

322 em_Q = energy_Q + m * self.om_Q

323 wBLS_Q += f0mmQ1_Q / (em_Q - omL)

324 wBSL_Q += f0Qmm1_Q / (em_Q - omL)

325 wCLS_Q += f0mmQ1_Q / (em_Q + omS_Q)

326 wCSL_Q += f0Qmm1_Q / (em_Q + omS_Q)

327

328 # unit e^2 Angstrom / sqrt(amu)

329 mdmdQ_Qcc = np.einsum('a,bc->bac', m_c, dmdQ_Qc.conj())

330 dmdQm_Qcc = np.einsum('ab,c->abc', dmdQ_Qc, m_c.conj())

331 if 'B' in term:

332 me_Qcc += np.multiply(wBLS_Q, mdmdQ_Qcc.T).T

333 me_Qcc += np.multiply(wBSL_Q, dmdQm_Qcc.T).T

334 if 'C' in term:

335 me_Qcc += np.multiply(wCLS_Q, mdmdQ_Qcc.T).T

336 me_Qcc += np.multiply(wCSL_Q, dmdQm_Qcc.T).T

337

338 self.comm.sum(me_Qcc)

339 return me_Qcc # unit e^2 Angstrom / eV / sqrt(amu)

340

341 def electronic_me_Qcc(self, omega, gamma):

342 self.calculate_energies_and_modes()

343

344 approx = self.approximation.lower()

345 assert self.combinations == 1

346 Vel_Qcc = np.zeros((len(self.om_Q), 3, 3), dtype=complex)

347 if approx == 'albrecht a' or approx == 'albrecht':

348 Vel_Qcc += self.meA(omega, gamma) # e^2 Angstrom^2 / eV

349 # divide through pre-factor

350 with np.errstate(divide='ignore'):

351 Vel_Qcc *= np.where(self.vib01_Q > 0,

352 1. / self.vib01_Q, 0)[:, None, None]

353 # -> e^2 Angstrom / eV / sqrt(amu)

354 if approx == 'albrecht bc' or approx == 'albrecht':

355 Vel_Qcc += self.meBC(omega, gamma) # e^2 Angstrom / eV / sqrt(amu)

356 if approx == 'albrecht b':

357 Vel_Qcc += self.meBC(omega, gamma, term='B')

358 if approx == 'albrecht c':

359 Vel_Qcc = self.meBC(omega, gamma, term='C')

360

361 Vel_Qcc *= u.Hartree * u.Bohr # e^2 Angstrom^2 / eV -> Angstrom^3

362

363 return Vel_Qcc # Angstrom^2 / sqrt(amu)

364

365 def me_Qcc(self, omega, gamma):

366 """Full matrix element"""

367 self.read()

368 approx = self.approximation.lower()

369 nv = len(self.om_v)

370 V_vcc = np.zeros((nv, 3, 3), dtype=complex)

371 if approx == 'albrecht a' or approx == 'albrecht':

372 if self.combinations == 1:

373 # e^2 Angstrom^2 / eV

374 V_vcc += self.meA(omega, gamma)[self.skip:]

375 else:

376 V_vcc += self.meAmult(omega, gamma)

377 if approx == 'albrecht bc' or approx == 'albrecht':

378 if self.combinations == 1:

379 vel_vcc = self.meBC(omega, gamma)

380 V_vcc += vel_vcc * self.vib01_Q[:, None, None]

381 else:

382 vel_vcc = self.meBCmult(omega, gamma)

383 V_vcc = 0

384 elif approx == 'albrecht b':

385 assert self.combinations == 1

386 vel_vcc = self.meBC(omega, gamma, term='B')

387 V_vcc = vel_vcc * self.vib01_Q[:, None, None]

388 if approx == 'albrecht c':

389 assert self.combinations == 1

390 vel_vcc = self.meBC(omega, gamma, term='C')

391 V_vcc = vel_vcc * self.vib01_Q[:, None, None]

392

393 return V_vcc # e^2 Angstrom^2 / eV

394

395 def summary(self, omega=0, gamma=0,

396 method='standard', direction='central',

397 log=sys.stdout):

398 """Print summary for given omega [eV]"""

399 if self.combinations > 1:

400 return self.extended_summary()

401

402 om_v = self.get_energies()

403 intensities = self.get_absolute_intensities(omega, gamma)[self.skip:]

404

405 if isinstance(log, str):

406 log = paropen(log, 'a')

407

408 parprint('-------------------------------------', file=log)

409 parprint(' excitation at ' + str(omega) + ' eV', file=log)

410 parprint(' gamma ' + str(gamma) + ' eV', file=log)

411 parprint(' approximation:', self.approximation, file=log)

412 parprint(' Mode Frequency Intensity', file=log)

413 parprint(' # meV cm^-1 [A^4/amu]', file=log)

414 parprint('-------------------------------------', file=log)

415 for n, e in enumerate(om_v):

416 if e.imag != 0:

417 c = 'i'

418 e = e.imag

419 else:

420 c = ' '

421 e = e.real

422 parprint('%3d %6.1f %7.1f%s %9.1f' %

423 (n, 1000 * e, e / u.invcm, c, intensities[n]),

424 file=log)

425 parprint('-------------------------------------', file=log)

426 parprint('Zero-point energy: %.3f eV' %

427 self.vibrations.get_zero_point_energy(),

428 file=log)

429

430 def extended_summary(self, omega=0, gamma=0,

431 method='standard', direction='central',

432 log=sys.stdout):

433 """Print summary for given omega [eV]"""

434 self.read(method, direction)

435 om_v = self.get_energies()

436 intens_v = self.intensity(omega, gamma)

437

438 if isinstance(log, str):

439 log = paropen(log, 'a')

440

441 parprint('-------------------------------------', file=log)

442 parprint(' excitation at ' + str(omega) + ' eV', file=log)

443 parprint(' gamma ' + str(gamma) + ' eV', file=log)

444 parprint(' approximation:', self.approximation, file=log)

445 parprint(' observation:', self.observation, file=log)

446 parprint(' Mode Frequency Intensity', file=log)

447 parprint(' # meV cm^-1 [e^4A^4/eV^2]', file=log)

448 parprint('-------------------------------------', file=log)

449 for v, e in enumerate(om_v):

450 parprint(self.ind_v[v], '{:6.1f} {:7.1f} {:9.1f}'.format(

451 1000 * e, e / u.invcm, 1e9 * intens_v[v]),

452 file=log)

453 parprint('-------------------------------------', file=log)

454 parprint('Zero-point energy: %.3f eV' %

455 self.vibrations.get_zero_point_energy(),

456 file=log)