Coverage for /builds/kinetik161/ase/ase/vibrations/franck_condon.py: 96.28%
242 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1from functools import reduce
2from itertools import chain, combinations
3from math import factorial
4from operator import mul
6import numpy as np
8from ase.units import C, _hbar, kB, kg
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