r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1from functools import reduce

2from itertools import combinations, chain

3from math import factorial

4from operator import mul

6import numpy as np

8from ase.units import kg, C, _hbar, kB

9from ase.vibrations import Vibrations

12class Factorial:

13 def __init__(self):

14 self._fac = [1]

15 self._inv = [1.]

17 def __call__(self, n):

18 try:

19 return self._fac[n]

20 except IndexError:

21 for i in range(len(self._fac), n + 1):

22 self._fac.append(i * self._fac[i - 1])

23 try:

24 self._inv.append(float(1. / self._fac[-1]))

25 except OverflowError:

26 self._inv.append(0.)

27 return self._fac[n]

29 def inv(self, n):

30 self(n)

31 return self._inv[n]

34class FranckCondonOverlap:

35 """Evaluate squared overlaps depending on the Huang-Rhys parameter."""

37 def __init__(self):

38 self.factorial = Factorial()

40 def directT0(self, n, S):

41 """|<0|n>|^2

43 Direct squared Franck-Condon overlap corresponding to T=0.

44 """

45 return np.exp(-S) * S**n * self.factorial.inv(n)

47 def direct(self, n, m, S_in):

48 """|<n|m>|^2

50 Direct squared Franck-Condon overlap.

51 """

52 if n > m:

53 # use symmetry

54 return self.direct(m, n, S_in)

56 S = np.array([S_in])

57 mask = np.where(S == 0)

58 S[mask] = 1 # hide zeros

59 s = 0

60 for k in range(n + 1):

61 s += (-1)**(n - k) * S**float(-k) / (

62 self.factorial(k) *

63 self.factorial(n - k) * self.factorial(m - k))

64 res = np.exp(-S) * S**(n + m) * s**2 * (

65 self.factorial(n) * self.factorial(m))

66 # use othogonality

67 res[mask] = int(n == m)

68 return res[0]

70 def direct0mm1(self, m, S):

71 """<0|m><m|1>"""

72 sum = S**m

73 if m:

74 sum -= m * S**(m - 1)

75 return np.exp(-S) * np.sqrt(S) * sum * self.factorial.inv(m)

77 def direct0mm2(self, m, S):

78 """<0|m><m|2>"""

79 sum = S**(m + 1)

80 if m >= 1:

81 sum -= 2 * m * S**m

82 if m >= 2:

83 sum += m * (m - 1) * S**(m - 1)

84 return np.exp(-S) / np.sqrt(2) * sum * self.factorial.inv(m)

87class FranckCondonRecursive:

88 """Recursive implementation of Franck-Condon overlaps

90 Notes

91 -----

92 The overlaps are signed according to the sign of the displacements.

94 Reference

95 ---------

96 Julien Guthmuller

97 The Journal of Chemical Physics 144, 064106 (2016); doi: 10.1063/1.4941449

98 """

100 def __init__(self):

101 self.factorial = Factorial()

103 def ov0m(self, m, delta):

104 if m == 0:

105 return np.exp(-0.25 * delta**2)

106 else:

107 assert m > 0

108 return - delta / np.sqrt(2 * m) * self.ov0m(m - 1, delta)

110 def ov1m(self, m, delta):

111 sum = delta * self.ov0m(m, delta) / np.sqrt(2.)

112 if m == 0:

113 return sum

114 else:

115 assert m > 0

116 return sum + np.sqrt(m) * self.ov0m(m - 1, delta)

118 def ov2m(self, m, delta):

119 sum = delta * self.ov1m(m, delta) / 2

120 if m == 0:

121 return sum

122 else:

123 assert m > 0

124 return sum + np.sqrt(m / 2.) * self.ov1m(m - 1, delta)

126 def ov3m(self, m, delta):

127 sum = delta * self.ov2m(m, delta) / np.sqrt(6.)

128 if m == 0:

129 return sum

130 else:

131 assert m > 0

132 return sum + np.sqrt(m / 3.) * self.ov2m(m - 1, delta)

134 def ov0mm1(self, m, delta):

135 if m == 0:

136 return delta / np.sqrt(2) * self.ov0m(m, delta)**2

137 else:

138 return delta / np.sqrt(2) * (

139 self.ov0m(m, delta)**2 - self.ov0m(m - 1, delta)**2)

141 def direct0mm1(self, m, delta):

142 """direct and fast <0|m><m|1>"""

143 S = delta**2 / 2.

144 sum = S**m

145 if m:

146 sum -= m * S**(m - 1)

147 return np.where(S == 0, 0,

148 (np.exp(-S) * delta / np.sqrt(2) * sum *

149 self.factorial.inv(m)))

151 def ov0mm2(self, m, delta):

152 if m == 0:

153 return delta**2 / np.sqrt(8) * self.ov0m(m, delta)**2

154 elif m == 1:

155 return delta**2 / np.sqrt(8) * (

156 self.ov0m(m, delta)**2 - 2 * self.ov0m(m - 1, delta)**2)

157 else:

158 return delta**2 / np.sqrt(8) * (

159 self.ov0m(m, delta)**2 - 2 * self.ov0m(m - 1, delta)**2 +

160 self.ov0m(m - 2, delta)**2)

162 def direct0mm2(self, m, delta):

163 """direct and fast <0|m><m|2>"""

164 S = delta**2 / 2.

165 sum = S**(m + 1)

166 if m >= 1:

167 sum -= 2 * m * S**m

168 if m >= 2:

169 sum += m * (m - 1) * S**(m - 1)

170 return np.exp(-S) / np.sqrt(2) * sum * self.factorial.inv(m)

172 def ov1mm2(self, m, delta):

173 p1 = delta**3 / 4.

174 sum = p1 * self.ov0m(m, delta)**2

175 if m == 0:

176 return sum

177 p2 = delta - 3. * delta**3 / 4

178 sum += p2 * self.ov0m(m - 1, delta)**2

179 if m == 1:

180 return sum

181 sum -= p2 * self.ov0m(m - 2, delta)**2

182 if m == 2:

183 return sum

184 return sum - p1 * self.ov0m(m - 3, delta)**2

186 def direct1mm2(self, m, delta):

187 S = delta**2 / 2.

188 sum = S**2

189 if m > 0:

190 sum -= 2 * m * S

191 if m > 1:

192 sum += m * (m - 1)

193 with np.errstate(divide='ignore', invalid='ignore'):

194 return np.where(S == 0, 0,

195 (np.exp(-S) * S**(m - 1) / delta

196 * (S - m) * sum * self.factorial.inv(m)))

198 def direct0mm3(self, m, delta):

199 S = delta**2 / 2.

200 with np.errstate(divide='ignore', invalid='ignore'):

201 return np.where(

202 S == 0, 0,

203 (np.exp(-S) * S**(m - 1) / delta * np.sqrt(12.) *

204 (S**3 / 6. - m * S**2 / 2

205 + m * (m - 1) * S / 2. - m * (m - 1) * (m - 2) / 6)

206 * self.factorial.inv(m)))

209class FranckCondon:

210 def __init__(self, atoms, vibname, minfreq=-np.inf, maxfreq=np.inf):

211 """Input is a atoms object and the corresponding vibrations.

212 With minfreq and maxfreq frequencies can

213 be excluded from the calculation"""

215 self.atoms = atoms

216 # V = a * v is the combined atom and xyz-index

217 self.mm05_V = np.repeat(1. / np.sqrt(atoms.get_masses()), 3)

218 self.minfreq = minfreq

219 self.maxfreq = maxfreq

220 self.shape = (len(self.atoms), 3)

222 vib = Vibrations(atoms, name=vibname)

223 self.energies = np.real(vib.get_energies(method='frederiksen')) # [eV]

224 self.frequencies = np.real(

225 vib.get_frequencies(method='frederiksen')) # [cm^-1]

226 self.modes = vib.modes

227 self.H = vib.H

229 def get_Huang_Rhys_factors(self, forces):

230 """Evaluate Huang-Rhys factors and corresponding frequencies

231 from forces on atoms in the exited electronic state.

232 The double harmonic approximation is used. HR factors are

233 the first approximation of FC factors,

234 no combinations or higher quanta (>1) exitations are considered"""

236 assert forces.shape == self.shape

238 # Hesse matrix

239 H_VV = self.H

240 # sqrt of inverse mass matrix

241 mm05_V = self.mm05_V

242 # mass weighted Hesse matrix

243 Hm_VV = mm05_V[:, None] * H_VV * mm05_V

244 # mass weighted displacements

245 Fm_V = forces.flat * mm05_V

246 X_V = np.linalg.solve(Hm_VV, Fm_V)

247 # projection onto the modes

248 modes_VV = self.modes

249 d_V = np.dot(modes_VV, X_V)

250 # Huang-Rhys factors S

251 s = 1.e-20 / kg / C / _hbar**2 # SI units

252 S_V = s * d_V**2 * self.energies / 2

254 # reshape for minfreq

255 indices = np.where(self.frequencies <= self.minfreq)

256 np.append(indices, np.where(self.frequencies >= self.maxfreq))

257 S_V = np.delete(S_V, indices)

258 frequencies = np.delete(self.frequencies, indices)

260 return S_V, frequencies

262 def get_Franck_Condon_factors(self, temperature, forces, order=1):

263 """Return FC factors and corresponding frequencies up to given order.

265 Parameters

266 ----------

267 temperature: float

268 Temperature in K. Vibronic levels are occupied by a

269 Boltzman distribution.

270 forces: array

271 Forces on atoms in the exited electronic state

272 order: int

273 number of quanta taken into account, default

275 Returns

276 --------

277 FC: 3 entry list

278 FC[0] = FC factors for 0-0 and +-1 vibrational quantum

279 FC[1] = FC factors for +-2 vibrational quanta

280 FC[2] = FC factors for combinations

281 frequencies: 3 entry list

282 frequencies[0] correspond to FC[0]

283 frequencies[1] correspond to FC[1]

284 frequencies[2] correspond to FC[2]

285 """

286 S, f = self.get_Huang_Rhys_factors(forces)

287 assert order > 0

288 n = order + 1

289 T = temperature

290 freq = np.array(f)

292 # frequencies and their multiples

293 freq_n = [[] * i for i in range(n - 1)]

294 freq_neg = [[] * i for i in range(n - 1)]

296 for i in range(1, n):

297 freq_n[i - 1] = freq * i

298 freq_neg[i - 1] = freq * (-i)

300 # combinations

301 freq_nn = [x for x in combinations(chain(*freq_n), 2)]

302 for i in range(len(freq_nn)):

303 freq_nn[i] = freq_nn[i][0] + freq_nn[i][1]

305 indices2 = []

306 for i, y in enumerate(freq):

307 ind = [j for j, x in enumerate(freq_nn) if y == 0 or x % y == 0]

308 indices2.append(ind)

309 indices2 = [x for x in chain(*indices2)]

310 freq_nn = np.delete(freq_nn, indices2)

312 frequencies = [[] * x for x in range(3)]

313 frequencies[0].append(freq_neg[0])

314 frequencies[0].append([0])

315 frequencies[0].append(freq_n[0])

316 frequencies[0] = [x for x in chain(*frequencies[0])]

318 for i in range(1, n - 1):

319 frequencies[1].append(freq_neg[i])

320 frequencies[1].append(freq_n[i])

321 frequencies[1] = [x for x in chain(*frequencies[1])]

323 frequencies[2] = freq_nn

325 # Franck-Condon factors

326 E = freq / 8065.5

327 f_n = [[] * i for i in range(n)]

329 for j in range(0, n):

330 f_n[j] = np.exp(-E * j / (kB * T))

332 # partition function

333 Z = np.empty(len(S))

334 Z = np.sum(f_n, 0)

336 # occupation probability

337 w_n = [[] * k for k in range(n)]

338 for lval in range(n):

339 w_n[lval] = f_n[lval] / Z

341 # overlap wavefunctions

342 O_n = [[] * m for m in range(n)]

343 O_neg = [[] * m for m in range(n)]

344 for o in range(n):

345 O_n[o] = [[] * p for p in range(n)]

346 O_neg[o] = [[] * p for p in range(n - 1)]

347 for q in range(o, n + o):

348 a = np.minimum(o, q)

349 summe = []

350 for k in range(a + 1):

351 s = ((-1)**(q - k) * np.sqrt(S)**(o + q - 2 * k) *

352 factorial(o) * factorial(q) /

353 (factorial(k) * factorial(o - k) * factorial(q - k)))

354 summe.append(s)

355 summe = np.sum(summe, 0)

356 O_n[o][q - o] = (np.exp(-S / 2) /

357 (factorial(o) * factorial(q))**(0.5) *

358 summe)**2 * w_n[o]

359 for q in range(n - 1):

360 O_neg[o][q] = [0 * b for b in range(len(S))]

361 for q in range(o - 1, -1, -1):

362 a = np.minimum(o, q)

363 summe = []

364 for k in range(a + 1):

365 s = ((-1)**(q - k) * np.sqrt(S)**(o + q - 2 * k) *

366 factorial(o) * factorial(q) /

367 (factorial(k) * factorial(o - k) * factorial(q - k)))

368 summe.append(s)

369 summe = np.sum(summe, 0)

370 O_neg[o][q] = (np.exp(-S / 2) /

371 (factorial(o) * factorial(q))**(0.5) *

372 summe)**2 * w_n[o]

373 O_neg = np.delete(O_neg, 0, 0)

375 # Franck-Condon factors

376 FC_n = [[] * i for i in range(n)]

377 FC_n = np.sum(O_n, 0)

378 zero = reduce(mul, FC_n[0])

379 FC_neg = [[] * i for i in range(n - 2)]

380 FC_neg = np.sum(O_neg, 0)

381 FC_n = np.delete(FC_n, 0, 0)

383 # combination FC factors

384 FC_nn = [x for x in combinations(chain(*FC_n), 2)]

385 for i in range(len(FC_nn)):

386 FC_nn[i] = FC_nn[i][0] * FC_nn[i][1]

388 FC_nn = np.delete(FC_nn, indices2)

390 FC = [[] * x for x in range(3)]

391 FC[0].append(FC_neg[0])

392 FC[0].append([zero])

393 FC[0].append(FC_n[0])

394 FC[0] = [x for x in chain(*FC[0])]

396 for i in range(1, n - 1):

397 FC[1].append(FC_neg[i])

398 FC[1].append(FC_n[i])

399 FC[1] = [x for x in chain(*FC[1])]

401 FC[2] = FC_nn

403 return FC, frequencies