Coverage for dynasor/post_processing/x_ray_form_factors.py: 98%
63 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-10 08:27 +0000
1import json
2from importlib.resources import files
3from typing import Optional
4from warnings import warn
6from numpy import exp, abs, pi
7from pandas import DataFrame, concat
8from .weights import Weights
11class XRayFormFactors(Weights):
12 r"""This class generates sample weights corresponding to X-ray form factors,
13 specifically the non-dispersion corrected parametrized form factors.
14 In general, the form factor may be written as
16 .. math::
17 f(q, \omega) = f_0(q) + f'(q, \omega) + if''(q, \omega)
19 where :math:`q` is a scalar, corresponding to the norm
20 of the desired reciprocal lattice point.
22 The weights generated by this class corresponds to :math:`f_0(q)`.
23 There are two possible parametrizations of :math:`f_0(q)`
24 to choose from, both based on a sum of exponentials of the form
26 .. math::
28 f_0(s) = \sum_{i=1}^k a_i \exp(-b_i s^2) + c,
30 where :math:`s = \sin\theta / \lambda = q / 4\pi`, since dynasor defines
31 :math:`q = 4\pi\sin\theta / \lambda`.
33 Two parametrizations are available:
35 * ``'waasmaier-1995'`` corresponds to Table 1 from D. Waasmaier, A. Kirfel,
36 Acta Crystallographica Section A **51**, 416 (1995);
37 `doi: 10.1107/S0108767394013292 <https://doi.org/10.1107/S0108767394013292>`_.
38 This parametrization uses five exponentials (:math:`k=5`) and extends up to
39 :math:`s=6.0\,\mathrm{Å}^{-1}`, i.e., :math:`q=75.4\,\mathrm{rad/Å}`.
40 * ``'itc-2006'`` corresponds to Table 6.1.1.4. from
41 *International Tables for Crystallography, Volume C: Mathematical, physical
42 and chemical tables* (2006);
43 `doi: 10.1107/97809553602060000103 <https://doi.org/10.1107/97809553602060000103>`_.
44 This parametrization uses four exponentials (:math:`k=4`) and extends up to
45 :math:`s=2.0\,\mathrm{Å}^{-1}`, i.e., :math:`q=25.1\,\mathrm{rad/Å}`.
47 In practice differences are expected to be insignificant. It is unlikely that
48 you have to deviate from the default, which is ``'waasmaier-1995'``.
50 Parameters
51 ----------
52 atom_types
53 List of atomic species for which to retrieve scattering lengths.
54 source
55 Source to use for parametrization of the form factors :math:`f_0(q)`.
56 Allowed values are ``'waasmaier-1995'`` and ``'itc-2006'``
57 (see above).
58 """
60 def __init__(
61 self,
62 atom_types: list[str],
63 source: Optional[str] = 'waasmaier-1995'
64 ):
65 self._source = source
66 form_factors = self._read_form_factors(source)
67 # Select the relevant species
68 form_factors = form_factors[form_factors.index.isin(atom_types)] # Atom species is index
69 self._form_factors = form_factors
71 # Check if any of the fetched form factors are missing,
72 # indicating that it is missing in the experimental database.
73 for s in atom_types:
74 row = form_factors[form_factors.index == s]
75 if row.empty:
76 if s == 'H':
77 # Manually insert values for H such that the form factor is
78 # zero and raise a warning, since it is such a common element.
79 warn('No parametrization for H. Setting form factor for H to zero.')
80 values = [['DUMMY', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
81 form_factors = concat([DataFrame(data=values,
82 index=['H'],
83 columns=form_factors.columns),
84 form_factors])
85 continue
86 raise ValueError('Missing tabulated values '
87 f'for requested species {s}.')
89 weights_coh = form_factors.to_dict(orient='index')
90 supports_currents = False
91 super().__init__(weights_coh, None, supports_currents=supports_currents)
93 def get_weight_coh(self, atom_type: str, q_norm: float) -> float:
94 """Get the coherent weight for a given atom type and q-vector norm."""
95 return self._compute_f0(self._weights_coh[atom_type], q_norm, self._source)
97 @property
98 def parameters(self) -> DataFrame:
99 """Parametrization used to compute the coherent form factors
100 f(q) for the selected species.
101 """
102 return self._form_factors
104 def _compute_f0(
105 self,
106 coefficients: dict,
107 q_norm: float,
108 source: str):
109 r"""Compute :math:`f_0(q)` based on the chosen parametrization.
110 There are two possible parametrizations of :math:`f_0(q)` to choose from,
111 both based on a sum of exponentials of the form
113 .. math::
115 f_0(q) = \sum_{i=1}^k a_i * exp(-b_i * s**2) + c
117 where :math:`s = sin(theta) / lambda`.
118 Parameters
119 ----------
120 coefficients
121 Parametrization parameters, read from the corresponding source file.
122 q_norm
123 The |q|-value at which to evaluate the form factor.
124 source
125 Either 'waasmaier-1995' or 'itc-2006',
126 corresponding to the two available sources for the
127 parametrization of the :math:`f_0(q)` term of the form factors.
128 """
129 s = q_norm / (4 * pi) # q in dynasor is q = 4 pi sin(theta) / lambda.
130 s_squared = s*s # s = sin(theta) / lambda
131 if source == 'waasmaier-1995':
132 if abs(s) > 6.0:
133 warn('Waasmaier parametrization is not reliable for q '
134 'above 75.398 rad/Å (corresponding to s=6.0 1/Å)')
135 return self._get_f0(coefficients, s_squared, nmax=5)
136 elif source == 'itc-2006': 136 ↛ 142line 136 didn't jump to line 142 because the condition on line 136 was always true
137 if abs(s) > 2.0:
138 warn('ITC.C parametrization is not reliable for q '
139 'above 25.132 rad/Å (corresponding to s=2.0 1/Å)')
140 return self._get_f0(coefficients, s_squared, nmax=4)
141 else:
142 raise ValueError(f'Unknown source {source}')
144 def _get_f0(self, coefficients: dict, s_squared: float, nmax: Optional[int] = 5):
145 r"""
146 Compute parametrizations of :math:`f_0(q)` on the form
148 .. math::
150 f_0(q) = \sum_{i=1}^k a_i * exp(-b_i * s**2) + c
152 where :math:`s = sin(theta) / lambda`.
153 """
154 f0 = coefficients['c']
155 for i in range(1, nmax+1):
156 f0 += coefficients[f'a{i}'] * exp(-coefficients[f'b{i}'] * s_squared)
157 return f0
159 def _read_form_factors(self, source: str) -> DataFrame:
160 r"""
161 Extracts the parametrization for the form factors :math:`f_0(q)`,
162 based on either of two sources.
164 Parameters
165 ----------
166 source
167 Either 'waasmaier-1995' or 'itc-2006',
168 corresponding to the two available sources for the
169 parametrization of the :math:`f_0(q)` term of the form factors.
170 """
171 if source == 'waasmaier-1995':
172 data_file = files(__package__) / \
173 'form-factors/x-ray-parameters-waasmaier-kirfel-1995.json'
174 with open(data_file) as fp:
175 coefficients = json.load(fp)
177 form_factors = DataFrame.from_dict(coefficients)
178 form_factors.index.names = ['species']
179 elif source == 'itc-2006':
180 data_file = files(__package__) / \
181 'form-factors/x-ray-parameters-itcc-2006.json'
182 with open(data_file) as fp:
183 coefficients = json.load(fp)
185 form_factors = DataFrame.from_dict(coefficients)
186 # There are two entries for H, one calculated using
187 # Hartree-Fock and one calculated with SDS.
188 # Both parametrizations are similar. We'll
189 # use HF since that has been used for the majority
190 # of species.
191 form_factors.drop(index='0', inplace=True) # H SDS is the first row
192 form_factors = form_factors.set_index('species')
193 else:
194 raise ValueError(f'Unknown source {source}')
196 return form_factors