/
spectrum.py
267 lines (232 loc) · 10.7 KB
/
spectrum.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
"""This module defines classes to represent all xas and stitching methods."""
from __future__ import annotations
import math
import warnings
from typing import Literal
import numpy as np
from scipy.interpolate import interp1d
from pymatgen.analysis.structure_matcher import StructureMatcher
from pymatgen.core.spectrum import Spectrum
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
__author__ = "Chen Zheng, Yiming Chen"
__copyright__ = "Copyright 2012, The Materials Project"
__version__ = "3.0"
__maintainer__ = "Yiming Chen"
__email__ = "chz022@ucsd.edu, yic111@ucsd.edu"
__date__ = "July 17, 2020"
class XAS(Spectrum):
"""
Basic XAS object.
Args:
x: A sequence of x-ray energies in eV
y: A sequence of mu(E)
structure (Structure): Structure associated with the spectrum
absorbing_element (Element): Element associated with the spectrum
edge (str): Absorption edge associated with the spectrum
spectrum_type (str): 'XANES' or 'EXAFS'
absorbing_index (None or int): If None, the spectrum is assumed to be a
site-weighted spectrum, which is comparable to experimental one.
Otherwise, it indicates that the absorbing_index for a site-wise spectrum.
Attributes:
x (Sequence[float]): The sequence of energies.
y (Sequence[float]): The sequence of mu(E).
absorbing_element (str): The absorbing element of the spectrum.
edge (str): The edge of the spectrum.
spectrum_type (str): The type of the spectrum (XANES or EXAFS).
absorbing_index (int): The absorbing index of the spectrum.
"""
XLABEL = "Energy"
YLABEL = "Intensity"
def __init__(
self,
x,
y,
structure,
absorbing_element,
edge="K",
spectrum_type="XANES",
absorbing_index=None,
):
"""Initializes a spectrum object."""
super().__init__(x, y, structure, absorbing_element, edge)
self.structure = structure
self.absorbing_element = absorbing_element
self.edge = edge
self.spectrum_type = spectrum_type
self.e0 = self.x[np.argmax(np.gradient(self.y) / np.gradient(self.x))]
# Wavenumber, k is calculated based on equation
# k^2=2*m/(hbar)^2*(E-E0)
self.k = [np.sqrt((i - self.e0) / 3.8537) if i > self.e0 else -np.sqrt((self.e0 - i) / 3.8537) for i in self.x]
self.absorbing_index = absorbing_index
# check for empty spectra and negative intensities
if sum(1 for i in self.y if i <= 0) / len(self.y) > 0.05:
raise ValueError("Double check the intensities. Most of them are non-positive.")
def __str__(self):
return (
f"{self.absorbing_element} {self.edge} Edge {self.spectrum_type} "
f"for {self.structure.reduced_formula}: {super()}"
)
def stitch(self, other: XAS, num_samples: int = 500, mode: Literal["XAFS", "L23"] = "XAFS") -> XAS:
"""
Stitch XAS objects to get the full XAFS spectrum or L23 edge XANES
spectrum depending on the mode.
1. Use XAFS mode for stitching XANES and EXAFS with same absorption edge.
The stitching will be performed based on wavenumber, k.
for k <= 3, XAS(k) = XAS[XANES(k)]
for 3 < k < max(xanes_k), will interpolate according to
XAS(k)=f(k)*mu[XANES(k)]+(1-f(k))*mu[EXAFS(k)]
where f(k)=cos^2((pi/2) (k-3)/(max(xanes_k)-3)
for k > max(xanes_k), XAS(k) = XAS[EXAFS(k)]
2. Use L23 mode for stitching L2 and L3 edge XANES for elements with
atomic number <=30.
Args:
other: Another XAS object.
num_samples (int): Number of samples for interpolation.
mode("XAFS" | "L23"): Either XAFS mode for stitching XANES and EXAFS
or L23 mode for stitching L2 and L3.
Returns:
XAS object: The stitched spectrum.
"""
matcher = StructureMatcher()
if not matcher.fit(self.structure, other.structure):
raise ValueError("The input structures for spectra mismatch")
if not self.absorbing_element == other.absorbing_element:
raise ValueError("The absorbing elements for spectra are different")
if not self.absorbing_index == other.absorbing_index:
raise ValueError("The absorbing indexes for spectra are different")
if mode == "XAFS":
if not self.edge == other.edge:
raise ValueError("Only spectrum with the same absorption edge can be stitched in XAFS mode.")
if self.spectrum_type == other.spectrum_type:
raise ValueError("Need one XANES and one EXAFS spectrum to stitch in XAFS mode")
xanes = self if self.spectrum_type == "XANES" else other
exafs = self if self.spectrum_type == "EXAFS" else other
if max(xanes.x) < min(exafs.x):
raise ValueError("Energy overlap between XANES and EXAFS is needed for stitching")
# for k <= 3
wavenumber: list[float] = []
mu: list[float] = []
idx = xanes.k.index(min(self.k, key=lambda x: (abs(x - 3), x)))
mu.extend(xanes.y[:idx])
wavenumber.extend(xanes.k[:idx])
# for 3 < k < max(xanes.k)
fs: list[float] = []
ks = np.linspace(3, max(xanes.k), 50)
for k in ks:
f = np.cos((math.pi / 2) * (k - 3) / (max(xanes.k) - 3)) ** 2
fs.append(f)
f_xanes = interp1d(
np.asarray(xanes.k),
np.asarray(xanes.y),
bounds_error=False,
fill_value=0,
)
f_exafs = interp1d(
np.asarray(exafs.k),
np.asarray(exafs.y),
bounds_error=False,
fill_value=0,
)
mu_xanes = f_xanes(ks)
mu_exafs = f_exafs(ks)
mus = [fs[i] * mu_xanes[i] + (1 - fs[i]) * mu_exafs[i] for i in np.arange(len(ks))]
mu.extend(mus)
wavenumber.extend(ks)
# for k > max(xanes.k)
idx = exafs.k.index(min(exafs.k, key=lambda x: (abs(x - max(xanes.k)))))
mu.extend(exafs.y[idx:])
wavenumber.extend(exafs.k[idx:])
# interpolation
f_final = interp1d(np.asarray(wavenumber), np.asarray(mu), bounds_error=False, fill_value=0)
wavenumber_final = np.linspace(min(wavenumber), max(wavenumber), num=num_samples)
mu_final = f_final(wavenumber_final)
energy_final = [3.8537 * i**2 + xanes.e0 if i > 0 else -3.8537 * i**2 + xanes.e0 for i in wavenumber_final]
return XAS(
energy_final,
mu_final,
self.structure,
self.absorbing_element,
xanes.edge,
"XAFS",
)
if mode == "L23":
if self.spectrum_type != "XANES" or other.spectrum_type != "XANES":
raise ValueError("Only XANES spectrum can be stitched in L23 mode.")
if self.edge not in ["L2", "L3"] or other.edge not in ["L2", "L3"] or self.edge == other.edge:
raise ValueError("Need one L2 and one L3 edge spectrum to stitch in L23 mode.")
l2_xanes = self if self.edge == "L2" else other
l3_xanes = self if self.edge == "L3" else other
if l2_xanes.absorbing_element.number > 30:
raise ValueError(f"Does not support L2,3-edge XANES for {l2_xanes.absorbing_element} element")
l2_f = interp1d(
l2_xanes.x,
l2_xanes.y,
bounds_error=False,
fill_value="extrapolate",
kind="cubic",
)
l3_f = interp1d(l3_xanes.x, l3_xanes.y, bounds_error=True, fill_value=0, kind="cubic")
energy = list(np.linspace(min(l3_xanes.x), max(l3_xanes.x), num=num_samples))
mu = [i + j for i, j in zip([0 if i < 0 else i for i in l2_f(energy)], l3_f(energy))]
# check for jumps at the onset of L2-edge XANES
idx = energy.index(min(energy, key=lambda x: (abs(x - l2_xanes.x[0]))))
if abs(mu[idx] - mu[idx - 1]) / (mu[idx - 1]) > 0.1:
warnings.warn(
"There might exist a jump at the L2 and L3-edge junction.",
UserWarning,
)
return XAS(energy, mu, self.structure, self.absorbing_element, "L23", "XANES")
raise ValueError("Invalid mode. Only XAFS and L23 are supported.")
def site_weighted_spectrum(xas_list: list[XAS], num_samples: int = 500) -> XAS:
"""
Obtain site-weighted XAS object based on site multiplicity for each
absorbing index and its corresponding site-wise spectrum.
Args:
xas_list([XAS]): List of XAS object to be weighted
num_samples(int): Number of samples for interpolation
Returns:
XAS object: The site-weighted spectrum
"""
matcher = StructureMatcher()
groups = matcher.group_structures([i.structure for i in xas_list])
if len(groups) > 1:
raise ValueError("The input structures mismatch")
if not len({i.absorbing_element for i in xas_list}) == len({i.edge for i in xas_list}) == 1:
raise ValueError(
"Can only perform site-weighting for spectra with same absorbing element and same absorbing edge."
)
if len({i.absorbing_index for i in xas_list}) == 1 or None in {i.absorbing_index for i in xas_list}:
raise ValueError("Need at least two site-wise spectra to perform site-weighting")
sa = SpacegroupAnalyzer(groups[0][0])
ss = sa.get_symmetrized_structure()
maxes, mines = [], []
fs = []
multiplicities = []
for xas in xas_list:
multiplicity = len(ss.find_equivalent_sites(ss[xas.absorbing_index]))
multiplicities.append(multiplicity)
maxes.append(max(xas.x))
mines.append(min(xas.x))
# use 3rd-order spline interpolation for mu (idx 3) vs energy (idx 0).
f = interp1d(
np.asarray(xas.x),
np.asarray(xas.y),
bounds_error=False,
fill_value=0,
kind="cubic",
)
fs.append(f)
# Interpolation within the intersection of x-axis ranges.
x_axis = np.linspace(max(mines), min(maxes), num=num_samples)
weighted_spectrum = np.zeros(num_samples)
sum_multiplicities = sum(multiplicities)
for i, j in enumerate(multiplicities):
weighted_spectrum += (j * fs[i](x_axis)) / sum_multiplicities
return XAS(
x_axis,
weighted_spectrum,
ss,
xas.absorbing_element,
xas.edge,
xas.spectrum_type,
)