Merge pull request #613 from scipp/atomic-weights

Add atomic weights and masses

Author: jl-wynen

src/scippneutron/atoms/__init__.py

@@ -1,17 +1,98 @@
 # SPDX-License-Identifier: BSD-3-Clause
-# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
-"""Parameters for neutron interactions with atoms."""
+# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
+"""Parameters of atoms and neutron interactions with atoms."""
 
 from __future__ import annotations
 
 import dataclasses
 import importlib.resources
+import re
 from functools import lru_cache
 from typing import TextIO
 
 import scipp as sc
 
 
+@dataclasses.dataclass(frozen=True, eq=False)
+class Atom:
+    """Atomic parameters of a specific element / isotope.
+
+    Values have been retrieved at 2025-03-20T15:11 from the lists at
+    https://www.ciaaw.org/atomic-weights.htm
+
+    Atomic weights are properties of an *element* while atomic masses
+    are properties of a specific *isotope* (nuclide).
+    The reported atomic weights are the abridged standard
+    atomic weight for each element.
+    The reported atomic masses are the most recent value for each isotope
+    (at the above date).
+    """
+
+    isotope: str
+    z: int | None
+    _atomic_weight: sc.Variable | None
+    _atomic_mass: sc.Variable | None
+
+    @property
+    def atomic_weight(self) -> sc.Variable:
+        """Return the atomic weight is available."""
+        if self._atomic_weight is None:
+            raise ValueError(
+                f"Atomic weight for '{self.isotope}' is not defined ."
+                "This likely means that there is no standard atomic weight "
+                "for this element."
+            )
+        return self._atomic_weight.copy()
+
+    @property
+    def atomic_mass(self) -> sc.Variable:
+        """Return the atomic mass is available."""
+        if self._atomic_mass is None:
+            raise ValueError(
+                f"Atomic mass for '{self.isotope}' is not defined ."
+                "This likely means that you specified an element name, not a specific "
+                "isotope; atomic masses are only defined for isotopes / nuclides."
+            )
+        return self._atomic_mass.copy()
+
+    def __eq__(self, other: object) -> bool | type(NotImplemented):
+        if not isinstance(other, Atom):
+            return NotImplemented
+        return all(
+            _eq_or_identical(getattr(self, field.name), getattr(other, field.name))
+            for field in dataclasses.fields(self)
+        )
+
+    @staticmethod
+    @lru_cache
+    def for_isotope(isotope: str) -> Atom:
+        """Return the atom parameters for the given element / isotope.
+
+        Parameters
+        ----------
+        isotope:
+            Name of the element or isotope.
+            For example, 'H', '3He', 'V', '50V'.
+
+        Returns
+        -------
+        :
+            Atom parameters.
+        """
+        element = _parse_isotope_name(isotope)
+        z, weight = _load_atomic_weight(element)
+        if element == isotope:
+            mass = None  # masses are only defined for specific isotopes
+        else:
+            mass = _load_atomic_mass(isotope)
+        return Atom(
+            isotope=isotope,
+            z=z,
+            _atomic_weight=weight,
+            _atomic_mass=mass,
+        )
+
+
 def reference_wavelength() -> sc.Variable:
     """Return the reference wavelength for absorption cross-sections.
 
@@ -62,9 +143,7 @@ def __eq__(self, other: object) -> bool | type(NotImplemented):
         if not isinstance(other, ScatteringParams):
             return NotImplemented
         return all(
-            self.isotope == other.isotope
-            if field.name == 'isotope'
-            else _eq_or_identical(getattr(self, field.name), getattr(other, field.name))
+            _eq_or_identical(getattr(self, field.name), getattr(other, field.name))
             for field in dataclasses.fields(self)
         )
 
@@ -84,37 +163,65 @@ def for_isotope(isotope: str) -> ScatteringParams:
         :
             Neutron scattering parameters.
         """
-        with _open_scattering_parameters_file() as f:
-            while line := f.readline():
-                name, rest = line.split(',', 1)
-                if name == isotope:
-                    return _parse_line(isotope, rest)
+        with _open_bundled_parameters_file('scattering_parameters.csv') as f:
+            if line_remainder := _find_line_with_isotope(isotope, f):
+                return ScatteringParams._parse_line(isotope, line_remainder)
         raise ValueError(f"No entry for element / isotope '{isotope}'")
 
+    @staticmethod
+    def _parse_line(isotope: str, line: str) -> ScatteringParams:
+        line = line.rstrip().split(',')
+        return ScatteringParams(
+            isotope=isotope,
+            coherent_scattering_length_re=_assemble_scalar(line[0], line[1], 'fm'),
+            coherent_scattering_length_im=_assemble_scalar(line[2], line[3], 'fm'),
+            incoherent_scattering_length_re=_assemble_scalar(line[4], line[5], 'fm'),
+            incoherent_scattering_length_im=_assemble_scalar(line[6], line[7], 'fm'),
+            coherent_scattering_cross_section=_assemble_scalar(
+                line[8], line[9], 'barn'
+            ),
+            incoherent_scattering_cross_section=_assemble_scalar(
+                line[10], line[11], 'barn'
+            ),
+            total_scattering_cross_section=_assemble_scalar(line[12], line[13], 'barn'),
+            absorption_cross_section=_assemble_scalar(line[14], line[15], 'barn'),
+        )
+
+
+def _open_bundled_parameters_file(name: str) -> TextIO:
+    return importlib.resources.files('scippneutron.atoms').joinpath(name).open('r')
+
 
-def _open_scattering_parameters_file() -> TextIO:
-    return (
-        importlib.resources.files('scippneutron.atoms')
-        .joinpath('scattering_parameters.csv')
-        .open('r')
-    )
+def _find_line_with_isotope(isotope: str, io: TextIO) -> str | None:
+    while line := io.readline():
+        name, rest = line.split(',', 1)
+        if name == isotope:
+            return rest
+    return None
 
 
-def _parse_line(isotope: str, line: str) -> ScatteringParams:
-    line = line.rstrip().split(',')
-    return ScatteringParams(
-        isotope=isotope,
-        coherent_scattering_length_re=_assemble_scalar(line[0], line[1], 'fm'),
-        coherent_scattering_length_im=_assemble_scalar(line[2], line[3], 'fm'),
-        incoherent_scattering_length_re=_assemble_scalar(line[4], line[5], 'fm'),
-        incoherent_scattering_length_im=_assemble_scalar(line[6], line[7], 'fm'),
-        coherent_scattering_cross_section=_assemble_scalar(line[8], line[9], 'barn'),
-        incoherent_scattering_cross_section=_assemble_scalar(
-            line[10], line[11], 'barn'
-        ),
-        total_scattering_cross_section=_assemble_scalar(line[12], line[13], 'barn'),
-        absorption_cross_section=_assemble_scalar(line[14], line[15], 'barn'),
-    )
+def _load_atomic_weight(element: str) -> tuple[int, sc.Variable | None]:
+    # The CSV file was extracted from https://www.ciaaw.org/abridged-atomic-weights.htm
+    # using the notebook in tools/atomic_weights.ipynb (in the ScippNeutron repo).
+    with _open_bundled_parameters_file('atomic_weights.csv') as f:
+        f.readline()  # skip copyright
+        f.readline()  # skip header
+        if line_remainder := _find_line_with_isotope(element, f):
+            z, weight, error = line_remainder.rstrip().split(',')
+            return int(z), _assemble_scalar(weight, error, 'Da')
+    raise ValueError(f"No entry for element '{element}'")
+
+
+def _load_atomic_mass(isotope: str) -> sc.Variable | None:
+    # The CSV file was extracted from https://www.ciaaw.org/atomic-masses.htm
+    # using the notebook in tools/atomic_weights.ipynb (in the ScippNeutron repo).
+    with _open_bundled_parameters_file('atomic_masses.csv') as f:
+        f.readline()  # skip copyright
+        f.readline()  # skip header
+        if line_remainder := _find_line_with_isotope(isotope, f):
+            weight, error = line_remainder.rstrip().split(',')
+            return _assemble_scalar(weight, error, 'Da')
+    raise ValueError(f"No entry for element / isotope '{isotope}'")
 
 
 def _assemble_scalar(value: str, std: str, unit: str) -> sc.Variable | None:
@@ -125,7 +232,14 @@ def _assemble_scalar(value: str, std: str, unit: str) -> sc.Variable | None:
     return sc.scalar(value, variance=variance, unit=unit)
 
 
-def _eq_or_identical(a: sc.Variable | None, b: sc.Variable | None) -> bool:
-    if a is None:
-        return b is None
-    return sc.identical(a, b)
+def _eq_or_identical(a: object, b: object) -> bool:
+    if isinstance(a, sc.Variable) or isinstance(b, sc.Variable):
+        return sc.identical(a, b)
+    return a == b
+
+
+def _parse_isotope_name(name: str) -> str:
+    # Extract the element name from an isotope name.
+    # 'H' -> 'H'
+    # '2H' -> 'H'
+    return re.match(r'(?:\d+)?([a-zA-Z]+)', name)[1]