Merge pull request #34 from DanPorter/32-feature-request-custom-atom-form-factors

32 feature request custom atom form factors

Author: DanPorter

Dans_Diffraction/__init__.py

@@ -31,8 +31,8 @@
 Diamond
 2017-2025
 
-Version 3.3.4
-Last updated: 12/04/2025
+Version 3.4.0
+Last updated: 05/08/2025
 
 Version History:
 02/03/18 1.0    Version History started.
@@ -85,6 +85,7 @@
 20/11/24 3.3.2  Added alternate option for neutron scattering lengths
 06/02/25 3.3.3  Added scattering options for polarised neutron and x-ray scattering. Thanks dragonyanglong!
 12/04/25 3.3.4  Improved superstructure calculations by fixing scale parameter
+05/08/25 3.4.0  Added custom atomic form factors and dispersion corrections
 
 Acknoledgements:
     2018        Thanks to Hepesu for help with Python3 support and ideas about breaking up calculations
@@ -112,6 +113,7 @@
     Sep 2024    Thanks to thamnos for suggestion to add complex neutron scattering lengths
     Oct 2024    Thanks to Lee Richter for pointing out the error in triclinic basis definition
     Dec 2024    Thanks to dragonyanglong for pointing out the error with magnetic neutron scattering
+    May 2025    Thanks to vbhartiya for suggestions about magnetic neutron scattering
 
 -----------------------------------------------------------------------------
    Copyright 2018-2025 Diamond Light Source Ltd.
@@ -166,8 +168,8 @@
            'Structures', 'Fdmnes', 'FdmnesAnalysis']
 
 
-__version__ = '3.3.3'
-__date__ = '2025/02/06'
+__version__ = '3.4.0'
+__date__ = '2025/08/05'
 
 
 # Build

Dans_Diffraction/classes_crystal.py

@@ -51,6 +51,7 @@
 22/05/23 3.2.4  Added Symmetry.wyckoff_label(), Symmetry.spacegroup_dict
 06/05/24 3.3.0  Symmetry.from_cif now loads operations from find_spacegroup if not already loaded
 06/04/25 3.3.1  scale parameter of superlattice improved
+15/09/25 3.3.2  Atoms.type changed to always be array type
 
 @author: DGPorter
 """
@@ -1088,6 +1089,7 @@ class Atoms:
         "_atom_site_fract_y",
         "_atom_site_fract_z",
     ]
+    _type_str_fmt = '<U8'
 
     def __init__(self, u=[0], v=[0], w=[0], type=None,
                  label=None, occupancy=None, uiso=None, mxmymz=None):
@@ -1100,9 +1102,9 @@ def __init__(self, u=[0], v=[0], w=[0], type=None,
         # ---Defaults---
         # type
         if type is None:
-            self.type = np.asarray([self._default_atom] * Natoms)
+            self.type = np.asarray([self._default_atom] * Natoms, dtype=self._type_str_fmt)
         else:
-            self.type = np.asarray(type, dtype=str).reshape(-1)
+            self.type = np.asarray(type, dtype=self._type_str_fmt).reshape(-1)
         # label
         if label is None:
             self.label = self.type.copy()
@@ -1136,7 +1138,7 @@ def __call__(self, u=[0], v=[0], w=[0], type=None,
 
     def __getitem__(self, idx):
         if isinstance(idx, str):
-            idx = self.label.index(idx)
+            idx = list(self.label).index(idx)
         return self.atom(idx)
 
     def fromcif(self, cifvals):
@@ -1221,8 +1223,8 @@ def fromcif(self, cifvals):
         self.u = u
         self.v = v
         self.w = w
-        self.type = element
-        self.label = label
+        self.type = np.array(element, dtype=self._type_str_fmt)
+        self.label = np.array(label, dtype=self._type_str_fmt)
         self.occupancy = occ
         self.uiso = uiso
         self.mx = mx
@@ -1276,19 +1278,21 @@ def atom(self, idx):
             return atoms[0]
         return atoms
 
-    def changeatom(self, idx=None, u=None, v=None, w=None, type=None,
+    def changeatom(self, idx, u=None, v=None, w=None, type=None,
                    label=None, occupancy=None, uiso=None, mxmymz=None):
         """
-        Change an atoms properties
-        :param idx:
-        :param u:
-        :param v:
-        :param w:
-        :param type:
-        :param label:
-        :param occupancy:
-        :param uiso:
-        :param mxmymz:
+        Change an atom site's properties.
+        If properties are given as None, they are not changed.
+
+        :param idx: atom array index
+        :param u: atomic position u in relative coordinates along basis vector a
+        :param v: atomic position u in relative coordinates along basis vector b
+        :param w: atomic position u in relative coordinates along basis vector c
+        :param type: atomic element type
+        :param label: atomic site label
+        :param occupancy: atom site occupancy
+        :param uiso: atom site isotropic thermal parameter
+        :param mxmymz: atom site magnetic vector (mx, my, mz)
         :return: None
         """
 
@@ -1307,7 +1311,7 @@ def changeatom(self, idx=None, u=None, v=None, w=None, type=None,
         if label is not None:
             old_labels = list(self.label)
             old_labels[idx] = label
-            self.label = np.array(old_labels)
+            self.label = np.array(old_labels, dtype=self._type_str_fmt)
 
         if occupancy is not None:
             self.occupancy[idx] = occupancy
@@ -1433,11 +1437,11 @@ def remove_duplicates(self, min_distance=0.01, all_types=False):
         """
 
         uvw = self.uvw()
-        type = self.type
+        atom_type = self.type
         rem_atom_idx = []
-        for n in range(0, len(type) - 1):
+        for n in range(0, len(atom_type) - 1):
             match_position = fg.mag(uvw[n, :] - uvw[n + 1:, :]) < min_distance
-            match_type = type[n] == type[n + 1:]
+            match_type = atom_type[n] == atom_type[n + 1:]
             if all_types:
                 rem_atom_idx += list(1 + n + np.where(match_position)[0])
             else:
@@ -1588,6 +1592,10 @@ def mass_fraction(self):
 
         return weights / total_weight
 
+    def scattering_factor_coefficients(self, table='itc'):
+        """Return scattering factor coefficients for the elements"""
+        return fc.scattering_factor_coefficients(*self.type, table=table)
+
     def info(self, idx=None, type=None):
         """
         Prints properties of all atoms

Dans_Diffraction/classes_plotting.py

@@ -1457,6 +1457,45 @@ def plot_ms_azimuth(self, hkl, energy_kev, azir=[0, 0, 1], pv=[1, 0], numsteps=3
         fp.labels(ttl, r'$\psi$ (deg)', 'Intensity')
         #plt.subplots_adjust(bottom=0.2)
 
+    def plot_scattering_factors(self, q_max=4, energy_range=None, q_range=None):
+        """
+        Plot atomic scattering factors across wavevector or energy
+
+        if q_range has more values than energy_range, figures will plot scattering factor vs Q
+        for each energy.
+        if energy_range has more values than q_range, figures will plot scattering factor vs energy
+        for each value of Q.
+
+        :param q_max: use a q_range of 0-q_max
+        :param energy_range: energy range in keV
+        :param q_range: range of wavecectors, or None to use q_max
+        :return: None
+        """
+
+        if q_range is None:
+            q_range = np.arange(0, q_max, 0.01)
+        atom_types, atom_idx = np.unique(self.xtl.Structure.type, return_index=True)
+        atom_scattering_factors = self.xtl.Scatter.scattering_factors(qmag=q_range, energy_kev=energy_range)
+        ttl = f"{self.xtl.Scatter._scattering_type}"
+        # plot multiple figures of lower dimension
+        if atom_scattering_factors.shape[2] > atom_scattering_factors.shape[0]: # energy_range > q_range
+            # different figures for different values of Q
+            for n in range(atom_scattering_factors.shape[0]):
+                plt.figure(figsize=self._figure_size, dpi=self._figure_dpi)
+                nttl = ttl + f" Q={q_range[n]:.2g}" + r" $\AA^{-1}$"
+                for atom, idx in zip(atom_types, atom_idx):
+                    plt.plot(energy_range, atom_scattering_factors[n, idx, :], label=atom)
+                fp.labels(nttl, 'Energy [keV]', 'Atomic scattering factor', legend=True)
+        else:
+            # different figures for different values of E
+            for n in range(atom_scattering_factors.shape[2]):
+                plt.figure(figsize=self._figure_size, dpi=self._figure_dpi)
+                nttl = ttl + (f" E = {energy_range[n]} keV" if atom_scattering_factors.shape[2] > 1 else "")
+                for atom, idx in zip(atom_types, atom_idx):
+                    plt.plot(q_range, np.abs(atom_scattering_factors[:, idx, n]), label=atom)
+                fp.labels(nttl, r'Q $\AA^{-1}$', 'Atomic scattering factor', legend=True)
+
+
     r''' Remove tensor_scattering 26/05/20
     def tensor_scattering_azimuth(self, atom_label, hkl, energy_kev, azir=[0, 0, 1], process='E1E1',
                                   rank=2, time=+1, parity=+1, mk=None, lk=None, sk=None):

Dans_Diffraction/classes_properties.py

@@ -31,6 +31,7 @@
 
 from . import functions_general as fg
 from . import functions_crystallography as fc
+from . import functions_scattering as fs
 from .classes_orbitals import CrystalOrbitals
 
 __version__ = '2.0'
@@ -119,7 +120,30 @@ def magnetic_form_factor(self, hkl):
         :return: [nxm] array of scattering factors for each atom and reflection
         """
         qmag = self.xtl.Cell.Qmag(hkl)
-        return fc.magnetic_form_factor(self.xtl.Structure.type, qmag)
+        return fc.magnetic_form_factor(*self.xtl.Structure.type, qmag=qmag)
+
+    def scattering_factors(self, scattering_type, hkl, energy_kev=None,
+                           use_sears=False, use_wasskirf=False):
+        """
+        Return an array of scattering factors based on the radiation
+        :param scattering_type: str radiation, see "get_scattering_function()"
+        :param hkl: [mx1] or None, float array of wavevector magnitudes for reflections
+        :param energy_kev: [ox1] or None, float array of energies in keV
+        :param use_sears: if True, use neutron scattering lengths from ITC Vol. C, By V. F. Sears
+        :param use_wasskirf: if True, use x-ray scattering factors from Waasmaier and Kirfel
+        :return: [nxmxo] array of scattering factors
+        """
+        qmag = self.xtl.Cell.Qmag(hkl)
+        # Scattering factors
+        ff = fs.scattering_factors(
+            scattering_type=scattering_type,
+            atom_type=self.xtl.Structure.type,
+            qmag=qmag,
+            enval=energy_kev,
+            use_sears=use_sears,
+            use_wasskirf=use_wasskirf,
+        )
+        return np.squeeze(ff)
 
     def xray_edges(self):
         """