Add nuclear zone data

app/celer-dump-data.cc

@@ -123,8 +123,8 @@ void print_isotopes(std::vector<ImportIsotope>& isotopes)
     cout << R"gfm(
 # Isotopes
 
-| Isotope ID | Name   | Atomic number | Atomic mass number | Binding energy (MeV) | Nuclear mass (MeV) |
-| ---------- | ------ | ------------- | ------------------ | -------------------- | ------------------ |
+| Isotope ID | Name   | Atomic number | Atomic mass number | Binding energy (MeV) | Proton loss energy (MeV) | Neutron loss energy (MeV) | Nuclear mass (MeV) |
+| ---------- | ------ | ------------- | ------------------ | -------------------- | ------------------------ | ------------------------- | ------------------ |
 )gfm";
 
     for (unsigned int isotope_id : range(isotopes.size()))
@@ -137,6 +137,8 @@ void print_isotopes(std::vector<ImportIsotope>& isotopes)
              << setw(13) << isotope.atomic_number << " | "
              << setw(18) << isotope.atomic_mass_number << " | "
              << setw(20) << isotope.binding_energy << " | "
+             << setw(24) << isotope.proton_loss_energy << " | "
+             << setw(25) << isotope.neutron_loss_energy << " | "
              << setw(18) << isotope.nuclear_mass << " |\n";
         // clang-format on
     }

src/celeritas/CMakeLists.txt

@@ -62,6 +62,7 @@ list(APPEND SOURCES
   mat/MaterialParams.cc
   mat/MaterialParamsOutput.cc
   mat/detail/Utils.cc
+  neutron/model/CascadeOptions.cc
   neutron/process/NeutronElasticProcess.cc
   neutron/process/NeutronInelasticProcess.cc
   optical/CerenkovParams.cc
@@ -151,6 +152,7 @@ if(CELERITAS_USE_JSON)
     ext/GeantPhysicsOptionsIO.json.cc
     field/FieldDriverOptionsIO.json.cc
     field/RZMapFieldInputIO.json.cc
+    neutron/model/CascadeOptionsIO.json.cc
     phys/PrimaryGeneratorOptionsIO.json.cc
     user/RootStepWriterIO.json.cc
   )

src/celeritas/ext/GeantImporter.cc

@@ -425,6 +425,29 @@ std::vector<ImportIsotope> import_isotopes()
         isotope.atomic_mass_number = g4isotope.GetN();
         isotope.binding_energy = G4NucleiProperties::GetBindingEnergy(
             isotope.atomic_mass_number, isotope.atomic_number);
+
+        // Binding energy difference for losing a nucleon
+        if (isotope.atomic_mass_number > 1 && isotope.atomic_number > 1
+            && isotope.atomic_mass_number >= isotope.atomic_number)
+        {
+            isotope.proton_loss_energy
+                = G4NucleiProperties::GetBindingEnergy(
+                      isotope.atomic_mass_number, isotope.atomic_number)
+                  - G4NucleiProperties::GetBindingEnergy(
+                      isotope.atomic_mass_number - 1,
+                      isotope.atomic_number - 1);
+            isotope.neutron_loss_energy
+                = G4NucleiProperties::GetBindingEnergy(
+                      isotope.atomic_mass_number, isotope.atomic_number)
+                  - G4NucleiProperties::GetBindingEnergy(
+                      isotope.atomic_mass_number - 1, isotope.atomic_number);
+        }
+        else
+        {
+            isotope.proton_loss_energy = 0;
+            isotope.neutron_loss_energy = 0;
+        }
+
         isotope.nuclear_mass = G4NucleiProperties::GetNuclearMass(
             isotope.atomic_mass_number, isotope.atomic_number);
     }

src/celeritas/io/ImportElement.hh

@@ -22,10 +22,12 @@ namespace celeritas
 struct ImportIsotope
 {
     std::string name;  //!< Isotope label
-    int atomic_number;  //!< Atomic number Z
-    int atomic_mass_number;  //!< Atomic number A
-    double binding_energy;  //!< Nuclear binding energy [MeV]
-    double nuclear_mass;  //!< Sum of nucleons' mass + binding energy [MeV]
+    int atomic_number{0};  //!< Atomic number Z
+    int atomic_mass_number{0};  //!< Atomic number A
+    double binding_energy{0};  //!< Nuclear binding energy (BE) [MeV]
+    double proton_loss_energy{0};  //!< BE(A, Z) - BE(A-1, Z-1) [MeV]
+    double neutron_loss_energy{0};  //!< BE(A, Z) - BE(A-1, Z) [MeV]
+    double nuclear_mass{0};  //!< Sum of nucleons' mass + binding energy [MeV]
 };
 
 //---------------------------------------------------------------------------//

src/celeritas/mat/IsotopeView.hh

@@ -50,6 +50,12 @@ class IsotopeView
     // Nuclear binding energy
     CELER_FORCEINLINE_FUNCTION MevEnergy binding_energy() const;
 
+    // Nuclear binding energy difference for a proton loss
+    CELER_FORCEINLINE_FUNCTION MevEnergy proton_loss_energy() const;
+
+    // Nuclear binding energy difference for a neutron loss
+    CELER_FORCEINLINE_FUNCTION MevEnergy neutron_loss_energy() const;
+
     // Sum of nucleons + binding energy
     CELER_FORCEINLINE_FUNCTION MevMass nuclear_mass() const;
 
@@ -115,6 +121,24 @@ CELER_FUNCTION units::MevEnergy IsotopeView::binding_energy() const
     return isotope_def().binding_energy;
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Nuclear binding energy difference for a proton loss.
+ */
+CELER_FUNCTION units::MevEnergy IsotopeView::proton_loss_energy() const
+{
+    return isotope_def().proton_loss_energy;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Nuclear binding energy difference for a neutron loss.
+ */
+CELER_FUNCTION units::MevEnergy IsotopeView::neutron_loss_energy() const
+{
+    return isotope_def().neutron_loss_energy;
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Nuclear mass, the sum of the nucleons' mass and their binding energy.

src/celeritas/mat/MaterialData.hh

@@ -33,7 +33,9 @@ struct IsotopeRecord
 
     AtomicNumber atomic_number;  //!< Atomic number Z
     AtomicMassNumber atomic_mass_number;  //!< Atomic number A
-    units::MevEnergy binding_energy;  //!< Nuclear binding energy
+    units::MevEnergy binding_energy;  //!< Nuclear binding energy (BE)
+    units::MevEnergy proton_loss_energy;  //!< BE(A, Z) - BE(A-1, Z-1)
+    units::MevEnergy neutron_loss_energy;  //!< BE(A, Z) - BE(A-1, Z)
     units::MevMass nuclear_mass;  //!< Nucleons' mass + binding energy
 };
 

src/celeritas/mat/MaterialParams.cc

@@ -78,6 +78,10 @@ MaterialParams::from_import(ImportData const& data)
             = AtomicNumber{isotope.atomic_mass_number};
         isotope_params.binding_energy
             = units::MevEnergy(isotope.binding_energy);
+        isotope_params.proton_loss_energy
+            = units::MevEnergy(isotope.proton_loss_energy);
+        isotope_params.neutron_loss_energy
+            = units::MevEnergy(isotope.neutron_loss_energy);
         // Convert from MeV (Geant4) to MeV/c^2 (Celeritas)
         isotope_params.nuclear_mass = units::MevMass(isotope.nuclear_mass);
 
@@ -373,6 +377,8 @@ void MaterialParams::append_isotope_def(IsotopeInput const& inp,
     CELER_EXPECT(inp.atomic_number);
     CELER_EXPECT(inp.atomic_mass_number);
     CELER_EXPECT(inp.binding_energy >= zero_quantity());
+    CELER_EXPECT(inp.proton_loss_energy >= zero_quantity());
+    CELER_EXPECT(inp.neutron_loss_energy >= zero_quantity());
     CELER_EXPECT(inp.nuclear_mass > zero_quantity());
 
     IsotopeRecord result;
@@ -381,6 +387,8 @@ void MaterialParams::append_isotope_def(IsotopeInput const& inp,
     result.atomic_number = inp.atomic_number;
     result.atomic_mass_number = inp.atomic_mass_number;
     result.binding_energy = inp.binding_energy;
+    result.proton_loss_energy = inp.proton_loss_energy;
+    result.neutron_loss_energy = inp.neutron_loss_energy;
     result.nuclear_mass = inp.nuclear_mass;
 
     // Add to host vector

src/celeritas/mat/MaterialParams.hh

@@ -53,11 +53,14 @@ class MaterialParams final : public ParamsDataInterface<MaterialParamsData>
         //!@{
         //! \name Type aliases
         using AtomicMassNumber = AtomicNumber;
+        using MevEnergy = units::MevEnergy;
         //!@}
 
         AtomicNumber atomic_number;  //!< Atomic number Z
         AtomicMassNumber atomic_mass_number;  //!< Atomic number A
-        units::MevEnergy binding_energy;  //!< Nuclear binding energy
+        MevEnergy binding_energy;  //!< Nuclear binding energy (BE)
+        MevEnergy proton_loss_energy;  //!< BE(A, Z) - BE(A-1, Z-1)
+        MevEnergy neutron_loss_energy;  //!< BE(A, Z) - BE(A-1, Z)
         units::MevMass nuclear_mass;  //!< Nucleons' mass + binding energy
         Label label;  //!< Isotope name
     };

src/celeritas/mat/MaterialParamsOutput.cc

@@ -54,6 +54,8 @@ void MaterialParamsOutput::output(JsonPimpl* j) const
         auto atomic_number = json::array();
         auto atomic_mass_number = json::array();
         auto binding_energy = json::array();
+        auto proton_loss_energy = json::array();
+        auto neutron_loss_energy = json::array();
         auto nuclear_mass = json::array();
 
         for (auto id : range(IsotopeId{material_->num_isotopes()}))
@@ -64,17 +66,26 @@ void MaterialParamsOutput::output(JsonPimpl* j) const
             atomic_mass_number.push_back(
                 iso_view.atomic_mass_number().unchecked_get());
             binding_energy.push_back(iso_view.binding_energy().value());
+            proton_loss_energy.push_back(iso_view.proton_loss_energy().value());
+            neutron_loss_energy.push_back(
+                iso_view.neutron_loss_energy().value());
             nuclear_mass.push_back(iso_view.nuclear_mass().value());
         }
         obj["isotopes"] = {
             {"label", std::move(label)},
             {"atomic_number", std::move(atomic_number)},
             {"atomic_mass_number", std::move(atomic_mass_number)},
             {"binding_energy", std::move(binding_energy)},
+            {"proton_loss_energy", std::move(proton_loss_energy)},
+            {"neutron_loss_energy", std::move(neutron_loss_energy)},
             {"nuclear_mass", std::move(nuclear_mass)},
         };
         units["binding_energy"]
             = accessor_unit_label<decltype(&IsotopeView::binding_energy)>();
+        units["proton_loss_energy"]
+            = accessor_unit_label<decltype(&IsotopeView::proton_loss_energy)>();
+        units["neutron_loss_energy"]
+            = accessor_unit_label<decltype(&IsotopeView::neutron_loss_energy)>();
         units["nuclear_mass"]
             = accessor_unit_label<decltype(&IsotopeView::nuclear_mass)>();
     }

src/celeritas/neutron/data/NeutronInelasticData.hh

@@ -26,9 +26,11 @@ struct NeutronInelasticScalars
     // Action and particle IDs
     ActionId action_id;
     ParticleId neutron_id;
+    ParticleId proton_id;
 
     // Particle mass * c^2 [MeV]
     units::MevMass neutron_mass;
+    units::MevMass proton_mass;
 
     //! Number of nucleon-nucleon channels
     static CELER_CONSTEXPR_FUNCTION size_type num_channels() { return 3; }
@@ -43,7 +45,9 @@ struct NeutronInelasticScalars
     //! Whether data are assigned
     explicit CELER_FUNCTION operator bool() const
     {
-        return action_id && neutron_id && neutron_mass > zero_quantity();
+        return action_id && neutron_id && proton_id
+               && neutron_mass > zero_quantity()
+               && neutron_mass > zero_quantity();
     }
 };
 
@@ -59,6 +63,69 @@ struct StepanovParameters
     Real3 coeffs;  //!< coefficients of a second order Stepanov's function
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * Components of nuclear zone properties of the Bertini cascade model.
+ */
+struct ZoneComponent
+{
+    using NucleonArray = Array<real_type, 2>;  //!< [proton, neutron]
+
+    real_type radius{};  //!< radius of zones in [femtometer]
+    real_type volume{};  //!< volume of zones in [femtometer^3]
+    NucleonArray density{0, 0};  //!< nucleon densities [1/femtometer^3]
+    NucleonArray fermi_mom{0, 0};  //!< fermi momenta in [MeV/c]
+    NucleonArray potential{0, 0};  //!< nucleon potentials [MeV]
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Data characterizing the nuclear zones.
+ */
+struct NuclearZones
+{
+    ItemRange<ZoneComponent> zones;
+
+    //! Whether all data are assigned and valid
+    explicit CELER_FUNCTION operator bool() const { return !zones.empty(); }
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Device data for nuclear zone properties
+ */
+template<Ownership W, MemSpace M>
+struct NuclearZoneData
+{
+    template<class T>
+    using Items = Collection<T, W, M>;
+    template<class T>
+    using IsotopeItems = Collection<T, W, M, IsotopeId>;
+
+    //// MEMBER DATA ////
+
+    // Nuclear zone data
+    Items<ZoneComponent> components;
+    IsotopeItems<NuclearZones> zones;
+
+    //! Whether the data are assigned
+    explicit CELER_FUNCTION operator bool() const
+    {
+        return !components.empty() && !zones.empty();
+    }
+
+    //! Assign from another set of data
+    template<Ownership W2, MemSpace M2>
+    NuclearZoneData& operator=(NuclearZoneData<W2, M2> const& other)
+    {
+        CELER_EXPECT(other);
+        components = other.components;
+        zones = other.zones;
+
+        return *this;
+    }
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * Device data for creating an interactor.
@@ -90,11 +157,14 @@ struct NeutronInelasticData
     // Backend data
     Items<real_type> reals;
 
+    // Nuclear zone data
+    NuclearZoneData<W, M> nuclear_zones;
+
     //! Whether the data are assigned
     explicit CELER_FUNCTION operator bool() const
     {
         return scalars && !micro_xs.empty() && !nucleon_xs.empty()
-               && !xs_params.empty() && !reals.empty();
+               && !xs_params.empty() && !reals.empty() && nuclear_zones;
     }
 
     //! Assign from another set of data
@@ -107,6 +177,7 @@ struct NeutronInelasticData
         nucleon_xs = other.nucleon_xs;
         xs_params = other.xs_params;
         reals = other.reals;
+        nuclear_zones = other.nuclear_zones;
         return *this;
     }
 };

src/celeritas/neutron/model/CascadeOptions.cc

@@ -0,0 +1,47 @@
+//----------------------------------*-C++-*----------------------------------//
+// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/neutron/model/CascadeOptions.cc
+//---------------------------------------------------------------------------//
+#include "CascadeOptions.hh"
+
+#include "corecel/Assert.hh"
+
+namespace celeritas
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Throw a runtime assertion if any of the input is invalid.
+ */
+void validate_input(CascadeOptions const& opts)
+{
+    CELER_VALIDATE(opts.prob_pion_absorption >= 0,
+                   << "invalid prob_pion_absorption "
+                   << opts.prob_pion_absorption);
+    CELER_VALIDATE(opts.radius_scale > 0,
+                   << "invalid radius_scale " << opts.radius_scale);
+    CELER_VALIDATE(opts.radius_small > 0,
+                   << "invalid radius_small " << opts.radius_small);
+    CELER_VALIDATE(opts.radius_alpha > 0,
+                   << "invalid radius_alpha " << opts.radius_alpha);
+    CELER_VALIDATE(opts.radius_trailing >= 0,
+                   << "invalid radius_trailing " << opts.radius_trailing);
+    CELER_VALIDATE(opts.fermi_scale > 0,
+                   << "invalid fermi_scale " << opts.fermi_scale);
+    CELER_VALIDATE(opts.xsec_scale > 0,
+                   << "invalid xsec_scale " << opts.xsec_scale);
+    CELER_VALIDATE(opts.gamma_qd_scale > 0,
+                   << "invalid gamma_qd_scale " << opts.gamma_qd_scale);
+    CELER_VALIDATE(opts.dp_max_doublet > 0,
+                   << "invalid dp_max_doublet " << opts.dp_max_doublet);
+    CELER_VALIDATE(opts.dp_max_triplet > 0,
+                   << "invalid dp_max_triplet " << opts.dp_max_triplet);
+    CELER_VALIDATE(opts.dp_max_alpha > 0,
+                   << "invalid dp_max_alpha " << opts.dp_max_alpha);
+    CELER_ENSURE(opts);
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace celeritas