Abins: improve handling of isotopic cross-sections

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,13 @@
+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,15 @@
 
 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.abinsalgorithm import AbinsAlgorithm, AtomInfo
 from abins.logging import get_logger, Logger
 
 
@@ -170,14 +172,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 +188,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 +204,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