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

1# fmt: off

3import sys

4from itertools import combinations_with_replacement

6import numpy as np

8import ase.units as u

9from ase.parallel import paropen, parprint

10from ase.vibrations.franck_condon import (

11 FranckCondonOverlap,

12 FranckCondonRecursive,

13)

14from ase.vibrations.resonant_raman import ResonantRaman

17class Albrecht(ResonantRaman):

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

19 """

20 Parameters

21 ----------

22 all from ResonantRaman.__init__

23 combinations: int

24 Combinations to consider for multiple excitations.

25 Default is 1, possible 2

26 skip: int

27 Number of first transitions to exclude. Default 0,

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

29 nm: int

30 Number of intermediate m levels to consider, default 20

31 """

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

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

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

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

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

39 self.set_approximation(approximation)

41 def set_approximation(self, value):

42 approx = value.lower()

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

44 if not self.overlap:

45 raise ValueError('Overlaps are needed')

46 elif approx != 'albrecht a':

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

48 self._approx = value

50 def calculate_energies_and_modes(self):

51 if hasattr(self, 'im_r'):

52 return

54 ResonantRaman.calculate_energies_and_modes(self)

56 # single transitions and their occupation

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

58 om_v = om_Q

59 ndof = len(om_Q)

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

62 l_Q = range(ndof)

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

65 if self.combinations > 1:

66 if self.combinations != 2:

67 raise NotImplementedError

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

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

72 nv = len(ind_v)

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

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

75 for j, wt in enumerate(ind_v):

76 for i in wt:

77 n_vQ[j, i] += 1

78 om_v = n_vQ.dot(om_Q)

80 self.ind_v = ind_v

81 self.om_v = om_v

82 self.n_vQ = n_vQ # how many of each

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

85 def get_energies(self):

86 self.calculate_energies_and_modes()

87 return self.om_v

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

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

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

92 return arr_ro

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

94 if len(arr_ro):

95 data_ro[self.slize] = arr_ro

96 self.comm.sum(data_ro)

97 return data_ro

99 def Huang_Rhys_factors(self, forces_r):

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

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

102

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

104 """Evaluate unitless displacements from forces

105

106 Parameters

107 ----------

108 forces_r: array

109 Forces in cartesian coordinates

110 mineigv: float

111 Minimal Eigenvalue to consider in matrix inversion to handle

112 numerical noise. Default 1e-12

113

114 Returns

115 -------

116 Unitless displacements in Eigenmode coordinates

117 """

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

119

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

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

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

123 # there might be zero or nearly zero eigenvalues

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

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

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

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

128 forces_r.flat * self.im_r)

129

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

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

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

133

134 return d_Q

135

136 def omegaLS(self, omega, gamma):

137 omL = omega + 1j * gamma

138 omS_Q = omL - self.om_Q

139 return omL, omS_Q

140

141 def init_parallel_excitations(self):

142 """Init for paralellization over excitations."""

143 n_p = len(self.ex0E_p)

144

145 # collect excited state forces

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

147

148 # select your work load

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

150 rank = self.comm.rank

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

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

153

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

155 """Evaluate Albrecht A term.

156

157 Returns

158 -------

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

160 """

161 self.read()

162

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

164 self.fcr = FranckCondonRecursive()

165

166 omL = omega + 1j * gamma

167 omS_Q = omL - self.om_Q

168

169 _n_p, myp, exF_pr = self.init_parallel_excitations()

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

171

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

173 for p in myp:

174 energy = self.ex0E_p[p]

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

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

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

178

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

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

181 for m in range(self.nm):

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

183 e_Q = energy_Q + m * self.om_Q

184 wm_Q += fco_Q / (e_Q - omL)

185 wp_Q += fco_Q / (e_Q + omS_Q)

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

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

188 self.comm.sum(m_Qcc)

189

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

191

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

193 """Evaluate Albrecht A term.

194

195 Returns

196 -------

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

198 """

199 self.read()

200

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

202 self.fcr = FranckCondonRecursive()

203

204 omL = omega + 1j * gamma

205 omS_v = omL - self.om_v

206 nv = len(self.om_v)

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

208 nQ = len(om_Q)

209

210 # n_v:

211 # how many FC factors are involved

212 # nvib_ov:

213 # delta functions to switch contributions depending on order o

214 # ind_ov:

215 # Q indicees

216 # n_ov:

217 # # of vibrational excitations

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

219

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

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

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

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

224 for o in range(self.combinations):

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

226 for v in range(nv):

227 try:

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

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

230 except IndexError:

231 pass

232 # XXXX change ????

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

234 n_ov[1] = nvib_ov[1]

235

236 _n_p, myp, exF_pr = self.init_parallel_excitations()

237

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

239 for p in myp:

240 energy = self.ex0E_p[p]

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

242 S_Q = d_Q**2 / 2.

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

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

245

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

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

248 for m in range(self.nm):

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

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

251

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

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

254 for m in range(self.nm):

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

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

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

258

259 em_v = energy_v + m * om_ov[0]

260 # multiples of same kind

261 fco_v = nvib_ov[0] * fco1_v

262 wm_v += fco_v / (em_v - omL)

263 wp_v += fco_v / (em_v + omS_v)

264 if nvib_ov[1].any():

265 # multiples of mixed type

266 for n in range(self.nm):

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

268 e_v = em_v + n * om_ov[1]

269 fco_v = nvib_ov[1] * fco1_v * fco2_v

270 wm_v += fco_v / (e_v - omL)

271 wp_v += fco_v / (e_v + omS_v)

272

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

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

275 self.comm.sum(m_vcc)

276

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

278

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

280 term='BC'):

281 """Evaluate Albrecht BC term.

282

283 Returns

284 -------

285 Full Albrecht BC matrix element.

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

287 """

288 self.read()

289

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

291 self.fco = FranckCondonOverlap()

292

293 omL = omega + 1j * gamma

294 omS_Q = omL - self.om_Q

295

296 # excited state forces

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

298 # derivatives after normal coordinates

299 exdmdr_rpc = self._collect_r(

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

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

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

303

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

305 for p in myp:

306 energy = self.ex0E_p[p]

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

308 # relaxed excited state energy

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

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

311 energy_Q = energy - self.om_Q * S_Q

312

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

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

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

316

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

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

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

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

321 for m in range(self.nm):

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

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

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

325

326 em_Q = energy_Q + m * self.om_Q

327 wBLS_Q += f0mmQ1_Q / (em_Q - omL)

328 wBSL_Q += f0Qmm1_Q / (em_Q - omL)

329 wCLS_Q += f0mmQ1_Q / (em_Q + omS_Q)

330 wCSL_Q += f0Qmm1_Q / (em_Q + omS_Q)

331

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

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

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

335 if 'B' in term:

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

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

338 if 'C' in term:

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

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

341

342 self.comm.sum(me_Qcc)

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

344

345 def electronic_me_Qcc(self, omega, gamma):

346 self.calculate_energies_and_modes()

347

348 approx = self.approximation.lower()

349 assert self.combinations == 1

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

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

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

353 # divide through pre-factor

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

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

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

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

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

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

360 if approx == 'albrecht b':

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

362 if approx == 'albrecht c':

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

364

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

366

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

368

369 def me_Qcc(self, omega, gamma):

370 """Full matrix element"""

371 self.read()

372 approx = self.approximation.lower()

373 nv = len(self.om_v)

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

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

376 if self.combinations == 1:

377 # e^2 Angstrom^2 / eV

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

379 else:

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

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

382 if self.combinations == 1:

383 vel_vcc = self.meBC(omega, gamma)

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

385 else:

386 vel_vcc = self.meBCmult(omega, gamma)

387 V_vcc = 0

388 elif approx == 'albrecht b':

389 assert self.combinations == 1

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

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

392 if approx == 'albrecht c':

393 assert self.combinations == 1

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

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

396

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

398

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

400 method='standard', direction='central',

401 log=sys.stdout):

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

403 if self.combinations > 1:

404 return self.extended_summary()

405

406 om_v = self.get_energies()

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

408

409 if isinstance(log, str):

410 log = paropen(log, 'a')

411

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

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

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

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

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

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

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

419 for n, e in enumerate(om_v):

420 if e.imag != 0:

421 c = 'i'

422 e = e.imag

423 else:

424 c = ' '

425 e = e.real

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

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

428 file=log)

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

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

431 self.vibrations.get_zero_point_energy(),

432 file=log)

433

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

435 method='standard', direction='central',

436 log=sys.stdout):

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

438 self.read(method, direction)

439 om_v = self.get_energies()

440 intens_v = self.intensity(omega, gamma)

441

442 if isinstance(log, str):

443 log = paropen(log, 'a')

444

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

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

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

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

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

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

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

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

453 for v, e in enumerate(om_v):

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

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

456 file=log)

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

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

459 self.vibrations.get_zero_point_energy(),

460 file=log)