Coverage for /builds/ase/ase/ase/vibrations/franck_condon.py: 96.28%
242 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3from functools import reduce
4from itertools import chain, combinations
5from math import factorial
6from operator import mul
8import numpy as np
10from ase.units import C, _hbar, kB, kg
11from ase.vibrations import Vibrations
14class Factorial:
15 def __init__(self):
16 self._fac = [1]
17 self._inv = [1.]
19 def __call__(self, n):
20 try:
21 return self._fac[n]
22 except IndexError:
23 for i in range(len(self._fac), n + 1):
24 self._fac.append(i * self._fac[i - 1])
25 try:
26 self._inv.append(float(1. / self._fac[-1]))
27 except OverflowError:
28 self._inv.append(0.)
29 return self._fac[n]
31 def inv(self, n):
32 self(n)
33 return self._inv[n]
36class FranckCondonOverlap:
37 """Evaluate squared overlaps depending on the Huang-Rhys parameter."""
39 def __init__(self):
40 self.factorial = Factorial()
42 def directT0(self, n, S):
43 """|<0|n>|^2
45 Direct squared Franck-Condon overlap corresponding to T=0.
46 """
47 return np.exp(-S) * S**n * self.factorial.inv(n)
49 def direct(self, n, m, S_in):
50 """|<n|m>|^2
52 Direct squared Franck-Condon overlap.
53 """
54 if n > m:
55 # use symmetry
56 return self.direct(m, n, S_in)
58 S = np.array([S_in])
59 mask = np.where(S == 0)
60 S[mask] = 1 # hide zeros
61 s = 0
62 for k in range(n + 1):
63 s += (-1)**(n - k) * S**float(-k) / (
64 self.factorial(k) *
65 self.factorial(n - k) * self.factorial(m - k))
66 res = np.exp(-S) * S**(n + m) * s**2 * (
67 self.factorial(n) * self.factorial(m))
68 # use othogonality
69 res[mask] = int(n == m)
70 return res[0]
72 def direct0mm1(self, m, S):
73 """<0|m><m|1>"""
74 sum = S**m
75 if m:
76 sum -= m * S**(m - 1)
77 return np.exp(-S) * np.sqrt(S) * sum * self.factorial.inv(m)
79 def direct0mm2(self, m, S):
80 """<0|m><m|2>"""
81 sum = S**(m + 1)
82 if m >= 1:
83 sum -= 2 * m * S**m
84 if m >= 2:
85 sum += m * (m - 1) * S**(m - 1)
86 return np.exp(-S) / np.sqrt(2) * sum * self.factorial.inv(m)
89class FranckCondonRecursive:
90 """Recursive implementation of Franck-Condon overlaps
92 Notes
93 -----
94 The overlaps are signed according to the sign of the displacements.
96 Reference
97 ---------
98 Julien Guthmuller
99 The Journal of Chemical Physics 144, 064106 (2016); doi: 10.1063/1.4941449
100 """
102 def __init__(self):
103 self.factorial = Factorial()
105 def ov0m(self, m, delta):
106 if m == 0:
107 return np.exp(-0.25 * delta**2)
108 else:
109 assert m > 0
110 return - delta / np.sqrt(2 * m) * self.ov0m(m - 1, delta)
112 def ov1m(self, m, delta):
113 sum = delta * self.ov0m(m, delta) / np.sqrt(2.)
114 if m == 0:
115 return sum
116 else:
117 assert m > 0
118 return sum + np.sqrt(m) * self.ov0m(m - 1, delta)
120 def ov2m(self, m, delta):
121 sum = delta * self.ov1m(m, delta) / 2
122 if m == 0:
123 return sum
124 else:
125 assert m > 0
126 return sum + np.sqrt(m / 2.) * self.ov1m(m - 1, delta)
128 def ov3m(self, m, delta):
129 sum = delta * self.ov2m(m, delta) / np.sqrt(6.)
130 if m == 0:
131 return sum
132 else:
133 assert m > 0
134 return sum + np.sqrt(m / 3.) * self.ov2m(m - 1, delta)
136 def ov0mm1(self, m, delta):
137 if m == 0:
138 return delta / np.sqrt(2) * self.ov0m(m, delta)**2
139 else:
140 return delta / np.sqrt(2) * (
141 self.ov0m(m, delta)**2 - self.ov0m(m - 1, delta)**2)
143 def direct0mm1(self, m, delta):
144 """direct and fast <0|m><m|1>"""
145 S = delta**2 / 2.
146 sum = S**m
147 if m:
148 sum -= m * S**(m - 1)
149 return np.where(S == 0, 0,
150 (np.exp(-S) * delta / np.sqrt(2) * sum *
151 self.factorial.inv(m)))
153 def ov0mm2(self, m, delta):
154 if m == 0:
155 return delta**2 / np.sqrt(8) * self.ov0m(m, delta)**2
156 elif m == 1:
157 return delta**2 / np.sqrt(8) * (
158 self.ov0m(m, delta)**2 - 2 * self.ov0m(m - 1, delta)**2)
159 else:
160 return delta**2 / np.sqrt(8) * (
161 self.ov0m(m, delta)**2 - 2 * self.ov0m(m - 1, delta)**2 +
162 self.ov0m(m - 2, delta)**2)
164 def direct0mm2(self, m, delta):
165 """direct and fast <0|m><m|2>"""
166 S = delta**2 / 2.
167 sum = S**(m + 1)
168 if m >= 1:
169 sum -= 2 * m * S**m
170 if m >= 2:
171 sum += m * (m - 1) * S**(m - 1)
172 return np.exp(-S) / np.sqrt(2) * sum * self.factorial.inv(m)
174 def ov1mm2(self, m, delta):
175 p1 = delta**3 / 4.
176 sum = p1 * self.ov0m(m, delta)**2
177 if m == 0:
178 return sum
179 p2 = delta - 3. * delta**3 / 4
180 sum += p2 * self.ov0m(m - 1, delta)**2
181 if m == 1:
182 return sum
183 sum -= p2 * self.ov0m(m - 2, delta)**2
184 if m == 2:
185 return sum
186 return sum - p1 * self.ov0m(m - 3, delta)**2
188 def direct1mm2(self, m, delta):
189 S = delta**2 / 2.
190 sum = S**2
191 if m > 0:
192 sum -= 2 * m * S
193 if m > 1:
194 sum += m * (m - 1)
195 with np.errstate(divide='ignore', invalid='ignore'):
196 return np.where(S == 0, 0,
197 (np.exp(-S) * S**(m - 1) / delta
198 * (S - m) * sum * self.factorial.inv(m)))
200 def direct0mm3(self, m, delta):
201 S = delta**2 / 2.
202 with np.errstate(divide='ignore', invalid='ignore'):
203 return np.where(
204 S == 0, 0,
205 (np.exp(-S) * S**(m - 1) / delta * np.sqrt(12.) *
206 (S**3 / 6. - m * S**2 / 2
207 + m * (m - 1) * S / 2. - m * (m - 1) * (m - 2) / 6)
208 * self.factorial.inv(m)))
211class FranckCondon:
212 def __init__(self, atoms, vibname, minfreq=-np.inf, maxfreq=np.inf):
213 """Input is a atoms object and the corresponding vibrations.
214 With minfreq and maxfreq frequencies can
215 be excluded from the calculation"""
217 self.atoms = atoms
218 # V = a * v is the combined atom and xyz-index
219 self.mm05_V = np.repeat(1. / np.sqrt(atoms.get_masses()), 3)
220 self.minfreq = minfreq
221 self.maxfreq = maxfreq
222 self.shape = (len(self.atoms), 3)
224 vib = Vibrations(atoms, name=vibname)
225 self.energies = np.real(vib.get_energies(method='frederiksen')) # [eV]
226 self.frequencies = np.real(
227 vib.get_frequencies(method='frederiksen')) # [cm^-1]
228 self.modes = vib.modes
229 self.H = vib.H
231 def get_Huang_Rhys_factors(self, forces):
232 """Evaluate Huang-Rhys factors and corresponding frequencies
233 from forces on atoms in the exited electronic state.
234 The double harmonic approximation is used. HR factors are
235 the first approximation of FC factors,
236 no combinations or higher quanta (>1) exitations are considered"""
238 assert forces.shape == self.shape
240 # Hesse matrix
241 H_VV = self.H
242 # sqrt of inverse mass matrix
243 mm05_V = self.mm05_V
244 # mass weighted Hesse matrix
245 Hm_VV = mm05_V[:, None] * H_VV * mm05_V
246 # mass weighted displacements
247 Fm_V = forces.flat * mm05_V
248 X_V = np.linalg.solve(Hm_VV, Fm_V)
249 # projection onto the modes
250 modes_VV = self.modes
251 d_V = np.dot(modes_VV, X_V)
252 # Huang-Rhys factors S
253 s = 1.e-20 / kg / C / _hbar**2 # SI units
254 S_V = s * d_V**2 * self.energies / 2
256 # reshape for minfreq
257 indices = np.where(self.frequencies <= self.minfreq)
258 np.append(indices, np.where(self.frequencies >= self.maxfreq))
259 S_V = np.delete(S_V, indices)
260 frequencies = np.delete(self.frequencies, indices)
262 return S_V, frequencies
264 def get_Franck_Condon_factors(self, temperature, forces, order=1):
265 """Return FC factors and corresponding frequencies up to given order.
267 Parameters
268 ----------
269 temperature: float
270 Temperature in K. Vibronic levels are occupied by a
271 Boltzman distribution.
272 forces: array
273 Forces on atoms in the exited electronic state
274 order: int
275 number of quanta taken into account, default
277 Returns
278 --------
279 FC: 3 entry list
280 FC[0] = FC factors for 0-0 and +-1 vibrational quantum
281 FC[1] = FC factors for +-2 vibrational quanta
282 FC[2] = FC factors for combinations
283 frequencies: 3 entry list
284 frequencies[0] correspond to FC[0]
285 frequencies[1] correspond to FC[1]
286 frequencies[2] correspond to FC[2]
287 """
288 S, f = self.get_Huang_Rhys_factors(forces)
289 assert order > 0
290 n = order + 1
291 T = temperature
292 freq = np.array(f)
294 # frequencies and their multiples
295 freq_n = [[] * i for i in range(n - 1)]
296 freq_neg = [[] * i for i in range(n - 1)]
298 for i in range(1, n):
299 freq_n[i - 1] = freq * i
300 freq_neg[i - 1] = freq * (-i)
302 # combinations
303 freq_nn = [x for x in combinations(chain(*freq_n), 2)]
304 for i in range(len(freq_nn)):
305 freq_nn[i] = freq_nn[i][0] + freq_nn[i][1]
307 indices2 = []
308 for i, y in enumerate(freq):
309 ind = [j for j, x in enumerate(freq_nn) if y == 0 or x % y == 0]
310 indices2.append(ind)
311 indices2 = [x for x in chain(*indices2)]
312 freq_nn = np.delete(freq_nn, indices2)
314 frequencies = [[] * x for x in range(3)]
315 frequencies[0].append(freq_neg[0])
316 frequencies[0].append([0])
317 frequencies[0].append(freq_n[0])
318 frequencies[0] = [x for x in chain(*frequencies[0])]
320 for i in range(1, n - 1):
321 frequencies[1].append(freq_neg[i])
322 frequencies[1].append(freq_n[i])
323 frequencies[1] = [x for x in chain(*frequencies[1])]
325 frequencies[2] = freq_nn
327 # Franck-Condon factors
328 E = freq / 8065.5
329 f_n = [[] * i for i in range(n)]
331 for j in range(n):
332 f_n[j] = np.exp(-E * j / (kB * T))
334 # partition function
335 Z = np.empty(len(S))
336 Z = np.sum(f_n, 0)
338 # occupation probability
339 w_n = [[] * k for k in range(n)]
340 for lval in range(n):
341 w_n[lval] = f_n[lval] / Z
343 # overlap wavefunctions
344 O_n = [[] * m for m in range(n)]
345 O_neg = [[] * m for m in range(n)]
346 for o in range(n):
347 O_n[o] = [[] * p for p in range(n)]
348 O_neg[o] = [[] * p for p in range(n - 1)]
349 for q in range(o, n + o):
350 a = np.minimum(o, q)
351 summe = []
352 for k in range(a + 1):
353 s = ((-1)**(q - k) * np.sqrt(S)**(o + q - 2 * k) *
354 factorial(o) * factorial(q) /
355 (factorial(k) * factorial(o - k) * factorial(q - k)))
356 summe.append(s)
357 summe = np.sum(summe, 0)
358 O_n[o][q - o] = (np.exp(-S / 2) /
359 (factorial(o) * factorial(q))**(0.5) *
360 summe)**2 * w_n[o]
361 for q in range(n - 1):
362 O_neg[o][q] = [0 * b for b in range(len(S))]
363 for q in range(o - 1, -1, -1):
364 a = np.minimum(o, q)
365 summe = []
366 for k in range(a + 1):
367 s = ((-1)**(q - k) * np.sqrt(S)**(o + q - 2 * k) *
368 factorial(o) * factorial(q) /
369 (factorial(k) * factorial(o - k) * factorial(q - k)))
370 summe.append(s)
371 summe = np.sum(summe, 0)
372 O_neg[o][q] = (np.exp(-S / 2) /
373 (factorial(o) * factorial(q))**(0.5) *
374 summe)**2 * w_n[o]
375 O_neg = np.delete(O_neg, 0, 0)
377 # Franck-Condon factors
378 FC_n = [[] * i for i in range(n)]
379 FC_n = np.sum(O_n, 0)
380 zero = reduce(mul, FC_n[0])
381 FC_neg = [[] * i for i in range(n - 2)]
382 FC_neg = np.sum(O_neg, 0)
383 FC_n = np.delete(FC_n, 0, 0)
385 # combination FC factors
386 FC_nn = [x for x in combinations(chain(*FC_n), 2)]
387 for i in range(len(FC_nn)):
388 FC_nn[i] = FC_nn[i][0] * FC_nn[i][1]
390 FC_nn = np.delete(FC_nn, indices2)
392 FC = [[] * x for x in range(3)]
393 FC[0].append(FC_neg[0])
394 FC[0].append([zero])
395 FC[0].append(FC_n[0])
396 FC[0] = [x for x in chain(*FC[0])]
398 for i in range(1, n - 1):
399 FC[1].append(FC_neg[i])
400 FC[1].append(FC_n[i])
401 FC[1] = [x for x in chain(*FC[1])]
403 FC[2] = FC_nn
405 return FC, frequencies