Coverage for /builds/debichem-team/python-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 (

9 FranckCondonOverlap,

10 FranckCondonRecursive,

11)

12from ase.vibrations.resonant_raman import ResonantRaman

15class Albrecht(ResonantRaman):

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

17 """

18 Parameters

19 ----------

20 all from ResonantRaman.__init__

21 combinations: int

22 Combinations to consider for multiple excitations.

23 Default is 1, possible 2

24 skip: int

25 Number of first transitions to exclude. Default 0,

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

27 nm: int

28 Number of intermediate m levels to consider, default 20

29 """

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

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

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

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

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

37 self.set_approximation(approximation)

39 def set_approximation(self, value):

40 approx = value.lower()

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

42 if not self.overlap:

43 raise ValueError('Overlaps are needed')

44 elif approx != 'albrecht a':

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

46 self._approx = value

48 def calculate_energies_and_modes(self):

49 if hasattr(self, 'im_r'):

50 return

52 ResonantRaman.calculate_energies_and_modes(self)

54 # single transitions and their occupation

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

56 om_v = om_Q

57 ndof = len(om_Q)

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

60 l_Q = range(ndof)

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

63 if self.combinations > 1:

64 if self.combinations != 2:

65 raise NotImplementedError

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

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

70 nv = len(ind_v)

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

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

73 for j, wt in enumerate(ind_v):

74 for i in wt:

75 n_vQ[j, i] += 1

76 om_v = n_vQ.dot(om_Q)

78 self.ind_v = ind_v

79 self.om_v = om_v

80 self.n_vQ = n_vQ # how many of each

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

83 def get_energies(self):

84 self.calculate_energies_and_modes()

85 return self.om_v

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

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

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

90 return arr_ro

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

92 if len(arr_ro):

93 data_ro[self.slize] = arr_ro

94 self.comm.sum(data_ro)

95 return data_ro

97 def Huang_Rhys_factors(self, forces_r):

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

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

100

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

102 """Evaluate unitless displacements from forces

103

104 Parameters

105 ----------

106 forces_r: array

107 Forces in cartesian coordinates

108 mineigv: float

109 Minimal Eigenvalue to consider in matrix inversion to handle

110 numerical noise. Default 1e-12

111

112 Returns

113 -------

114 Unitless displacements in Eigenmode coordinates

115 """

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

117

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

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

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

121 # there might be zero or nearly zero eigenvalues

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

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

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

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

126 forces_r.flat * self.im_r)

127

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

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

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

131

132 return d_Q

133

134 def omegaLS(self, omega, gamma):

135 omL = omega + 1j * gamma

136 omS_Q = omL - self.om_Q

137 return omL, omS_Q

138

139 def init_parallel_excitations(self):

140 """Init for paralellization over excitations."""

141 n_p = len(self.ex0E_p)

142

143 # collect excited state forces

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

145

146 # select your work load

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

148 rank = self.comm.rank

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

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

151

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

153 """Evaluate Albrecht A term.

154

155 Returns

156 -------

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

158 """

159 self.read()

160

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

162 self.fcr = FranckCondonRecursive()

163

164 omL = omega + 1j * gamma

165 omS_Q = omL - self.om_Q

166

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

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

169

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

171 for p in myp:

172 energy = self.ex0E_p[p]

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

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

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

176

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

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

179 for m in range(self.nm):

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

181 e_Q = energy_Q + m * self.om_Q

182 wm_Q += fco_Q / (e_Q - omL)

183 wp_Q += fco_Q / (e_Q + omS_Q)

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

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

186 self.comm.sum(m_Qcc)

187

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

189

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

191 """Evaluate Albrecht A term.

192

193 Returns

194 -------

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

196 """

197 self.read()

198

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

200 self.fcr = FranckCondonRecursive()

201

202 omL = omega + 1j * gamma

203 omS_v = omL - self.om_v

204 nv = len(self.om_v)

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

206 nQ = len(om_Q)

207

208 # n_v:

209 # how many FC factors are involved

210 # nvib_ov:

211 # delta functions to switch contributions depending on order o

212 # ind_ov:

213 # Q indicees

214 # n_ov:

215 # # of vibrational excitations

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

217

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

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

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

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

222 for o in range(self.combinations):

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

224 for v in range(nv):

225 try:

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

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

228 except IndexError:

229 pass

230 # XXXX change ????

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

232 n_ov[1] = nvib_ov[1]

233

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

235

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

237 for p in myp:

238 energy = self.ex0E_p[p]

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

240 S_Q = d_Q**2 / 2.

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

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

243

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

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

246 for m in range(self.nm):

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

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

249

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

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

252 for m in range(self.nm):

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

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

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

256

257 em_v = energy_v + m * om_ov[0]

258 # multiples of same kind

259 fco_v = nvib_ov[0] * fco1_v

260 wm_v += fco_v / (em_v - omL)

261 wp_v += fco_v / (em_v + omS_v)

262 if nvib_ov[1].any():

263 # multiples of mixed type

264 for n in range(self.nm):

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

266 e_v = em_v + n * om_ov[1]

267 fco_v = nvib_ov[1] * fco1_v * fco2_v

268 wm_v += fco_v / (e_v - omL)

269 wp_v += fco_v / (e_v + omS_v)

270

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

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

273 self.comm.sum(m_vcc)

274

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

276

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

278 term='BC'):

279 """Evaluate Albrecht BC term.

280

281 Returns

282 -------

283 Full Albrecht BC matrix element.

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

285 """

286 self.read()

287

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

289 self.fco = FranckCondonOverlap()

290

291 omL = omega + 1j * gamma

292 omS_Q = omL - self.om_Q

293

294 # excited state forces

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

296 # derivatives after normal coordinates

297 exdmdr_rpc = self._collect_r(

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

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

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

301

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

303 for p in myp:

304 energy = self.ex0E_p[p]

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

306 # relaxed excited state energy

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

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

309 energy_Q = energy - self.om_Q * S_Q

310

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

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

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

314

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

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

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

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

319 for m in range(self.nm):

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

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

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

323

324 em_Q = energy_Q + m * self.om_Q

325 wBLS_Q += f0mmQ1_Q / (em_Q - omL)

326 wBSL_Q += f0Qmm1_Q / (em_Q - omL)

327 wCLS_Q += f0mmQ1_Q / (em_Q + omS_Q)

328 wCSL_Q += f0Qmm1_Q / (em_Q + omS_Q)

329

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

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

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

333 if 'B' in term:

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

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

336 if 'C' in term:

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

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

339

340 self.comm.sum(me_Qcc)

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

342

343 def electronic_me_Qcc(self, omega, gamma):

344 self.calculate_energies_and_modes()

345

346 approx = self.approximation.lower()

347 assert self.combinations == 1

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

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

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

351 # divide through pre-factor

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

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

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

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

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

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

358 if approx == 'albrecht b':

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

360 if approx == 'albrecht c':

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

362

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

364

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

366

367 def me_Qcc(self, omega, gamma):

368 """Full matrix element"""

369 self.read()

370 approx = self.approximation.lower()

371 nv = len(self.om_v)

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

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

374 if self.combinations == 1:

375 # e^2 Angstrom^2 / eV

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

377 else:

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

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

380 if self.combinations == 1:

381 vel_vcc = self.meBC(omega, gamma)

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

383 else:

384 vel_vcc = self.meBCmult(omega, gamma)

385 V_vcc = 0

386 elif approx == 'albrecht b':

387 assert self.combinations == 1

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

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

390 if approx == 'albrecht c':

391 assert self.combinations == 1

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

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

394

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

396

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

398 method='standard', direction='central',

399 log=sys.stdout):

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

401 if self.combinations > 1:

402 return self.extended_summary()

403

404 om_v = self.get_energies()

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

406

407 if isinstance(log, str):

408 log = paropen(log, 'a')

409

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

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

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

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

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

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

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

417 for n, e in enumerate(om_v):

418 if e.imag != 0:

419 c = 'i'

420 e = e.imag

421 else:

422 c = ' '

423 e = e.real

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

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

426 file=log)

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

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

429 self.vibrations.get_zero_point_energy(),

430 file=log)

431

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

433 method='standard', direction='central',

434 log=sys.stdout):

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

436 self.read(method, direction)

437 om_v = self.get_energies()

438 intens_v = self.intensity(omega, gamma)

439

440 if isinstance(log, str):

441 log = paropen(log, 'a')

442

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

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

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

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

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

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

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

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

451 for v, e in enumerate(om_v):

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

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

454 file=log)

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

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

457 self.vibrations.get_zero_point_energy(),

458 file=log)