feat: added a gaussian t_r function with x-ray weighing

Author: CPrescher

glassure/fitting.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from .utility import calculate_weighting_factor
+from .utility import calculate_weighting_factor, calculate_weighting_factor_wkm
 from .transform import calculate_fr
 from .pattern import Pattern
 
@@ -24,10 +24,28 @@ def i_q_peak(q, n, position, sigma, composition, element_1, element_2):
     num_atoms = sum([val for _, val in composition.items()])
     c_2 = composition[element_2] / num_atoms
     w = calculate_weighting_factor(composition, element_1, element_2, q)
-    return n * w / c_2 * np.sin(q * position) / (q * position) * np.exp(-q ** 2 * sigma ** 2 / 2)
+    return (
+        n
+        * w
+        / c_2
+        * np.sin(q * position)
+        / (q * position)
+        * np.exp(-(q**2) * sigma**2 / 2)
+    )
 
 
-def t_r_peak(r, n, position, sigma, composition, element_1, element_2, q, use_modification_fcn=False, method='fft'):
+def t_r_peak(
+    r,
+    n,
+    position,
+    sigma,
+    composition,
+    element_1,
+    element_2,
+    q,
+    use_modification_fcn=False,
+    method="fft",
+):
     """
     Calculates the contribution of one element 1 - element 2 peak in real space to t(r). We assume a gaussian
     broadening. The math is explained in the paper about NXFit (Pickup et al. 2014, J. Appl. Cryst. 47, 1790-1796).
@@ -54,4 +72,35 @@ def t_r_peak(r, n, position, sigma, composition, element_1, element_2, q, use_mo
                     - 'fft' solves the Fourier integral by using fast fourier transformation
     """
     q_peak = i_q_peak(q, n, position, sigma, composition, element_1, element_2)
-    return calculate_fr(Pattern(q, q_peak + 1), r, use_modification_fcn=use_modification_fcn, method=method).y
+    return calculate_fr(
+        Pattern(q, q_peak + 1),
+        r,
+        use_modification_fcn=use_modification_fcn,
+        method=method,
+    ).y
+
+
+def t_r_peak_gaussian(r, n, position, sigma, composition, element_1, element_2):
+    """
+    Calculates the contribution of one element 1 - element 2 peak in real space to t(r). We assume a gaussian
+    distribution. The math is well explained in the LiquidDiffract paper (Heinen and Drewitt 2022. Physics and
+    Chemistry of Minerals 49, no. 5: 9. https://doi.org/10.1007/s00269-022-01186-6).
+
+    The gaussian is weighted based on the X-ray form factors of the two elements.
+
+    :param r: numpy array giving the r-values for which the peak will be calculated
+    :param n: coordination number of element 2 to element 1
+    :param position: average distance between the two elements
+    :param sigma: measure for broadness of distances distribution
+    :param composition: composition: dictionary with elements as key and abundances as relative numbers
+    :param element_1: string giving element 1
+    :param element_2: string giving element 1
+
+    :return: T(r) pattern
+    """
+
+    wf = 1 / calculate_weighting_factor_wkm(composition, element_1, element_2, 0)
+
+    part1 = n * wf / (composition[element_2] * sigma * r * np.sqrt(2 * np.pi))
+    part2 = np.exp(-((r - position) ** 2) / (2 * sigma**2))
+    return part1 * part2

glassure/utility.py

@@ -9,6 +9,7 @@
 import scipy.constants as const
 
 import lmfit
+import xraylib as xr
 
 from .scattering_factors import (
     calculate_coherent_scattering_factor,
@@ -252,6 +253,129 @@ def calculate_weighting_factor(
     )
 
 
+def calculate_weighting_factor_wkm(
+    composition: Composition,
+    element_1: str,
+    element_2: str,
+    q: np.ndarray | float = 0.0,
+    sf_source="hajdu",
+) -> float:
+    """
+    Calculates the weighting factor for an element-element contribution in a given composition (e.g. for Si-O in SiO2)
+    using the Warren-Krutter-Morningstar approximation for the form factors or effective atomic numbers (Warren et al.
+    1936).
+
+    This means we calculate:
+    .. math::
+        W_{\alpha \beta} = \frac{c_\alpha c_\beta K_\alpha(q) K_\beta(q) (2-\delta_{\alpha \beta})}{\sum_{\alpha} c_\alpha K_\alpha(q)^2}
+
+    :param composition: dictionary with elements as key and abundances as relative numbers
+    :param element_1: string giving element 1
+    :param element_2: string giving element 2
+    :param q: Q value or numpy array with a unit of A^-1. Default is 0.0.
+    :param sf_source: source of the scattering factors. Possible sources are 'hajdu' and 'brown_hubbell'.
+    :return: weighting factor
+    """
+    if element_1 == element_2:
+        factor = 1
+    else:
+        factor = 2
+
+    effective_atomic_numbers = {
+        element: calculate_wkm_effective_atomic_number(
+            composition, element, q, sf_source
+        )
+        for element in composition.keys()
+    }
+
+    denominator = (
+        np.sum(
+            [
+                effective_atomic_numbers[element] * c
+                for element, c in composition.items()
+            ]
+        )
+        ** 2
+    )
+    numerator = (
+        composition[element_1]
+        * composition[element_2]
+        * effective_atomic_numbers[element_1]
+        * effective_atomic_numbers[element_2]
+        * factor
+    )
+
+    return numerator / denominator
+
+
+def calculate_wkm_effective_atomic_number(
+    composition: Composition,
+    element: str,
+    q: np.ndarray,
+    sf_source="hajdu",
+) -> float:
+    """
+    Calculates the effective atomic number for a given element using the Warren-Krutter-Morningstar approximation
+    (Warren et al. 1936).
+
+    This means we calculate:
+    .. math::
+       K_\alpha(q) = \left\langle\frac{f_\alpha( q)}{f_e(q)}\right\rangle
+    where :math:`f_\alpha(q)` is the form factor of the element :math:`\alpha` and :math:`f_e(q)` is the effective
+    form factor of the composition.
+
+    .. math::
+       f_e(q) = \frac{\sum_{\alpha} f_\alpha(q)}{Z_{tot}}
+    where :math:`Z_{tot}` is the total atomic number of the compositional unit
+
+    :param composition: dictionary with elements as key and abundances as relative numbers
+    :param element: string giving element for which the effective atomic number will be calculated
+    :param q: Q value or numpy array with a unit of A^-1
+    :param sf_source: source of the scattering factors. Possible sources are 'hajdu', 'brown_hubbell' and 'xraylib'.
+    """
+
+    effective_form_factor = calculate_effective_form_factor(composition, q, sf_source)
+    form_factor = calculate_coherent_scattering_factor(element, q, sf_source)
+    return np.mean(form_factor / effective_form_factor)
+
+
+def calculate_total_atomic_number(composition: Composition) -> float:
+    """
+    Calculates the total atomic number of a given composition.
+    """
+    return sum(
+        [xr.SymbolToAtomicNumber(element) * c for element, c in composition.items()]
+    )
+
+
+def calculate_effective_form_factor(
+    composition: Composition, q: np.ndarray, sf_source="hajdu"
+) -> np.ndarray:
+    """
+    Calculates the effective form factor for a given composition, which is given by
+
+    .. math::
+        f_e(q) = \frac{\sum_{\alpha} f_\alpha(q)}{Z_{tot}}
+    where :math:`f_\alpha(q)` is the form factor of the element :math:`\alpha` and :math:`Z_{tot}` is the total atomic
+    number of the compositional unit
+
+    :param composition: dictionary with elements as key and abundances as relative numbers
+    :param q: Q value or numpy array with a unit of A^-1
+    :param sf_source: source of the scattering factors. Possible sources are 'hajdu', 'brown_hubbell' and 'xraylib'.
+
+    :return: effective form factor array
+    """
+    total_atomic_number = calculate_total_atomic_number(composition)
+
+    form_factors = [
+        calculate_coherent_scattering_factor(element, q, sf_source) * c
+        for element, c in composition.items()
+    ]
+    form_factor_sum = np.sum(form_factors, axis=0)
+
+    return form_factor_sum / total_atomic_number
+
+
 def normalize_composition(composition: Composition) -> dict[str, float]:
     """
     normalizes elemental abundances to 1

tests/test_fitting.py

@@ -3,7 +3,7 @@
 import unittest
 import numpy as np
 
-from glassure.fitting import i_q_peak, t_r_peak
+from glassure.fitting import i_q_peak, t_r_peak, t_r_peak_gaussian
 
 
 class FittingTest(unittest.TestCase):
@@ -22,3 +22,12 @@ def test_t_r_peak(self):
         q = np.linspace(0.01, 10, 1001)
         peak = t_r_peak(self.r, 4, 1.6, 0.15, self.composition, 'Si', 'O', q)
         self.assertEqual(len(peak), len(self.r))
+
+    def test_t_r_peak_gaussian(self):
+        r = np.linspace(0.1, 10, 1001)
+        peak = t_r_peak_gaussian(r, 4, 1.6, 0.15, self.composition, 'Si', 'O')
+        self.assertEqual(len(peak), len(r))
+
+        import matplotlib.pyplot as plt
+        plt.plot(r, peak)
+        plt.show()

tests/test_utility.py

@@ -18,6 +18,10 @@
     calculate_s0,
     calculate_kn_correction,
     parse_str_to_composition,
+    calculate_total_atomic_number,
+    calculate_effective_form_factor,
+    calculate_wkm_effective_atomic_number,
+    calculate_weighting_factor_wkm,
 )
 from glassure import Pattern
 
@@ -39,14 +43,12 @@ def test_convert_density_to_atoms_per_cubic_angstrom(self):
 
         self.assertAlmostEqual(density_au, 0.0662, places=4)
 
-
     def test_convert_density_to_grams_per_cubic_centimeter(self):
         composition = {"Si": 1, "O": 2}
         density_au = convert_density_to_atoms_per_cubic_angstrom(composition, 2.2)
         density = convert_density_to_grams_per_cubic_centimeter(composition, density_au)
         self.assertAlmostEqual(density, 2.2, places=4)
 
-
     def test_calculate_f_mean_squared(self):
         q = np.linspace(0, 10)
         composition = {"Si": 1, "O": 2}
@@ -238,15 +240,14 @@ def test_convert_two_theta_to_q_space(self):
             np.max(pattern_q.x), 4 * np.pi * np.sin(25.0 / 360 * np.pi) / wavelength
         )
 
-
     def test_parse_string_to_composition(self):
         inputs = [
             "Si O 0.5",
             "(Mg)(Si O2)2.5",
             "[Al2O3]3",
             "\\{Fe3O4\\}2.5",
             "(H2O)2[NaCl]\\{KOH\\}0.5",
-            "\\{[Co(NH3)4(OH)2]3[Co(CN)6]\\}2"
+            "\\{[Co(NH3)4(OH)2]3[Co(CN)6]\\}2",
         ]
 
         outputs = [
@@ -255,9 +256,54 @@ def test_parse_string_to_composition(self):
             {"Al": 6, "O": 9},
             {"Fe": 7.5, "O": 10},
             {"H": 4.5, "O": 2.5, "Na": 1, "Cl": 1, "K": 0.5},
-            {"Co": 8, "N": 36, "H": 84, "O": 12, "C": 12}
+            {"Co": 8, "N": 36, "H": 84, "O": 12, "C": 12},
         ]
 
         for input, output in zip(inputs, outputs):
             parsed_formula = parse_str_to_composition(input)
             self.assertEqual(parsed_formula, output)
+
+    def test_calculate_total_atomic_number(self):
+        def test_composition(composition, total_atomic_number):
+            total_atomic_number = calculate_total_atomic_number(composition)
+            self.assertEqual(total_atomic_number, total_atomic_number)
+
+        test_composition({"Si": 1, "O": 2}, 30)
+        test_composition({"Ge": 1, "O": 2}, 48)
+        test_composition({"H": 2, "O": 1}, 10)
+        test_composition({"C": 1}, 6)
+
+    def test_calculate_effective_form_factor(self):
+        q = np.linspace(0, 10, 234)
+        composition = {"Si": 1, "O": 2}
+        effective_form_factor = calculate_effective_form_factor(composition, q)
+        self.assertEqual(len(q), len(effective_form_factor))
+
+        effective_form_factor = calculate_effective_form_factor(
+            composition, q, sf_source="brown_hubbell"
+        )
+        self.assertEqual(len(q), len(effective_form_factor))
+
+        effective_form_factor = calculate_effective_form_factor(
+            composition, q, sf_source="xraylib"
+        )
+        self.assertEqual(len(q), len(effective_form_factor))
+
+        effective_form_factor = calculate_effective_form_factor(
+            composition, 0, sf_source="hajdu"
+        )
+        self.assertAlmostEqual(effective_form_factor, 1.0, places=3)
+
+    def test_calculate_wkm_form_factor(self):
+        q = 0.0
+        composition = {"Si": 1, "O": 2}
+        wkm_form_factor = calculate_wkm_effective_atomic_number(composition, "Si", q)
+        self.assertAlmostEqual(wkm_form_factor, 14.005, places=2)
+        wkm_form_factor = calculate_wkm_effective_atomic_number(composition, "O", q)
+        self.assertAlmostEqual(wkm_form_factor, 7.9975, places=2)
+
+    def test_calculate_wkm_weighting_factor(self):
+        q = 0.0
+        composition = {"Si": 1, "O": 2}
+        wkm_weighting_factor = calculate_weighting_factor_wkm(composition, "Si", "O", q)
+        print(wkm_weighting_factor)