Coverage for larch/xafs/sigma2_models.py: 24%
139 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1#!/usr/bin/env python
2# models for debye-waller factors for xafs
4import ctypes
5import numpy as np
6from larch.larchlib import get_dll
8import scipy.constants as consts
9from scipy.special import gamma
11# EINS_FACTOR = hbarc*hbarc/(2 * k_boltz * amu) = 24.254360157751783
12# k_boltz = 8.6173324e-5 # [eV / K]
13# amu = 931.494061e6 # [eV / (c*c)]
14# hbarc = 1973.26938 # [eV * A]
15EINS_FACTOR = 1.e20*consts.hbar**2/(2*consts.k*consts.atomic_mass)
17FEFF6LIB = None
19def gnxas(r0, sigma, beta, path=None):
20 """calculate GNXAS amplitude for values of r0, sigma, beta for a feffpath
22 amplitude = gnxas(r0, sigma, beta, path)
24 Parameters:
25 -----------
26 r0 peak R0 value (in Ang)
27 sigma variance (approximately sqrt(sigma^2) Ang
28 beta asymmetry factor (unitless)
29 path feff path or None [None]
31 Notes:
32 ------
33 For reff = reff of the provided path:
34 q = 4. / beta**2
35 alpha = q + (2*(reff-r0))/(beta*sigma)
36 gr = max(0, 2*n*exp(-alpha)*alpha**(q-1)/(sigma*abs(beta)*gamma(q)))
38 """
39 feffpath = path._feffdat
40 if feffpath is None:
41 return 0.
43 reff = feffpath.reff
44 if abs(beta) < 1.e-15:
45 beta = 1.e-15
46 q = 4.0/beta**2
47 x = (reff-r0) * beta / (2*sigma)
48 alpha = q *(1+x)
49 try:
50 amp = np.exp(-alpha) * alpha**(q-1)
51 except:
52 amp = 0.0
54 out = 2*amp/(sigma*abs(beta)*gamma(q))
55 if isinstance(out, np.ndarray):
56 out[np.where(out<0)] = 0
57 out[np.where(np.isnan(out))] = 0
58 out[np.where(np.isinf(out))] = 0
59 else:
60 out = max(0, out)
61 return out
65def sigma2_eins(t, theta, path):
66 """calculate sigma2 for a Feff Path wih the einstein model
68 sigma2 = sigma2_eins(t, theta, path)
70 Parameters:
71 -----------
72 t sample temperature (in K)
73 theta Einstein temperature (in K)
74 path FeffPath to calculate sigma2 for
76 Notes:
77 sigma2 = FACTOR*coth(2*t/theta)/(theta * mass_red)
79 mass_red = reduced mass of Path (in amu)
80 FACTOR = hbarc*hbarc/(2*k_boltz*amu) ~= 24.25 Ang^2 * K * amu
81 """
82 feffpath = path._feffdat
83 if feffpath is None:
84 return 0.
85 theta = max(float(theta), 1.e-5)
86 t = max(float(t), 1.e-5)
87 rmass = 0.
88 for sym, iz, ipot, amass, x, y, z in feffpath.geom:
89 rmass = rmass + 1.0/max(0.1, amass)
90 rmass = 1.0/max(1.e-12, rmass)
91 return EINS_FACTOR/(theta * rmass * np.tanh(theta/(2.0*t)))
93def sigma2_debye(t, theta, path):
94 """calculate sigma2 for a Feff Path wih the correlated Debye model
96 sigma2 = sigma2_debye(t, theta, path)
98 Parameters:
99 -----------
100 t sample temperature (in K)
101 theta Debye temperature (in K)
102 path FeffPath to calculate sigma2 for
103 """
104 feffpath = path._feffdat
105 if feffpath is None:
106 return 0.
107 thetad = max(float(theta), 1.e-5)
108 tempk = max(float(t), 1.e-5)
109 natoms = len(feffpath.geom)
110 rnorm = feffpath.rnorman
111 atomx, atomy, atomz, atomm = [], [], [], []
112 for sym, iz, ipot, am, x, y, z in feffpath.geom:
113 atomx.append(x)
114 atomy.append(y)
115 atomz.append(z)
116 atomm.append(am)
118 return sigma2_correldebye(natoms, tempk, thetad, rnorm,
119 atomx, atomy, atomz, atomm)
121def sigma2_correldebye(natoms, tk, theta, rnorm, x, y, z, atwt):
122 """
123 internal sigma2 calc for a Feff Path wih the correlated Debye model
125 these routines come courtesy of jj rehr and si zabinsky.
127 Arguments:
128 natoms *int, lengths for x, y, z, atwt [in]
129 tk *double, sample temperature (K) [in]
130 theta *double, Debye temperature (K) [in]
131 rnorm *double, Norman radius (Ang) [in]
132 x *double, array of x coord (Ang) [in]
133 y *double, array of y coord (Ang) [in]
134 x *double, array of z coord (Ang) [in]
135 atwt *double, array of atomic_weight (amu) [in]
137 Returns:
138 sig2_cordby double, calculated sigma2
139 """
140 global FEFF6LIB
141 if FEFF6LIB is None:
142 FEFF6LIB = get_dll('feff6')
143 FEFF6LIB.sigma2_debye.restype = ctypes.c_double
145 na = ctypes.pointer(ctypes.c_int(natoms))
146 t = ctypes.pointer(ctypes.c_double(tk))
147 th = ctypes.pointer(ctypes.c_double(theta))
148 rs = ctypes.pointer(ctypes.c_double(rnorm))
150 ax = (natoms*ctypes.c_double)()
151 ay = (natoms*ctypes.c_double)()
152 az = (natoms*ctypes.c_double)()
153 am = (natoms*ctypes.c_double)()
155 for i in range(natoms):
156 ax[i], ay[i], az[i], am[i] = x[i], y[i], z[i], atwt[i]
158 return FEFF6LIB.sigma2_debye(na, t, th, rs, ax, ay, az, am)
161def sigma2_correldebye_py(natoms, tk, theta, rnorm, x, y, z, atwt):
162 """calculate the XAFS debye-waller factor for a path based
163 on the temperature, debye temperature, average norman radius,
164 atoms in the path, and their positions.
166 these routines come courtesy of jj rehr and si zabinsky.
168 Arguments:
169 natoms *int, lengths for x, y, z, atwt [in]
170 tk *double, sample temperature (K) [in]
171 theta *double, Debye temperature (K) [in]
172 rnorm *double, Norman radius (Ang) [in]
173 x *double, array of x coord (Ang) [in]
174 y *double, array of y coord (Ang) [in]
175 x *double, array of z coord (Ang) [in]
176 atwt *double, array of atomic_weight (amu) [in]
178 Returns:
179 sig2_cordby double, calculated sigma2
181 Notes:
182 1. natoms must be >= 2.
183 2. rnorman is the wigner-seitz or norman radius,
184 averaged over entire problem:
185 (4pi/3)*rs**3 = sum( (4pi/3)rnrm**3 ) / n
186 (sum is over all atoms in the problem)
187 3. all distances are in Angstroms
189 moved from Feff6 sigms.f, original copyright:
190 copyright 1993 university of washington
191 john rehr, steve zabinsky, matt newville
192 """
193 sig2 = 0.0
194 for i0 in range(natoms):
195 i1 = (i0 + 1) % natoms
196 for j0 in range(i0, natoms):
197 j1 = (j0 + 1) % natoms
198 # calculate r_i-r_i-1 and r_j-r_j-1 and the rest of the
199 # distances, and get the partial cosine term:
200 # cosine(i,j) = r_i.r_j / ((r_i0- r_i-1) * (r_j - r_j-1))
201 ri0j0 = dist(x[i0], y[i0], z[i0], x[j0], y[j0], z[j0])
202 ri1j1 = dist(x[i1], y[i1], z[i1], x[j1], y[j1], z[j1])
203 ri0j1 = dist(x[i0], y[i0], z[i0], x[j1], y[j1], z[j1])
204 ri1j0 = dist(x[i1], y[i1], z[i1], x[j0], y[j0], z[j0])
205 ri0i1 = dist(x[i0], y[i0], z[i0], x[i1], y[i1], z[i1])
206 rj0j1 = dist(x[j0], y[j0], z[j0], x[j1], y[j1], z[j1])
207 ridotj = ( (x[i0] - x[i1]) * (x[j0] - x[j1]) +
208 (y[i0] - y[i1]) * (y[j0] - y[j1]) +
209 (z[i0] - z[i1]) * (z[j0] - z[j1]) )
211 # call corrfn to get the correlations between atom pairs
212 ci0j0 = corrfn(ri0j0, theta, tk, atwt[i0], atwt[j0], rnorm)
213 ci1j1 = corrfn(ri1j1, theta, tk, atwt[i1], atwt[j1], rnorm)
214 ci0j1 = corrfn(ri0j1, theta, tk, atwt[i0], atwt[j1], rnorm)
215 ci1j0 = corrfn(ri1j0, theta, tk, atwt[i1], atwt[j0], rnorm)
217 # combine outputs of corrfn to give the debye-waller factor for
218 # this atom pair. !! note: don't double count (i.eq.j) terms !!!
219 sig2ij = ridotj*(ci0j0 + ci1j1 - ci0j1 - ci1j0)/(ri0i1*rj0j1)
220 if j0 == i0:
221 sig2ij /= 2.0
222 sig2 += sig2ij
224 return sig2/2.0
228def dist(x0, y0, z0, x1, y1, z1):
229 """find distance between cartesian points
230 (x, y, z)0 and (x, y, z)1
231 port of Fortran from feff6 sigms.f
232 """
233 return np.sqrt( (x0-x1)**2 + (y0-y1)**2 + (z0-z1)**2 )
235def corrfn(rij, theta, tk, am1, am2, rs):
236 """calculate correlation function
237 c(ri, rj) = <xi xj> in the debye approximation
239 ported from feff6 sigms.f
241 copyright 1993 university of washington
242 john rehr, steve zabinsky, matt newville
244 subroutine calculates correlation function
245 c(ri, rj) = <xi xj> in the debye approximation
246 = (1/n)sum_k exp(ik.(ri-rj)) (1/sqrt(mi*mj))*
247 (hbar/2w_k)*coth(beta hbar w_k/2)
249 = (3kt/mu w_d**2) * sqrt(mu**2/mi*mj) * int
250 where :
251 x k_d*r (distance parameter) r distance in angstroms
252 theta debye temp in degrees k
253 tk temperature in degrees k
254 temper theta / tk = hbar omegad/kt
255 k_d debye wave number = (6*pi**2 n/v)
256 n/v free electron number density = 1/(4pi/3rs**3)
257 rs wigner seitz or norman radius in bohr
258 ami atomic mass at sites i in amu
259 amj atomic mass at sites j in amu
260 int int_0^1 (temper/x) dw sin(wx)coth(w*temper/2)
262 solution by numerical integration, with parameters pi, bohr, con:
263 con=hbar**2/kb*amu)*10**20 in ang**2 units
264 k_boltz = 8.6173324e-5 # [eV / K]
265 amu = 931.494061e6 # [eV / (c*c)]
266 hbarc = 1973.26938 # [eV * A]
267 bohr = 0.52917721 # [A]
268 conh = (3/2.)* hbar**2 / (kb*amu) ~= 72.76
269 conr = (9*pi/2)**(1/3.0) / bohr ~= 4.57
271 NOTE: for backward compatibility, the constants used by feff6 are
272 retained, even though some have been refined later.
273 """
274 conh = 72.7630804732553
275 conr = 4.5693349700844
277 # theta in degrees k, t temperature in degrees k
278 rx = conr * rij / rs
279 tx = theta / tk
280 rmass = theta * np.sqrt(am1 * am2)
281 return conh * debint(rx, tx) / rmass
283def debfun(w, rx, tx):
284 """ debye function, ported from feff6 sigms.f
286 copyright 1993 university of washington
287 john rehr, steve zabinsky, matt newville
289 debfun = (sin(w*rx)/rx) * coth(w*tx/2)
290 """
291 # print(" debfun ", w, rx, tx)
292 wmin = 1.e-20
293 argmax = 50.0
294 result = 2.0 / tx
295 # allow t = 0 without bombing
296 if w > wmin:
297 result = w
298 if rx > 0:
299 result = np.sin(w*rx) / rx
300 emwt = np.exp( -min(w*tx, argmax))
301 result *= (1 + emwt) / (1 - emwt)
302 return result
304def debint(rx, tx):
305 """ calculates integrals between [0,1] b = int_0^1 f(z) dz
306 by trapezoidal rule and binary refinement (romberg integration)
307 ported from feff6 sigms.f:
309 copyright 1993 university of washington
310 john rehr, steve zabinsky, matt newville
312 subroutine calculates integrals between [0,1] b = int_0^1 f(z) dz
313 by trapezoidal rule and binary refinement (romberg integration)
314 coded by j rehr (10 feb 92) see, e.g., numerical recipes
315 for discussion and a much fancier version
316 """
317 MAXITER = 12
318 tol = 1.e-9
319 itn = 1
320 step = 1.0
321 result = 0.0
322 bo = bn = (debfun(0.0, rx, tx) + debfun(1.0, rx, tx))/2.0
323 for iter in range(MAXITER):
324 # nth iteration
325 # b_n+1=(b_n)/2+deln*sum_0^2**n f([2n-1]deln)
326 step = step / 2.
327 sum = 0
328 for i in range(itn):
329 sum += debfun(step*(2*i + 1), rx, tx)
330 itn = 2*itn
331 # bnp1=b_n+1 is current value of integral
332 # cancel leading error terms b=[4b-bn]/3
333 # note: this is the first term in the neville table - remaining
334 # errors were found too small to justify the added code
335 bnp1 = step * sum + (bn / 2.0)
336 result = (4 * bnp1 - bn) / 3.0
337 if (abs( (result - bo) / result) < tol):
338 break
339 bn = bnp1
340 bo = result
341 return result
344####################################################
345## sigma2_eins and sigma2_debye are defined here to
346## be injected as Procedures within lmfit's asteval
347## for calculating XAFS sigma2 for a scattering path
348## these use `reff` or `feffpath.geom` which will be updated
349## for each path during an XAFS path calculation
350##
351_sigma2_funcs = """
352def sigma2_eins(t, theta):
353 if feffpath is None:
354 return 0.
355 theta = max(float(theta), 1.e-5)
356 t = max(float(t), 1.e-5)
357 rmass = 0.
358 for sym, iz, ipot, amass, x, y, z in feffpath.geom:
359 rmass = rmass + 1.0/max(0.1, amass)
360 rmass = 1.0/max(1.e-12, rmass)
361 return EINS_FACTOR/(theta * rmass * tanh(theta/(2.0*t)))
363def sigma2_debye(t, theta):
364 if feffpath is None:
365 return 0.
366 thetad = max(float(theta), 1.e-5)
367 tempk = max(float(t), 1.e-5)
368 natoms = len(feffpath.geom)
369 rnorm = feffpath.rnorman
370 atomx, atomy, atomz, atomm = [], [], [], []
371 for sym, iz, ipot, am, x, y, z in feffpath.geom:
372 atomx.append(x)
373 atomy.append(y)
374 atomz.append(z)
375 atomm.append(am)
377 return sigma2_correldebye(natoms, tempk, thetad, rnorm,
378 atomx, atomy, atomz, atomm)
381def gnxas(r0, sigma, beta):
382 if feffpath is None:
383 return 0.
384 r = feffpath.reff
385 if abs(beta) < 1.e-15:
386 beta = 1.e-15
387 q = 4. / beta**2
388 alpha = q + (2*(r-r0)) / (beta*sigma)
389 try:
390 a2q = alpha**(q-1)
391 except:
392 a2q = 1.e100
394 print('> ', reff, alpha, q)
396 return max(0, 2*exp(-alpha)*a2q/(sigma*abs(beta)*gamma(q)))
398"""
399def add_sigma2funcs(params):
400 """set sigma2funcs into Parameters' asteval"""
401 f_eval = params._asteval
402 f_eval.symtable['EINS_FACTOR'] = EINS_FACTOR
403 f_eval.symtable['sigma2_correldebye'] = sigma2_correldebye
404 f_eval.symtable['feffpath'] = None
405 f_eval.symtable['gamma'] = gamma
406 f_eval(_sigma2_funcs)