From cbe9cc8c290d675b515fa40507fd4179c74a9b5b Mon Sep 17 00:00:00 2001
From: Joffrey Dorville <joffreydorville@mines.edu>
Date: Fri, 1 Apr 2022 17:31:16 -0600
Subject: [PATCH 01/12] Kalbach-Mann slope calculation for ENDF files when the
 projectile is a neutron

---
 openmc/data/kalbach_mann.py                | 528 ++++++++++++++++++++-
 openmc/data/njoy.py                        |  11 +-
 openmc/data/reaction.py                    |  29 +-
 tests/unit_tests/test_data_kalbach_mann.py | 280 +++++++++++
 4 files changed, 834 insertions(+), 14 deletions(-)
 create mode 100644 tests/unit_tests/test_data_kalbach_mann.py

diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py
index f4c914d90..11e3acf40 100644
--- a/openmc/data/kalbach_mann.py
+++ b/openmc/data/kalbach_mann.py
@@ -5,6 +5,7 @@ from warnings import warn
 import numpy as np
 
 import openmc.checkvalue as cv
+from openmc.mixin import EqualityMixin
 from openmc.stats import Tabular, Univariate, Discrete, Mixture
 from .function import Tabulated1D, INTERPOLATION_SCHEME
 from .angle_energy import AngleEnergy
@@ -12,6 +13,468 @@ from .data import EV_PER_MEV
 from .endf import get_list_record, get_tab2_record
 
 
+# Kalbach-Mann constants as defined in ENDF-6 manual BNL-203218-2018-INRE,
+# Revision 215, File 6 description for LAW=1 and LANG=2.
+_C1 = 0.04  # [1/MeV]
+_C2 = 1.8E-6  # [1/MeV^3]
+_C3 = 6.7E-7  # [1/MeV^4]
+_ET1 = 130.  # [MeV]
+_ET3 = 41.  # [MeV]
+_M_NEUTRON = 1.
+_M_PROTON = 1.
+_M_DEUTERON = 1.
+_M_TRITON = None
+_M_3HE = None
+_M_ALPHA = 0.
+_SM_NEUTRON = 1/2.
+_SM_PROTON = 1.
+_SM_DEUTERON = 1.
+_SM_TRITON = 1.
+_SM_3HE = 1.
+_SM_ALPHA = 2.
+
+# Kalbach-Mann M coefficients
+_TABULATED_PARTICLE_M = {
+    1: _M_NEUTRON,
+    1001: _M_PROTON,
+    1002: _M_DEUTERON,
+    1003: _M_TRITON,
+    2003: _M_3HE,
+    2004: _M_ALPHA
+}
+
+# Kalbach-Mann m coefficients
+_TABULATED_PARTICLE_SM = {
+    1: _SM_NEUTRON,
+    1001: _SM_PROTON,
+    1002: _SM_DEUTERON,
+    1003: _SM_TRITON,
+    2003: _SM_3HE,
+    2004: _SM_ALPHA
+}
+
+# Breaking energy as defined in ENDF-6 manual BNL-203218-2018-INRE,
+# Revision 215, Appendix H, Table 3.
+_BREAKING_ENERGY_NEUTRON = 0.
+_BREAKING_ENERGY_PROTON = 0.
+_BREAKING_ENERGY_DEUTERON = 2.224566  # [MeV]
+_BREAKING_ENERGY_TRITON = 8.481798  # [MeV]
+_BREAKING_ENERGY_3HE = 7.718043  # [MeV]
+_BREAKING_ENERGY_ALPHA = 28.29566  # [MeV]
+
+_TABULATED_BREAKING_ENERGY = {
+    1: _BREAKING_ENERGY_NEUTRON,
+    1001: _BREAKING_ENERGY_PROTON,
+    1002: _BREAKING_ENERGY_DEUTERON,
+    1003: _BREAKING_ENERGY_TRITON,
+    2003: _BREAKING_ENERGY_3HE,
+    2004: _BREAKING_ENERGY_ALPHA
+}
+
+# Abundant IZA translation in merged library
+_IZA_TRANSLATION = {
+    6000: 6012,
+}
+
+
+class AtomicRepresentation(EqualityMixin):
+    """Atomic representation of an isotope or a particle.
+
+    Parameters
+    ----------
+    z: int
+        Number of protons (atomic number)
+    a: int
+        Number of nucleons (mass number)
+
+    Raises
+    ------
+    IOError:
+        When the number of protons (z) declared is higher than the number
+        of nucleons (a)
+
+    Attributes
+    ----------
+    z: int
+        Number of protons (atomic number)
+    a: int
+        Number of nucleons (mass number)
+    n: int
+        Number of neutrons
+    breaking_energy: float
+        Energy required to break the isotope or particle into their
+        constituent nucleons from tabulated values
+    M: float
+        Kalbach-Mann M coefficient
+    m: float
+        Kalbach-Mann m coefficient
+    iza: int
+        ZA identifier defined as:
+        iza = Z x 1000 + A,
+        where Z is the number of protons and A the number of nucleons
+
+    """
+    def __init__(self, z, a):
+        self._consistency_check(z, a)
+        self._z = z
+        self._a = a
+        self.z = self._z
+        self.a = self._a
+
+    def __add__(self, other):
+        """Adds two AtomicRepresentations.
+
+        """
+        z = self.z + other.z
+        a = self.a + other.a
+        return AtomicRepresentation(z=z, a=a)
+
+    def __sub__(self, other):
+        """Substracts two AtomicRepresentations.
+
+        """
+        z = self.z - other.z
+        a = self.a - other.a
+        return AtomicRepresentation(z=z, a=a)
+
+    @property
+    def a(self):
+        self._consistency_check(self._z, self._a)
+        return self._a
+
+    @property
+    def z(self):
+        self._consistency_check(self._z, self._a)
+        return self._z
+
+    @property
+    def n(self):
+        return self.a - self.z
+
+    @property
+    def breaking_energy(self):
+        breaking_energy = None
+        if self.iza in _TABULATED_BREAKING_ENERGY:
+            breaking_energy = _TABULATED_BREAKING_ENERGY[self.iza]
+        return breaking_energy
+
+    @property
+    def M(self):
+        M = None
+        if self.iza in _TABULATED_PARTICLE_M:
+            M = _TABULATED_PARTICLE_M[self.iza]
+        return M
+
+    @property
+    def m(self):
+        m = None
+        if self.iza in _TABULATED_PARTICLE_SM:
+            m = _TABULATED_PARTICLE_SM[self.iza]
+        return m
+
+    @property
+    def iza(self):
+        iza = self.z * 1000 + self.a
+        return iza
+
+    @a.setter
+    def a(self, an):
+        cv.check_type('a', an, Integral)
+        cv.check_greater_than('a', an, 0, equality=True)
+        self._consistency_check(self._z, an)
+        self._a = an
+
+    @z.setter
+    def z(self, zn):
+        cv.check_type('z', zn, Integral)
+        cv.check_greater_than('z', zn, 0, equality=True)
+        self._consistency_check(zn, self._a)
+        self._z = zn
+
+    @staticmethod
+    def _consistency_check(z, a):
+        """Simple consistency check.
+
+        Parameters
+        ----------
+        z: int
+            Number of protons (atomic number)
+        a: int
+            Number of nucleons (mass number)
+
+        Raises
+        ------
+        IOError:
+            When the number of protons (z) declared is higher than the number
+            of nucleons (a)
+
+        """
+        if z > a:
+            raise IOError(
+                "Number of protons (%i) incompatible with number of "
+                "nucleons (%i)" % (z, a)
+            )
+
+    @classmethod
+    def from_iza(cls, iza):
+        """Instantiates an AtomicRepresentation from a ZA identifier.
+
+        The ZA identifier is defined as:
+        iza = Z x 1000 + A,
+        where Z is the number of protons and A the number of nucleons.
+
+        Parameters
+        ----------
+        iza: int
+            ZA identifier
+
+        Returns
+        -------
+        AtomicRepresentation
+            Atomic representation of the isotope/particle
+
+        """
+        if iza in _IZA_TRANSLATION.keys():
+            iza = _IZA_TRANSLATION[iza]
+
+        z = int(iza/1000)
+        a = np.mod(iza, 1000)
+
+        return cls(z, a)
+
+
+def _calculate_separation_energy(compound, nucleus, particle):
+    """Calculates the separation energy as defined in ENDF-6 manual
+    BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1
+    and LANG=2. This function can be used for the incident or emitted
+    particle of the following reaction: A + a -> C -> B + b
+
+    Parameters
+    ----------
+    compound: AtomicRepresentation
+        Atomic representation of the compound (C)
+    nucleus: AtomicRepresentation
+        Atomic representation of the nucleus (A or B)
+    particle: AtomicRepresentation
+        Atomic representation of the particle (a or b)
+
+    Returns
+    -------
+    separation_energy: float
+        Separation energy in MeV
+
+    """
+    coef_1 = 15.68 * (compound.a - nucleus.a)
+    coef_2 = 28.07 * ((compound.n - compound.z)**2 / float(compound.a) - \
+        (nucleus.n - nucleus.z)**2 / float(nucleus.a))
+    coef_3 = 18.56 * (compound.a**(2./3.) - nucleus.a**(2./3.))
+    coef_4 = 33.22 * ((compound.n - compound.z)**2 / float(compound.a)**(4./3.) - \
+        (nucleus.n - nucleus.z)**2 / float(nucleus.a)**(4./3.))
+    coef_5 = 0.717 * (compound.z**2 / float(compound.a)**(1./3.) - \
+        nucleus.z**2 / float(nucleus.a)**(1./3.))
+    coef_6 = 1.211 * (compound.z**2 / float(compound.a) - \
+        nucleus.z**2 / float(nucleus.a))
+
+    separation_energy = coef_1 - coef_2 - coef_3 + coef_4 \
+        - coef_5 + coef_6 - particle.breaking_energy
+
+    return separation_energy
+
+
+def _return_entrance_channel_energy(e_p, awr_t, awr_p):
+    """Returns the entrance channel energy as defined in ENDF-6 manual
+    BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1
+    and LANG=2.
+
+    Parameters
+    ----------
+    e_p: float
+        Energy of the incident projectile in the laboratory system in eV
+    awr_t: float
+        Atomic weight ratio of the target
+    awr_p: float
+        Atomic weight ratio of the projectile
+
+    Returns
+    -------
+    epsilon_p: float
+        Entrance channel energy in eV
+
+    """
+    epsilon_p = e_p * awr_t / (awr_t + awr_p)
+    return epsilon_p
+
+
+def _return_emission_channel_energy(e_e, awr_r, awr_e):
+    """Returns the emission channel energy as defined in ENDF-6 manual
+    BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1
+    and LANG=2.
+
+    Parameters
+    ----------
+    e_e: float
+        Energy of the emitted particle in the center of mass system in eV
+    awr_r: float
+        Atomic weight ratio of the residual nucleus
+    awr_e: float
+        Atomic weight ratio of the emitted particle
+
+    Returns
+    -------
+    epsilon_e: float
+        Emission channel energy in eV
+
+    """
+    epsilon_e = e_e * (awr_r + awr_e) / awr_r
+    return epsilon_e
+
+
+def _calculate_kalbach_slope(energy_projectile,
+                             energy_emitted,
+                             projectile,
+                             target,
+                             compound,
+                             emitted,
+                             residual):
+    """Calculate the Kalbach slope for projectiles other than photons
+    as defined in ENDF-6 manual BNL-203218-2018-INRE, Revision 215,
+    File 6 description for LAW=1 and LANG=2.
+
+    The entrance and emission channel energies are not calculated with
+    the AWR number, but approximated with the number of mass instead.
+
+    Parameters
+    ----------
+    energy_projectile: float
+        Energy of the projectile in the laboratory system in eV
+    energy_emitted: float
+        Energy of the emitted particle in the center of mass system in eV
+    projectile: AtomicRepresentation
+        Atomic representation of the projectile
+    target: AtomicRepresentation
+        Atomic representation of the target
+    compound: AtomicRepresentation
+        Atomic representation of the compound
+    emitted: AtomicRepresentation
+        Atomic representation of the emitted particle
+    residual: AtomicRepresentation
+        Atomic representation of the residual nucleus
+
+    Returns
+    -------
+    slope: float
+        Kalbach-Mann slope
+
+    """
+    epsilon_a = _return_entrance_channel_energy(
+        energy_projectile,
+        target.a,
+        projectile.a
+    ) / EV_PER_MEV
+    epsilon_b = _return_emission_channel_energy(
+        energy_emitted,
+        residual.a,
+        emitted.a
+    ) / EV_PER_MEV
+
+    s_a = _calculate_separation_energy(compound, target, projectile)
+    s_b = _calculate_separation_energy(compound, residual, emitted)
+
+    e_a = epsilon_a + s_a
+    e_b = epsilon_b + s_b
+
+    r_1 = min(e_a, _ET1)
+    r_3 = min(e_a, _ET3)
+
+    x_1 = r_1 * e_b / e_a
+    x_3 = r_3 * e_b / e_a
+
+    slope = _C1 * x_1 \
+        + _C2 * x_1**3 \
+        + _C3 * projectile.M * emitted.m * x_3**4
+
+    return slope
+
+
+def return_kalbach_slope(energy_projectile,
+                         energy_emitted,
+                         iza_projectile,
+                         iza_emitted,
+                         iza_target):
+    """Returns Kalbach-Mann slope from calculations.
+    
+    The associated reaction is defined as:
+    A + a -> C -> B + b
+
+    Where:
+
+    - A is the targeted nucleus,
+    - a is the projectile,
+    - C is the compound,
+    - B is the residual nucleus,
+    - b is the emitted particle.
+
+    This function uses the concept of ZA identifier defined as:
+    iza = Z x 1000 + A,
+    where Z is the number of protons and A the number of nucleons.
+
+    The Kalbach-Mann slope calculation is done as defined in ENDF-6 manual
+    BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 and
+    LANG=2. One exception to this, is that the entrance and emission channel
+    energies are not calculated with the AWR number, but approximated with
+    the number of mass instead.
+
+    Parameters
+    ----------
+    energy_projectile: float
+        Energy of the projectile in the laboratory system in eV
+    energy_emitted: float
+        Energy of the emitted particle in the center of mass system in eV
+    iza_projectile: int
+        ZA identifier of the projectile
+    iza_emitted: int
+        ZA identifier of the emitted particle
+    iza_target: int
+        ZA identifier of the targeted nucleus
+
+    Raises
+    ------
+    NotImplementedError:
+        When the ZA identifier of the projectile is not equal to 1
+        (ie. other than a neutron).
+
+    Returns
+    -------
+    slope: float
+        Kalbach-Mann slope given with the same format as ACE file.
+
+    """
+    # TODO: develop for photons as projectile
+    # TODO: test for other particles than neutron
+    if iza_projectile != 1:
+        raise NotImplementedError(
+            "Developed and tested for neutron projectile only."
+        )
+
+    projectile = AtomicRepresentation.from_iza(iza_projectile)
+    emitted = AtomicRepresentation.from_iza(iza_emitted)
+    target = AtomicRepresentation.from_iza(iza_target)
+    compound = projectile + target
+    residual = compound - emitted
+
+    slope = _calculate_kalbach_slope(
+        energy_projectile,
+        energy_emitted,
+        projectile,
+        target,
+        compound,
+        emitted,
+        residual
+    )
+
+    return float("%7e" % slope)
+
+
 class KalbachMann(AngleEnergy):
     """Kalbach-Mann distribution
 
@@ -319,7 +782,7 @@ class KalbachMann(AngleEnergy):
             n_energy_out = int(ace.xss[idx + 1])
             data = ace.xss[idx + 2:idx + 2 + 5*n_energy_out].copy()
             data.shape = (5, n_energy_out)
-            data[0,:] *= EV_PER_MEV
+            data[0, :] *= EV_PER_MEV
 
             # Create continuous distribution
             eout_continuous = Tabular(data[0][n_discrete_lines:],
@@ -352,13 +815,28 @@ class KalbachMann(AngleEnergy):
         return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a)
 
     @classmethod
-    def from_endf(cls, file_obj):
-        """Generate Kalbach-Mann distribution from an ENDF evaluation
+    def from_endf(cls, file_obj, iza_emitted, iza_target, projectile_mass):
+        """Generate Kalbach-Mann distribution from an ENDF evaluation.
+
+        If the projectile is a neutron, the slope is calculated when it is
+        not given explicitly.
 
         Parameters
         ----------
         file_obj : file-like object
             ENDF file positioned at the start of the Kalbach-Mann distribution
+        iza_emitted : int
+            ZA identifier of the emitted particle
+        iza_target : int
+            ZA identifier of the target
+        projectile_mass: float
+            Mass of the projectile
+
+        Warns
+        -----
+        UserWarning
+            If the mass of the projectile is not equal to 1 (other than
+            a neutron), the slope is not calculated and set to 0 if missing.
 
         Returns
         -------
@@ -374,6 +852,7 @@ class KalbachMann(AngleEnergy):
         energy_out = []
         precompound = []
         slope = []
+        calculated_slope = []
         for i in range(ne):
             items, values = get_list_record(file_obj)
             energy[i] = items[1]
@@ -385,19 +864,46 @@ class KalbachMann(AngleEnergy):
             values.shape = (n_energy_out, n_angle + 2)
 
             # Outgoing energy distribution at the i-th incoming energy
-            eout_i = values[:,0]
-            eout_p_i = values[:,1]
+            eout_i = values[:, 0]
+            eout_p_i = values[:, 1]
             energy_out_i = Tabular(eout_i, eout_p_i, INTERPOLATION_SCHEME[lep])
             energy_out.append(energy_out_i)
 
-            # Precompound and slope factors for Kalbach-Mann
-            r_i = values[:,2]
+            # Precompound factors for Kalbach-Mann
+            r_i = values[:, 2]
+
+            # Slope factors for Kalbach-Mann
             if n_angle == 2:
-                a_i = values[:,3]
+                a_i = values[:, 3]
+                calculated_slope.append(False)
             else:
-                a_i = np.zeros_like(r_i)
+                # Check if the projectile is not a neutron
+                if not np.isclose(projectile_mass, 1.0, atol=1.0e-12, rtol=0.):
+                    warn(
+                        "Kalbach-Mann slope calculation is only available with "
+                        "neutrons as projectile. Slope coefficients are set to 0."
+                    )
+                    a_i = np.zeros_like(r_i)
+                    calculated_slope.append(False)
+
+                else:
+                    # TODO: retrieve IZA of the projectile
+                    iza_projectile = 1
+                    a_i = [return_kalbach_slope(energy_projectile=energy[i],
+                                                energy_emitted=e,
+                                                iza_projectile=iza_projectile,
+                                                iza_emitted=iza_emitted,
+                                                iza_target=iza_target)
+                           for e in eout_i]
+                    calculated_slope.append(True)
+
             precompound.append(Tabulated1D(eout_i, r_i))
             slope.append(Tabulated1D(eout_i, a_i))
 
-        return cls(tab2.breakpoints, tab2.interpolation, energy,
-                   energy_out, precompound, slope)
+        km_distribution = cls(tab2.breakpoints, tab2.interpolation, energy,
+                              energy_out, precompound, slope)
+
+        # List of bool to indicate slope calculation by OpenMC
+        km_distribution._calculated_slope = calculated_slope
+
+        return km_distribution
diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py
index 305edaf4f..fbb26be05 100644
--- a/openmc/data/njoy.py
+++ b/openmc/data/njoy.py
@@ -127,7 +127,7 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
 1 0 1 .{ext} /
 '{library}: {zsymam} at {temperature}'/
 {mat} {temperature}
-1 1/
+1 1 {ismoothing}/
 /
 """
 
@@ -248,7 +248,8 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
 
 def make_ace(filename, temperatures=None, acer=True, xsdir=None,
              output_dir=None, pendf=False, error=0.001, broadr=True,
-             heatr=True, gaspr=True, purr=True, evaluation=None, **kwargs):
+             heatr=True, gaspr=True, purr=True, evaluation=None,
+             smoothing=True, **kwargs):
     """Generate incident neutron ACE file from an ENDF file
 
     File names can be passed to
@@ -298,6 +299,8 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
     evaluation : openmc.data.endf.Evaluation, optional
         If the ENDF file contains multiple material evaluations, this argument
         indicates which evaluation should be used.
+    smoothing : bool, optional
+        If the smoothing option in ACER is on (1) or off (0) in the card 6.
     **kwargs
         Keyword arguments passed to :func:`openmc.data.njoy.run`
 
@@ -380,6 +383,10 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
 
     # acer
     if acer:
+        if smoothing is True:
+            ismoothing = 1
+        else:
+            ismoothing = 0
         nacer_in = nlast
         for i, temperature in enumerate(temperatures):
             # Extend input with an ACER run for each temperature
diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py
index 5e4287f16..42fc576da 100644
--- a/openmc/data/reaction.py
+++ b/openmc/data/reaction.py
@@ -80,6 +80,14 @@ def _get_products(ev, mt):
     mt : int
         The MT value of the reaction to get products for
 
+    Raises
+    ------
+    IOError:
+        When the Kalbach-Mann systematics is used, but the product
+        is not defined in the 'center-of-mass' system. The breakup logic
+        is not implemented which can lead to this error being raised while
+        the definition of the product is correct.
+
     Returns
     -------
     products : list of openmc.data.Product
@@ -141,7 +149,26 @@ def _get_products(ev, mt):
             if lang == 1:
                 p.distribution = [CorrelatedAngleEnergy.from_endf(file_obj)]
             elif lang == 2:
-                p.distribution = [KalbachMann.from_endf(file_obj)]
+                # Products need to be described in the center-of-mass system
+                product_center_of_mass = False
+                if reference_frame == 'center-of-mass':
+                    product_center_of_mass = True
+                elif reference_frame == 'light-heavy':
+                    product_center_of_mass = (awr <= 4.0)
+                # TODO: 'breakup' logic not implemented
+
+                if product_center_of_mass is False:
+                    raise IOError(
+                        "Kalbach-Mann representation must be defined in the "
+                        "'center-of-mass' system"
+                    )
+
+                zat = ev.target["atomic_number"] * 1000 + ev.target["mass_number"]
+                projectile_mass = ev.projectile["mass"]
+                p.distribution = [KalbachMann.from_endf(file_obj,
+                                                        za,
+                                                        zat,
+                                                        projectile_mass)]
 
         elif law == 2:
             # Discrete two-body scattering
diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py
new file mode 100644
index 000000000..99869808d
--- /dev/null
+++ b/tests/unit_tests/test_data_kalbach_mann.py
@@ -0,0 +1,280 @@
+"""Test of the Kalbach-Mann slope calculation when data are
+retrieved from ENDF files."""
+
+import os
+import pytest
+import numpy as np
+
+from openmc.data import IncidentNeutron
+from openmc.data import AtomicRepresentation
+from openmc.data.kalbach_mann import _BREAKING_ENERGY_TRITON
+from openmc.data.kalbach_mann import _M_TRITON
+from openmc.data.kalbach_mann import _SM_TRITON
+from openmc.data.kalbach_mann import _calculate_separation_energy
+from openmc.data.kalbach_mann import _return_emission_channel_energy
+from openmc.data.kalbach_mann import _return_entrance_channel_energy
+from openmc.data.kalbach_mann import _calculate_kalbach_slope
+from openmc.data import return_kalbach_slope
+from openmc.data import KalbachMann
+
+from . import needs_njoy
+
+
+@pytest.fixture(scope='module')
+def neutron():
+    """Neutron AtomicRepresentation."""
+    return AtomicRepresentation(z=0, a=1)
+
+
+@pytest.fixture(scope='module')
+def triton():
+    """Triton AtomicRepresentation."""
+    return AtomicRepresentation(z=1, a=3)
+
+
+@pytest.fixture(scope='module')
+def b10():
+    """B10 AtomicRepresentation."""
+    return AtomicRepresentation(z=5, a=10)
+
+
+@pytest.fixture(scope='module')
+def c12():
+    """C12 AtomicRepresentation."""
+    return AtomicRepresentation(z=6, a=12)
+
+
+@pytest.fixture(scope='module')
+def c13():
+    """C13 AtomicRepresentation."""
+    return AtomicRepresentation(z=6, a=13)
+
+
+@pytest.fixture(scope='module')
+def na23():
+    """Na23 AtomicRepresentation."""
+    return AtomicRepresentation(z=11, a=23)
+
+
+def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
+    """Test the AtomicRepresentation class."""
+    # Test instanciation from_iza
+    assert b10 == AtomicRepresentation.from_iza(5010)
+
+    # Test instanciation from_iza using IZA translation
+    assert c12 == AtomicRepresentation.from_iza(6000)
+
+    # Test addition
+    assert c13 + b10 == na23
+
+    # Test substraction
+    assert c13 - c12 == neutron
+    assert c13 - b10 == triton
+
+    # Test properties when no information for Kalbach-Mann are given
+    assert c13.a == 13
+    assert c13.z == 6
+    assert c13.n == 7
+    assert c13.breaking_energy is None
+    assert c13.M is None
+    assert c13.m is None
+    assert c13.iza == 6013
+
+    # Test properties when information for Kalbach-Mann are given
+    assert triton.a == 3
+    assert triton.z == 1
+    assert triton.n == 2
+    assert triton.breaking_energy == _BREAKING_ENERGY_TRITON
+    assert triton.M == _M_TRITON
+    assert triton.m == _SM_TRITON
+    assert triton.iza == 1003
+
+    # Test instanciation errors
+    with pytest.raises(IOError):
+        AtomicRepresentation(z=5, a=1)
+    with pytest.raises(ValueError):
+        AtomicRepresentation(z=-1, a=1)
+    with pytest.raises(IOError):
+        AtomicRepresentation(z=5, a=0)
+    with pytest.raises(IOError):
+        AtomicRepresentation(z=5, a=-2)
+    with pytest.raises(OSError):
+        neutron - triton
+
+
+def test__calculate_separation_energy(triton, b10, c13):
+    """Comparison to hand-calculations on a simple example."""
+    assert _calculate_separation_energy(
+        compound=c13,
+        nucleus=b10,
+        particle=triton
+    ) == pytest.approx(18.6880713)
+
+
+def test__return_entrance_channel_energy():
+    """Comparison to hand-calculations on a simple example."""
+    assert _return_entrance_channel_energy(
+        e_p=5.2,
+        awr_t=13.7,
+        awr_p=7.2
+    ) == pytest.approx(3.4086124)
+
+
+def test__return_emission_channel_energy():
+    """Comparison to hand-calculations on a simple example."""
+    assert _return_emission_channel_energy(
+        e_e=6.8,
+        awr_r=49.2,
+        awr_e=5.4
+    ) == pytest.approx(7.5463415)
+
+
+def test__calculate_kalbach_slope(neutron, triton, b10, c12, c13):
+    """Comparison to hand-calculations for n + c12 -> c13 -> triton + b10."""
+    energy_projectile = 10.2  # [eV]
+    energy_emitted = 5.4  # [eV]
+
+    assert _calculate_kalbach_slope(
+        energy_projectile=energy_projectile,
+        energy_emitted=energy_emitted,
+        projectile=neutron,
+        target=c12,
+        compound=c13,
+        emitted=triton,
+        residual=b10
+    ) == pytest.approx(0.8409921475)
+
+
+def test_return_kalbach_slope():
+    """Comparison to hand-calculations for n + c12 -> c13 -> triton + b10."""
+    energy_projectile = 10.2  # [eV]
+    energy_emitted = 5.4  # [eV]
+
+    # Check that NotImplementedError is raised if the projectile is not
+    # a neutron
+    with pytest.raises(NotImplementedError):
+        return_kalbach_slope(
+            energy_projectile=energy_projectile,
+            energy_emitted=energy_emitted,
+            iza_projectile=1000,
+            iza_emitted=1,
+            iza_target=6012
+        )
+
+    assert return_kalbach_slope(
+        energy_projectile=energy_projectile,
+        energy_emitted=energy_emitted,
+        iza_projectile=1,
+        iza_emitted=1003,
+        iza_target=6012
+    ) == pytest.approx(0.8409921475)
+
+
+@pytest.mark.parametrize(
+    "hdf5_filename, endf_type, endf_filename", [
+        ('O16.h5', 'neutrons', 'n-008_O_016.endf'),
+        ('Ca46.h5', 'neutrons', 'n-020_Ca_046.endf'),
+        ('Hg204.h5', 'neutrons', 'n-080_Hg_204.endf')
+    ]
+)
+def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
+    """Test the calculation of the Kalbach-Mann slope done by OpenMC
+    by comparing it to HDF5 data. The test is based on the first product
+    of MT=5 (neutron). The isotopes tested have been selected because the
+    corresponding products in ENDF/B-VII.1 are described using MF=6, LAW=1,
+    LANG=2 (ie. Kalbach-Mann systematics) and the slope is not given
+    explicitly.
+
+    If an error occurs during the "validity check", this means that
+    the nuclear data evaluation has evolved and the distribution might
+    no longer be described using Kalbach-Mann systematics. Another
+    isotope needs to be identified and tested.
+
+    Warning: This test is valid as long as ENDF files are not directly
+    used to generate the HDF5 files used in the tests.
+
+    """
+    # HDF5 data
+    hdf5_directory = os.path.dirname(os.environ['OPENMC_CROSS_SECTIONS'])
+    hdf5_path = os.path.join(hdf5_directory, hdf5_filename)
+    hdf5_data = IncidentNeutron.from_hdf5(hdf5_path)
+    hdf5_product = hdf5_data[5].products[0]
+    hdf5_distribution = hdf5_product.distribution[0]
+
+    # ENDF data
+    endf_directory = os.environ['OPENMC_ENDF_DATA']
+    endf_path = os.path.join(endf_directory, endf_type, endf_filename)
+    endf_data = IncidentNeutron.from_endf(endf_path)
+    endf_product = endf_data[5].products[0]
+    endf_distribution = endf_product.distribution[0]
+
+    # Validity check
+    assert isinstance(endf_distribution, KalbachMann)
+    assert isinstance(hdf5_distribution, KalbachMann)
+    assert endf_product.particle == hdf5_product.particle
+    assert len(endf_distribution.slope) == len(hdf5_distribution.slope)
+
+    # Results check
+    for i, hdf5_slope in enumerate(hdf5_distribution.slope):
+
+        assert endf_distribution._calculated_slope[i] is True
+
+        np.testing.assert_array_almost_equal(
+            endf_distribution.slope[i].y,
+            hdf5_slope.y,
+            decimal=6
+        )
+
+
+@needs_njoy
+@pytest.mark.parametrize(
+    "endf_type, endf_filename", [
+        ('neutrons', 'n-008_O_016.endf'),
+        ('neutrons', 'n-020_Ca_046.endf'),
+        ('neutrons', 'n-080_Hg_204.endf')
+    ]
+)
+def test_comparison_slope_njoy(endf_type, endf_filename):
+    """Test the calculation of the Kalbach-Mann slope done by OpenMC
+    by comparing it to an NJOY calculation. The test is based on
+    the first product of MT=5 (neutron). The isotopes tested have
+    been selected because the corresponding products in ENDF/B-VII.1
+    are described using MF=6, LAW=1, LANG=2 (ie. Kalbach-Mann
+    systematics) and the slope is not given explicitly.
+
+    If an error occurs during the "validity check", this means that
+    the nuclear data evaluation has evolved and the distribution might
+    no longer be described using Kalbach-Mann systematics. Another
+    isotope needs to be identified and tested.
+
+    """
+    endf_directory = os.environ['OPENMC_ENDF_DATA']
+    endf_path = os.path.join(endf_directory, endf_type, endf_filename)
+
+    # ENDF data
+    endf_data = IncidentNeutron.from_endf(endf_path)
+    endf_product = endf_data[5].products[0]
+    endf_distribution = endf_product.distribution[0]
+
+    # NJOY data
+    njoy_data = IncidentNeutron.from_njoy(endf_path, heatr=False, gaspr=False,
+                                          purr=False, smoothing=False)
+    njoy_product = njoy_data[5].products[0]
+    njoy_distribution = njoy_product.distribution[0]
+
+    # Validity check
+    assert isinstance(endf_distribution, KalbachMann)
+    assert isinstance(njoy_distribution, KalbachMann)
+    assert endf_product.particle == njoy_product.particle
+    assert len(endf_distribution.slope) == len(njoy_distribution.slope)
+
+    # Results check
+    for i, njoy_slope in enumerate(njoy_distribution.slope):
+
+        assert endf_distribution._calculated_slope[i] is True
+
+        np.testing.assert_array_almost_equal(
+            endf_distribution.slope[i].y,
+            njoy_slope.y,
+            decimal=6
+        )

From 8f564052e5e6f818bd4d2b0e9ea20116c497622c Mon Sep 17 00:00:00 2001
From: Joffrey Dorville <joffreydorville@mines.edu>
Date: Sat, 2 Apr 2022 00:05:32 -0600
Subject: [PATCH 02/12] Add the Kalbach-Mann slope documentation

---
 docs/source/pythonapi/data.rst | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst
index 95fdfeca9..471be43d8 100644
--- a/docs/source/pythonapi/data.rst
+++ b/docs/source/pythonapi/data.rst
@@ -24,6 +24,7 @@ and product yields.
     FissionProductYields
     WindowedMultipole
     ProbabilityTables
+    AtomicRepresentation
 
 The following classes are used for storing atomic data (incident photon cross
 sections, atomic relaxation):
@@ -68,6 +69,7 @@ Core Functions
     thin
     water_density
     zam
+    return_kalbach_slope
 
 One-dimensional Functions
 -------------------------

From af579895dd5933dc42c07bbd2c025e233c6abd8d Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 7 Apr 2022 10:55:33 -0500
Subject: [PATCH 03/12] Refactor K-M slope calculation

---
 openmc/data/kalbach_mann.py                | 317 +++++----------------
 tests/unit_tests/test_data_kalbach_mann.py |  63 +---
 2 files changed, 80 insertions(+), 300 deletions(-)

diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py
index 11e3acf40..5c0423f30 100644
--- a/openmc/data/kalbach_mann.py
+++ b/openmc/data/kalbach_mann.py
@@ -13,70 +13,6 @@ from .data import EV_PER_MEV
 from .endf import get_list_record, get_tab2_record
 
 
-# Kalbach-Mann constants as defined in ENDF-6 manual BNL-203218-2018-INRE,
-# Revision 215, File 6 description for LAW=1 and LANG=2.
-_C1 = 0.04  # [1/MeV]
-_C2 = 1.8E-6  # [1/MeV^3]
-_C3 = 6.7E-7  # [1/MeV^4]
-_ET1 = 130.  # [MeV]
-_ET3 = 41.  # [MeV]
-_M_NEUTRON = 1.
-_M_PROTON = 1.
-_M_DEUTERON = 1.
-_M_TRITON = None
-_M_3HE = None
-_M_ALPHA = 0.
-_SM_NEUTRON = 1/2.
-_SM_PROTON = 1.
-_SM_DEUTERON = 1.
-_SM_TRITON = 1.
-_SM_3HE = 1.
-_SM_ALPHA = 2.
-
-# Kalbach-Mann M coefficients
-_TABULATED_PARTICLE_M = {
-    1: _M_NEUTRON,
-    1001: _M_PROTON,
-    1002: _M_DEUTERON,
-    1003: _M_TRITON,
-    2003: _M_3HE,
-    2004: _M_ALPHA
-}
-
-# Kalbach-Mann m coefficients
-_TABULATED_PARTICLE_SM = {
-    1: _SM_NEUTRON,
-    1001: _SM_PROTON,
-    1002: _SM_DEUTERON,
-    1003: _SM_TRITON,
-    2003: _SM_3HE,
-    2004: _SM_ALPHA
-}
-
-# Breaking energy as defined in ENDF-6 manual BNL-203218-2018-INRE,
-# Revision 215, Appendix H, Table 3.
-_BREAKING_ENERGY_NEUTRON = 0.
-_BREAKING_ENERGY_PROTON = 0.
-_BREAKING_ENERGY_DEUTERON = 2.224566  # [MeV]
-_BREAKING_ENERGY_TRITON = 8.481798  # [MeV]
-_BREAKING_ENERGY_3HE = 7.718043  # [MeV]
-_BREAKING_ENERGY_ALPHA = 28.29566  # [MeV]
-
-_TABULATED_BREAKING_ENERGY = {
-    1: _BREAKING_ENERGY_NEUTRON,
-    1001: _BREAKING_ENERGY_PROTON,
-    1002: _BREAKING_ENERGY_DEUTERON,
-    1003: _BREAKING_ENERGY_TRITON,
-    2003: _BREAKING_ENERGY_3HE,
-    2004: _BREAKING_ENERGY_ALPHA
-}
-
-# Abundant IZA translation in merged library
-_IZA_TRANSLATION = {
-    6000: 6012,
-}
-
-
 class AtomicRepresentation(EqualityMixin):
     """Atomic representation of an isotope or a particle.
 
@@ -151,27 +87,6 @@ class AtomicRepresentation(EqualityMixin):
     def n(self):
         return self.a - self.z
 
-    @property
-    def breaking_energy(self):
-        breaking_energy = None
-        if self.iza in _TABULATED_BREAKING_ENERGY:
-            breaking_energy = _TABULATED_BREAKING_ENERGY[self.iza]
-        return breaking_energy
-
-    @property
-    def M(self):
-        M = None
-        if self.iza in _TABULATED_PARTICLE_M:
-            M = _TABULATED_PARTICLE_M[self.iza]
-        return M
-
-    @property
-    def m(self):
-        m = None
-        if self.iza in _TABULATED_PARTICLE_SM:
-            m = _TABULATED_PARTICLE_SM[self.iza]
-        return m
-
     @property
     def iza(self):
         iza = self.z * 1000 + self.a
@@ -234,16 +149,11 @@ class AtomicRepresentation(EqualityMixin):
             Atomic representation of the isotope/particle
 
         """
-        if iza in _IZA_TRANSLATION.keys():
-            iza = _IZA_TRANSLATION[iza]
-
-        z = int(iza/1000)
-        a = np.mod(iza, 1000)
-
+        z, a = divmod(iza, 1000)
         return cls(z, a)
 
 
-def _calculate_separation_energy(compound, nucleus, particle):
+def _separation_energy(compound, nucleus, particle):
     """Calculates the separation energy as defined in ENDF-6 manual
     BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1
     and LANG=2. This function can be used for the incident or emitted
@@ -251,158 +161,56 @@ def _calculate_separation_energy(compound, nucleus, particle):
 
     Parameters
     ----------
-    compound: AtomicRepresentation
+    compound : AtomicRepresentation
         Atomic representation of the compound (C)
-    nucleus: AtomicRepresentation
+    nucleus : AtomicRepresentation
         Atomic representation of the nucleus (A or B)
-    particle: AtomicRepresentation
+    particle : AtomicRepresentation
         Atomic representation of the particle (a or b)
 
     Returns
     -------
-    separation_energy: float
+    separation_energy : float
         Separation energy in MeV
 
     """
-    coef_1 = 15.68 * (compound.a - nucleus.a)
-    coef_2 = 28.07 * ((compound.n - compound.z)**2 / float(compound.a) - \
-        (nucleus.n - nucleus.z)**2 / float(nucleus.a))
-    coef_3 = 18.56 * (compound.a**(2./3.) - nucleus.a**(2./3.))
-    coef_4 = 33.22 * ((compound.n - compound.z)**2 / float(compound.a)**(4./3.) - \
-        (nucleus.n - nucleus.z)**2 / float(nucleus.a)**(4./3.))
-    coef_5 = 0.717 * (compound.z**2 / float(compound.a)**(1./3.) - \
-        nucleus.z**2 / float(nucleus.a)**(1./3.))
-    coef_6 = 1.211 * (compound.z**2 / float(compound.a) - \
-        nucleus.z**2 / float(nucleus.a))
+    # Determine A, Z, and N for compound and nucleus
+    A_c = compound.a
+    Z_c = compound.z
+    N_c = compound.n
+    A_a = nucleus.a
+    Z_a = nucleus.z
+    N_a = nucleus.n
 
-    separation_energy = coef_1 - coef_2 - coef_3 + coef_4 \
-        - coef_5 + coef_6 - particle.breaking_energy
+    # Determine breakup energy of incident particle (ENDF-6 Formats Manual,
+    # Appendix H, Table 3)
+    iza_to_breaking_energy = {
+        1: 0.0,
+        1001: 0.0,
+        1002: 2.224566,
+        1003: 8.481798,
+        2003: 7.718043,
+        2004: 28.29566
+    }
+    I_b = iza_to_breaking_energy[particle.iza]
 
-    return separation_energy
+    # Eq. 4 in in doi:10.1103/PhysRevC.37.2350 or ENDF-6 Formats Manual section
+    # 6.2.3.2
+    return (
+        15.68 * (A_c - A_a) -
+        28.07 * ((N_c - Z_c)**2 / A_c - (N_a - Z_a)**2 / A_a) -
+        18.56 * (A_c**(2./3.) - A_a**(2./3.)) +
+        33.22 * ((N_c - Z_c)**2 / A_c**(4./3.) - (N_a - Z_a)**2 / A_a**(4./3.)) -
+        0.717 * (Z_c**2 / A_c**(1./3.) - Z_a**2 / A_a**(1./3.)) +
+        1.211 * (Z_c**2 / A_c - Z_a**2 / A_a) -
+        I_b
+    )
 
 
-def _return_entrance_channel_energy(e_p, awr_t, awr_p):
-    """Returns the entrance channel energy as defined in ENDF-6 manual
-    BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1
-    and LANG=2.
-
-    Parameters
-    ----------
-    e_p: float
-        Energy of the incident projectile in the laboratory system in eV
-    awr_t: float
-        Atomic weight ratio of the target
-    awr_p: float
-        Atomic weight ratio of the projectile
-
-    Returns
-    -------
-    epsilon_p: float
-        Entrance channel energy in eV
-
-    """
-    epsilon_p = e_p * awr_t / (awr_t + awr_p)
-    return epsilon_p
-
-
-def _return_emission_channel_energy(e_e, awr_r, awr_e):
-    """Returns the emission channel energy as defined in ENDF-6 manual
-    BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1
-    and LANG=2.
-
-    Parameters
-    ----------
-    e_e: float
-        Energy of the emitted particle in the center of mass system in eV
-    awr_r: float
-        Atomic weight ratio of the residual nucleus
-    awr_e: float
-        Atomic weight ratio of the emitted particle
-
-    Returns
-    -------
-    epsilon_e: float
-        Emission channel energy in eV
-
-    """
-    epsilon_e = e_e * (awr_r + awr_e) / awr_r
-    return epsilon_e
-
-
-def _calculate_kalbach_slope(energy_projectile,
-                             energy_emitted,
-                             projectile,
-                             target,
-                             compound,
-                             emitted,
-                             residual):
-    """Calculate the Kalbach slope for projectiles other than photons
-    as defined in ENDF-6 manual BNL-203218-2018-INRE, Revision 215,
-    File 6 description for LAW=1 and LANG=2.
-
-    The entrance and emission channel energies are not calculated with
-    the AWR number, but approximated with the number of mass instead.
-
-    Parameters
-    ----------
-    energy_projectile: float
-        Energy of the projectile in the laboratory system in eV
-    energy_emitted: float
-        Energy of the emitted particle in the center of mass system in eV
-    projectile: AtomicRepresentation
-        Atomic representation of the projectile
-    target: AtomicRepresentation
-        Atomic representation of the target
-    compound: AtomicRepresentation
-        Atomic representation of the compound
-    emitted: AtomicRepresentation
-        Atomic representation of the emitted particle
-    residual: AtomicRepresentation
-        Atomic representation of the residual nucleus
-
-    Returns
-    -------
-    slope: float
-        Kalbach-Mann slope
-
-    """
-    epsilon_a = _return_entrance_channel_energy(
-        energy_projectile,
-        target.a,
-        projectile.a
-    ) / EV_PER_MEV
-    epsilon_b = _return_emission_channel_energy(
-        energy_emitted,
-        residual.a,
-        emitted.a
-    ) / EV_PER_MEV
-
-    s_a = _calculate_separation_energy(compound, target, projectile)
-    s_b = _calculate_separation_energy(compound, residual, emitted)
-
-    e_a = epsilon_a + s_a
-    e_b = epsilon_b + s_b
-
-    r_1 = min(e_a, _ET1)
-    r_3 = min(e_a, _ET3)
-
-    x_1 = r_1 * e_b / e_a
-    x_3 = r_3 * e_b / e_a
-
-    slope = _C1 * x_1 \
-        + _C2 * x_1**3 \
-        + _C3 * projectile.M * emitted.m * x_3**4
-
-    return slope
-
-
-def return_kalbach_slope(energy_projectile,
-                         energy_emitted,
-                         iza_projectile,
-                         iza_emitted,
-                         iza_target):
+def kalbach_slope(energy_projectile, energy_emitted, iza_projectile,
+                  iza_emitted, iza_target):
     """Returns Kalbach-Mann slope from calculations.
-    
+
     The associated reaction is defined as:
     A + a -> C -> B + b
 
@@ -456,21 +264,42 @@ def return_kalbach_slope(energy_projectile,
             "Developed and tested for neutron projectile only."
         )
 
+    # Special handling of elemental carbon
+    if iza_emitted == 6000:
+        iza_emitted = 6012
+    if iza_target == 6000:
+        iza_target = 6012
+
     projectile = AtomicRepresentation.from_iza(iza_projectile)
     emitted = AtomicRepresentation.from_iza(iza_emitted)
     target = AtomicRepresentation.from_iza(iza_target)
     compound = projectile + target
     residual = compound - emitted
 
-    slope = _calculate_kalbach_slope(
-        energy_projectile,
-        energy_emitted,
-        projectile,
-        target,
-        compound,
-        emitted,
-        residual
-    )
+    # Calculate entrance and emission channel energy in MeV, defined in section
+    # 6.2.3.2 in the ENDF-6 Formats Manual
+    epsilon_a = energy_projectile * target.a / (target.a + projectile.a) / EV_PER_MEV
+    epsilon_b = energy_emitted * (residual.a + emitted.a) \
+        / (residual.a * EV_PER_MEV)
+
+    # Calculate separation energies using Eq. 4 in doi:10.1103/PhysRevC.37.2350
+    # or ENDF-6 Formats Manual section 6.2.3.2
+    s_a = _separation_energy(compound, target, projectile)
+    s_b = _separation_energy(compound, residual, emitted)
+
+    # See Eq. 10 in doi:10.1103/PhysRevC.37.2350 or section 6.2.3.2 in the
+    # ENDF-6 Formats Manual
+    iza_to_M = {1: 1.0, 1001: 1.0, 1002: 1.0, 2004: 0.0}
+    iza_to_m = {1: 0.5, 1001: 1.0, 1002: 1.0, 1003: 1.0, 2003: 1.0, 2004: 2.0}
+    M = iza_to_M[projectile.iza]
+    m = iza_to_m[emitted.iza]
+    e_a = epsilon_a + s_a
+    e_b = epsilon_b + s_b
+    r_1 = min(e_a, 130.)
+    r_3 = min(e_a, 41.)
+    x_1 = r_1 * e_b / e_a
+    x_3 = r_3 * e_b / e_a
+    slope = 0.04 * x_1 + 1.8e-6 * x_1**3 + 6.7e-7 * M * m * x_3**4
 
     return float("%7e" % slope)
 
@@ -889,11 +718,11 @@ class KalbachMann(AngleEnergy):
                 else:
                     # TODO: retrieve IZA of the projectile
                     iza_projectile = 1
-                    a_i = [return_kalbach_slope(energy_projectile=energy[i],
-                                                energy_emitted=e,
-                                                iza_projectile=iza_projectile,
-                                                iza_emitted=iza_emitted,
-                                                iza_target=iza_target)
+                    a_i = [kalbach_slope(energy_projectile=energy[i],
+                                         energy_emitted=e,
+                                         iza_projectile=iza_projectile,
+                                         iza_emitted=iza_emitted,
+                                         iza_target=iza_target)
                            for e in eout_i]
                     calculated_slope.append(True)
 
diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py
index 99869808d..643fffcfc 100644
--- a/tests/unit_tests/test_data_kalbach_mann.py
+++ b/tests/unit_tests/test_data_kalbach_mann.py
@@ -7,14 +7,8 @@ import numpy as np
 
 from openmc.data import IncidentNeutron
 from openmc.data import AtomicRepresentation
-from openmc.data.kalbach_mann import _BREAKING_ENERGY_TRITON
-from openmc.data.kalbach_mann import _M_TRITON
-from openmc.data.kalbach_mann import _SM_TRITON
-from openmc.data.kalbach_mann import _calculate_separation_energy
-from openmc.data.kalbach_mann import _return_emission_channel_energy
-from openmc.data.kalbach_mann import _return_entrance_channel_energy
-from openmc.data.kalbach_mann import _calculate_kalbach_slope
-from openmc.data import return_kalbach_slope
+from openmc.data.kalbach_mann import _separation_energy
+from openmc.data import kalbach_slope
 from openmc.data import KalbachMann
 
 from . import needs_njoy
@@ -61,9 +55,6 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
     # Test instanciation from_iza
     assert b10 == AtomicRepresentation.from_iza(5010)
 
-    # Test instanciation from_iza using IZA translation
-    assert c12 == AtomicRepresentation.from_iza(6000)
-
     # Test addition
     assert c13 + b10 == na23
 
@@ -75,18 +66,12 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
     assert c13.a == 13
     assert c13.z == 6
     assert c13.n == 7
-    assert c13.breaking_energy is None
-    assert c13.M is None
-    assert c13.m is None
     assert c13.iza == 6013
 
     # Test properties when information for Kalbach-Mann are given
     assert triton.a == 3
     assert triton.z == 1
     assert triton.n == 2
-    assert triton.breaking_energy == _BREAKING_ENERGY_TRITON
-    assert triton.M == _M_TRITON
-    assert triton.m == _SM_TRITON
     assert triton.iza == 1003
 
     # Test instanciation errors
@@ -102,50 +87,16 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
         neutron - triton
 
 
-def test__calculate_separation_energy(triton, b10, c13):
+def test_separation_energy(triton, b10, c13):
     """Comparison to hand-calculations on a simple example."""
-    assert _calculate_separation_energy(
+    assert _separation_energy(
         compound=c13,
         nucleus=b10,
         particle=triton
     ) == pytest.approx(18.6880713)
 
 
-def test__return_entrance_channel_energy():
-    """Comparison to hand-calculations on a simple example."""
-    assert _return_entrance_channel_energy(
-        e_p=5.2,
-        awr_t=13.7,
-        awr_p=7.2
-    ) == pytest.approx(3.4086124)
-
-
-def test__return_emission_channel_energy():
-    """Comparison to hand-calculations on a simple example."""
-    assert _return_emission_channel_energy(
-        e_e=6.8,
-        awr_r=49.2,
-        awr_e=5.4
-    ) == pytest.approx(7.5463415)
-
-
-def test__calculate_kalbach_slope(neutron, triton, b10, c12, c13):
-    """Comparison to hand-calculations for n + c12 -> c13 -> triton + b10."""
-    energy_projectile = 10.2  # [eV]
-    energy_emitted = 5.4  # [eV]
-
-    assert _calculate_kalbach_slope(
-        energy_projectile=energy_projectile,
-        energy_emitted=energy_emitted,
-        projectile=neutron,
-        target=c12,
-        compound=c13,
-        emitted=triton,
-        residual=b10
-    ) == pytest.approx(0.8409921475)
-
-
-def test_return_kalbach_slope():
+def test_kalbach_slope():
     """Comparison to hand-calculations for n + c12 -> c13 -> triton + b10."""
     energy_projectile = 10.2  # [eV]
     energy_emitted = 5.4  # [eV]
@@ -153,7 +104,7 @@ def test_return_kalbach_slope():
     # Check that NotImplementedError is raised if the projectile is not
     # a neutron
     with pytest.raises(NotImplementedError):
-        return_kalbach_slope(
+        kalbach_slope(
             energy_projectile=energy_projectile,
             energy_emitted=energy_emitted,
             iza_projectile=1000,
@@ -161,7 +112,7 @@ def test_return_kalbach_slope():
             iza_target=6012
         )
 
-    assert return_kalbach_slope(
+    assert kalbach_slope(
         energy_projectile=energy_projectile,
         energy_emitted=energy_emitted,
         iza_projectile=1,

From e7595c32015e73671a96a3b2cec6ab9b8ca1c491 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 7 Apr 2022 10:58:15 -0500
Subject: [PATCH 04/12] Don't throw away precision on K-M slope calculation

---
 openmc/data/kalbach_mann.py                | 4 +---
 tests/unit_tests/test_data_kalbach_mann.py | 4 ++--
 2 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py
index 5c0423f30..932deea87 100644
--- a/openmc/data/kalbach_mann.py
+++ b/openmc/data/kalbach_mann.py
@@ -299,9 +299,7 @@ def kalbach_slope(energy_projectile, energy_emitted, iza_projectile,
     r_3 = min(e_a, 41.)
     x_1 = r_1 * e_b / e_a
     x_3 = r_3 * e_b / e_a
-    slope = 0.04 * x_1 + 1.8e-6 * x_1**3 + 6.7e-7 * M * m * x_3**4
-
-    return float("%7e" % slope)
+    return 0.04 * x_1 + 1.8e-6 * x_1**3 + 6.7e-7 * M * m * x_3**4
 
 
 class KalbachMann(AngleEnergy):
diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py
index 643fffcfc..27890dacd 100644
--- a/tests/unit_tests/test_data_kalbach_mann.py
+++ b/tests/unit_tests/test_data_kalbach_mann.py
@@ -173,7 +173,7 @@ def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
         np.testing.assert_array_almost_equal(
             endf_distribution.slope[i].y,
             hdf5_slope.y,
-            decimal=6
+            decimal=5
         )
 
 
@@ -227,5 +227,5 @@ def test_comparison_slope_njoy(endf_type, endf_filename):
         np.testing.assert_array_almost_equal(
             endf_distribution.slope[i].y,
             njoy_slope.y,
-            decimal=6
+            decimal=5
         )

From a8fd1ddf7144b657d7050eba01a4bd4255aed406 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 7 Apr 2022 11:11:51 -0500
Subject: [PATCH 05/12] Change iza --> za

---
 openmc/data/kalbach_mann.py                | 119 +++++++++------------
 tests/unit_tests/test_data_kalbach_mann.py |  20 ++--
 2 files changed, 61 insertions(+), 78 deletions(-)

diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py
index 932deea87..8d1a5c2e1 100644
--- a/openmc/data/kalbach_mann.py
+++ b/openmc/data/kalbach_mann.py
@@ -18,36 +18,28 @@ class AtomicRepresentation(EqualityMixin):
 
     Parameters
     ----------
-    z: int
+    z : int
         Number of protons (atomic number)
-    a: int
+    a : int
         Number of nucleons (mass number)
 
     Raises
     ------
-    IOError:
+    IOError
         When the number of protons (z) declared is higher than the number
         of nucleons (a)
 
     Attributes
     ----------
-    z: int
+    z : int
         Number of protons (atomic number)
-    a: int
+    a : int
         Number of nucleons (mass number)
-    n: int
+    n : int
         Number of neutrons
-    breaking_energy: float
-        Energy required to break the isotope or particle into their
-        constituent nucleons from tabulated values
-    M: float
-        Kalbach-Mann M coefficient
-    m: float
-        Kalbach-Mann m coefficient
-    iza: int
-        ZA identifier defined as:
-        iza = Z x 1000 + A,
-        where Z is the number of protons and A the number of nucleons
+    za : int
+        ZA identifier, 1000*Z + A, where Z is the atomic number and A the mass
+        number
 
     """
     def __init__(self, z, a):
@@ -88,9 +80,8 @@ class AtomicRepresentation(EqualityMixin):
         return self.a - self.z
 
     @property
-    def iza(self):
-        iza = self.z * 1000 + self.a
-        return iza
+    def za(self):
+        return self.z * 1000 + self.a
 
     @a.setter
     def a(self, an):
@@ -112,9 +103,9 @@ class AtomicRepresentation(EqualityMixin):
 
         Parameters
         ----------
-        z: int
+        z : int
             Number of protons (atomic number)
-        a: int
+        a : int
             Number of nucleons (mass number)
 
         Raises
@@ -131,17 +122,14 @@ class AtomicRepresentation(EqualityMixin):
             )
 
     @classmethod
-    def from_iza(cls, iza):
+    def from_za(cls, za):
         """Instantiates an AtomicRepresentation from a ZA identifier.
 
-        The ZA identifier is defined as:
-        iza = Z x 1000 + A,
-        where Z is the number of protons and A the number of nucleons.
-
         Parameters
         ----------
-        iza: int
-            ZA identifier
+        za : int
+            ZA identifier, 1000*Z + A, where Z is the atomic number and A the
+            mass number
 
         Returns
         -------
@@ -149,7 +137,7 @@ class AtomicRepresentation(EqualityMixin):
             Atomic representation of the isotope/particle
 
         """
-        z, a = divmod(iza, 1000)
+        z, a = divmod(za, 1000)
         return cls(z, a)
 
 
@@ -184,7 +172,7 @@ def _separation_energy(compound, nucleus, particle):
 
     # Determine breakup energy of incident particle (ENDF-6 Formats Manual,
     # Appendix H, Table 3)
-    iza_to_breaking_energy = {
+    za_to_breaking_energy = {
         1: 0.0,
         1001: 0.0,
         1002: 2.224566,
@@ -192,7 +180,7 @@ def _separation_energy(compound, nucleus, particle):
         2003: 7.718043,
         2004: 28.29566
     }
-    I_b = iza_to_breaking_energy[particle.iza]
+    I_b = za_to_breaking_energy[particle.za]
 
     # Eq. 4 in in doi:10.1103/PhysRevC.37.2350 or ENDF-6 Formats Manual section
     # 6.2.3.2
@@ -207,8 +195,8 @@ def _separation_energy(compound, nucleus, particle):
     )
 
 
-def kalbach_slope(energy_projectile, energy_emitted, iza_projectile,
-                  iza_emitted, iza_target):
+def kalbach_slope(energy_projectile, energy_emitted, za_projectile,
+                  za_emitted, za_target):
     """Returns Kalbach-Mann slope from calculations.
 
     The associated reaction is defined as:
@@ -222,10 +210,6 @@ def kalbach_slope(energy_projectile, energy_emitted, iza_projectile,
     - B is the residual nucleus,
     - b is the emitted particle.
 
-    This function uses the concept of ZA identifier defined as:
-    iza = Z x 1000 + A,
-    where Z is the number of protons and A the number of nucleons.
-
     The Kalbach-Mann slope calculation is done as defined in ENDF-6 manual
     BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 and
     LANG=2. One exception to this, is that the entrance and emission channel
@@ -234,45 +218,44 @@ def kalbach_slope(energy_projectile, energy_emitted, iza_projectile,
 
     Parameters
     ----------
-    energy_projectile: float
+    energy_projectile : float
         Energy of the projectile in the laboratory system in eV
-    energy_emitted: float
+    energy_emitted : float
         Energy of the emitted particle in the center of mass system in eV
-    iza_projectile: int
+    za_projectile : int
         ZA identifier of the projectile
-    iza_emitted: int
+    za_emitted : int
         ZA identifier of the emitted particle
-    iza_target: int
+    za_target : int
         ZA identifier of the targeted nucleus
 
     Raises
     ------
-    NotImplementedError:
-        When the ZA identifier of the projectile is not equal to 1
-        (ie. other than a neutron).
+    NotImplementedError
+        When the projectile is not a neutron
 
     Returns
     -------
-    slope: float
+    slope : float
         Kalbach-Mann slope given with the same format as ACE file.
 
     """
     # TODO: develop for photons as projectile
     # TODO: test for other particles than neutron
-    if iza_projectile != 1:
+    if za_projectile != 1:
         raise NotImplementedError(
             "Developed and tested for neutron projectile only."
         )
 
     # Special handling of elemental carbon
-    if iza_emitted == 6000:
-        iza_emitted = 6012
-    if iza_target == 6000:
-        iza_target = 6012
+    if za_emitted == 6000:
+        za_emitted = 6012
+    if za_target == 6000:
+        za_target = 6012
 
-    projectile = AtomicRepresentation.from_iza(iza_projectile)
-    emitted = AtomicRepresentation.from_iza(iza_emitted)
-    target = AtomicRepresentation.from_iza(iza_target)
+    projectile = AtomicRepresentation.from_za(za_projectile)
+    emitted = AtomicRepresentation.from_za(za_emitted)
+    target = AtomicRepresentation.from_za(za_target)
     compound = projectile + target
     residual = compound - emitted
 
@@ -289,10 +272,10 @@ def kalbach_slope(energy_projectile, energy_emitted, iza_projectile,
 
     # See Eq. 10 in doi:10.1103/PhysRevC.37.2350 or section 6.2.3.2 in the
     # ENDF-6 Formats Manual
-    iza_to_M = {1: 1.0, 1001: 1.0, 1002: 1.0, 2004: 0.0}
-    iza_to_m = {1: 0.5, 1001: 1.0, 1002: 1.0, 1003: 1.0, 2003: 1.0, 2004: 2.0}
-    M = iza_to_M[projectile.iza]
-    m = iza_to_m[emitted.iza]
+    za_to_M = {1: 1.0, 1001: 1.0, 1002: 1.0, 2004: 0.0}
+    za_to_m = {1: 0.5, 1001: 1.0, 1002: 1.0, 1003: 1.0, 2003: 1.0, 2004: 2.0}
+    M = za_to_M[projectile.za]
+    m = za_to_m[emitted.za]
     e_a = epsilon_a + s_a
     e_b = epsilon_b + s_b
     r_1 = min(e_a, 130.)
@@ -642,7 +625,7 @@ class KalbachMann(AngleEnergy):
         return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a)
 
     @classmethod
-    def from_endf(cls, file_obj, iza_emitted, iza_target, projectile_mass):
+    def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass):
         """Generate Kalbach-Mann distribution from an ENDF evaluation.
 
         If the projectile is a neutron, the slope is calculated when it is
@@ -652,11 +635,11 @@ class KalbachMann(AngleEnergy):
         ----------
         file_obj : file-like object
             ENDF file positioned at the start of the Kalbach-Mann distribution
-        iza_emitted : int
+        za_emitted : int
             ZA identifier of the emitted particle
-        iza_target : int
+        za_target : int
             ZA identifier of the target
-        projectile_mass: float
+        projectile_mass : float
             Mass of the projectile
 
         Warns
@@ -714,13 +697,13 @@ class KalbachMann(AngleEnergy):
                     calculated_slope.append(False)
 
                 else:
-                    # TODO: retrieve IZA of the projectile
-                    iza_projectile = 1
+                    # TODO: retrieve ZA of the projectile
+                    za_projectile = 1
                     a_i = [kalbach_slope(energy_projectile=energy[i],
                                          energy_emitted=e,
-                                         iza_projectile=iza_projectile,
-                                         iza_emitted=iza_emitted,
-                                         iza_target=iza_target)
+                                         za_projectile=za_projectile,
+                                         za_emitted=za_emitted,
+                                         za_target=za_target)
                            for e in eout_i]
                     calculated_slope.append(True)
 
diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py
index 27890dacd..36282f7c5 100644
--- a/tests/unit_tests/test_data_kalbach_mann.py
+++ b/tests/unit_tests/test_data_kalbach_mann.py
@@ -52,8 +52,8 @@ def na23():
 
 def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
     """Test the AtomicRepresentation class."""
-    # Test instanciation from_iza
-    assert b10 == AtomicRepresentation.from_iza(5010)
+    # Test instantiation from_za
+    assert b10 == AtomicRepresentation.from_za(5010)
 
     # Test addition
     assert c13 + b10 == na23
@@ -66,13 +66,13 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
     assert c13.a == 13
     assert c13.z == 6
     assert c13.n == 7
-    assert c13.iza == 6013
+    assert c13.za == 6013
 
     # Test properties when information for Kalbach-Mann are given
     assert triton.a == 3
     assert triton.z == 1
     assert triton.n == 2
-    assert triton.iza == 1003
+    assert triton.za == 1003
 
     # Test instanciation errors
     with pytest.raises(IOError):
@@ -107,17 +107,17 @@ def test_kalbach_slope():
         kalbach_slope(
             energy_projectile=energy_projectile,
             energy_emitted=energy_emitted,
-            iza_projectile=1000,
-            iza_emitted=1,
-            iza_target=6012
+            za_projectile=1000,
+            za_emitted=1,
+            za_target=6012
         )
 
     assert kalbach_slope(
         energy_projectile=energy_projectile,
         energy_emitted=energy_emitted,
-        iza_projectile=1,
-        iza_emitted=1003,
-        iza_target=6012
+        za_projectile=1,
+        za_emitted=1003,
+        za_target=6012
     ) == pytest.approx(0.8409921475)
 
 

From 751d09836c5d0710552073f265694e14c8599a32 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 7 Apr 2022 13:20:56 -0500
Subject: [PATCH 06/12] Make 'a' and 'z' properties read-only

---
 openmc/data/kalbach_mann.py                | 52 ++++------------------
 tests/unit_tests/test_data_kalbach_mann.py |  8 ++--
 2 files changed, 13 insertions(+), 47 deletions(-)

diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py
index 8d1a5c2e1..4c6b89e41 100644
--- a/openmc/data/kalbach_mann.py
+++ b/openmc/data/kalbach_mann.py
@@ -43,11 +43,17 @@ class AtomicRepresentation(EqualityMixin):
 
     """
     def __init__(self, z, a):
-        self._consistency_check(z, a)
+        # Sanity checks on values
+        cv.check_type('z', z, Integral)
+        cv.check_greater_than('z', z, 0, equality=True)
+        cv.check_type('a', a, Integral)
+        cv.check_greater_than('a', a, 0, equality=True)
+        if z > a:
+            raise ValueError(f"Number of protons ({z}) must be less than or "
+                             f"equal to number of nucleons ({a}).")
+
         self._z = z
         self._a = a
-        self.z = self._z
-        self.a = self._a
 
     def __add__(self, other):
         """Adds two AtomicRepresentations.
@@ -67,12 +73,10 @@ class AtomicRepresentation(EqualityMixin):
 
     @property
     def a(self):
-        self._consistency_check(self._z, self._a)
         return self._a
 
     @property
     def z(self):
-        self._consistency_check(self._z, self._a)
         return self._z
 
     @property
@@ -83,44 +87,6 @@ class AtomicRepresentation(EqualityMixin):
     def za(self):
         return self.z * 1000 + self.a
 
-    @a.setter
-    def a(self, an):
-        cv.check_type('a', an, Integral)
-        cv.check_greater_than('a', an, 0, equality=True)
-        self._consistency_check(self._z, an)
-        self._a = an
-
-    @z.setter
-    def z(self, zn):
-        cv.check_type('z', zn, Integral)
-        cv.check_greater_than('z', zn, 0, equality=True)
-        self._consistency_check(zn, self._a)
-        self._z = zn
-
-    @staticmethod
-    def _consistency_check(z, a):
-        """Simple consistency check.
-
-        Parameters
-        ----------
-        z : int
-            Number of protons (atomic number)
-        a : int
-            Number of nucleons (mass number)
-
-        Raises
-        ------
-        IOError:
-            When the number of protons (z) declared is higher than the number
-            of nucleons (a)
-
-        """
-        if z > a:
-            raise IOError(
-                "Number of protons (%i) incompatible with number of "
-                "nucleons (%i)" % (z, a)
-            )
-
     @classmethod
     def from_za(cls, za):
         """Instantiates an AtomicRepresentation from a ZA identifier.
diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py
index 36282f7c5..bd540addb 100644
--- a/tests/unit_tests/test_data_kalbach_mann.py
+++ b/tests/unit_tests/test_data_kalbach_mann.py
@@ -75,15 +75,15 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
     assert triton.za == 1003
 
     # Test instanciation errors
-    with pytest.raises(IOError):
+    with pytest.raises(ValueError):
         AtomicRepresentation(z=5, a=1)
     with pytest.raises(ValueError):
         AtomicRepresentation(z=-1, a=1)
-    with pytest.raises(IOError):
+    with pytest.raises(ValueError):
         AtomicRepresentation(z=5, a=0)
-    with pytest.raises(IOError):
+    with pytest.raises(ValueError):
         AtomicRepresentation(z=5, a=-2)
-    with pytest.raises(OSError):
+    with pytest.raises(ValueError):
         neutron - triton
 
 

From f6858279123e7ed413df405f22e0cefb1274462c Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 26 Apr 2022 16:41:49 -0500
Subject: [PATCH 07/12] Make AtomicRepresentation private (leading underscore)

---
 docs/source/pythonapi/data.rst             |  3 +--
 openmc/data/kalbach_mann.py                | 26 +++++++++++-----------
 tests/unit_tests/test_data_kalbach_mann.py | 25 ++++++++++-----------
 3 files changed, 26 insertions(+), 28 deletions(-)

diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst
index 471be43d8..a07fa9bc3 100644
--- a/docs/source/pythonapi/data.rst
+++ b/docs/source/pythonapi/data.rst
@@ -24,7 +24,6 @@ and product yields.
     FissionProductYields
     WindowedMultipole
     ProbabilityTables
-    AtomicRepresentation
 
 The following classes are used for storing atomic data (incident photon cross
 sections, atomic relaxation):
@@ -65,11 +64,11 @@ Core Functions
     dose_coefficients
     gnd_name
     isotopes
+    kalbach_slope
     linearize
     thin
     water_density
     zam
-    return_kalbach_slope
 
 One-dimensional Functions
 -------------------------
diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py
index 4c6b89e41..10c002234 100644
--- a/openmc/data/kalbach_mann.py
+++ b/openmc/data/kalbach_mann.py
@@ -13,7 +13,7 @@ from .data import EV_PER_MEV
 from .endf import get_list_record, get_tab2_record
 
 
-class AtomicRepresentation(EqualityMixin):
+class _AtomicRepresentation(EqualityMixin):
     """Atomic representation of an isotope or a particle.
 
     Parameters
@@ -56,20 +56,20 @@ class AtomicRepresentation(EqualityMixin):
         self._a = a
 
     def __add__(self, other):
-        """Adds two AtomicRepresentations.
+        """Adds two _AtomicRepresentations.
 
         """
         z = self.z + other.z
         a = self.a + other.a
-        return AtomicRepresentation(z=z, a=a)
+        return _AtomicRepresentation(z=z, a=a)
 
     def __sub__(self, other):
-        """Substracts two AtomicRepresentations.
+        """Substracts two _AtomicRepresentations.
 
         """
         z = self.z - other.z
         a = self.a - other.a
-        return AtomicRepresentation(z=z, a=a)
+        return _AtomicRepresentation(z=z, a=a)
 
     @property
     def a(self):
@@ -89,7 +89,7 @@ class AtomicRepresentation(EqualityMixin):
 
     @classmethod
     def from_za(cls, za):
-        """Instantiates an AtomicRepresentation from a ZA identifier.
+        """Instantiate an _AtomicRepresentation from a ZA identifier.
 
         Parameters
         ----------
@@ -99,7 +99,7 @@ class AtomicRepresentation(EqualityMixin):
 
         Returns
         -------
-        AtomicRepresentation
+        _AtomicRepresentation
             Atomic representation of the isotope/particle
 
         """
@@ -115,11 +115,11 @@ def _separation_energy(compound, nucleus, particle):
 
     Parameters
     ----------
-    compound : AtomicRepresentation
+    compound : _AtomicRepresentation
         Atomic representation of the compound (C)
-    nucleus : AtomicRepresentation
+    nucleus : _AtomicRepresentation
         Atomic representation of the nucleus (A or B)
-    particle : AtomicRepresentation
+    particle : _AtomicRepresentation
         Atomic representation of the particle (a or b)
 
     Returns
@@ -219,9 +219,9 @@ def kalbach_slope(energy_projectile, energy_emitted, za_projectile,
     if za_target == 6000:
         za_target = 6012
 
-    projectile = AtomicRepresentation.from_za(za_projectile)
-    emitted = AtomicRepresentation.from_za(za_emitted)
-    target = AtomicRepresentation.from_za(za_target)
+    projectile = _AtomicRepresentation.from_za(za_projectile)
+    emitted = _AtomicRepresentation.from_za(za_emitted)
+    target = _AtomicRepresentation.from_za(za_target)
     compound = projectile + target
     residual = compound - emitted
 
diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py
index bd540addb..2a579b703 100644
--- a/tests/unit_tests/test_data_kalbach_mann.py
+++ b/tests/unit_tests/test_data_kalbach_mann.py
@@ -6,8 +6,7 @@ import pytest
 import numpy as np
 
 from openmc.data import IncidentNeutron
-from openmc.data import AtomicRepresentation
-from openmc.data.kalbach_mann import _separation_energy
+from openmc.data.kalbach_mann import _separation_energy, _AtomicRepresentation
 from openmc.data import kalbach_slope
 from openmc.data import KalbachMann
 
@@ -17,43 +16,43 @@ from . import needs_njoy
 @pytest.fixture(scope='module')
 def neutron():
     """Neutron AtomicRepresentation."""
-    return AtomicRepresentation(z=0, a=1)
+    return _AtomicRepresentation(z=0, a=1)
 
 
 @pytest.fixture(scope='module')
 def triton():
     """Triton AtomicRepresentation."""
-    return AtomicRepresentation(z=1, a=3)
+    return _AtomicRepresentation(z=1, a=3)
 
 
 @pytest.fixture(scope='module')
 def b10():
     """B10 AtomicRepresentation."""
-    return AtomicRepresentation(z=5, a=10)
+    return _AtomicRepresentation(z=5, a=10)
 
 
 @pytest.fixture(scope='module')
 def c12():
     """C12 AtomicRepresentation."""
-    return AtomicRepresentation(z=6, a=12)
+    return _AtomicRepresentation(z=6, a=12)
 
 
 @pytest.fixture(scope='module')
 def c13():
     """C13 AtomicRepresentation."""
-    return AtomicRepresentation(z=6, a=13)
+    return _AtomicRepresentation(z=6, a=13)
 
 
 @pytest.fixture(scope='module')
 def na23():
     """Na23 AtomicRepresentation."""
-    return AtomicRepresentation(z=11, a=23)
+    return _AtomicRepresentation(z=11, a=23)
 
 
 def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
     """Test the AtomicRepresentation class."""
     # Test instantiation from_za
-    assert b10 == AtomicRepresentation.from_za(5010)
+    assert b10 == _AtomicRepresentation.from_za(5010)
 
     # Test addition
     assert c13 + b10 == na23
@@ -76,13 +75,13 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
 
     # Test instanciation errors
     with pytest.raises(ValueError):
-        AtomicRepresentation(z=5, a=1)
+        _AtomicRepresentation(z=5, a=1)
     with pytest.raises(ValueError):
-        AtomicRepresentation(z=-1, a=1)
+        _AtomicRepresentation(z=-1, a=1)
     with pytest.raises(ValueError):
-        AtomicRepresentation(z=5, a=0)
+        _AtomicRepresentation(z=5, a=0)
     with pytest.raises(ValueError):
-        AtomicRepresentation(z=5, a=-2)
+        _AtomicRepresentation(z=5, a=-2)
     with pytest.raises(ValueError):
         neutron - triton
 

From 7ca7b1c7cffd2c7469fc2c869857f6c14bcc133a Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 26 Apr 2022 16:53:17 -0500
Subject: [PATCH 08/12] Fix some docstrings

---
 openmc/data/kalbach_mann.py | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py
index 10c002234..02403b3cc 100644
--- a/openmc/data/kalbach_mann.py
+++ b/openmc/data/kalbach_mann.py
@@ -25,7 +25,7 @@ class _AtomicRepresentation(EqualityMixin):
 
     Raises
     ------
-    IOError
+    ValueError
         When the number of protons (z) declared is higher than the number
         of nucleons (a)
 
@@ -56,17 +56,13 @@ class _AtomicRepresentation(EqualityMixin):
         self._a = a
 
     def __add__(self, other):
-        """Adds two _AtomicRepresentations.
-
-        """
+        """Add two _AtomicRepresentations"""
         z = self.z + other.z
         a = self.a + other.a
         return _AtomicRepresentation(z=z, a=a)
 
     def __sub__(self, other):
-        """Substracts two _AtomicRepresentations.
-
-        """
+        """Substract two _AtomicRepresentations"""
         z = self.z - other.z
         a = self.a - other.a
         return _AtomicRepresentation(z=z, a=a)

From cdc1bd02179f9df7379225cf80af801ed6d5c679 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 26 Apr 2022 17:00:11 -0500
Subject: [PATCH 09/12] Changes in njoy.py

---
 openmc/data/njoy.py | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py
index fbb26be05..70351a772 100644
--- a/openmc/data/njoy.py
+++ b/openmc/data/njoy.py
@@ -300,7 +300,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
         If the ENDF file contains multiple material evaluations, this argument
         indicates which evaluation should be used.
     smoothing : bool, optional
-        If the smoothing option in ACER is on (1) or off (0) in the card 6.
+        If the smoothing option (ACER card 6) is on (True) or off (False).
     **kwargs
         Keyword arguments passed to :func:`openmc.data.njoy.run`
 
@@ -383,10 +383,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
 
     # acer
     if acer:
-        if smoothing is True:
-            ismoothing = 1
-        else:
-            ismoothing = 0
+        ismoothing = int(smoothing)
         nacer_in = nlast
         for i, temperature in enumerate(temperatures):
             # Extend input with an ACER run for each temperature

From f4d0440d4f72e911ee59b82665de88336fef07f0 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 26 Apr 2022 17:07:25 -0500
Subject: [PATCH 10/12] Remove KM slope NJOY test, other small changes

---
 tests/unit_tests/test_data_kalbach_mann.py | 80 ++++------------------
 1 file changed, 13 insertions(+), 67 deletions(-)

diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py
index 2a579b703..931e236bd 100644
--- a/tests/unit_tests/test_data_kalbach_mann.py
+++ b/tests/unit_tests/test_data_kalbach_mann.py
@@ -2,7 +2,9 @@
 retrieved from ENDF files."""
 
 import os
+from pathlib import Path
 import pytest
+
 import numpy as np
 
 from openmc.data import IncidentNeutron
@@ -50,7 +52,7 @@ def na23():
 
 
 def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
-    """Test the AtomicRepresentation class."""
+    """Test the _AtomicRepresentation class."""
     # Test instantiation from_za
     assert b10 == _AtomicRepresentation.from_za(5010)
 
@@ -121,10 +123,10 @@ def test_kalbach_slope():
 
 
 @pytest.mark.parametrize(
-    "hdf5_filename, endf_type, endf_filename", [
-        ('O16.h5', 'neutrons', 'n-008_O_016.endf'),
-        ('Ca46.h5', 'neutrons', 'n-020_Ca_046.endf'),
-        ('Hg204.h5', 'neutrons', 'n-080_Hg_204.endf')
+    "hdf5_filename, endf_filename", [
+        ('O16.h5', 'n-008_O_016.endf'),
+        ('Ca46.h5', 'n-020_Ca_046.endf'),
+        ('Hg204.h5', 'n-080_Hg_204.endf')
     ]
 )
 def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
@@ -132,7 +134,7 @@ def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
     by comparing it to HDF5 data. The test is based on the first product
     of MT=5 (neutron). The isotopes tested have been selected because the
     corresponding products in ENDF/B-VII.1 are described using MF=6, LAW=1,
-    LANG=2 (ie. Kalbach-Mann systematics) and the slope is not given
+    LANG=2 (i.e., Kalbach-Mann systematics) and the slope is not given
     explicitly.
 
     If an error occurs during the "validity check", this means that
@@ -145,15 +147,14 @@ def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
 
     """
     # HDF5 data
-    hdf5_directory = os.path.dirname(os.environ['OPENMC_CROSS_SECTIONS'])
-    hdf5_path = os.path.join(hdf5_directory, hdf5_filename)
-    hdf5_data = IncidentNeutron.from_hdf5(hdf5_path)
+    hdf5_directory = Path(os.environ['OPENMC_CROSS_SECTIONS']).parent
+    hdf5_data = IncidentNeutron.from_hdf5(hdf5_directory / hdf5_filename)
     hdf5_product = hdf5_data[5].products[0]
     hdf5_distribution = hdf5_product.distribution[0]
 
     # ENDF data
-    endf_directory = os.environ['OPENMC_ENDF_DATA']
-    endf_path = os.path.join(endf_directory, endf_type, endf_filename)
+    endf_directory = Path(os.environ['OPENMC_ENDF_DATA'])
+    endf_path = endf_directory / 'neutrons' / endf_filename
     endf_data = IncidentNeutron.from_endf(endf_path)
     endf_product = endf_data[5].products[0]
     endf_distribution = endf_product.distribution[0]
@@ -166,65 +167,10 @@ def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
 
     # Results check
     for i, hdf5_slope in enumerate(hdf5_distribution.slope):
-
-        assert endf_distribution._calculated_slope[i] is True
+        assert endf_distribution._calculated_slope[i]
 
         np.testing.assert_array_almost_equal(
             endf_distribution.slope[i].y,
             hdf5_slope.y,
             decimal=5
         )
-
-
-@needs_njoy
-@pytest.mark.parametrize(
-    "endf_type, endf_filename", [
-        ('neutrons', 'n-008_O_016.endf'),
-        ('neutrons', 'n-020_Ca_046.endf'),
-        ('neutrons', 'n-080_Hg_204.endf')
-    ]
-)
-def test_comparison_slope_njoy(endf_type, endf_filename):
-    """Test the calculation of the Kalbach-Mann slope done by OpenMC
-    by comparing it to an NJOY calculation. The test is based on
-    the first product of MT=5 (neutron). The isotopes tested have
-    been selected because the corresponding products in ENDF/B-VII.1
-    are described using MF=6, LAW=1, LANG=2 (ie. Kalbach-Mann
-    systematics) and the slope is not given explicitly.
-
-    If an error occurs during the "validity check", this means that
-    the nuclear data evaluation has evolved and the distribution might
-    no longer be described using Kalbach-Mann systematics. Another
-    isotope needs to be identified and tested.
-
-    """
-    endf_directory = os.environ['OPENMC_ENDF_DATA']
-    endf_path = os.path.join(endf_directory, endf_type, endf_filename)
-
-    # ENDF data
-    endf_data = IncidentNeutron.from_endf(endf_path)
-    endf_product = endf_data[5].products[0]
-    endf_distribution = endf_product.distribution[0]
-
-    # NJOY data
-    njoy_data = IncidentNeutron.from_njoy(endf_path, heatr=False, gaspr=False,
-                                          purr=False, smoothing=False)
-    njoy_product = njoy_data[5].products[0]
-    njoy_distribution = njoy_product.distribution[0]
-
-    # Validity check
-    assert isinstance(endf_distribution, KalbachMann)
-    assert isinstance(njoy_distribution, KalbachMann)
-    assert endf_product.particle == njoy_product.particle
-    assert len(endf_distribution.slope) == len(njoy_distribution.slope)
-
-    # Results check
-    for i, njoy_slope in enumerate(njoy_distribution.slope):
-
-        assert endf_distribution._calculated_slope[i] is True
-
-        np.testing.assert_array_almost_equal(
-            endf_distribution.slope[i].y,
-            njoy_slope.y,
-            decimal=5
-        )

From 1e849fe16239883531c1bb27ce80d826f4583dad Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 13 Jul 2022 11:09:35 -0500
Subject: [PATCH 11/12] Apply @JoffreyDorville suggestions from code review

Co-authored-by: Joffrey Dorville <54550047+JoffreyDorville@users.noreply.github.com>
---
 openmc/data/kalbach_mann.py                | 6 +++---
 tests/unit_tests/test_data_kalbach_mann.py | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py
index 02403b3cc..a036da694 100644
--- a/openmc/data/kalbach_mann.py
+++ b/openmc/data/kalbach_mann.py
@@ -133,7 +133,7 @@ def _separation_energy(compound, nucleus, particle):
     N_a = nucleus.n
 
     # Determine breakup energy of incident particle (ENDF-6 Formats Manual,
-    # Appendix H, Table 3)
+    # Appendix H, Table 3) in MeV
     za_to_breaking_energy = {
         1: 0.0,
         1001: 0.0,
@@ -142,7 +142,7 @@ def _separation_energy(compound, nucleus, particle):
         2003: 7.718043,
         2004: 28.29566
     }
-    I_b = za_to_breaking_energy[particle.za]
+    I_a = za_to_breaking_energy[particle.za]
 
     # Eq. 4 in in doi:10.1103/PhysRevC.37.2350 or ENDF-6 Formats Manual section
     # 6.2.3.2
@@ -153,7 +153,7 @@ def _separation_energy(compound, nucleus, particle):
         33.22 * ((N_c - Z_c)**2 / A_c**(4./3.) - (N_a - Z_a)**2 / A_a**(4./3.)) -
         0.717 * (Z_c**2 / A_c**(1./3.) - Z_a**2 / A_a**(1./3.)) +
         1.211 * (Z_c**2 / A_c - Z_a**2 / A_a) -
-        I_b
+        I_a
     )
 
 
diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py
index 931e236bd..3b837af7d 100644
--- a/tests/unit_tests/test_data_kalbach_mann.py
+++ b/tests/unit_tests/test_data_kalbach_mann.py
@@ -75,7 +75,7 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23):
     assert triton.n == 2
     assert triton.za == 1003
 
-    # Test instanciation errors
+    # Test instantiation errors
     with pytest.raises(ValueError):
         _AtomicRepresentation(z=5, a=1)
     with pytest.raises(ValueError):
@@ -129,7 +129,7 @@ def test_kalbach_slope():
         ('Hg204.h5', 'n-080_Hg_204.endf')
     ]
 )
-def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
+def test_comparison_slope_hdf5(hdf5_filename, endf_filename):
     """Test the calculation of the Kalbach-Mann slope done by OpenMC
     by comparing it to HDF5 data. The test is based on the first product
     of MT=5 (neutron). The isotopes tested have been selected because the

From fc3d94b781b26e9ad4a885b77f847d329040f856 Mon Sep 17 00:00:00 2001
From: Joffrey Dorville <joffreydorville@mines.edu>
Date: Thu, 4 Aug 2022 11:24:46 -0600
Subject: [PATCH 12/12] Apply @paulromano suggestions

---
 openmc/data/njoy.py     | 4 ++--
 openmc/data/reaction.py | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py
index 70351a772..d8ecaae18 100644
--- a/openmc/data/njoy.py
+++ b/openmc/data/njoy.py
@@ -127,7 +127,7 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
 1 0 1 .{ext} /
 '{library}: {zsymam} at {temperature}'/
 {mat} {temperature}
-1 1 {ismoothing}/
+1 1 {ismooth}/
 /
 """
 
@@ -383,7 +383,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
 
     # acer
     if acer:
-        ismoothing = int(smoothing)
+        ismooth = int(smoothing)
         nacer_in = nlast
         for i, temperature in enumerate(temperatures):
             # Extend input with an ACER run for each temperature
diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py
index 42fc576da..a2431cf1c 100644
--- a/openmc/data/reaction.py
+++ b/openmc/data/reaction.py
@@ -82,7 +82,7 @@ def _get_products(ev, mt):
 
     Raises
     ------
-    IOError:
+    IOError
         When the Kalbach-Mann systematics is used, but the product
         is not defined in the 'center-of-mass' system. The breakup logic
         is not implemented which can lead to this error being raised while