Implement D+D and D+He fusion (#3257)

* Implement D+D -> n+He3 fusion

* Fix logic for fusion reaction

* Check products in a different place

* Correct compilation error

* Implement D+D -> T+p cross-section

* Update example

* Use clearer naming convention for fusion types

* Revert changes to example input script

* Add analysis script

* Progress on tests

* Use 2 species in test

* Correct momentum of colliding species

* Update test

* Update test

* Generalize species names in fusion tests

* Update benchmarks

* Correct typo

* Updated scripts

* Update script so that it works for D+T and D+D

* Update CI

* Add benchmark file

* Correct typo

* Fix compilation on GPU

* Update RZ CI test

* Implement Deuterium-Helium reaction

* Apply suggestions from code review

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

8 files changed, 215 insertions, 111 deletions

diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H
index 2333e1f41..09213ba03 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H
@@ -12,9 +12,20 @@
 
 #include "Particles/MultiParticleContainer.H"
 
-enum struct CollisionType { DeuteriumTritiumFusion, ProtonBoronFusion, Undefined };
-
-enum struct NuclearFusionType { DeuteriumTritium, ProtonBoron, Undefined };
+enum struct CollisionType { DeuteriumTritiumToNeutronHeliumFusion,
+                            DeuteriumDeuteriumToProtonTritiumFusion,
+                            DeuteriumDeuteriumToNeutronHeliumFusion,
+                            DeuteriumHeliumToProtonHeliumFusion,
+                            ProtonBoronToAlphasFusion,
+                            Undefined };
+
+enum struct NuclearFusionType {
+                                DeuteriumTritiumToNeutronHelium,
+                                DeuteriumDeuteriumToProtonTritium,
+                                DeuteriumDeuteriumToNeutronHelium,
+                                DeuteriumHeliumToProtonHelium,
+                                ProtonBoronToAlphas,
+                                Undefined };
 
 namespace BinaryCollisionUtils{
 
diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
index 4bbdb4bc4..81acecd3e 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
@@ -25,20 +25,75 @@ namespace BinaryCollisionUtils{
             pp_collision_name.getarr("species", species_names);
             auto& species1 = mypc->GetParticleContainerFromName(species_names[0]);
             auto& species2 = mypc->GetParticleContainerFromName(species_names[1]);
+            amrex::Vector<std::string> product_species_name;
+            pp_collision_name.getarr("product_species", product_species_name);
 
             if ((species1.AmIA<PhysicalSpecies::hydrogen2>() && species2.AmIA<PhysicalSpecies::hydrogen3>())
                 ||
                 (species1.AmIA<PhysicalSpecies::hydrogen3>() && species2.AmIA<PhysicalSpecies::hydrogen2>())
                 )
             {
-                return NuclearFusionType::DeuteriumTritium;
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species_name.size() == 2u,
+                    "ERROR: Deuterium-tritium fusion must contain exactly two product species");
+                auto& product_species1 = mypc->GetParticleContainerFromName(product_species_name[0]);
+                auto& product_species2 = mypc->GetParticleContainerFromName(product_species_name[1]);
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    (product_species1.AmIA<PhysicalSpecies::helium4>() && product_species2.AmIA<PhysicalSpecies::neutron>())
+                    ||
+                    (product_species1.AmIA<PhysicalSpecies::neutron>() && product_species2.AmIA<PhysicalSpecies::helium4>()),
+                    "ERROR: Product species of deuterium-tritium fusion must be of type neutron and helium4");
+                return NuclearFusionType::DeuteriumTritiumToNeutronHelium;
+            }
+            else if (species1.AmIA<PhysicalSpecies::hydrogen2>() && species2.AmIA<PhysicalSpecies::hydrogen2>())
+            {
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species_name.size() == 2u,
+                    "ERROR: Deuterium-deuterium fusion must contain exactly two product species");
+                auto& product_species1 = mypc->GetParticleContainerFromName(product_species_name[0]);
+                auto& product_species2 = mypc->GetParticleContainerFromName(product_species_name[1]);
+                if (
+                    (product_species1.AmIA<PhysicalSpecies::helium3>() && product_species2.AmIA<PhysicalSpecies::neutron>())
+                  ||(product_species1.AmIA<PhysicalSpecies::neutron>() && product_species2.AmIA<PhysicalSpecies::helium3>())){
+                    return NuclearFusionType::DeuteriumDeuteriumToNeutronHelium;
+                } else if (
+                    (product_species1.AmIA<PhysicalSpecies::hydrogen3>() && product_species2.AmIA<PhysicalSpecies::proton>())
+                  ||(product_species1.AmIA<PhysicalSpecies::proton>() && product_species2.AmIA<PhysicalSpecies::hydrogen3>())){
+                    return NuclearFusionType::DeuteriumDeuteriumToProtonTritium;
+                } else {
+                    amrex::Abort("ERROR: Product species of proton-boron fusion must be of type helium3 and neutron, or tritium and proton");
+                }
+            }
+            else if ((species1.AmIA<PhysicalSpecies::hydrogen2>() && species2.AmIA<PhysicalSpecies::helium3>())
+                ||
+                (species1.AmIA<PhysicalSpecies::helium3>() && species2.AmIA<PhysicalSpecies::hydrogen2>())
+                )
+            {
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species_name.size() == 2u,
+                    "ERROR: Deuterium-helium fusion must contain exactly two product species");
+                auto& product_species1 = mypc->GetParticleContainerFromName(product_species_name[0]);
+                auto& product_species2 = mypc->GetParticleContainerFromName(product_species_name[1]);
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    (product_species1.AmIA<PhysicalSpecies::helium4>() && product_species2.AmIA<PhysicalSpecies::proton>())
+                    ||
+                    (product_species1.AmIA<PhysicalSpecies::proton>() && product_species2.AmIA<PhysicalSpecies::helium4>()),
+                    "ERROR: Product species of deuterium-helium fusion must be of type proton and helium4");
+                return NuclearFusionType::DeuteriumHeliumToProtonHelium;
             }
             else if ((species1.AmIA<PhysicalSpecies::proton>() && species2.AmIA<PhysicalSpecies::boron11>())
                 ||
                 (species1.AmIA<PhysicalSpecies::boron11>() && species2.AmIA<PhysicalSpecies::proton>())
                 )
             {
-                return NuclearFusionType::ProtonBoron;
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species_name.size() == 1,
+                    "ERROR: Proton-boron must contain exactly one product species");
+                auto& product_species = mypc->GetParticleContainerFromName(product_species_name[0]);
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species.AmIA<PhysicalSpecies::helium>(),
+                    "ERROR: Product species of proton-boron fusion must be of type alpha");
+                return NuclearFusionType::ProtonBoronToAlphas;
             }
             amrex::Abort("Binary nuclear fusion not implemented between species " +
                          species_names[0] + " of type " + species1.getSpeciesTypeName() +
@@ -63,10 +118,16 @@ namespace BinaryCollisionUtils{
 
     CollisionType nuclear_fusion_type_to_collision_type (const NuclearFusionType fusion_type)
         {
-            if (fusion_type == NuclearFusionType::DeuteriumTritium)
-                return CollisionType::DeuteriumTritiumFusion;
-            if (fusion_type == NuclearFusionType::ProtonBoron)
-                return CollisionType::ProtonBoronFusion;
+            if (fusion_type == NuclearFusionType::DeuteriumTritiumToNeutronHelium)
+                return CollisionType::DeuteriumTritiumToNeutronHeliumFusion;
+            if (fusion_type == NuclearFusionType::DeuteriumDeuteriumToProtonTritium)
+                return CollisionType::DeuteriumDeuteriumToProtonTritiumFusion;
+            if (fusion_type == NuclearFusionType::DeuteriumDeuteriumToNeutronHelium)
+                return CollisionType::DeuteriumDeuteriumToNeutronHeliumFusion;
+            if (fusion_type == NuclearFusionType::DeuteriumHeliumToProtonHelium)
+                return CollisionType::DeuteriumHeliumToProtonHeliumFusion;
+            if (fusion_type == NuclearFusionType::ProtonBoronToAlphas)
+                return CollisionType::ProtonBoronToAlphasFusion;
             amrex::Abort("Invalid nuclear fusion type");
             return CollisionType::Undefined;
         }
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/BoschHaleFusionCrossSection.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/BoschHaleFusionCrossSection.H
new file mode 100644
index 000000000..08af65ca6
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/BoschHaleFusionCrossSection.H
@@ -0,0 +1,111 @@
+/* Copyright 2022 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef BOSCH_HALE_FUSION_CROSS_SECTION_H
+#define BOSCH_HALE_FUSION_CROSS_SECTION_H
+
+#include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
+#include "Utils/WarpXConst.H"
+
+#include <AMReX_REAL.H>
+
+#include <cmath>
+
+/**
+ * \brief Computes the fusion cross section, using the analytical fits given in
+ * H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611
+ *
+ * @param[in] E_kin_star the kinetic energy of the reactants in their center of mass frame, in SI units.
+ * @param[in] fusion_type indicates which fusion reaction to calculate the cross-section for
+ * @return The total cross section in SI units (square meters).
+ */
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+amrex::ParticleReal BoschHaleFusionCrossSection (
+    const amrex::ParticleReal& E_kin_star,
+    const NuclearFusionType& fusion_type,
+    const amrex::ParticleReal& m1,
+    const amrex::ParticleReal& m2 )
+{
+    using namespace amrex::literals;
+
+    constexpr amrex::ParticleReal joule_to_keV = 1.e-3_prt/PhysConst::q_e;
+    const amrex::ParticleReal E_keV = E_kin_star*joule_to_keV;
+
+    // If kinetic energy is 0, return a 0 cross section and avoid later division by 0.
+    if (E_keV == 0._prt) {return 0._prt;}
+
+    // Compute the Gamow constant B_G (in keV^{1/2})
+    // (See Eq. 3 in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
+    const amrex::ParticleReal m_reduced = m1 / (1._prt + m1/m2);
+    // The formula for `B_G` below assumes that both reactants have Z=1
+    // When one of the reactants is helium (Z=2), this formula is corrected further below.
+    amrex::ParticleReal B_G = MathConst::pi * PhysConst::alpha
+         * std::sqrt( 2._prt*m_reduced*PhysConst::c*PhysConst::c * joule_to_keV );
+    if (fusion_type == NuclearFusionType::DeuteriumHeliumToProtonHelium) {
+        // Take into account the fact that Z=2 for one of the reactant (helium) in this case
+        B_G *= 2;
+    }
+
+    // Compute astrophysical_factor
+    // (See Eq. 9 and Table IV in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
+    amrex::ParticleReal A1=0_prt, A2=0_prt, A3=0_prt, A4=0_prt, A5=0_prt, B1=0_prt, B2=0_prt, B3=0_prt, B4=0_prt;
+    if (fusion_type == NuclearFusionType::DeuteriumTritiumToNeutronHelium) {
+        A1 = 6.927e4_prt;
+        A2 = 7.454e8_prt;
+        A3 = 2.050e6_prt;
+        A4 = 5.2002e4_prt;
+        A5 = 0_prt;
+        B1 = 6.38e1_prt;
+        B2 = -9.95e-1_prt;
+        B3 = 6.981e-5_prt;
+        B4 = 1.728e-4_prt;
+    }
+    else if (fusion_type == NuclearFusionType::DeuteriumDeuteriumToProtonTritium) {
+        A1 = 5.5576e4_prt;
+        A2 = 2.1054e2_prt;
+        A3 = -3.2638e-2_prt;
+        A4 = 1.4987e-6_prt;
+        A5 = 1.8181e-10_prt;
+        B1 = 0_prt;
+        B2 = 0_prt;
+        B3 = 0_prt;
+        B4 = 0_prt;
+    }
+    else if (fusion_type == NuclearFusionType::DeuteriumDeuteriumToNeutronHelium) {
+        A1 = 5.3701e4_prt;
+        A2 = 3.3027e2_prt;
+        A3 = -1.2706e-1_prt;
+        A4 = 2.9327e-5_prt;
+        A5 = -2.5151e-9_prt;
+        B1 = 0_prt;
+        B2 = 0_prt;
+        B3 = 0_prt;
+        B4 = 0_prt;
+    }
+    else if (fusion_type == NuclearFusionType::DeuteriumHeliumToProtonHelium) {
+        A1 = 5.7501e6_prt;
+        A2 = 2.5226e3_prt;
+        A3 = 4.5566e1_prt;
+        A4 = 0_prt;
+        A5 = 0_prt;
+        B1 = -3.1995e-3_prt;
+        B2 = -8.5530e-6_prt;
+        B3 = 5.9014e-8_prt;
+        B4 = 0_prt;
+    }
+
+    amrex::ParticleReal astrophysical_factor =
+        (A1 + E_keV*(A2 + E_keV*(A3 + E_keV*(A4 + E_keV*A5)))) /
+        (1_prt + E_keV*(B1 + E_keV*(B2 + E_keV*(B3 + E_keV*B4))));
+
+    // Compute cross-section in SI units
+    // (See Eq. 8 in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
+    constexpr amrex::ParticleReal millibarn_to_sqm = 1.e-31_prt;
+    return millibarn_to_sqm * astrophysical_factor/E_keV * std::exp(-B_G/std::sqrt(E_keV));
+}
+
+#endif // BOSCH_HALE_FUSION_CROSS_SECTION_H
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumFusionCrossSection.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumFusionCrossSection.H
deleted file mode 100644
index 58aaa986b..000000000
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumFusionCrossSection.H
+++ /dev/null
@@ -1,64 +0,0 @@
-/* Copyright 2022 Remi Lehe
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-
-#ifndef DEUTERIUM_TRITIUM_FUSION_CROSS_SECTION_H
-#define DEUTERIUM_TRITIUM_FUSION_CROSS_SECTION_H
-
-#include "Utils/WarpXConst.H"
-
-#include <AMReX_REAL.H>
-
-#include <cmath>
-
-/**
- * \brief Computes the total deuterium-tritium fusion cross section, using
- * the analytical fits given in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611
- *
- * @param[in] E_kin_star the kinetic energy of the deuterium-tritium pair in its center of mass frame, in SI units.
- * @return The total cross section in SI units (square meters).
- */
-AMREX_GPU_HOST_DEVICE AMREX_INLINE
-amrex::ParticleReal DeuteriumTritiumFusionCrossSection (const amrex::ParticleReal& E_kin_star)
-{
-    using namespace amrex::literals;
-
-    constexpr amrex::ParticleReal joule_to_keV = 1.e-3_prt/PhysConst::q_e;
-    const amrex::ParticleReal E_keV = E_kin_star*joule_to_keV;
-
-    // If kinetic energy is 0, return a 0 cross section and avoid later division by 0.
-    if (E_keV == 0._prt) {return 0._prt;}
-
-    // Compute the Gamow constant B_G (in keV^{1/2})
-    // (See Eq. 3 in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
-    constexpr amrex::ParticleReal m_D = 2.01410177812 * PhysConst::m_u;
-    constexpr amrex::ParticleReal m_T = 3.0160492779 * PhysConst::m_u;
-    constexpr amrex::ParticleReal m_reduced = m_D / (1._prt + m_D/m_T);
-    amrex::ParticleReal B_G = MathConst::pi * PhysConst::alpha * 1._prt * 1._prt
-         * std::sqrt( 2._prt*m_reduced*PhysConst::c*PhysConst::c * joule_to_keV );
-
-    // Compute astrophysical_factor
-    // (See Eq. 9 and Table IV in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
-    constexpr amrex::ParticleReal A1 = 6.927e4;
-    constexpr amrex::ParticleReal A2 = 7.454e8;
-    constexpr amrex::ParticleReal A3 = 2.050e6;
-    constexpr amrex::ParticleReal A4 = 5.2002e4;
-    constexpr amrex::ParticleReal B1 = 6.38e1;
-    constexpr amrex::ParticleReal B2 = -9.95e-1;
-    constexpr amrex::ParticleReal B3 = 6.981e-5;
-    constexpr amrex::ParticleReal B4 = 1.728e-4;
-
-    amrex::ParticleReal astrophysical_factor =
-        (A1 + E_keV*(A2 + E_keV*(A3 + E_keV*A4))) /
-        (1_prt + E_keV*(B1 + E_keV*(B2 + E_keV*(B3 + E_keV*B4))));
-
-    // Compute cross-section in SI units
-    // (See Eq. 8 in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
-    constexpr amrex::ParticleReal millibarn_to_sqm = 1.e-31_prt;
-    return millibarn_to_sqm * astrophysical_factor/E_keV * std::exp(-B_G/std::sqrt(E_keV));
-}
-
-#endif // DEUTERIUM_TRITIUM_TRITIUM_FUSION_CROSS_SECTION_H
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
index 25624469d..2dcbe22e9 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
@@ -65,33 +65,6 @@ public:
         m_fusion_type = BinaryCollisionUtils::get_nuclear_fusion_type(collision_name, mypc);
 
         amrex::ParmParse pp_collision_name(collision_name);
-        amrex::Vector<std::string> product_species_name;
-        pp_collision_name.getarr("product_species", product_species_name);
-
-        if (m_fusion_type == NuclearFusionType::DeuteriumTritium)
-        {
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                product_species_name.size() == 2u,
-                "ERROR: Deuterium-tritium fusion must contain exactly two product species");
-            auto& product_species1 = mypc->GetParticleContainerFromName(product_species_name[0]);
-            auto& product_species2 = mypc->GetParticleContainerFromName(product_species_name[1]);
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                (product_species1.AmIA<PhysicalSpecies::helium4>() && product_species2.AmIA<PhysicalSpecies::neutron>())
-                ||
-                (product_species1.AmIA<PhysicalSpecies::neutron>() && product_species2.AmIA<PhysicalSpecies::helium4>()),
-                "ERROR: Product species of deuterium-tritium fusion must be of type neutron and helium4");
-        }
-        if (m_fusion_type == NuclearFusionType::ProtonBoron)
-        {
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                product_species_name.size() == 1,
-                "ERROR: Proton-boron must contain exactly one product species");
-            auto& product_species = mypc->GetParticleContainerFromName(product_species_name[0]);
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                product_species.AmIA<PhysicalSpecies::helium>(),
-                "ERROR: Product species of proton-boron fusion must be of type alpha");
-        }
-
         // default fusion multiplier
         m_fusion_multiplier = 1.0_prt;
         utils::parser::queryWithParser(
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
index 0bab4bef4..e4704f746 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
@@ -8,7 +8,7 @@
 #ifndef SINGLE_NUCLEAR_FUSION_EVENT_H_
 #define SINGLE_NUCLEAR_FUSION_EVENT_H_
 
-#include "DeuteriumTritiumFusionCrossSection.H"
+#include "BoschHaleFusionCrossSection.H"
 #include "ProtonBoronFusionCrossSection.H"
 
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
@@ -112,13 +112,15 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part
 
     // Compute fusion cross section as a function of kinetic energy in the center of mass frame
     auto fusion_cross_section = amrex::ParticleReal(0.);
-    if (fusion_type == NuclearFusionType::DeuteriumTritium)
+    if (fusion_type == NuclearFusionType::ProtonBoronToAlphas)
     {
-        fusion_cross_section = DeuteriumTritiumFusionCrossSection(E_kin_star);
+        fusion_cross_section = ProtonBoronFusionCrossSection(E_kin_star);
     }
-    else if (fusion_type == NuclearFusionType::ProtonBoron)
+    else if ((fusion_type == NuclearFusionType::DeuteriumTritiumToNeutronHelium)
+          || (fusion_type == NuclearFusionType::DeuteriumDeuteriumToProtonTritium)
+          || (fusion_type == NuclearFusionType::DeuteriumDeuteriumToNeutronHelium))
     {
-        fusion_cross_section = ProtonBoronFusionCrossSection(E_kin_star);
+        fusion_cross_section = BoschHaleFusionCrossSection(E_kin_star, fusion_type, m1, m2);
     }
 
     // Square of the norm of the momentum of one of the particles in the center of mass frame
diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
index be77ae9e0..085adcc88 100644
--- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
@@ -216,7 +216,7 @@ public:
 
                 // Initialize the product particles' momentum, using a function depending on the
                 // specific collision type
-                if (t_collision_type == CollisionType::ProtonBoronFusion)
+                if (t_collision_type == CollisionType::ProtonBoronToAlphasFusion)
                 {
                     const index_type product_start_index = products_np_data[0] + 2*p_offsets[i]*
                                                            p_num_products_device[0];
@@ -224,11 +224,19 @@ public:
                                                         p_pair_indices_1[i], p_pair_indices_2[i],
                                                         product_start_index, m1, m2, engine);
                 }
-                else if (t_collision_type == CollisionType::DeuteriumTritiumFusion)
+                else if ((t_collision_type == CollisionType::DeuteriumTritiumToNeutronHeliumFusion)
+                      || (t_collision_type == CollisionType::DeuteriumDeuteriumToProtonTritiumFusion)
+                      || (t_collision_type == CollisionType::DeuteriumDeuteriumToNeutronHeliumFusion))
                 {
                     amrex::ParticleReal fusion_energy = 0.0_prt;
-                    if (t_collision_type == CollisionType::DeuteriumTritiumFusion) {
-                        fusion_energy = 17.5893e6_prt * PhysConst::q_e; // 17.6 MeV
+                    if (t_collision_type == CollisionType::DeuteriumTritiumToNeutronHeliumFusion) {
+                        fusion_energy = 17.5893e6_prt * PhysConst::q_e;
+                    }
+                    else if (t_collision_type == CollisionType::DeuteriumDeuteriumToProtonTritiumFusion) {
+                        fusion_energy = 4.032667e6_prt * PhysConst::q_e;
+                    }
+                    else if (t_collision_type == CollisionType::DeuteriumDeuteriumToNeutronHeliumFusion) {
+                        fusion_energy = 3.268911e6_prt * PhysConst::q_e;
                     }
                     TwoProductFusionInitializeMomentum(soa_1, soa_2,
                         soa_products_data[0], soa_products_data[1],
diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp
index bab92e486..3b1951107 100644
--- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp
@@ -23,7 +23,7 @@ ParticleCreationFunc::ParticleCreationFunc (const std::string collision_name,
 
         m_collision_type = BinaryCollisionUtils::get_collision_type(collision_name, mypc);
 
-        if (m_collision_type == CollisionType::ProtonBoronFusion)
+        if (m_collision_type == CollisionType::ProtonBoronToAlphasFusion)
             {
                 // Proton-Boron fusion only produces alpha particles
                 m_num_product_species = 1;
@@ -34,7 +34,9 @@ ParticleCreationFunc::ParticleCreationFunc (const std::string collision_name,
                 m_num_products_device.push_back(3);
 #endif
             }
-        else if (m_collision_type == CollisionType::DeuteriumTritiumFusion)
+        else if ((m_collision_type == CollisionType::DeuteriumTritiumToNeutronHeliumFusion)
+              || (m_collision_type == CollisionType::DeuteriumDeuteriumToProtonTritiumFusion)
+              || (m_collision_type == CollisionType::DeuteriumDeuteriumToNeutronHeliumFusion))
         {
             m_num_product_species = 2;
             m_num_products_host.push_back(1);