Abins: improve handling of isotopic cross-sections

This is a squashed version of mantidproject#38448

Create failing test: Abins get_cross_section for Zn65

This isotope currently gives a NaN cross-section. We should ensure
that unsupported isotopes give a noisy error instead.

Implement ValueError for NaN cross-section data

Fix missing return value, test good cases

Begin refactoring Abins isotope handling

- remove "substitution" logic branch
- instead delegate isotope detection and workspace naming to a new class
- in this commit a new Atom is still (unnecessarily) spawned to get cross-section

Continue refactor: streamline _create_workspace

We can now pass on the AtomInfo object rather than its consitituent data.

Continue refactor: move AtomInfo into _fill_s_workspace methods

Simplify get_cross_section from an established AtomInfo

Handle neutron data for atoms where one isotope = mixture

In cases such as F, one isotope is so dominant that essentially we
would always use the properties of the standard mixture. Oddly, in
this case Mantid will not find neutron data if the atomic mass was
provided.

i.e. Atom('F').neutron() works, but Atom('F', 19).neutron() returns a
dict full of NaN, despite them having exactly the same mass.

Update Abins2D to use AtomInfo class

Handle some more Abins isotope challenging cases

Abins input files usually give the element symbol and mass. These
masses often do not match exactly to those in Mantid's data.

- In the case that the nearest isotope has no data, but the "mixture"
  is _quite_ near the target mass (0.01 amu, for now), log a warning and
  use the mixture.

  This helps handle nasty cases like 127I where Mantid has _slightly_
  different masses for the pure isotope and mixture, but essentially
  it is a 100% mixture and there is no neutron-scattering data for the
  isolated isotope.

- If the nearest isotope has no data and the standard mixture is too
  far away, raise an error.

Direct unit test for AtomInfo

Test refactor: move error check closer to source

The error is now raised by the data class when we try to access a
property that uses Atom.

As such it is not really a get_cross_section test any more and can be
moved out.

Add release note

Tweak release note format

Update scripts/abins/abinsalgorithm.py

Co-authored-by: Adri Diaz <146007827+adriazalvarez@users.noreply.github.com>

Move AtomInfo to separate file

Tighten epsilon for "same" mass values, re-use in atominfo

This gets rid of a hard-coded epsilon. Strictly these two comparisons
are not quite doing the same thing, but this precision should be good
for both.

Author: ajjackson

Framework/PythonInterface/plugins/algorithms/Abins.py

@@ -7,14 +7,16 @@
 
 from typing import Dict
 
+import numpy as np
+
 from mantid.api import AlgorithmFactory, FileAction, FileProperty, PythonAlgorithm, Progress
 from mantid.api import WorkspaceFactory, AnalysisDataService
 
 # noinspection PyProtectedMember
 from mantid.simpleapi import ConvertUnits, GroupWorkspaces, Load
 from mantid.kernel import Direction, StringListValidator
 import abins
-from abins.abinsalgorithm import AbinsAlgorithm
+from abins.abinsalgorithm import AbinsAlgorithm, AtomInfo
 
 
 # noinspection PyPep8Naming,PyMethodMayBeStatic
@@ -184,14 +186,13 @@ def PyExec(self):
         self.setProperty("OutputWorkspace", self._out_ws_name)
         prog_reporter.report("Group workspace with all required  dynamical structure factors has been constructed.")
 
-    def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None):
+    def _fill_s_workspace(self, *, s_points: np.ndarray, workspace: str, species: AtomInfo | None = None):
         """
         Puts S into workspace(s).
 
         :param s_points: dynamical factor for the given atom
         :param workspace:  workspace to be filled with S
-        :param protons_number: number of protons in the given type fo atom
-        :param nucleons_number: number of nucleons in the given type of atom
+        :param species: atom/isotope identity and data
         """
 
         from abins.constants import FUNDAMENTALS, ONE_DIMENSIONAL_INSTRUMENTS, ONE_DIMENSIONAL_SPECTRUM
@@ -201,15 +202,11 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None,
 
         # only FUNDAMENTALS [data is 2d with one row]
         if s_points.shape[0] == FUNDAMENTALS:
-            self._fill_s_1d_workspace(
-                s_points=s_points[0], workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number
-            )
+            self._fill_s_1d_workspace(s_points=s_points[0], workspace=workspace, species=species)
 
         # total workspaces [data is 1d vector]
         elif len(s_points.shape) == ONE_DIMENSIONAL_SPECTRUM:
-            self._fill_s_1d_workspace(
-                s_points=s_points, workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number
-            )
+            self._fill_s_1d_workspace(s_points=s_points, workspace=workspace, species=species)
 
         # quantum order events (fundamentals  or  overtones + combinations for the given order)
         # [data is 2d table of S with a row for each quantum order]
@@ -221,28 +218,19 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None,
                 wrk_name = f"{workspace}_quantum_event_{n + 1}"
                 partial_wrk_names.append(wrk_name)
 
-                self._fill_s_1d_workspace(
-                    s_points=s_points[n], workspace=wrk_name, protons_number=protons_number, nucleons_number=nucleons_number
-                )
+                self._fill_s_1d_workspace(s_points=s_points[n], workspace=wrk_name, species=species)
 
             GroupWorkspaces(InputWorkspaces=partial_wrk_names, OutputWorkspace=workspace)
 
-    def _fill_s_1d_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None):
+    def _fill_s_1d_workspace(self, *, s_points: np.ndarray, workspace: str, species: AtomInfo | None = None):
         """
         Puts 1D S into workspace.
-        :param protons_number: number of protons in the given type of atom
-        :param nucleons_number: number of nucleons in the given type of atom
         :param s_points: dynamical factor for the given atom
         :param workspace: workspace to be filled with S
+        :param species: atom/isotope identity and data
         """
-        if protons_number is not None:
-            s_points = (
-                s_points
-                * self._scale
-                * self.get_cross_section(
-                    scattering=self._scale_by_cross_section, protons_number=protons_number, nucleons_number=nucleons_number
-                )
-            )
+        if species is not None:
+            s_points = s_points * self._scale * self.get_cross_section(scattering=self._scale_by_cross_section, species=species)
         dim = 1
         length = s_points.size
 

Framework/PythonInterface/plugins/algorithms/Abins2D.py

@@ -8,14 +8,16 @@
 import numbers
 from typing import Dict
 
+import numpy as np
+
 from mantid.api import AlgorithmFactory, PythonAlgorithm, Progress
 from mantid.api import WorkspaceFactory, AnalysisDataService
 from mantid.kernel import logger
 
 # noinspection PyProtectedMember
 from mantid.simpleapi import GroupWorkspaces
 import abins
-from abins.abinsalgorithm import AbinsAlgorithm
+from abins.abinsalgorithm import AbinsAlgorithm, AtomInfo
 
 
 # noinspection PyPep8Naming,PyMethodMayBeStatic
@@ -184,14 +186,13 @@ def PyExec(self):
         self.setProperty("OutputWorkspace", self._out_ws_name)
         prog_reporter.report("Group workspace with all required  dynamical structure factors has been constructed.")
 
-    def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None):
+    def _fill_s_workspace(self, *, s_points=None, workspace: str, species: AtomInfo | None = None):
         """
         Puts S into workspace(s).
 
         :param s_points: dynamical factor for the given atom
         :param workspace:  workspace to be filled with S
-        :param protons_number: number of protons in the given type fo atom
-        :param nucleons_number: number of nucleons in the given type of atom
+        :param species: Atom/isotope data
         """
         from abins.constants import FUNDAMENTALS, TWO_DIMENSIONAL_INSTRUMENTS
 
@@ -200,15 +201,11 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None,
 
         # only FUNDAMENTALS [data is 3d with length 1 in axis 0]
         if s_points.shape[0] == FUNDAMENTALS:
-            self._fill_s_2d_workspace(
-                s_points=s_points[0], workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number
-            )
+            self._fill_s_2d_workspace(s_points=s_points[0], workspace=workspace, species=species)
 
         # total workspaces [data is 2d array of S]
         elif s_points.shape[0] == abins.parameters.instruments[self._instrument.get_name()]["q_size"]:
-            self._fill_s_2d_workspace(
-                s_points=s_points, workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number
-            )
+            self._fill_s_2d_workspace(s_points=s_points, workspace=workspace, species=species)
 
         # Multiple quantum order events [data is 3d table of S using axis 0 for quantum orders]
         else:
@@ -219,9 +216,7 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None,
                 wrk_name = f"{workspace}_quantum_event_{n + 1}"
                 partial_wrk_names.append(wrk_name)
 
-                self._fill_s_2d_workspace(
-                    s_points=s_points[n], workspace=wrk_name, protons_number=protons_number, nucleons_number=nucleons_number
-                )
+                self._fill_s_2d_workspace(s_points=s_points[n], workspace=wrk_name, species=species)
 
                 GroupWorkspaces(InputWorkspaces=partial_wrk_names, OutputWorkspace=workspace)
 
@@ -232,14 +227,12 @@ def _create_dummy_workspace(self, name):
         AnalysisDataService.addOrReplace(name, wrk)
         return wrk
 
-    def _fill_s_2d_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None):
+    def _fill_s_2d_workspace(self, *, s_points: np.ndarray, workspace: str, species: AtomInfo | None = None):
         from mantid.api import NumericAxis
         from abins.constants import MILLI_EV_TO_WAVENUMBER
 
-        if protons_number is not None:
-            s_points = s_points * self.get_cross_section(
-                scattering=self._scale_by_cross_section, protons_number=protons_number, nucleons_number=nucleons_number
-            )
+        if species is not None:
+            s_points = s_points * self.get_cross_section(scattering=self._scale_by_cross_section, species=species)
 
         n_q_values, n_freq_bins = s_points.shape
         n_q_bins = self._q_bins.size

docs/source/release/v6.12.0/Indirect/Algorithms/Bugfixes/38448.rst

@@ -0,0 +1,12 @@
+- The identification of isotopes has been improved in Abins/Abins2D.
+
+  - Previously, a species would only be identified as the standard
+    isotopic mixture if the mass is very close to the Mantid reference
+    data. In some cases the values used by external phonon calculators
+    are significantly different and this could lead to misassignment.
+    (e.g. Zn in CASTEP.) Now, Abins will initially choose the _nearest_
+    mass option between an isotope and the standard isotopic mixture.
+  - Many isotopes lack cross-section data and would lead to NaN
+    intensities. Now, if a NaN cross-section is identified Abins
+    will either use the standard mixture data (if the mass is within
+    0.01) or raise an error.

scripts/abins/abins2.py

@@ -7,13 +7,16 @@
 
 from typing import Dict, Optional
 
+import numpy as np
+
 from mantid.api import AlgorithmFactory, PythonAlgorithm, Progress
 from mantid.api import WorkspaceFactory, AnalysisDataService
 
 # noinspection PyProtectedMember
 from mantid.simpleapi import ConvertUnits, GroupWorkspaces
 import abins
 from abins.abinsalgorithm import AbinsAlgorithm
+from abins.atominfo import AtomInfo
 from abins.logging import get_logger, Logger
 
 
@@ -170,14 +173,13 @@ def PyExec(self):
         self.setProperty("OutputWorkspace", self._out_ws_name)
         prog_reporter.report("Group workspace with all required  dynamical structure factors has been constructed.")
 
-    def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None):
+    def _fill_s_workspace(self, *, s_points: np.ndarray, workspace: str, species: AtomInfo | None = None) -> None:
         """
         Puts S into workspace(s).
 
         :param s_points: dynamical factor for the given atom
-        :param workspace:  workspace to be filled with S
-        :param protons_number: number of protons in the given type fo atom
-        :param nucleons_number: number of nucleons in the given type of atom
+        :param workspace: workspace to be filled with S
+        :param species: Element/isotope data
         """
 
         from abins.constants import FUNDAMENTALS, ONE_DIMENSIONAL_INSTRUMENTS, ONE_DIMENSIONAL_SPECTRUM
@@ -187,15 +189,11 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None,
 
         # only FUNDAMENTALS [data is 2d with one row]
         if s_points.shape[0] == FUNDAMENTALS:
-            self._fill_s_1d_workspace(
-                s_points=s_points[0], workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number
-            )
+            self._fill_s_1d_workspace(s_points=s_points[0], workspace=workspace, species=species)
 
         # total workspaces [data is 1d vector]
         elif len(s_points.shape) == ONE_DIMENSIONAL_SPECTRUM:
-            self._fill_s_1d_workspace(
-                s_points=s_points, workspace=workspace, protons_number=protons_number, nucleons_number=nucleons_number
-            )
+            self._fill_s_1d_workspace(s_points=s_points, workspace=workspace, species=species)
 
         # quantum order events (fundamentals  or  overtones + combinations for the given order)
         # [data is 2d table of S with a row for each quantum order]
@@ -207,24 +205,20 @@ def _fill_s_workspace(self, s_points=None, workspace=None, protons_number=None,
                 wrk_name = f"{workspace}_quantum_event_{n + 1}"
                 partial_wrk_names.append(wrk_name)
 
-                self._fill_s_1d_workspace(
-                    s_points=s_points[n], workspace=wrk_name, protons_number=protons_number, nucleons_number=nucleons_number
-                )
+                self._fill_s_1d_workspace(s_points=s_points[n], workspace=wrk_name, species=species)
 
             GroupWorkspaces(InputWorkspaces=partial_wrk_names, OutputWorkspace=workspace)
 
-    def _fill_s_1d_workspace(self, s_points=None, workspace=None, protons_number=None, nucleons_number=None):
+    def _fill_s_1d_workspace(self, s_points=None, workspace=None, species: AtomInfo | None = None):
         """
         Puts 1D S into workspace.
         :param protons_number: number of protons in the given type of atom
         :param nucleons_number: number of nucleons in the given type of atom
         :param s_points: dynamical factor for the given atom
         :param workspace: workspace to be filled with S
         """
-        if protons_number is not None:
-            s_points = s_points * self.get_cross_section(
-                scattering=self._scale_by_cross_section, protons_number=protons_number, nucleons_number=nucleons_number
-            )
+        if species is not None:
+            s_points = s_points * self.get_cross_section(scattering=self._scale_by_cross_section, species=species)
         dim = 1
         length = s_points.size