Add structure for proton-boron fusion module (#2217)

* Add structure for proton-boron fusion module

* Fix clang compile error

* Generalize structure for other nuclear fusion types

* Workaround: try to use std::vectors to fix clang compile error

* Simplify workaround for clang

* Rebase on 2245 and add empty particle creation functor

* Apply suggestions from code review

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Add couple of consts and runtime assert

* Add another const

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

13 files changed, 1495 insertions, 0 deletions

diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
new file mode 100644
index 000000000..338ada68f
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
@@ -0,0 +1,554 @@
+/* Copyright 2020-2021 Yinjian Zhao, David Grote, Neil Zaim
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_BINARYCOLLISION_H_
+#define WARPX_PARTICLES_COLLISION_BINARYCOLLISION_H_
+
+#include "Particles/Collision/BinaryCollision/NuclearFusionFunc.H"
+#include "Particles/Collision/BinaryCollision/PairWiseCoulombCollisionFunc.H"
+#include "Particles/Collision/BinaryCollision/ParticleCreationFunc.H"
+#include "Particles/Collision/BinaryCollision/ShuffleFisherYates.H"
+#include "Particles/Collision/CollisionBase.H"
+#include "Particles/ParticleCreation/SmartCopy.H"
+#include "Particles/ParticleCreation/SmartUtils.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
+#include "Particles/MultiParticleContainer.H"
+#include "Particles/WarpXParticleContainer.H"
+#include "Utils/ParticleUtils.H"
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "WarpX.H"
+
+#include "Particles/MultiParticleContainer_fwd.H"
+#include "Particles/WarpXParticleContainer_fwd.H"
+
+#include <AMReX.H>
+#include <AMReX_Algorithm.H>
+#include <AMReX_BLassert.H>
+#include <AMReX_Config.H>
+#include <AMReX_DenseBins.H>
+#include <AMReX_Extension.H>
+#include <AMReX_Geometry.H>
+#include <AMReX_GpuAtomic.H>
+#include <AMReX_GpuContainers.H>
+#include <AMReX_GpuControl.H>
+#include <AMReX_GpuDevice.H>
+#include <AMReX_GpuLaunch.H>
+#include <AMReX_GpuQualifiers.H>
+#include <AMReX_LayoutData.H>
+#include <AMReX_MFIter.H>
+#include <AMReX_PODVector.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_Particles.H>
+#include <AMReX_ParticleTile.H>
+#include <AMReX_Random.H>
+#include <AMReX_REAL.H>
+#include <AMReX_Scan.H>
+#include <AMReX_Utility.H>
+#include <AMReX_Vector.H>
+
+#include <AMReX_BaseFwd.H>
+
+#include <cmath>
+#include <string>
+
+/**
+ * \brief This class performs generic binary collisions.
+ *
+ * \tparam CollisionFunctorType the type of the specific binary collision functor that acts on a
+ *         single cell
+ * \tparam CopyTransformFunctorType the type of the second functor used in the case of
+ *         particle creation
+ *
+ */
+template <typename CollisionFunctorType,
+          typename CopyTransformFunctorType = NoParticleCreationFunc>
+class BinaryCollision final
+    : public CollisionBase
+{
+    // Define shortcuts for frequently-used type names
+    using ParticleType = WarpXParticleContainer::ParticleType;
+    using ParticleTileType = WarpXParticleContainer::ParticleTileType;
+    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+    using index_type = ParticleBins::index_type;
+
+public:
+    /**
+     * \brief Constructor of the BinaryCollision class.
+     *
+     * @param collision_name the name of the collision
+     *
+     */
+    BinaryCollision (std::string collision_name, MultiParticleContainer const * const mypc)
+        : CollisionBase(collision_name)
+    {
+        if(m_species_names.size() != 2)
+            amrex::Abort("Binary collision " + collision_name + " must have exactly two species.");
+
+        if (m_species_names[0] == m_species_names[1])
+            m_isSameSpecies = true;
+        else
+            m_isSameSpecies = false;
+
+        m_binary_collision_functor = CollisionFunctorType(collision_name, mypc);
+
+        amrex::ParmParse pp_collision_name(collision_name);
+        pp_collision_name.queryarr("product_species", m_product_species);
+        m_have_product_species = m_product_species.size() > 0;
+        m_copy_transform_functor = CopyTransformFunctorType(collision_name, mypc);
+    }
+
+    virtual ~BinaryCollision () = default;
+
+    /** Perform the collisions
+     *
+     * @param lev AMR level of the tile
+     * @param cur_time Current time
+     * @param mypc Container of species involved
+     *
+     */
+    void doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc) override
+    {
+        const amrex::Real dt = WarpX::GetInstance().getdt(0);
+        if ( int(std::floor(cur_time/dt)) % m_ndt != 0 ) return;
+
+        auto& species1 = mypc->GetParticleContainerFromName(m_species_names[0]);
+        auto& species2 = mypc->GetParticleContainerFromName(m_species_names[1]);
+
+        // In case of particle creation, create the necessary vectors
+        const int n_product_species = m_product_species.size();
+        amrex::Vector<WarpXParticleContainer*> product_species_vector;
+        amrex::Vector<SmartCopy> copy_species1;
+        amrex::Vector<SmartCopy> copy_species2;
+        for (int i = 0; i < n_product_species; i++)
+        {
+            auto& product = mypc->GetParticleContainerFromName(m_product_species[i]);
+            product.defineAllParticleTiles();
+            product_species_vector.push_back(&product);
+            SmartCopyFactory copy_factory_species1(species1, product);
+            SmartCopyFactory copy_factory_species2(species2, product);
+            copy_species1.push_back(copy_factory_species1.getSmartCopy());
+            copy_species2.push_back(copy_factory_species2.getSmartCopy());
+        }
+#ifdef AMREX_USE_GPU
+        amrex::Gpu::DeviceVector<SmartCopy> device_copy_species1(n_product_species);
+        amrex::Gpu::DeviceVector<SmartCopy> device_copy_species2(n_product_species);
+        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, copy_species1.begin(),
+                              copy_species1.end(), device_copy_species1.begin());
+        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, copy_species2.begin(),
+                              copy_species2.end(), device_copy_species2.begin());
+        amrex::Gpu::streamSynchronize();
+        auto copy_species1_data = device_copy_species1.data();
+        auto copy_species2_data = device_copy_species2.data();
+#else
+        auto copy_species1_data = copy_species1.data();
+        auto copy_species2_data = copy_species2.data();
+#endif
+        if (m_have_product_species){
+            species1.defineAllParticleTiles();
+            species2.defineAllParticleTiles();
+        }
+
+        // Enable tiling
+        amrex::MFItInfo info;
+        if (amrex::Gpu::notInLaunchRegion()) info.EnableTiling(species1.tile_size);
+
+        // Loop over refinement levels
+        for (int lev = 0; lev <= species1.finestLevel(); ++lev){
+
+        amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
+
+            // Loop over all grids/tiles at this level
+#ifdef AMREX_USE_OMP
+            info.SetDynamic(true);
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+            for (amrex::MFIter mfi = species1.MakeMFIter(lev, info); mfi.isValid(); ++mfi){
+                if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+                {
+                    amrex::Gpu::synchronize();
+                }
+                amrex::Real wt = amrex::second();
+
+                doCollisionsWithinTile( lev, mfi, species1, species2, product_species_vector,
+                                         copy_species1_data, copy_species2_data);
+
+                if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+                {
+                    amrex::Gpu::synchronize();
+                    wt = amrex::second() - wt;
+                    amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+                }
+            }
+        }
+    }
+
+    /** Perform all binary collisions within a tile
+     *
+     * \param mfi iterator for multifab
+     * \param species_1 first species container
+     * \param species_2 second species container
+     * \param product_species_vector vector of pointers to product species containers
+     * \param copy_species1 vector of SmartCopy functors used to copy species 1 to product species
+     * \param copy_species2 vector of SmartCopy functors used to copy species 2 to product species
+     *
+     */
+    void doCollisionsWithinTile (
+        int const lev, amrex::MFIter const& mfi,
+        WarpXParticleContainer& species_1,
+        WarpXParticleContainer& species_2,
+        amrex::Vector<WarpXParticleContainer*> product_species_vector,
+        SmartCopy* copy_species1, SmartCopy* copy_species2)
+    {
+        using namespace ParticleUtils;
+        using namespace amrex::literals;
+
+        int const ndt = m_ndt;
+        CollisionFunctorType binary_collision_functor = m_binary_collision_functor;
+        const bool have_product_species = m_have_product_species;
+
+        // Store product species data in vectors
+        const int n_product_species = m_product_species.size();
+        amrex::Vector<SoaData_type> soa_products;
+        amrex::Vector<GetParticlePosition> get_position_products;
+        amrex::Vector<index_type> products_np;
+        constexpr int getpos_offset = 0;
+        for (int i = 0; i < n_product_species; i++)
+        {
+            ParticleTileType& ptile_product = product_species_vector[i]->ParticlesAt(lev, mfi);
+            soa_products.push_back(ptile_product.getParticleTileData());
+            get_position_products.push_back(GetParticlePosition(ptile_product,
+                                                                      getpos_offset));
+            products_np.push_back(ptile_product.numParticles());
+        }
+#ifdef AMREX_USE_GPU
+        amrex::Gpu::DeviceVector<SoaData_type> device_soa_products(n_product_species);
+        amrex::Gpu::DeviceVector<GetParticlePosition> device_get_position_products(n_product_species);
+        amrex::Gpu::DeviceVector<index_type> device_products_np(n_product_species);
+        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, soa_products.begin(),
+                              soa_products.end(),
+                              device_soa_products.begin());
+        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, get_position_products.begin(),
+                              get_position_products.end(),
+                              device_get_position_products.begin());
+        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, products_np.begin(),
+                              products_np.end(),
+                              device_products_np.begin());
+        amrex::Gpu::streamSynchronize();
+        auto soa_products_data = device_soa_products.data();
+        auto get_position_products_data = device_get_position_products.data();
+        auto products_np_data = device_products_np.data();
+#else
+        auto soa_products_data = soa_products.data();
+        auto get_position_products_data = get_position_products.data();
+        auto products_np_data = products_np.data();
+#endif
+
+        if ( m_isSameSpecies ) // species_1 == species_2
+        {
+            // Extract particles in the tile that `mfi` points to
+            ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
+
+            // Find the particles that are in each cell of this tile
+            ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
+
+            // Loop over cells, and collide the particles in each cell
+
+            // Extract low-level data
+            int const n_cells = bins_1.numBins();
+            // - Species 1
+            const auto soa_1 = ptile_1.getParticleTileData();
+            index_type* indices_1 = bins_1.permutationPtr();
+            index_type const* cell_offsets_1 = bins_1.offsetsPtr();
+            amrex::Real q1 = species_1.getCharge();
+            amrex::Real m1 = species_1.getMass();
+            auto get_position_1  = GetParticlePosition(ptile_1, getpos_offset);
+
+            const amrex::Real dt = WarpX::GetInstance().getdt(lev);
+            amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
+#if defined WARPX_DIM_XZ
+            auto dV = geom.CellSize(0) * geom.CellSize(1);
+#elif defined WARPX_DIM_RZ
+            amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box
+            const auto lo = lbound(cbx);
+            const auto hi = ubound(cbx);
+            int const nz = hi.y-lo.y+1;
+            auto dr = geom.CellSize(0);
+            auto dz = geom.CellSize(1);
+#elif (AMREX_SPACEDIM == 3)
+            auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
+#endif
+
+
+            /*
+              The following calculations are only required when creating product particles
+            */
+            const int n_cells_products = have_product_species ? n_cells: 0;
+            amrex::Gpu::DeviceVector<index_type> n_pairs_in_each_cell(n_cells_products);
+            auto p_n_pairs_in_each_cell = n_pairs_in_each_cell.dataPtr();
+
+            // Compute how many pairs in each cell and store in n_pairs_in_each_cell array
+            // For a single species, the number of pair in a cell is half the number of particles
+            // in that cell, rounded up to the next higher integer.
+            amrex::ParallelFor( n_cells_products,
+                [=] AMREX_GPU_DEVICE (int i_cell) noexcept
+                {
+                    const auto n_part_in_cell = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
+                    // Particular case: if there's only 1 particle in a cell, then there's no pair
+                    p_n_pairs_in_each_cell[i_cell] = (n_part_in_cell == 1)? 0: (n_part_in_cell+1)/2;
+                }
+            );
+
+            // Start indices of the pairs in a cell. Will be used for particle creation.
+            amrex::Gpu::DeviceVector<index_type> pair_offsets(n_cells_products);
+            const auto n_total_pairs = amrex::Scan::ExclusiveSum(n_cells_products,
+                                                     p_n_pairs_in_each_cell, pair_offsets.data());
+            auto p_pair_offsets = pair_offsets.dataPtr();
+
+            // mask: equal to 1 if particle creation occurs for a given pair, 0 otherwise
+            amrex::Gpu::DeviceVector<index_type> mask(n_total_pairs);
+            auto p_mask = mask.dataPtr();
+            // Will be filled with the index of the first particle of a given pair
+            amrex::Gpu::DeviceVector<index_type> pair_indices_1(n_total_pairs);
+            auto p_pair_indices_1 = pair_indices_1.dataPtr();
+            // Will be filled with the index of the second particle of a given pair
+            amrex::Gpu::DeviceVector<index_type> pair_indices_2(n_total_pairs);
+            auto p_pair_indices_2 = pair_indices_2.dataPtr();
+            // How much weight should be given to the produced particles (and removed from the
+            // reacting particles)
+            amrex::Gpu::DeviceVector<amrex::ParticleReal> pair_reaction_weight(n_total_pairs);
+            auto p_pair_reaction_weight = pair_reaction_weight.dataPtr();
+            /*
+              End of calculations only required when creating product particles
+            */
+
+
+            // Loop over cells
+            amrex::ParallelForRNG( n_cells,
+                [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
+                {
+                    // The particles from species1 that are in the cell `i_cell` are
+                    // given by the `indices_1[cell_start_1:cell_stop_1]`
+                    index_type const cell_start_1 = cell_offsets_1[i_cell];
+                    index_type const cell_stop_1  = cell_offsets_1[i_cell+1];
+                    index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2;
+
+                    // Same but for the pairs
+                    index_type const cell_start_pair = have_product_species?
+                                                           p_pair_offsets[i_cell] : 0;
+
+                    // Do not collide if there is only one particle in the cell
+                    if ( cell_stop_1 - cell_start_1 <= 1 ) return;
+
+                    // shuffle
+                    ShuffleFisherYates(
+                        indices_1, cell_start_1, cell_half_1, engine );
+#if defined WARPX_DIM_RZ
+                    int ri = (i_cell - i_cell%nz) / nz;
+                    auto dV = MathConst::pi*(2.0_rt*ri+1.0_rt)*dr*dr*dz;
+#endif
+                    // Call the function in order to perform collisions
+                    // If there are product species, mask, p_pair_indices_1/2, and
+                    // p_pair_reaction_weight are filled here
+                    binary_collision_functor(
+                        cell_start_1, cell_half_1,
+                        cell_half_1, cell_stop_1,
+                        indices_1, indices_1,
+                        soa_1, soa_1, get_position_1, get_position_1,
+                        q1, q1, m1, m1, dt*ndt, dV,
+                        cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
+                        p_pair_reaction_weight, engine );
+                }
+            );
+
+            // Create the new product particles and define their initial values
+            // num_added: how many particles of each product species have been created
+            const amrex::Vector<int> num_added = m_copy_transform_functor(soa_1, soa_1,
+                                                    soa_products_data,
+                                                    get_position_1, get_position_1,
+                                                    get_position_products_data,
+                                                    p_mask, products_np_data,
+                                                    copy_species1, copy_species2,
+                                                    p_pair_indices_1, p_pair_indices_2,
+                                                    p_pair_reaction_weight);
+
+            for (int i = 0; i < n_product_species; i++)
+            {
+                ParticleTileType& ptile_product = product_species_vector[i]->ParticlesAt(lev, mfi);
+                setNewParticleIDs(ptile_product, products_np[i], num_added[i]);
+            }
+        }
+        else // species_1 != species_2
+        {
+            // Extract particles in the tile that `mfi` points to
+            ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
+            ParticleTileType& ptile_2 = species_2.ParticlesAt(lev, mfi);
+
+            // Find the particles that are in each cell of this tile
+            ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
+            ParticleBins bins_2 = findParticlesInEachCell( lev, mfi, ptile_2 );
+
+            // Loop over cells, and collide the particles in each cell
+
+            // Extract low-level data
+            int const n_cells = bins_1.numBins();
+            // - Species 1
+            const auto soa_1 = ptile_1.getParticleTileData();
+            index_type* indices_1 = bins_1.permutationPtr();
+            index_type const* cell_offsets_1 = bins_1.offsetsPtr();
+            amrex::Real q1 = species_1.getCharge();
+            amrex::Real m1 = species_1.getMass();
+            auto get_position_1  = GetParticlePosition(ptile_1, getpos_offset);
+            // - Species 2
+            const auto soa_2 = ptile_2.getParticleTileData();
+            index_type* indices_2 = bins_2.permutationPtr();
+            index_type const* cell_offsets_2 = bins_2.offsetsPtr();
+            amrex::Real q2 = species_2.getCharge();
+            amrex::Real m2 = species_2.getMass();
+            auto get_position_2  = GetParticlePosition(ptile_2, getpos_offset);
+
+            const amrex::Real dt = WarpX::GetInstance().getdt(lev);
+            amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
+#if defined WARPX_DIM_XZ
+            auto dV = geom.CellSize(0) * geom.CellSize(1);
+#elif defined WARPX_DIM_RZ
+            amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box
+            const auto lo = lbound(cbx);
+            const auto hi = ubound(cbx);
+            int nz = hi.y-lo.y+1;
+            auto dr = geom.CellSize(0);
+            auto dz = geom.CellSize(1);
+#elif (AMREX_SPACEDIM == 3)
+            auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
+#endif
+
+
+            /*
+              The following calculations are only required when creating product particles
+            */
+            const int n_cells_products = have_product_species ? n_cells: 0;
+            amrex::Gpu::DeviceVector<index_type> n_pairs_in_each_cell(n_cells_products);
+            auto p_n_pairs_in_each_cell = n_pairs_in_each_cell.dataPtr();
+
+            // Compute how many pairs in each cell and store in n_pairs_in_each_cell array
+            // For different species, the number of pairs in a cell is the number of particles of
+            // the species that has the most particles in that cell
+            amrex::ParallelFor( n_cells_products,
+                [=] AMREX_GPU_DEVICE (int i_cell) noexcept
+                {
+                    const auto n_part_in_cell_1 = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
+                    const auto n_part_in_cell_2 = cell_offsets_2[i_cell+1] - cell_offsets_2[i_cell];
+                    // Particular case: no pair if a species has no particle in that cell
+                    if (n_part_in_cell_1 == 0 || n_part_in_cell_2 == 0)
+                        p_n_pairs_in_each_cell[i_cell] = 0;
+                    else
+                        p_n_pairs_in_each_cell[i_cell] =
+                                                      amrex::max(n_part_in_cell_1,n_part_in_cell_2);
+                }
+            );
+
+            // Start indices of the pairs in a cell. Will be used for particle creation
+            amrex::Gpu::DeviceVector<index_type> pair_offsets(n_cells_products);
+            auto n_total_pairs = amrex::Scan::ExclusiveSum(n_cells_products,
+                                                     p_n_pairs_in_each_cell, pair_offsets.data());
+            auto p_pair_offsets = pair_offsets.dataPtr();
+
+            // mask: equal to 1 if particle creation occurs for a given pair, 0 otherwise
+            amrex::Gpu::DeviceVector<index_type> mask(n_total_pairs);
+            auto p_mask = mask.dataPtr();
+            // Will be filled with the index of the first particle of a given pair
+            amrex::Gpu::DeviceVector<index_type> pair_indices_1(n_total_pairs);
+            auto p_pair_indices_1 = pair_indices_1.dataPtr();
+            // Will be filled with the index of the second particle of a given pair
+            amrex::Gpu::DeviceVector<index_type> pair_indices_2(n_total_pairs);
+            auto p_pair_indices_2 = pair_indices_2.dataPtr();
+            // How much weight should be given to the produced particles (and removed from the
+            // reacting particles)
+            amrex::Gpu::DeviceVector<amrex::ParticleReal> pair_reaction_weight(n_total_pairs);
+            auto p_pair_reaction_weight = pair_reaction_weight.dataPtr();
+            /*
+              End of calculations only required when creating product particles
+            */
+
+
+            // Loop over cells
+            amrex::ParallelForRNG( n_cells,
+                [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
+                {
+                    // The particles from species1 that are in the cell `i_cell` are
+                    // given by the `indices_1[cell_start_1:cell_stop_1]`
+                    index_type const cell_start_1 = cell_offsets_1[i_cell];
+                    index_type const cell_stop_1  = cell_offsets_1[i_cell+1];
+                    // Same for species 2
+                    index_type const cell_start_2 = cell_offsets_2[i_cell];
+                    index_type const cell_stop_2  = cell_offsets_2[i_cell+1];
+                    // Same but for the pairs
+                    index_type const cell_start_pair = have_product_species?
+                                                           p_pair_offsets[i_cell] : 0;
+
+                    // ux from species1 can be accessed like this:
+                    // ux_1[ indices_1[i] ], where i is between
+                    // cell_start_1 (inclusive) and cell_start_2 (exclusive)
+
+                    // Do not collide if one species is missing in the cell
+                    if ( cell_stop_1 - cell_start_1 < 1 ||
+                        cell_stop_2 - cell_start_2 < 1 ) return;
+
+                    // shuffle
+                    ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine);
+                    ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2, engine);
+#if defined WARPX_DIM_RZ
+                    int ri = (i_cell - i_cell%nz) / nz;
+                    auto dV = MathConst::pi*(2.0_rt*ri+1.0_rt)*dr*dr*dz;
+#endif
+                    // Call the function in order to perform collisions
+                    // If there are product species, p_mask, p_pair_indices_1/2, and
+                    // p_pair_reaction_weight are filled here
+                    binary_collision_functor(
+                        cell_start_1, cell_stop_1, cell_start_2, cell_stop_2,
+                        indices_1, indices_2,
+                        soa_1, soa_2, get_position_1, get_position_2,
+                        q1, q2, m1, m2, dt*ndt, dV,
+                        cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
+                        p_pair_reaction_weight, engine );
+                }
+            );
+
+
+            // Create the new product particles and define their initial values
+            // num_added: how many particles of each product species have been created
+            const amrex::Vector<int> num_added = m_copy_transform_functor(soa_1, soa_2,
+                                                    soa_products_data,
+                                                    get_position_1, get_position_2,
+                                                    get_position_products_data,
+                                                    p_mask, products_np_data,
+                                                    copy_species1, copy_species2,
+                                                    p_pair_indices_1, p_pair_indices_2,
+                                                    p_pair_reaction_weight);
+
+            for (int i = 0; i < n_product_species; i++)
+            {
+                ParticleTileType& ptile_product = product_species_vector[i]->ParticlesAt(lev, mfi);
+                setNewParticleIDs(ptile_product, products_np[i], num_added[i]);
+            }
+
+        } // end if ( m_isSameSpecies)
+
+    }
+
+private:
+
+    bool m_isSameSpecies;
+    bool m_have_product_species;
+    amrex::Vector<std::string> m_product_species;
+    // functor that performs collisions within a cell
+    CollisionFunctorType m_binary_collision_functor;
+    // functor that creates new particles and initializes their parameters
+    CopyTransformFunctorType m_copy_transform_functor;
+
+};
+
+#endif // WARPX_PARTICLES_COLLISION_BINARYCOLLISION_H_
diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H
new file mode 100644
index 000000000..c20e9716c
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H
@@ -0,0 +1,30 @@
+/* Copyright 2021 Neil Zaim
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef BINARY_COLLISION_UTILS_H_
+#define BINARY_COLLISION_UTILS_H_
+
+#include <string>
+
+#include "Particles/MultiParticleContainer.H"
+
+enum struct CollisionType { ProtonBoronFusion, Undefined };
+
+enum struct NuclearFusionType { ProtonBoron, Undefined };
+
+namespace BinaryCollisionUtils{
+
+    NuclearFusionType get_nuclear_fusion_type (const std::string collision_name,
+                                               MultiParticleContainer const * const mypc);
+
+    CollisionType get_collision_type (const std::string collision_name,
+                                      MultiParticleContainer const * const mypc);
+
+    CollisionType nuclear_fusion_type_to_collision_type (const NuclearFusionType fusion_type);
+}
+
+#endif // BINARY_COLLISION_UTILS_H_
diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
new file mode 100644
index 000000000..fc70b39ca
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
@@ -0,0 +1,64 @@
+/* Copyright 2021 Neil Zaim
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "BinaryCollisionUtils.H"
+
+#include "Particles/MultiParticleContainer.H"
+#include "Particles/WarpXParticleContainer.H"
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_Vector.H>
+
+#include <string>
+
+namespace BinaryCollisionUtils{
+
+    NuclearFusionType get_nuclear_fusion_type (const std::string collision_name,
+                                               MultiParticleContainer const * const mypc)
+        {
+            amrex::ParmParse pp_collision_name(collision_name);
+            amrex::Vector<std::string> species_names;
+            pp_collision_name.getarr("species", species_names);
+            auto& species1 = mypc->GetParticleContainerFromName(species_names[0]);
+            auto& species2 = mypc->GetParticleContainerFromName(species_names[1]);
+
+            if ((species1.AmIA<PhysicalSpecies::hydrogen>() && species2.AmIA<PhysicalSpecies::boron>())
+                ||
+                (species1.AmIA<PhysicalSpecies::boron>() && species2.AmIA<PhysicalSpecies::hydrogen>())
+                )
+            {
+                return NuclearFusionType::ProtonBoron;
+            }
+            amrex::Abort("Binary nuclear fusion not implemented between species " +
+                         species_names[0] + " of type " + species1.getSpeciesTypeName() +
+                         " and species " + species_names[1] + " of type " +
+                         species2.getSpeciesTypeName());
+            return NuclearFusionType::Undefined;
+        }
+
+    CollisionType get_collision_type (const std::string collision_name,
+                                      MultiParticleContainer const * const mypc)
+        {
+            amrex::ParmParse pp_collision_name(collision_name);
+            std::string type;
+            pp_collision_name.get("type", type);
+            if (type == "nuclearfusion") {
+                NuclearFusionType fusion_type = get_nuclear_fusion_type(collision_name, mypc);
+                return nuclear_fusion_type_to_collision_type(fusion_type);
+            }
+            amrex::Abort(type + " is not a valid type of collision that creates new particles");
+            return CollisionType::Undefined;
+        }
+
+    CollisionType nuclear_fusion_type_to_collision_type (const NuclearFusionType fusion_type)
+        {
+            if (fusion_type == NuclearFusionType::ProtonBoron)
+                return CollisionType::ProtonBoronFusion;
+            amrex::Abort("Invalid nuclear fusion type");
+            return CollisionType::Undefined;
+        }
+}
diff --git a/Source/Particles/Collision/BinaryCollision/CMakeLists.txt b/Source/Particles/Collision/BinaryCollision/CMakeLists.txt
new file mode 100644
index 000000000..c9a632f2e
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/CMakeLists.txt
@@ -0,0 +1,6 @@
+target_sources(WarpX
+  PRIVATE
+    BinaryCollisionUtils.cpp
+    ParticleCreationFunc.cpp
+)
+
diff --git a/Source/Particles/Collision/BinaryCollision/ComputeTemperature.H b/Source/Particles/Collision/BinaryCollision/ComputeTemperature.H
new file mode 100644
index 000000000..50e1f4a4b
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/ComputeTemperature.H
@@ -0,0 +1,47 @@
+/* Copyright 2019-2020 Andrew Myers, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_
+#define WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_
+
+#include "Utils/WarpXConst.H"
+
+
+template <typename T_index, typename T_R>
+AMREX_GPU_HOST_DEVICE
+T_R ComputeTemperature (
+    T_index const Is, T_index const Ie, T_index const *I,
+    T_R const *ux, T_R const *uy, T_R const *uz, T_R const m )
+{
+
+    T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c );
+
+    int N = Ie - Is;
+    if ( N == 0 ) { return T_R(0.0); }
+
+    T_R vx = T_R(0.0);    T_R vy = T_R(0.0);
+    T_R vz = T_R(0.0);    T_R vs = T_R(0.0);
+    T_R gm = T_R(0.0);    T_R us = T_R(0.0);
+
+    for (int i = Is; i < static_cast<int>(Ie); ++i)
+    {
+        us = ( ux[ I[i] ] * ux[ I[i] ] +
+               uy[ I[i] ] * uy[ I[i] ] +
+               uz[ I[i] ] * uz[ I[i] ] );
+        gm = std::sqrt( T_R(1.0) + us*inv_c2 );
+        vx += ux[ I[i] ] / gm;
+        vy += uy[ I[i] ] / gm;
+        vz += uz[ I[i] ] / gm;
+        vs += us / gm / gm;
+    }
+
+    vx = vx / N;    vy = vy / N;
+    vz = vz / N;    vs = vs / N;
+
+    return m/T_R(3.0)*(vs-(vx*vx+vy*vy+vz*vz));
+}
+
+#endif // WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_
diff --git a/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H b/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H
new file mode 100644
index 000000000..bee415860
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/ElasticCollisionPerez.H
@@ -0,0 +1,126 @@
+/* Copyright 2019 Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
+#define WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
+
+#include "ComputeTemperature.H"
+#include "UpdateMomentumPerezElastic.H"
+#include "Particles/WarpXParticleContainer.H"
+#include "Utils/WarpXConst.H"
+
+#include <AMReX_Random.H>
+
+
+/** \brief Prepare information for and call
+ *        UpdateMomentumPerezElastic().
+ * @param[in] I1s,I2s is the start index for I1,I2 (inclusive).
+ * @param[in] I1e,I2e is the stop index for I1,I2 (exclusive).
+ * @param[in] I1 and I2 are the index arrays. They determine all elements that will be used.
+ * @param[in,out] u1 and u2 are the velocity arrays (u=v*gamma),
+ *                they could be either different or the same,
+ *                their lengths are not needed,
+ * @param[in] w1 and w2 are arrays of weights.
+ * @param[in] q1 and q2 are charges. m1 and m2 are masses.
+ * @param[in] T1 and T2 are temperatures (Joule)
+ *            and will be used if greater than zero,
+ *            otherwise will be computed.
+ * @param[in] dt is the time step length between two collision calls.
+ * @param[in] L is the Coulomb log and will be used if greater than zero,
+ *            otherwise will be computed.
+ * @param[in] dV is the volume of the corresponding cell.
+*/
+
+template <typename T_index, typename T_R, typename SoaData_type>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void ElasticCollisionPerez (
+    T_index const I1s, T_index const I1e,
+    T_index const I2s, T_index const I2e,
+    T_index *I1,       T_index *I2,
+    SoaData_type soa_1, SoaData_type soa_2,
+    T_R const  q1, T_R const  q2,
+    T_R const  m1, T_R const  m2,
+    T_R const  T1, T_R const  T2,
+    T_R const  dt, T_R const   L, T_R const dV,
+    amrex::RandomEngine const& engine)
+{
+    int NI1 = I1e - I1s;
+    int NI2 = I2e - I2s;
+
+    T_R * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
+    T_R * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
+    T_R * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
+    T_R * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
+
+    T_R * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
+    T_R * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
+    T_R * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
+    T_R * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
+
+    // get local T1t and T2t
+    T_R T1t; T_R T2t;
+    if ( T1 <= T_R(0.0) && L <= T_R(0.0) )
+    {
+        T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,u1z,m1);
+    }
+    else { T1t = T1; }
+    if ( T2 <= T_R(0.0) && L <= T_R(0.0) )
+    {
+        T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,u2z,m2);
+    }
+    else { T2t = T2; }
+
+    // local density
+    T_R n1  = T_R(0.0);
+    T_R n2  = T_R(0.0);
+    T_R n12 = T_R(0.0);
+    for (int i1=I1s; i1<static_cast<int>(I1e); ++i1) { n1 += w1[ I1[i1] ]; }
+    for (int i2=I2s; i2<static_cast<int>(I2e); ++i2) { n2 += w2[ I2[i2] ]; }
+    n1 = n1 / dV; n2 = n2 / dV;
+    {
+      int i1 = I1s; int i2 = I2s;
+      for (int k = 0; k < amrex::max(NI1,NI2); ++k)
+      {
+        n12 += amrex::min( w1[ I1[i1] ], w2[ I2[i2] ] );
+        ++i1; if ( i1 == static_cast<int>(I1e) ) { i1 = I1s; }
+        ++i2; if ( i2 == static_cast<int>(I2e) ) { i2 = I2s; }
+      }
+      n12 = n12 / dV;
+    }
+
+    // compute Debye length lmdD
+    T_R lmdD;
+    if ( T1t < T_R(0.0) || T2t < T_R(0.0) ) {
+        lmdD = T_R(0.0);
+    }
+    else {
+        lmdD = T_R(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) +
+                                   n2*q2*q2/(T2t*PhysConst::ep0) );
+    }
+    T_R rmin = std::pow( T_R(4.0) * MathConst::pi / T_R(3.0) *
+               amrex::max(n1,n2), T_R(-1.0/3.0) );
+    lmdD = amrex::max(lmdD, rmin);
+
+    // call UpdateMomentumPerezElastic()
+    {
+      int i1 = I1s; int i2 = I2s;
+      for (int k = 0; k < amrex::max(NI1,NI2); ++k)
+      {
+          UpdateMomentumPerezElastic(
+              u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
+              u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
+              n1, n2, n12,
+              q1, m1, w1[ I1[i1] ], q2, m2, w2[ I2[i2] ],
+              dt, L, lmdD,
+              engine);
+          ++i1; if ( i1 == static_cast<int>(I1e) ) { i1 = I1s; }
+          ++i2; if ( i2 == static_cast<int>(I2e) ) { i2 = I2s; }
+      }
+    }
+
+}
+
+#endif // WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_
diff --git a/Source/Particles/Collision/BinaryCollision/Make.package b/Source/Particles/Collision/BinaryCollision/Make.package
new file mode 100644
index 000000000..92e2f0b5d
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/Make.package
@@ -0,0 +1,4 @@
+CEXE_sources += BinaryCollisionUtils.cpp
+CEXE_sources += ParticleCreationFunc.cpp
+
+VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/Collision/BinaryCollision
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusionFunc.H
new file mode 100644
index 000000000..6dfe718a0
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusionFunc.H
@@ -0,0 +1,104 @@
+/* Copyright 2021 Neil Zaim
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef NUCLEAR_FUSION_FUNC_H_
+#define NUCLEAR_FUSION_FUNC_H_
+
+#include "BinaryCollisionUtils.H"
+
+#include "Particles/Pusher/GetAndSetPosition.H"
+#include "Particles/MultiParticleContainer.H"
+#include "Particles/WarpXParticleContainer.H"
+#include "Utils/WarpXUtil.H"
+#include "WarpX.H"
+
+#include <AMReX_DenseBins.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_Random.H>
+#include <AMReX_REAL.H>
+#include <AMReX_Vector.H>
+
+/**
+ * \brief This functor does binary nuclear fusions on a single cell.
+ *  Particles of the two reacting species are paired with each other and for each pair we compute
+ *  if a fusion event occurs. If so, we fill a mask (input parameter p_mask) with true so that
+ *  product particles corresponding to a given pair can be effectively created in the particle
+ *  creation functor.
+ *  This functor also reads and contains the fusion multiplier.
+ */
+class NuclearFusionFunc{
+    // Define shortcuts for frequently-used type names
+    using ParticleType = WarpXParticleContainer::ParticleType;
+    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using index_type = ParticleBins::index_type;
+    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+
+public:
+    /**
+     * \brief Default constructor of the NuclearFusionFunc class.
+     */
+    NuclearFusionFunc () = default;
+
+    /**
+     * \brief Constructor of the NuclearFusionFunc class
+     *
+     * @param[in] collision_name the name of the collision
+     * @param[in] mypc pointer to the MultiParticleContainer
+     */
+    NuclearFusionFunc (const std::string collision_name, MultiParticleContainer const * const mypc)
+    {
+        using namespace amrex::literals;
+
+        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::ProtonBoron)
+        {
+            AMREX_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]);
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+                product_species.AmIA<PhysicalSpecies::helium>(),
+                "ERROR: Product species of proton-boron fusion must be of type helium");
+        }
+
+        // default fusion multiplier
+        m_fusion_multiplier = 1.0_rt;
+        queryWithParser(pp_collision_name, "fusion_multiplier", m_fusion_multiplier);
+    }
+
+    /**
+     * \brief operator() of the NuclearFusionFunc class. Needs to be implemented.
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_INLINE
+    void operator() (
+        index_type const /*I1s*/, index_type const /*I1e*/,
+        index_type const /*I2s*/, index_type const /*I2e*/,
+        index_type* /*I1*/,       index_type* /*I2*/,
+        SoaData_type /*soa_1*/, SoaData_type /*soa_2*/,
+        GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/,
+        amrex::Real const  /*q1*/, amrex::Real const  /*q2*/,
+        amrex::Real const  /*m1*/, amrex::Real const  /*m2*/,
+        amrex::Real const  /*dt*/, amrex::Real const /*dV*/,
+        index_type const /*cell_start_pair*/, index_type* /*p_mask*/,
+        index_type* /*p_pair_indices_1*/, index_type* /*p_pair_indices_2*/,
+        amrex::ParticleReal* /*p_pair_reaction_weight*/,
+        amrex::RandomEngine const& /*engine*/) const
+        {
+            amrex::Abort("Nuclear fusion module not yet implemented");
+        }
+
+private:
+    amrex::Real m_fusion_multiplier;
+    NuclearFusionType m_fusion_type;
+};
+
+#endif // NUCLEAR_FUSION_FUNC_H_
diff --git a/Source/Particles/Collision/BinaryCollision/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/PairWiseCoulombCollisionFunc.H
new file mode 100644
index 000000000..086803bba
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/PairWiseCoulombCollisionFunc.H
@@ -0,0 +1,96 @@
+/* Copyright 2021 Neil Zaim
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef PAIRWISE_COULOMB_COLLISION_FUNC_H_
+#define PAIRWISE_COULOMB_COLLISION_FUNC_H_
+
+#include "ElasticCollisionPerez.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
+#include "Particles/WarpXParticleContainer.H"
+#include "Utils/WarpXUtil.H"
+
+#include <AMReX_DenseBins.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_Random.H>
+#include <AMReX_REAL.H>
+
+
+/**
+ * \brief This functor performs pairwise Coulomb collision on a single cell by calling the function
+ *        ElasticCollisionPerez. It also reads and contains the Coulomb logarithm.
+ */
+class PairWiseCoulombCollisionFunc{
+    // Define shortcuts for frequently-used type names
+    using ParticleType = WarpXParticleContainer::ParticleType;
+    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using index_type = ParticleBins::index_type;
+    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+
+public:
+    /**
+     * \brief Default constructor of the PairWiseCoulombCollisionFunc class.
+     */
+    PairWiseCoulombCollisionFunc () = default;
+
+    /**
+     * \brief Constructor of the PairWiseCoulombCollisionFunc class
+     *
+     * @param[in] collision_name the name of the collision
+     */
+    PairWiseCoulombCollisionFunc (const std::string collision_name,
+                                  MultiParticleContainer const * const /*mypc*/)
+    {
+        using namespace amrex::literals;
+        amrex::ParmParse pp_collision_name(collision_name);
+        // default Coulomb log, if < 0, will be computed automatically
+        m_CoulombLog = -1.0_rt;
+        queryWithParser(pp_collision_name, "CoulombLog", m_CoulombLog);
+    }
+
+    /**
+     * \brief operator() of the PairWiseCoulombCollisionFunc class. Performs Coulomb collisions
+     * at the cell level by calling ElasticCollisionPerez.
+     *
+     * @param[in] I1s,I2s is the start index for I1,I2 (inclusive).
+     * @param[in] I1e,I2e is the stop index for I1,I2 (exclusive).
+     * @param[in] I1 and I2 are the index arrays. They determine all elements that will be used.
+     * @param[in,out] u1 and u2 are the velocity arrays (u=v*gamma),
+     *                they could be either different or the same,
+     *                their lengths are not needed,
+     * @param[in] w1 and w2 are arrays of weights.
+     * @param[in] q1 and q2 are charges. m1 and m2 are masses.
+     * @param[in] dt is the time step length between two collision calls.
+     * @param[in] dV is the volume of the corresponding cell.
+     * @param[in] engine the random engine.
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_INLINE
+    void operator() (
+        index_type const I1s, index_type const I1e,
+        index_type const I2s, index_type const I2e,
+        index_type* I1,       index_type* I2,
+        SoaData_type soa_1, SoaData_type soa_2,
+        GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/,
+        amrex::Real const  q1, amrex::Real const  q2,
+        amrex::Real const  m1, amrex::Real const  m2,
+        amrex::Real const  dt, amrex::Real const dV,
+        index_type const /*cell_start_pair*/, index_type* /*p_mask*/,
+        index_type* /*p_pair_indices_1*/, index_type* /*p_pair_indices_2*/,
+        amrex::ParticleReal* /*p_pair_reaction_weight*/,
+        amrex::RandomEngine const& engine) const
+        {
+            ElasticCollisionPerez(
+                    I1s, I1e, I2s, I2e, I1, I2,
+                    soa_1, soa_2,
+                    q1, q2, m1, m2, amrex::Real(-1.0), amrex::Real(-1.0),
+                    dt, m_CoulombLog, dV, engine );
+        }
+
+private:
+    amrex::Real m_CoulombLog;
+};
+
+#endif // PAIRWISE_COULOMB_COLLISION_FUNC_H_
diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
new file mode 100644
index 000000000..578f08741
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
@@ -0,0 +1,112 @@
+/* Copyright 2021 Neil Zaim
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef PARTICLE_CREATION_FUNC_H_
+#define PARTICLE_CREATION_FUNC_H_
+
+#include "BinaryCollisionUtils.H"
+
+#include "Particles/ParticleCreation/SmartCopy.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
+#include "Particles/MultiParticleContainer.H"
+#include "Particles/WarpXParticleContainer.H"
+#include "Utils/WarpXUtil.H"
+#include "WarpX.H"
+
+#include <AMReX_DenseBins.H>
+#include <AMReX_GpuContainers.H>
+#include <AMReX_REAL.H>
+#include <AMReX_Vector.H>
+
+/**
+ * \brief This functor creates particles produced from a binary collision and sets their initial
+ * properties (position, momentum, weight).
+ */
+class ParticleCreationFunc{
+    // Define shortcuts for frequently-used type names
+    using ParticleType = WarpXParticleContainer::ParticleType;
+    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using index_type = ParticleBins::index_type;
+    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+
+public:
+    /**
+     * \brief Default constructor of the ParticleCreationFunc class.
+     */
+    ParticleCreationFunc () = default;
+
+    /**
+     * \brief Constructor of the ParticleCreationFunc class
+     *
+     * @param[in] collision_name the name of the collision
+     * @param[in] mypc pointer to the MultiParticleContainer
+     */
+    ParticleCreationFunc (const std::string collision_name, MultiParticleContainer const * const mypc);
+
+    /**
+     * \brief operator() of the ParticleCreationFunc class. Needs to be implemented.
+     */
+    AMREX_INLINE
+    amrex::Vector<int> operator() (
+        const SoaData_type /*soa_1*/, const SoaData_type /*soa_2*/,
+                    const SoaData_type* /*soa_products*/,
+                    const GetParticlePosition /*get_position_1*/,
+                    const GetParticlePosition /*get_position_2*/,
+                    const GetParticlePosition* /*get_position_products*/,
+                    const index_type* /*p_mask*/, const index_type* /*products_np*/,
+                    const SmartCopy* /*copy_species1*/, const SmartCopy* /*copy_species2*/,
+                    const index_type* /*p_pair_indices_1*/, const index_type* /*p_pair_indices_2*/,
+                    const amrex::ParticleReal* /*p_pair_reaction_weight*/
+                    ) const
+        {
+            return amrex::Vector<int>(m_num_product_species, 0);
+        }
+
+private:
+    // How many different type of species the collision produces
+    int m_num_product_species;
+    // Vector of size m_num_product_species storing how many particles of a given species are
+    // produced by a collision event
+    amrex::Gpu::DeviceVector<int> m_num_products;
+    CollisionType m_collision_type;
+};
+
+
+/**
+ * \brief This class does nothing and is used as second template parameter for binary collisions
+ * that do not create particles.
+ */
+class NoParticleCreationFunc{
+    using ParticleType = WarpXParticleContainer::ParticleType;
+    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using index_type = ParticleBins::index_type;
+    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+
+public:
+    NoParticleCreationFunc () = default;
+
+    NoParticleCreationFunc (const std::string /*collision_name*/,
+                            MultiParticleContainer const * const /*mypc*/) {}
+
+    AMREX_INLINE
+    amrex::Vector<int> operator() (
+                    const SoaData_type /*soa_1*/, const SoaData_type /*soa_2*/,
+                    const SoaData_type* /*soa_products*/,
+                    const GetParticlePosition /*get_position_1*/,
+                    const GetParticlePosition /*get_position_2*/,
+                    const GetParticlePosition* /*get_position_products*/,
+                    const index_type* /*p_mask*/, const index_type* /*products_np*/,
+                    const SmartCopy* /*copy_species1*/, const SmartCopy* /*copy_species2*/,
+                    const index_type* /*p_pair_indices_1*/, const index_type* /*p_pair_indices_2*/,
+                    const amrex::ParticleReal* /*p_pair_reaction_weight*/
+                    ) const
+        {
+            return amrex::Vector<int>();
+        }
+};
+
+#endif // PARTICLE_CREATION_FUNC_H_
diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp
new file mode 100644
index 000000000..fff4150c1
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp
@@ -0,0 +1,45 @@
+/* Copyright 2021 Neil Zaim
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "ParticleCreationFunc.H"
+
+#include "BinaryCollisionUtils.H"
+#include "Particles/MultiParticleContainer.H"
+
+#include <AMReX_GpuContainers.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_Vector.H>
+
+#include <string>
+
+ParticleCreationFunc::ParticleCreationFunc (const std::string collision_name,
+                                            MultiParticleContainer const * const mypc)
+    {
+        amrex::ParmParse pp_collision_name(collision_name);
+
+        m_collision_type = BinaryCollisionUtils::get_collision_type(collision_name, mypc);
+
+#ifdef AMREX_USE_GPU
+        amrex::Vector<int> host_num_products;
+#else
+        amrex::Gpu::DeviceVector<int>& host_num_products = m_num_products;
+#endif
+        if (m_collision_type == CollisionType::ProtonBoronFusion)
+            {
+                // Proton-Boron fusion only produces alpha particles
+                m_num_product_species = 1;
+                // Proton-Boron fusion produces 3 alpha particles per fusion reaction
+                host_num_products.push_back(3);
+            }
+#ifdef AMREX_USE_GPU
+        m_num_products.resize(m_num_product_species);
+        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, host_num_products.begin(),
+                              host_num_products.end(),
+                              m_num_products.begin());
+        amrex::Gpu::streamSynchronize();
+#endif
+    }
diff --git a/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H b/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H
new file mode 100644
index 000000000..d7a302e3a
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H
@@ -0,0 +1,36 @@
+/* Copyright 2019 Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_
+#define WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_
+
+#include <AMReX_Random.H>
+
+/* \brief Shuffle array according to Fisher-Yates algorithm.
+ *        Only shuffle the part between is <= i < ie, n = ie-is.
+ *        T_index shall be
+ *        amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type
+*/
+
+template <typename T_index>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void ShuffleFisherYates (T_index *array, T_index const is, T_index const ie,
+                         amrex::RandomEngine const& engine)
+{
+    int     j;
+    T_index buf;
+    for (int i = ie-1; i >= static_cast<int>(is+1); --i)
+    {
+        // get random number j: is <= j <= i
+        j = amrex::Random_int(i-is+1, engine) + is;
+        // swop the ith array element with the jth
+        buf      = array[i];
+        array[i] = array[j];
+        array[j] = buf;
+    }
+}
+
+#endif // WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_
diff --git a/Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H
new file mode 100644
index 000000000..b5acd85a7
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/UpdateMomentumPerezElastic.H
@@ -0,0 +1,271 @@
+/* Copyright 2019 Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
+#define WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_
+
+#include "Utils/WarpXConst.H"
+
+#include <AMReX_Math.H>
+#include <AMReX_Random.H>
+
+#include <cmath>  // isnan() isinf()
+#include <limits> // numeric_limits<float>::min()
+
+/* \brief Update particle velocities according to
+ *        F. Perez et al., Phys.Plasmas.19.083104 (2012),
+ *        which is based on Nanbu's method, PhysRevE.55.4642 (1997).
+ *        @param[in] LmdD is max(Debye length, minimal interparticle distance).
+ *        @param[in] L is the Coulomb log. A fixed L will be used if L > 0,
+ *        otherwise L will be calculated based on the algorithm.
+ *        To see if there are nan or inf updated velocities,
+ *        compile with USE_ASSERTION=TRUE.
+*/
+
+template <typename T_R>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void UpdateMomentumPerezElastic (
+    T_R& u1x, T_R& u1y, T_R& u1z, T_R& u2x, T_R& u2y, T_R& u2z,
+    T_R const n1, T_R const n2, T_R const n12,
+    T_R const q1, T_R const m1, T_R const w1,
+    T_R const q2, T_R const m2, T_R const w2,
+    T_R const dt, T_R const L,  T_R const lmdD,
+    amrex::RandomEngine const& engine)
+{
+
+    // If g = u1 - u2 = 0, do not collide.
+    if ( amrex::Math::abs(u1x-u2x) < std::numeric_limits<T_R>::min() &&
+         amrex::Math::abs(u1y-u2y) < std::numeric_limits<T_R>::min() &&
+         amrex::Math::abs(u1z-u2z) < std::numeric_limits<T_R>::min() )
+    { return; }
+
+    T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c );
+
+    // Compute Lorentz factor gamma
+    T_R const g1 = std::sqrt( T_R(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 );
+    T_R const g2 = std::sqrt( T_R(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 );
+
+    // Compute momenta
+    T_R const p1x = u1x * m1;
+    T_R const p1y = u1y * m1;
+    T_R const p1z = u1z * m1;
+    T_R const p2x = u2x * m2;
+    T_R const p2y = u2y * m2;
+    T_R const p2z = u2z * m2;
+
+    // Compute center-of-mass (COM) velocity and gamma
+    T_R const mass_g = m1 * g1 + m2 * g2;
+    T_R const vcx    = (p1x+p2x) / mass_g;
+    T_R const vcy    = (p1y+p2y) / mass_g;
+    T_R const vcz    = (p1z+p2z) / mass_g;
+    T_R const vcms   = vcx*vcx + vcy*vcy + vcz*vcz;
+    T_R const gc     = T_R(1.0) / std::sqrt( T_R(1.0) - vcms*inv_c2 );
+
+    // Compute vc dot v1 and v2
+    T_R const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1;
+    T_R const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2;
+
+    // Compute p1 star
+    T_R p1sx;
+    T_R p1sy;
+    T_R p1sz;
+    if ( vcms > std::numeric_limits<T_R>::min() )
+    {
+        T_R const lorentz_tansform_factor =
+            ( (gc-T_R(1.0))/vcms*vcDv1 - gc )*m1*g1;
+        p1sx = p1x + vcx*lorentz_tansform_factor;
+        p1sy = p1y + vcy*lorentz_tansform_factor;
+        p1sz = p1z + vcz*lorentz_tansform_factor;
+    }
+    else // If vcms = 0, don't do Lorentz-transform.
+    {
+        p1sx = p1x;
+        p1sy = p1y;
+        p1sz = p1z;
+    }
+    T_R const p1sm = std::sqrt( p1sx*p1sx + p1sy*p1sy + p1sz*p1sz );
+
+    // Compute gamma star
+    T_R const g1s = ( T_R(1.0) - vcDv1*inv_c2 )*gc*g1;
+    T_R const g2s = ( T_R(1.0) - vcDv2*inv_c2 )*gc*g2;
+
+    // Compute the Coulomb log lnLmd
+    T_R lnLmd;
+    if ( L > T_R(0.0) ) { lnLmd = L; }
+    else
+    {
+        // Compute b0
+        T_R const b0 = amrex::Math::abs(q1*q2) * inv_c2 /
+               (T_R(4.0)*MathConst::pi*PhysConst::ep0) * gc/mass_g *
+               ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_R(1.0) );
+
+        // Compute the minimal impact parameter
+        T_R bmin = amrex::max(PhysConst::hbar*MathConst::pi/p1sm,b0);
+
+        // Compute the Coulomb log lnLmd
+        lnLmd = amrex::max( T_R(2.0),
+                T_R(0.5)*std::log(T_R(1.0)+lmdD*lmdD/(bmin*bmin)) );
+    }
+
+    // Compute s
+    const auto tts = m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_R(1.0);
+    const auto tts2 = tts*tts;
+    T_R s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 /
+          ( T_R(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 *
+            m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g * tts2;
+
+    // Compute s'
+    const auto cbrt_n1 = std::cbrt(n1);
+    const auto cbrt_n2 = std::cbrt(n2);
+    const auto coeff = static_cast<T_R>(
+        std::pow(4.0*MathConst::pi/3.0,1.0/3.0));
+    T_R const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc);
+    T_R const sp = coeff * n1*n2/n12 * dt * vrel * (m1+m2) /
+          amrex::max( m1*cbrt_n1*cbrt_n1,
+                      m2*cbrt_n2*cbrt_n2);
+
+    // Determine s
+    s = amrex::min(s,sp);
+
+    // Get random numbers
+    T_R r = amrex::Random(engine);
+
+    // Compute scattering angle
+    T_R cosXs;
+    T_R sinXs;
+    if ( s <= T_R(0.1) )
+    {
+        while ( true )
+        {
+            cosXs = T_R(1.0) + s * std::log(r);
+            // Avoid the bug when r is too small such that cosXs < -1
+            if ( cosXs >= T_R(-1.0) ) { break; }
+            r = amrex::Random(engine);
+        }
+    }
+    else if ( s > T_R(0.1) && s <= T_R(3.0) )
+    {
+        T_R const Ainv = static_cast<T_R>(
+            0.0056958 + 0.9560202*s - 0.508139*s*s +
+            0.47913906*s*s*s - 0.12788975*s*s*s*s + 0.02389567*s*s*s*s*s);
+        cosXs = Ainv * std::log( std::exp(T_R(-1.0)/Ainv) +
+                T_R(2.0) * r * std::sinh(T_R(1.0)/Ainv) );
+    }
+    else if ( s > T_R(3.0) && s <= T_R(6.0) )
+    {
+        T_R const A = T_R(3.0) * std::exp(-s);
+        cosXs = T_R(1.0)/A * std::log( std::exp(-A) +
+                T_R(2.0) * r * std::sinh(A) );
+    }
+    else
+    {
+        cosXs = T_R(2.0) * r - T_R(1.0);
+    }
+    sinXs = std::sqrt(T_R(1.0) - cosXs*cosXs);
+
+    // Get random azimuthal angle
+    T_R const phis = amrex::Random(engine) * T_R(2.0) * MathConst::pi;
+    T_R const cosphis = std::cos(phis);
+    T_R const sinphis = std::sin(phis);
+
+    // Compute post-collision momenta pfs in COM
+    T_R p1fsx;
+    T_R p1fsy;
+    T_R p1fsz;
+    // p1sp is the p1s perpendicular
+    T_R p1sp = std::sqrt( p1sx*p1sx + p1sy*p1sy );
+    // Make sure p1sp is not almost zero
+    if ( p1sp > std::numeric_limits<T_R>::min() )
+    {
+        p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis +
+                ( p1sy*p1sm/p1sp ) * sinXs*sinphis +
+                ( p1sx           ) * cosXs;
+        p1fsy = ( p1sy*p1sz/p1sp ) * sinXs*cosphis +
+                (-p1sx*p1sm/p1sp ) * sinXs*sinphis +
+                ( p1sy           ) * cosXs;
+        p1fsz = (-p1sp           ) * sinXs*cosphis +
+                ( T_R(0.0)       ) * sinXs*sinphis +
+                ( p1sz           ) * cosXs;
+        // Note a negative sign is different from
+        // Eq. (12) in Perez's paper,
+        // but they are the same due to the random nature of phis.
+    }
+    else
+    {
+        // If the previous p1sp is almost zero
+        // x->y  y->z  z->x
+        // This set is equivalent to the one in Nanbu's paper
+        p1sp = std::sqrt( p1sy*p1sy + p1sz*p1sz );
+        p1fsy = ( p1sy*p1sx/p1sp ) * sinXs*cosphis +
+                ( p1sz*p1sm/p1sp ) * sinXs*sinphis +
+                ( p1sy           ) * cosXs;
+        p1fsz = ( p1sz*p1sx/p1sp ) * sinXs*cosphis +
+                (-p1sy*p1sm/p1sp ) * sinXs*sinphis +
+                ( p1sz           ) * cosXs;
+        p1fsx = (-p1sp           ) * sinXs*cosphis +
+                ( T_R(0.0)       ) * sinXs*sinphis +
+                ( p1sx           ) * cosXs;
+    }
+
+    T_R const p2fsx = -p1fsx;
+    T_R const p2fsy = -p1fsy;
+    T_R const p2fsz = -p1fsz;
+
+    // Transform from COM to lab frame
+    T_R p1fx;    T_R p2fx;
+    T_R p1fy;    T_R p2fy;
+    T_R p1fz;    T_R p2fz;
+    if ( vcms > std::numeric_limits<T_R>::min() )
+    {
+        T_R const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz;
+        T_R const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz;
+        T_R const factor = (gc-T_R(1.0))/vcms;
+        T_R const factor1 = factor*vcDp1fs + m1*g1s*gc;
+        T_R const factor2 = factor*vcDp2fs + m2*g2s*gc;
+        p1fx = p1fsx + vcx * factor1;
+        p1fy = p1fsy + vcy * factor1;
+        p1fz = p1fsz + vcz * factor1;
+        p2fx = p2fsx + vcx * factor2;
+        p2fy = p2fsy + vcy * factor2;
+        p2fz = p2fsz + vcz * factor2;
+    }
+    else // If vcms = 0, don't do Lorentz-transform.
+    {
+        p1fx = p1fsx;
+        p1fy = p1fsy;
+        p1fz = p1fsz;
+        p2fx = p2fsx;
+        p2fy = p2fsy;
+        p2fz = p2fsz;
+    }
+
+    // Rejection method
+    r = amrex::Random(engine);
+    if ( w2 > r*amrex::max(w1, w2) )
+    {
+        u1x  = p1fx / m1;
+        u1y  = p1fy / m1;
+        u1z  = p1fz / m1;
+#ifndef AMREX_USE_DPCPP
+        AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z));
+        AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z));
+#endif
+    }
+    r = amrex::Random(engine);
+    if ( w1 > r*amrex::max(w1, w2) )
+    {
+        u2x  = p2fx / m2;
+        u2y  = p2fy / m2;
+        u2z  = p2fz / m2;
+#ifndef AMREX_USE_DPCPP
+        AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z));
+        AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z));
+#endif
+    }
+
+}
+
+#endif // WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_