Coverage for /builds/debichem-team/python-ase/ase/vibrations/franck_condon.py : 96.28%
Hot-keys on this page
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."""
36 def __init__(self):
37 self.factorial = Factorial()
39 def directT0(self, n, S):
40 """|<0|n>|^2
42 Direct squared Franck-Condon overlap corresponding to T=0.
43 """
44 return np.exp(-S) * S**n * self.factorial.inv(n)
46 def direct(self, n, m, S_in):
47 """|<n|m>|^2
49 Direct squared Franck-Condon overlap.
50 """
51 if n > m:
52 # use symmetry
53 return self.direct(m, n, S_in)
55 S = np.array([S_in])
56 mask = np.where(S == 0)
57 S[mask] = 1 # hide zeros
58 s = 0
59 for k in range(n + 1):
60 s += (-1)**(n - k) * S**float(-k) / (
61 self.factorial(k) *
62 self.factorial(n - k) * self.factorial(m - k))
63 res = np.exp(-S) * S**(n + m) * s**2 * (
64 self.factorial(n) * self.factorial(m))
65 # use othogonality
66 res[mask] = int(n == m)
67 return res[0]
69 def direct0mm1(self, m, S):
70 """<0|m><m|1>"""
71 sum = S**m
72 if m:
73 sum -= m * S**(m - 1)
74 return np.exp(-S) * np.sqrt(S) * sum * self.factorial.inv(m)
76 def direct0mm2(self, m, S):
77 """<0|m><m|2>"""
78 sum = S**(m + 1)
79 if m >= 1:
80 sum -= 2 * m * S**m
81 if m >= 2:
82 sum += m * (m - 1) * S**(m - 1)
83 return np.exp(-S) / np.sqrt(2) * sum * self.factorial.inv(m)
86class FranckCondonRecursive:
87 """Recursive implementation of Franck-Condon overlaps
89 Notes
90 -----
91 The ovelaps are signed according to the sign of the displacements.
93 Reference
94 ---------
95 Julien Guthmuller
96 The Journal of Chemical Physics 144, 064106 (2016); doi: 10.1063/1.4941449
97 """
98 def __init__(self):
99 self.factorial = Factorial()
101 def ov0m(self, m, delta):
102 if m == 0:
103 return np.exp(-0.25 * delta**2)
104 else:
105 assert(m > 0)
106 return - delta / np.sqrt(2 * m) * self.ov0m(m - 1, delta)
108 def ov1m(self, m, delta):
109 sum = delta * self.ov0m(m, delta) / np.sqrt(2.)
110 if m == 0:
111 return sum
112 else:
113 assert(m > 0)
114 return sum + np.sqrt(m) * self.ov0m(m - 1, delta)
116 def ov2m(self, m, delta):
117 sum = delta * self.ov1m(m, delta) / 2
118 if m == 0:
119 return sum
120 else:
121 assert(m > 0)
122 return sum + np.sqrt(m / 2.) * self.ov1m(m - 1, delta)
124 def ov3m(self, m, delta):
125 sum = delta * self.ov2m(m, delta) / np.sqrt(6.)
126 if m == 0:
127 return sum
128 else:
129 assert(m > 0)
130 return sum + np.sqrt(m / 3.) * self.ov2m(m - 1, delta)
132 def ov0mm1(self, m, delta):
133 if m == 0:
134 return delta / np.sqrt(2) * self.ov0m(m, delta)**2
135 else:
136 return delta / np.sqrt(2) * (
137 self.ov0m(m, delta)**2 - self.ov0m(m - 1, delta)**2)
139 def direct0mm1(self, m, delta):
140 """direct and fast <0|m><m|1>"""
141 S = delta**2 / 2.
142 sum = S**m
143 if m:
144 sum -= m * S**(m - 1)
145 return np.where(S == 0, 0,
146 (np.exp(-S) * delta / np.sqrt(2) * sum *
147 self.factorial.inv(m)))
149 def ov0mm2(self, m, delta):
150 if m == 0:
151 return delta**2 / np.sqrt(8) * self.ov0m(m, delta)**2
152 elif m == 1:
153 return delta**2 / np.sqrt(8) * (
154 self.ov0m(m, delta)**2 - 2 * self.ov0m(m - 1, delta)**2)
155 else:
156 return delta**2 / np.sqrt(8) * (
157 self.ov0m(m, delta)**2 - 2 * self.ov0m(m - 1, delta)**2 +
158 self.ov0m(m - 2, delta)**2)
160 def direct0mm2(self, m, delta):
161 """direct and fast <0|m><m|2>"""
162 S = delta**2 / 2.
163 sum = S**(m + 1)
164 if m >= 1:
165 sum -= 2 * m * S**m
166 if m >= 2:
167 sum += m * (m - 1) * S**(m - 1)
168 return np.exp(-S) / np.sqrt(2) * sum * self.factorial.inv(m)
170 def ov1mm2(self, m, delta):
171 p1 = delta**3 / 4.
172 sum = p1 * self.ov0m(m, delta)**2
173 if m == 0:
174 return sum
175 p2 = delta - 3. * delta**3 / 4
176 sum += p2 * self.ov0m(m - 1, delta)**2
177 if m == 1:
178 return sum
179 sum -= p2 * self.ov0m(m - 2, delta)**2
180 if m == 2:
181 return sum
182 return sum - p1 * self.ov0m(m - 3, delta)**2
184 def direct1mm2(self, m, delta):
185 S = delta**2 / 2.
186 sum = S**2
187 if m > 0:
188 sum -= 2 * m * S
189 if m > 1:
190 sum += m * (m - 1)
191 with np.errstate(divide='ignore', invalid='ignore'):
192 return np.where(S == 0, 0,
193 (np.exp(-S) * S**(m - 1) / delta
194 * (S - m) * sum * self.factorial.inv(m)))
196 def direct0mm3(self, m, delta):
197 S = delta**2 / 2.
198 with np.errstate(divide='ignore', invalid='ignore'):
199 return np.where(
200 S == 0, 0,
201 (np.exp(-S) * S**(m - 1) / delta * np.sqrt(12.) *
202 (S**3 / 6. - m * S**2 / 2
203 + m * (m - 1) * S / 2. - m * (m - 1) * (m - 2) / 6)
204 * self.factorial.inv(m)))
207class FranckCondon:
208 def __init__(self, atoms, vibname, minfreq=-np.inf, maxfreq=np.inf):
209 """Input is a atoms object and the corresponding vibrations.
210 With minfreq and maxfreq frequencies can
211 be excluded from the calculation"""
213 self.atoms = atoms
214 # V = a * v is the combined atom and xyz-index
215 self.mm05_V = np.repeat(1. / np.sqrt(atoms.get_masses()), 3)
216 self.minfreq = minfreq
217 self.maxfreq = maxfreq
218 self.shape = (len(self.atoms), 3)
220 vib = Vibrations(atoms, name=vibname)
221 self.energies = np.real(vib.get_energies(method='frederiksen')) # [eV]
222 self.frequencies = np.real(
223 vib.get_frequencies(method='frederiksen')) # [cm^-1]
224 self.modes = vib.modes
225 self.H = vib.H
227 def get_Huang_Rhys_factors(self, forces):
228 """Evaluate Huang-Rhys factors and corresponding frequencies
229 from forces on atoms in the exited electronic state.
230 The double harmonic approximation is used. HR factors are
231 the first approximation of FC factors,
232 no combinations or higher quanta (>1) exitations are considered"""
234 assert(forces.shape == self.shape)
236 # Hesse matrix
237 H_VV = self.H
238 # sqrt of inverse mass matrix
239 mm05_V = self.mm05_V
240 # mass weighted Hesse matrix
241 Hm_VV = mm05_V[:, None] * H_VV * mm05_V
242 # mass weighted displacements
243 Fm_V = forces.flat * mm05_V
244 X_V = np.linalg.solve(Hm_VV, Fm_V)
245 # projection onto the modes
246 modes_VV = self.modes
247 d_V = np.dot(modes_VV, X_V)
248 # Huang-Rhys factors S
249 s = 1.e-20 / kg / C / _hbar**2 # SI units
250 S_V = s * d_V**2 * self.energies / 2
252 # reshape for minfreq
253 indices = np.where(self.frequencies <= self.minfreq)
254 np.append(indices, np.where(self.frequencies >= self.maxfreq))
255 S_V = np.delete(S_V, indices)
256 frequencies = np.delete(self.frequencies, indices)
258 return S_V, frequencies
260 def get_Franck_Condon_factors(self, temperature, forces, order=1):
261 """Return FC factors and corresponding frequencies up to given order.
263 Parameters
264 ----------
265 temperature: float
266 Temperature in K. Vibronic levels are occupied by a
267 Boltzman distribution.
268 forces: array
269 Forces on atoms in the exited electronic state
270 order: int
271 number of quanta taken into account, default
273 Returns
274 --------
275 FC: 3 entry list
276 FC[0] = FC factors for 0-0 and +-1 vibrational quantum
277 FC[1] = FC factors for +-2 vibrational quanta
278 FC[2] = FC factors for combinations
279 frequencies: 3 entry list
280 frequencies[0] correspond to FC[0]
281 frequencies[1] correspond to FC[1]
282 frequencies[2] correspond to FC[2]
283 """
284 S, f = self.get_Huang_Rhys_factors(forces)
285 assert order > 0
286 n = order + 1
287 T = temperature
288 freq = np.array(f)
290 # frequencies and their multiples
291 freq_n = [[] * i for i in range(n - 1)]
292 freq_neg = [[] * i for i in range(n - 1)]
294 for i in range(1, n):
295 freq_n[i - 1] = freq * i
296 freq_neg[i - 1] = freq * (-i)
298 # combinations
299 freq_nn = [x for x in combinations(chain(*freq_n), 2)]
300 for i in range(len(freq_nn)):
301 freq_nn[i] = freq_nn[i][0] + freq_nn[i][1]
303 indices2 = []
304 for i, y in enumerate(freq):
305 ind = [j for j, x in enumerate(freq_nn) if y == 0 or x % y == 0]
306 indices2.append(ind)
307 indices2 = [x for x in chain(*indices2)]
308 freq_nn = np.delete(freq_nn, indices2)
310 frequencies = [[] * x for x in range(3)]
311 frequencies[0].append(freq_neg[0])
312 frequencies[0].append([0])
313 frequencies[0].append(freq_n[0])
314 frequencies[0] = [x for x in chain(*frequencies[0])]
316 for i in range(1, n - 1):
317 frequencies[1].append(freq_neg[i])
318 frequencies[1].append(freq_n[i])
319 frequencies[1] = [x for x in chain(*frequencies[1])]
321 frequencies[2] = freq_nn
323 # Franck-Condon factors
324 E = freq / 8065.5
325 f_n = [[] * i for i in range(n)]
327 for j in range(0, n):
328 f_n[j] = np.exp(-E * j / (kB * T))
330 # partition function
331 Z = np.empty(len(S))
332 Z = np.sum(f_n, 0)
334 # occupation probability
335 w_n = [[] * k for k in range(n)]
336 for l in range(n):
337 w_n[l] = f_n[l] / Z
339 # overlap wavefunctions
340 O_n = [[] * m for m in range(n)]
341 O_neg = [[] * m for m in range(n)]
342 for o in range(n):
343 O_n[o] = [[] * p for p in range(n)]
344 O_neg[o] = [[] * p for p in range(n - 1)]
345 for q in range(o, n + o):
346 a = np.minimum(o, q)
347 summe = []
348 for k in range(a + 1):
349 s = ((-1)**(q - k) * np.sqrt(S)**(o + q - 2 * k) *
350 factorial(o) * factorial(q) /
351 (factorial(k) * factorial(o - k) * factorial(q - k)))
352 summe.append(s)
353 summe = np.sum(summe, 0)
354 O_n[o][q - o] = (np.exp(-S / 2) /
355 (factorial(o) * factorial(q))**(0.5) *
356 summe)**2 * w_n[o]
357 for q in range(n - 1):
358 O_neg[o][q] = [0 * b for b in range(len(S))]
359 for q in range(o - 1, -1, -1):
360 a = np.minimum(o, q)
361 summe = []
362 for k in range(a + 1):
363 s = ((-1)**(q - k) * np.sqrt(S)**(o + q - 2 * k) *
364 factorial(o) * factorial(q) /
365 (factorial(k) * factorial(o - k) * factorial(q - k)))
366 summe.append(s)
367 summe = np.sum(summe, 0)
368 O_neg[o][q] = (np.exp(-S / 2) /
369 (factorial(o) * factorial(q))**(0.5) *
370 summe)**2 * w_n[o]
371 O_neg = np.delete(O_neg, 0, 0)
373 # Franck-Condon factors
374 FC_n = [[] * i for i in range(n)]
375 FC_n = np.sum(O_n, 0)
376 zero = reduce(mul, FC_n[0])
377 FC_neg = [[] * i for i in range(n - 2)]
378 FC_neg = np.sum(O_neg, 0)
379 FC_n = np.delete(FC_n, 0, 0)
381 # combination FC factors
382 FC_nn = [x for x in combinations(chain(*FC_n), 2)]
383 for i in range(len(FC_nn)):
384 FC_nn[i] = FC_nn[i][0] * FC_nn[i][1]
386 FC_nn = np.delete(FC_nn, indices2)
388 FC = [[] * x for x in range(3)]
389 FC[0].append(FC_neg[0])
390 FC[0].append([zero])
391 FC[0].append(FC_n[0])
392 FC[0] = [x for x in chain(*FC[0])]
394 for i in range(1, n - 1):
395 FC[1].append(FC_neg[i])
396 FC[1].append(FC_n[i])
397 FC[1] = [x for x in chain(*FC[1])]
399 FC[2] = FC_nn
401 return FC, frequencies