From e3c049a742db500b32a94d867674fb0f36704da8 Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Thu, 15 Jun 2023 15:38:58 +0200
Subject: [PATCH 01/14] prepare functions

---
 src/particles/collisions/CoulombCollision.H | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/particles/collisions/CoulombCollision.H b/src/particles/collisions/CoulombCollision.H
index 675bbca18b..a3799f790f 100644
--- a/src/particles/collisions/CoulombCollision.H
+++ b/src/particles/collisions/CoulombCollision.H
@@ -8,7 +8,8 @@
 #include <AMReX_ParmParse.H>
 
 /**
- * \brief This class handles Coulomb collisions between 2 plasma species (can be the same).
+ * \brief This class handles Coulomb collisions between 2 particle species
+ * (can be plasma-plasma or beam-plasma, the species can be the same)
  */
 class CoulombCollision
 {

From e67b4192f718113bea71e828df21ad0a0763babf Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Thu, 15 Jun 2023 18:18:37 +0200
Subject: [PATCH 02/14] add basic framework

---
 src/Hipace.H                                  |  12 ++
 src/Hipace.cpp                                |  49 +++++-
 src/particles/beam/MultiBeam.H                |   4 +-
 src/particles/collisions/CoulombCollision.H   |  25 ++-
 src/particles/collisions/CoulombCollision.cpp | 145 +++++++++++++++++-
 src/particles/plasma/MultiPlasma.H            |  27 ++--
 src/particles/plasma/MultiPlasma.cpp          |  36 +----
 7 files changed, 240 insertions(+), 58 deletions(-)

diff --git a/src/Hipace.H b/src/Hipace.H
index 29cfca7d45..d78f0796d9 100644
--- a/src/Hipace.H
+++ b/src/Hipace.H
@@ -162,6 +162,11 @@ public:
     /** \brief reset all laser slices to 0 */
     void ResetLaser ();
 
+    /**
+     * \brief does Coulomb collisions between plasmas and beams
+     */
+    void doCoulombCollision ();
+
     /**
        \brief Returns true on the head rank, otherwise false.
      */
@@ -310,6 +315,8 @@ public:
     inline static amrex::Real m_background_density_SI = 0.;
     /** Time step for the beam evolution */
     amrex::Real m_dt = 0.0;
+    /** Number of binary collisions */
+    inline static int m_ncollisions = 0;
     /** Adaptive time step instance */
     AdaptiveTimeStep m_adaptive_time_step;
     /** Laser instance (soon to be multi laser container) */
@@ -368,6 +375,11 @@ private:
     /** Diagnostics */
     Diagnostic m_diags;
 
+    /** User-input names of the binary collisions to be used */
+    std::vector<std::string> m_collision_names;
+    /** Vector of binary collisions */
+    amrex::Vector< CoulombCollision > m_all_collisions;
+
     /** \brief resizes the diagnostic fab to the correct box in a loop over boxes
      *
      * \param[in] it index of box to be resized to
diff --git a/src/Hipace.cpp b/src/Hipace.cpp
index 79e06e3649..1eb59482ea 100644
--- a/src/Hipace.cpp
+++ b/src/Hipace.cpp
@@ -170,6 +170,18 @@ Hipace::Hipace () :
     MakeGeometry();
 
     m_use_laser = m_multi_laser.m_use_laser;
+
+    queryWithParser(pph, "collisions", m_collision_names);
+    /** Initialize the collision objects */
+    m_ncollisions = m_collision_names.size();
+     for (int i = 0; i < m_ncollisions; ++i) {
+         m_all_collisions.emplace_back(CoulombCollision(m_multi_plasma.m_names, m_multi_beam.m_names, m_collision_names[i]));
+     }
+     if (m_normalized_units && m_ncollisions > 0) {
+         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_background_density_SI!=0,
+             "For collisions with normalized units, a background plasma density must "
+             "be specified via 'hipace.background_density_SI'");
+     }
 }
 
 Hipace::~Hipace ()
@@ -618,8 +630,8 @@ Hipace::SolveOneSlice (int islice, const int islice_local, int step)
                                            m_external_E_uniform, m_external_B_uniform,
                                            m_external_E_slope, m_external_B_slope);
 
-    // collisions for all particles calculated on level 0
-    m_multi_plasma.doCoulombCollision(0, m_slice_geom[0].Domain(), m_slice_geom[0]);
+    // collisions for plasmas and beams
+    doCoulombCollision();
 
     // Advance laser slice by 1 step and store result to 3D array
     // no MR for laser
@@ -1350,6 +1362,39 @@ Hipace::NotifyFinish (const int it, bool only_ghost, bool only_time)
 #endif
 }
 
+void
+Hipace::doCoulombCollision ()
+{
+
+    // collisions for all particles calculated on level 0
+    const int lev = 0;
+
+    for (int i = 0; i < m_ncollisions; ++i)
+    {
+        if (m_all_collisions[i].m_nbeams == 1) {
+            // do beam-plasma collisions
+            auto& species1 = m_multi_beam.m_all_beams[ m_all_collisions[i].m_species1_index ];
+            auto& species2 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species2_index ];
+
+            // TODO: enable tiling
+
+            CoulombCollision::doBeamPlasmaCoulombCollision(
+                lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2,
+                m_all_collisions[i].m_CoulombLog, m_background_density_SI);
+        } else {
+            // do plasma-plasma collisions
+            auto& species1 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species1_index ];
+            auto& species2 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species2_index ];
+
+            // TODO: enable tiling
+
+            CoulombCollision::doPlasmaPlasmaCoulombCollision(
+                lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2, m_all_collisions[i].m_isSameSpecies,
+                m_all_collisions[i].m_CoulombLog, m_background_density_SI);
+        }
+    }
+}
+
 void
 Hipace::ResizeFDiagFAB (const int it, const int step)
 {
diff --git a/src/particles/beam/MultiBeam.H b/src/particles/beam/MultiBeam.H
index ebd1b5a4f9..13d0d2b2cf 100644
--- a/src/particles/beam/MultiBeam.H
+++ b/src/particles/beam/MultiBeam.H
@@ -193,11 +193,11 @@ public:
     /** \brief returns if any beam uses the SALAME algorithm
      */
     bool AnySpeciesSalame ();
+    amrex::Vector<std::string> m_names {"no_beam"}; /**< names of all beam containers */
+    amrex::Vector<BeamParticleContainer> m_all_beams; /**< contains all beam containers */
 
 private:
 
-    amrex::Vector<BeamParticleContainer> m_all_beams; /**< contains all beam containers */
-    amrex::Vector<std::string> m_names {"no_beam"}; /**< names of all beam containers */
     int m_nbeams {0}; /**< number of beam containers */
     /** number of real particles per beam, as opposed to ghost particles */
     amrex::Vector<amrex::Long> m_n_real_particles;
diff --git a/src/particles/collisions/CoulombCollision.H b/src/particles/collisions/CoulombCollision.H
index a3799f790f..ab26435a52 100644
--- a/src/particles/collisions/CoulombCollision.H
+++ b/src/particles/collisions/CoulombCollision.H
@@ -2,6 +2,7 @@
 #define HIPACE_COULOMB_COLLISION_H_
 
 #include "particles/plasma/PlasmaParticleContainer.H"
+#include "particles/beam/BeamParticleContainer.H"
 
 #include <AMReX_DenseBins.H>
 #include <AMReX_REAL.H>
@@ -16,12 +17,14 @@ class CoulombCollision
 public:
     int  m_species1_index {-1};
     int  m_species2_index {-1};
+    int  m_nbeams {0};
     bool m_isSameSpecies {false};
     amrex::Real m_CoulombLog {-1.};
 
     /** Constructor */
     CoulombCollision(
-        const std::vector<std::string>& species_names,
+        const std::vector<std::string>& plasma_species_names,
+        const std::vector<std::string>& beam_species_names,
         std::string const collision_name);
 
     /**
@@ -38,11 +41,29 @@ public:
      *            Coulomb logarithm is deduced from the plasma temperature, measured in each cell.
      * \param[in] background_density_SI background plasma density (only needed for normalized units)
      **/
-    static void doCoulombCollision (
+    static void doPlasmaPlasmaCoulombCollision (
         int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1,
         PlasmaParticleContainer& species2, bool is_same_species, amrex::Real CoulombLog,
         amrex::Real background_density_SI);
 
+    /**
+     * \brief Perform Coulomb collisions of a beam with a plasma species over a push by one beam time step
+     *        Particles of both species are sorted per cell, paired, and collided pairwise.
+     *
+     * \param[in] lev MR level
+     * \param[in] bx transverse box (plasma particles will be sorted per-cell on this box)
+     * \param[in] geom corresponding geometry object
+     * \param[in,out] species1 beam species
+     * \param[in,out] species2 plasma species
+     * \param[in] CoulombLog Value of the Coulomb logarithm used for the collisions. If <0, the
+     *            Coulomb logarithm is deduced from the plasma temperature, measured in each cell.
+     * \param[in] background_density_SI background plasma density (only needed for normalized units)
+     **/
+    static void doBeamPlasmaCoulombCollision (
+        int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1,
+        PlasmaParticleContainer& species2, amrex::Real CoulombLog,
+        amrex::Real background_density_SI);
+
 };
 
 #endif // HIPACE_COULOMB_COLLISION_H_
diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp
index dd905672a1..20a7b50e62 100644
--- a/src/particles/collisions/CoulombCollision.cpp
+++ b/src/particles/collisions/CoulombCollision.cpp
@@ -6,7 +6,8 @@
 #include "utils/HipaceProfilerWrapper.H"
 
 CoulombCollision::CoulombCollision(
-    const std::vector<std::string>& species_names,
+    const std::vector<std::string>& plasma_species_names,
+    const std::vector<std::string>& beam_species_names,
     std::string const collision_name)
 {
     using namespace amrex::literals;
@@ -24,20 +25,39 @@ CoulombCollision::CoulombCollision(
     // default Coulomb log is -1, if < 0 (e.g. not specified), will be computed automatically
     pp.query("CoulombLog", m_CoulombLog);
 
-    for (int i=0; i<(int) species_names.size(); i++)
-    {
-        if (species_names[i] == collision_species[0]) m_species1_index = i;
-        if (species_names[i] == collision_species[1]) m_species2_index = i;
+    for (int i=0; i<(int) beam_species_names.size(); i++) {
+        if (beam_species_names[i] == collision_species[0]) m_nbeams += 1;
+        if (beam_species_names[i] == collision_species[1]) m_nbeams += 1;
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+                m_nbeams <= 1,
+                "<collision name>.species must contain at maximum one beam species name!");
     }
+    if (m_nbeams == 1) {
+        for (int i=0; i<(int) beam_species_names.size(); i++) {
+            if (beam_species_names[i] == collision_species[0] ||
+                beam_species_names[i] == collision_species[1]) m_species1_index = i;
+        }
+        for (int i=0; i<(int) plasma_species_names.size(); i++) {
+            if (plasma_species_names[i] == collision_species[0] ||
+                plasma_species_names[i] == collision_species[1]) m_species2_index = i;
+        }
+    } else {
+        for (int i=0; i<(int) plasma_species_names.size(); i++) {
+            if (plasma_species_names[i] == collision_species[0]) m_species1_index = i;
+            if (plasma_species_names[i] == collision_species[1]) m_species2_index = i;
+        }
+    }
+
     m_isSameSpecies = collision_species[0] == collision_species[1] ? true : false;
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         m_species1_index >= 0 && m_species2_index >= 0,
-        "<collision name>.species must contain exactly the name of 2 species in plasma.names"
+        "<collision name>.species must contain exactly the name of 2 species in plasmas.names "
+        "or one in plasmas.names and one in beams.names"
         );
 }
 
 void
-CoulombCollision::doCoulombCollision (
+CoulombCollision::doPlasmaPlasmaCoulombCollision (
     int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1,
     PlasmaParticleContainer& species2, bool is_same_species, amrex::Real CoulombLog,
     amrex::Real background_density_SI)
@@ -213,3 +233,114 @@ CoulombCollision::doCoulombCollision (
         AMREX_ALWAYS_ASSERT(count == 1);
     }
 }
+
+void
+CoulombCollision::doBeamPlasmaCoulombCollision (
+    int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1,
+    PlasmaParticleContainer& species2, amrex::Real CoulombLog,
+    amrex::Real background_density_SI)
+{
+    HIPACE_PROFILE("CoulombCollision::doBeamPlasmaCoulombCollision()");
+    AMREX_ALWAYS_ASSERT(lev == 0);
+
+    if (species1.TotalNumberOfParticles() == 0 || species2.TotalNumberOfParticles() == 0) return;
+
+    using namespace amrex::literals;
+    const PhysConst cst = get_phys_const();
+    bool normalized_units = Hipace::m_normalized_units;
+
+    // const amrex::Real clight = cst.c;
+    // const amrex::Real inv_c = 1.0_rt / cst.c;
+    // const amrex::Real inv_c2 = 1.0_rt / ( cst.c * cst.c );
+    // constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c;
+    // constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c );
+
+    // Logically particles per-cell, and return indices of particles in each cell
+    // PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom);
+    PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom);
+
+    int const n_cells = bins2.numBins();
+
+    // Counter to check there is only 1 box
+    int count = 0;
+    for (PlasmaParticleIterator pti(species2); pti.isValid(); ++pti) {
+
+        // // Get particles SoA data for species 1
+        // auto& soa1 = pti.GetStructOfArrays();
+        // amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux_half_step).data();
+        // amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy_half_step).data();
+        // amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi_half_step).data();
+        // const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data();
+        // const int* const ion_lev1 = soa1.GetIntData(PlasmaIdx::ion_lev).data();
+        // PlasmaBins::index_type * const indices1 = bins1.permutationPtr();
+        // PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr();
+        // amrex::Real q1 = species1.GetCharge();
+        // amrex::Real m1 = species1.GetMass();
+        // const bool can_ionize1 = species1.m_can_ionize;
+
+        // Get particles SoA data for species 2
+        auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex());
+        auto& soa2 = ptile2.GetStructOfArrays();
+        amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux_half_step).data();
+        amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy_half_step).data();
+        amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi_half_step).data();
+        const amrex::Real* const w2 = soa2.GetRealData(PlasmaIdx::w).data();
+        const int* const ion_lev2 = soa2.GetIntData(PlasmaIdx::ion_lev).data();
+        PlasmaBins::index_type * const indices2 = bins2.permutationPtr();
+        PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr();
+        amrex::Real q2 = species2.GetCharge();
+        amrex::Real m2 = species2.GetMass();
+        const bool can_ionize2 = species2.m_can_ionize;
+
+        // volume is used to calculate density, but weights already represent density in normalized units
+        const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2);
+        // static_cast<double> to avoid precision problems in FP32
+        const amrex::Real wp = std::sqrt(static_cast<double>(background_density_SI) *
+                                         PhysConstSI::q_e*PhysConstSI::q_e /
+                                         (PhysConstSI::ep0*PhysConstSI::m_e));
+        const amrex::Real dt = normalized_units ? geom.CellSize(2)/wp
+                                                : geom.CellSize(2)/PhysConstSI::c;
+        // 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);
+        // Loop over cells, and collide the particles in each cell
+
+        // 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]`
+                // PlasmaBins::index_type const cell_start1 = offsets1[i_cell];
+                // PlasmaBins::index_type const cell_stop1  = offsets1[i_cell+1];
+                // // Same for species 2
+                // PlasmaBins::index_type const cell_start2 = offsets2[i_cell];
+                // PlasmaBins::index_type const cell_stop2  = offsets2[i_cell+1];
+                //
+                // // 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_stop1 - cell_start1 < 1 ||
+                //      cell_stop2 - cell_start2 < 1 ) return;
+                // // shuffle
+                // ShuffleFisherYates(indices1, cell_start1, cell_stop1, engine);
+                // ShuffleFisherYates(indices2, cell_start2, cell_stop2, engine);
+                //
+                // // TODO: FIX DT.
+                // // Call the function in order to perform collisions
+                // ElasticCollisionPerez(
+                //     cell_start1, cell_stop1, cell_start2, cell_stop2,
+                //     indices1, indices2,
+                //     ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev1, ion_lev2,
+                //     q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2,
+                //     dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI,
+                //     normalized_units, background_density_SI, is_same_species, engine );
+            }
+            );
+        count++;
+    }
+    AMREX_ALWAYS_ASSERT(count == 1);
+}
diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H
index 2ac504ea80..f22c8ede07 100644
--- a/src/particles/plasma/MultiPlasma.H
+++ b/src/particles/plasma/MultiPlasma.H
@@ -122,15 +122,15 @@ public:
     /** returns number of plasma species */
     int GetNPlasmas() const {return m_nplasmas;}
 
-    /** Perform binary elastic Coulomb collision, inter- and/or intra-species.
-     *
-     * The algorithm implemented is that of [Perez et al., Phys. Plasmas 19, 083104 (2012)]
-     *
-     * \param[in] lev MR level
-     * \param[in] bx box on which plasma particles are sorted per-cell and collide together
-     * \param[in] geom Corresponding gemetry
-     */
-    void doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom);
+    // /** Perform binary elastic Coulomb collision, inter- and/or intra-species.
+    //  *
+    //  * The algorithm implemented is that of [Perez et al., Phys. Plasmas 19, 083104 (2012)]
+    //  *
+    //  * \param[in] lev MR level
+    //  * \param[in] bx box on which plasma particles are sorted per-cell and collide together
+    //  * \param[in] geom Corresponding gemetry
+    //  */
+    // void doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom);
 
     /** Reorder particles to speed-up current deposition
      * \param[in] islice zeta slice index
@@ -166,22 +166,15 @@ public:
 
     int m_sort_bin_size {32}; /**< Tile size to sort plasma particles */
 
-    /** Number of binary collisions */
-    int m_ncollisions = 0;
-
     amrex::Vector<PlasmaParticleContainer> m_all_plasmas; /**< contains all plasma containers */
     int m_nplasmas = 0; /**< number of plasma containers */
 
     amrex::Vector<PlasmaBins> m_all_bins; /**< Logical tile bins for all plasma containers */
-private:
     amrex::Vector<std::string> m_names {"no_plasma"}; /**< names of all plasma containers */
+private:
     /** Background (hypothetical) density, used to compute the adaptive time step */
     amrex::Real m_adaptive_density = 0.;
 
-    /** User-input names of the binary collisions to be used */
-    std::vector<std::string> m_collision_names;
-    /** Vector of binary collisions */
-    amrex::Vector< CoulombCollision > m_all_collisions;
 };
 
 #endif // MULTIPLASMA_H_
diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp
index 46233015a2..2984a13bb6 100644
--- a/src/particles/plasma/MultiPlasma.cpp
+++ b/src/particles/plasma/MultiPlasma.cpp
@@ -21,8 +21,9 @@ MultiPlasma::MultiPlasma ()
     queryWithParser(pp, "names", m_names);
     queryWithParser(pp, "adaptive_density", m_adaptive_density);
     queryWithParser(pp, "sort_bin_size", m_sort_bin_size);
-    queryWithParser(pp, "collisions", m_collision_names);
 
+    DeprecatedInput ("plasmas", "collisions",
+                    "hipace.collisions", "", true);
     DeprecatedInput ("plasmas", "background_density_SI",
                     "hipace.background_density_SI", "", true);
 
@@ -33,16 +34,6 @@ MultiPlasma::MultiPlasma ()
         m_all_plasmas.emplace_back(PlasmaParticleContainer(m_names[i]));
     }
 
-    /** Initialize the collision objects */
-    m_ncollisions = m_collision_names.size();
-     for (int i = 0; i < m_ncollisions; ++i) {
-         m_all_collisions.emplace_back(CoulombCollision(m_names, m_collision_names[i]));
-     }
-     if (Hipace::m_normalized_units && m_ncollisions > 0) {
-         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(Hipace::m_background_density_SI!=0,
-             "For collisions with normalized units, a background plasma density must "
-             "be specified via 'hipace.background_density_SI'");
-     }
 }
 
 void
@@ -168,23 +159,12 @@ MultiPlasma::TileSort (amrex::Box bx, amrex::Geometry geom)
     }
 }
 
-void
-MultiPlasma::doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom)
-{
-    HIPACE_PROFILE("MultiPlasma::doCoulombCollision()");
-    for (int i = 0; i < m_ncollisions; ++i)
-    {
-        AMREX_ALWAYS_ASSERT(lev == 0);
-        auto& species1 = m_all_plasmas[ m_all_collisions[i].m_species1_index ];
-        auto& species2 = m_all_plasmas[ m_all_collisions[i].m_species2_index ];
-
-        // TODO: enable tiling
-
-        CoulombCollision::doCoulombCollision(
-            lev, bx, geom, species1, species2, m_all_collisions[i].m_isSameSpecies,
-            m_all_collisions[i].m_CoulombLog, Hipace::m_background_density_SI);
-    }
-}
+// void
+// MultiPlasma::doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom)
+// {
+//     HIPACE_PROFILE("MultiPlasma::doCoulombCollision()");
+//
+// }
 
 void
 MultiPlasma::ReorderParticles (const int islice)

From d82ed72a1657d5e374c39b287de30fd881af5783 Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Thu, 15 Jun 2023 21:08:18 +0200
Subject: [PATCH 03/14] add transverse beam particle binning per slice

---
 src/Hipace.H                                  |  3 +-
 src/Hipace.cpp                                |  6 +--
 src/particles/collisions/CoulombCollision.H   |  5 +-
 src/particles/collisions/CoulombCollision.cpp |  6 +--
 src/particles/sorting/TileSort.H              | 16 ++++++
 src/particles/sorting/TileSort.cpp            | 52 +++++++++++++++++++
 6 files changed, 79 insertions(+), 9 deletions(-)

diff --git a/src/Hipace.H b/src/Hipace.H
index d78f0796d9..367f3d1bb0 100644
--- a/src/Hipace.H
+++ b/src/Hipace.H
@@ -164,8 +164,9 @@ public:
 
     /**
      * \brief does Coulomb collisions between plasmas and beams
+     * \param[in] islice_local slice on which collisions are calculated, needed for beam binning
      */
-    void doCoulombCollision ();
+    void doCoulombCollision (const int islice_local);
 
     /**
        \brief Returns true on the head rank, otherwise false.
diff --git a/src/Hipace.cpp b/src/Hipace.cpp
index 1eb59482ea..552288211e 100644
--- a/src/Hipace.cpp
+++ b/src/Hipace.cpp
@@ -631,7 +631,7 @@ Hipace::SolveOneSlice (int islice, const int islice_local, int step)
                                            m_external_E_slope, m_external_B_slope);
 
     // collisions for plasmas and beams
-    doCoulombCollision();
+    doCoulombCollision(islice_local);
 
     // Advance laser slice by 1 step and store result to 3D array
     // no MR for laser
@@ -1363,7 +1363,7 @@ Hipace::NotifyFinish (const int it, bool only_ghost, bool only_time)
 }
 
 void
-Hipace::doCoulombCollision ()
+Hipace::doCoulombCollision (const int islice_local)
 {
 
     // collisions for all particles calculated on level 0
@@ -1379,7 +1379,7 @@ Hipace::doCoulombCollision ()
             // TODO: enable tiling
 
             CoulombCollision::doBeamPlasmaCoulombCollision(
-                lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2,
+                lev, m_slice_geom[0].Domain(), m_slice_geom[0], islice_local, species1, species2,
                 m_all_collisions[i].m_CoulombLog, m_background_density_SI);
         } else {
             // do plasma-plasma collisions
diff --git a/src/particles/collisions/CoulombCollision.H b/src/particles/collisions/CoulombCollision.H
index ab26435a52..4c38d85415 100644
--- a/src/particles/collisions/CoulombCollision.H
+++ b/src/particles/collisions/CoulombCollision.H
@@ -53,6 +53,7 @@ public:
      * \param[in] lev MR level
      * \param[in] bx transverse box (plasma particles will be sorted per-cell on this box)
      * \param[in] geom corresponding geometry object
+     * \param[in] islice_local local slice on which collisions are calculated
      * \param[in,out] species1 beam species
      * \param[in,out] species2 plasma species
      * \param[in] CoulombLog Value of the Coulomb logarithm used for the collisions. If <0, the
@@ -60,8 +61,8 @@ public:
      * \param[in] background_density_SI background plasma density (only needed for normalized units)
      **/
     static void doBeamPlasmaCoulombCollision (
-        int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1,
-        PlasmaParticleContainer& species2, amrex::Real CoulombLog,
+        int lev, const amrex::Box& bx, const amrex::Geometry& geom, int islice_local,
+        BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog,
         amrex::Real background_density_SI);
 
 };
diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp
index 20a7b50e62..6a6492d68a 100644
--- a/src/particles/collisions/CoulombCollision.cpp
+++ b/src/particles/collisions/CoulombCollision.cpp
@@ -236,8 +236,8 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
 
 void
 CoulombCollision::doBeamPlasmaCoulombCollision (
-    int lev, const amrex::Box& bx, const amrex::Geometry& geom, BeamParticleContainer& species1,
-    PlasmaParticleContainer& species2, amrex::Real CoulombLog,
+    int lev, const amrex::Box& bx, const amrex::Geometry& geom, int islice_local,
+    BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog,
     amrex::Real background_density_SI)
 {
     HIPACE_PROFILE("CoulombCollision::doBeamPlasmaCoulombCollision()");
@@ -256,7 +256,7 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
     // constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c );
 
     // Logically particles per-cell, and return indices of particles in each cell
-    // PlasmaBins bins1 = findParticlesInEachTile(bx, 1, species1, geom);
+    BeamBins bins1 = findBeamParticlesInEachTile(bx, 1, islice_local, species1, geom);
     PlasmaBins bins2 = findParticlesInEachTile(bx, 1, species2, geom);
 
     int const n_cells = bins2.numBins();
diff --git a/src/particles/sorting/TileSort.H b/src/particles/sorting/TileSort.H
index 9137806881..e7254d2467 100644
--- a/src/particles/sorting/TileSort.H
+++ b/src/particles/sorting/TileSort.H
@@ -9,6 +9,7 @@
 #define HIPACE_TILESORT_H_
 
 #include "particles/plasma/PlasmaParticleContainer.H"
+#include "particles/beam/BeamParticleContainer.H"
 
 #include <AMReX_MultiFab.H>
 
@@ -28,4 +29,19 @@ findParticlesInEachTile (
     amrex::Box bx, int bin_size,
     PlasmaParticleContainer& plasma, const amrex::Geometry& geom);
 
+/** \brief Find beam particles in each bin, and return collections of indices per bin (tile).
+ *
+ * Note that this does *not* rearrange particle arrays
+ *
+ * \param[in] bx 3d box in which particles are sorted per slice
+ * \param[in] bin_size number of cells per tile (square)
+ * \param[in] islice_local slice of which the beam particles should be sorted
+ * \param[in] beam beam particle container
+ * \param[in] geom Geometry
+ */
+BeamBins
+findBeamParticlesInEachTile (
+    amrex::Box bx, int bin_size, int islice_local,
+    BeamParticleContainer& beam, const amrex::Geometry& geom);
+
 #endif // HIPACE_TILESORT_H_
diff --git a/src/particles/sorting/TileSort.cpp b/src/particles/sorting/TileSort.cpp
index f1abb2bf29..85d16343f4 100644
--- a/src/particles/sorting/TileSort.cpp
+++ b/src/particles/sorting/TileSort.cpp
@@ -53,3 +53,55 @@ findParticlesInEachTile (
     AMREX_ALWAYS_ASSERT(count <= 1);
     return bins;
 }
+
+BeamBins
+findBeamParticlesInEachTile (
+    amrex::Box bx, int bin_size, int islice_local,
+    BeamParticleContainer& beam, const amrex::Geometry& geom)
+{
+    HIPACE_PROFILE("findBeamParticlesInEachTile()");
+
+    // Tile box: only 1 cell longitudinally, same as bx transversally.
+    const amrex::Box tcbx = bx.coarsen(bin_size);
+    const amrex::Box cbx = {{tcbx.smallEnd(0),tcbx.smallEnd(1),0}, {tcbx.bigEnd(0),tcbx.bigEnd(1),0}};
+
+    BeamBins bins;
+
+    const int box_offset = beam.m_box_sorter.boxOffsetsPtr()[beam.m_ibox];
+
+    BeamBins::index_type const * const indices = beam.m_slice_bins.permutationPtr();
+    BeamBins::index_type const * const slice_offsets = beam.m_slice_bins.offsetsPtrCpu();
+    BeamBins::index_type const
+        cell_start = slice_offsets[islice_local], cell_stop = slice_offsets[islice_local+1];
+    // The particles that are in slice islice_local are
+    // given by the indices[cell_start:cell_stop]
+
+    int const num_particles = cell_stop-cell_start;
+
+    // Extract box properties
+    const auto lo = lbound(cbx);
+    const auto dxi = amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>({
+            geom.InvCellSizeArray()[0]/bin_size,
+            geom.InvCellSizeArray()[1]/bin_size,
+            1.});
+    const auto plo = geom.ProbLoArray();
+
+    // Find the particles that are in each slice and return collections of indices per slice.
+    bins.build(
+        num_particles, beam.getParticleTileData(), cbx,
+        // Pass lambda function that returns the slice index
+        [=] AMREX_GPU_DEVICE (const BeamParticleContainer::ParticleTileDataType& pdt,
+                              unsigned int idx)
+        noexcept -> amrex::IntVect
+        {
+            auto p = pdt[indices[cell_start+idx] + box_offset];
+
+            return amrex::IntVect(
+                AMREX_D_DECL(
+                    static_cast<int>((p.pos(0)-plo[0])*dxi[0]-lo.x),
+                    static_cast<int>((p.pos(1)-plo[1])*dxi[1]-lo.y),
+                    0));
+        });
+
+    return bins;
+}

From 296c6208925ab9958857753d5865249f013f8f6f Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Thu, 15 Jun 2023 22:19:51 +0200
Subject: [PATCH 04/14] add beam collisions

---
 src/particles/beam/BeamParticleContainer.H    |   4 +
 src/particles/collisions/CoulombCollision.cpp | 101 +++++++++---------
 .../collisions/ElasticCollisionPerez.H        |  22 ++--
 3 files changed, 65 insertions(+), 62 deletions(-)

diff --git a/src/particles/beam/BeamParticleContainer.H b/src/particles/beam/BeamParticleContainer.H
index 0e87f192ed..82888fb7c7 100644
--- a/src/particles/beam/BeamParticleContainer.H
+++ b/src/particles/beam/BeamParticleContainer.H
@@ -157,6 +157,10 @@ public:
     std::string get_name () const {return m_name;}
     amrex::Real m_charge; /**< charge of each particle of this species */
     amrex::Real m_mass; /**< mass of each particle of this species */
+    /** Returns elementary charge q_e (or -q_e for electrons). */
+    amrex::Real GetCharge () const {return m_charge;}
+    /** Returns mass of physical species */
+    amrex::Real GetMass () const {return m_mass;}
     bool m_do_z_push {true}; /**< Pushing beam particles in z direction */
     int m_n_subcycles {10}; /**< Number of sub-cycles in the beam pusher */
     bool m_do_radiation_reaction {false}; /**< whether to calculate radiation losses */
diff --git a/src/particles/collisions/CoulombCollision.cpp b/src/particles/collisions/CoulombCollision.cpp
index 6a6492d68a..b71a4a52eb 100644
--- a/src/particles/collisions/CoulombCollision.cpp
+++ b/src/particles/collisions/CoulombCollision.cpp
@@ -134,7 +134,7 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
                         ux1, uy1, psi1, ux1, uy1, psi1, w1, w1, ion_lev1, ion_lev1,
                         q1, q1, m1, m1, -1.0_rt, -1.0_rt, can_ionize1, can_ionize1,
                         dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI,
-                        normalized_units, background_density_SI, is_same_species, engine );
+                        normalized_units, background_density_SI, is_same_species, false, engine );
                 }
                 );
             count++;
@@ -225,7 +225,7 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
                         ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev1, ion_lev2,
                         q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2,
                         dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI,
-                        normalized_units, background_density_SI, is_same_species, engine );
+                        normalized_units, background_density_SI, is_same_species, false, engine );
                 }
                 );
             count++;
@@ -249,11 +249,11 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
     const PhysConst cst = get_phys_const();
     bool normalized_units = Hipace::m_normalized_units;
 
-    // const amrex::Real clight = cst.c;
-    // const amrex::Real inv_c = 1.0_rt / cst.c;
-    // const amrex::Real inv_c2 = 1.0_rt / ( cst.c * cst.c );
-    // constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c;
-    // constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c );
+    const amrex::Real clight = cst.c;
+    const amrex::Real inv_c = 1.0_rt / cst.c;
+    const amrex::Real inv_c2 = 1.0_rt / ( cst.c * cst.c );
+    constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c;
+    constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c );
 
     // Logically particles per-cell, and return indices of particles in each cell
     BeamBins bins1 = findBeamParticlesInEachTile(bx, 1, islice_local, species1, geom);
@@ -266,21 +266,21 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
     for (PlasmaParticleIterator pti(species2); pti.isValid(); ++pti) {
 
         // // Get particles SoA data for species 1
-        // auto& soa1 = pti.GetStructOfArrays();
-        // amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux_half_step).data();
-        // amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy_half_step).data();
-        // amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi_half_step).data();
-        // const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data();
-        // const int* const ion_lev1 = soa1.GetIntData(PlasmaIdx::ion_lev).data();
-        // PlasmaBins::index_type * const indices1 = bins1.permutationPtr();
-        // PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr();
-        // amrex::Real q1 = species1.GetCharge();
-        // amrex::Real m1 = species1.GetMass();
-        // const bool can_ionize1 = species1.m_can_ionize;
+        const int box_offset = species1.m_box_sorter.boxOffsetsPtr()[species1.m_ibox];
+        auto& soa1 = species1.GetStructOfArrays();
+        amrex::Real* const ux1 = soa1.GetRealData(BeamIdx::ux).data() + box_offset;
+        amrex::Real* const uy1 = soa1.GetRealData(BeamIdx::uy).data() + box_offset;
+        amrex::Real* const psi1 = soa1.GetRealData(BeamIdx::uz).data() + box_offset;
+        const amrex::Real* const w1 = soa1.GetRealData(BeamIdx::w).data() + box_offset;
+        BeamBins::index_type * const indices1 = bins1.permutationPtr();
+        BeamBins::index_type const * const offsets1 = bins1.offsetsPtr();
+        amrex::Real q1 = species1.GetCharge();
+        amrex::Real m1 = species1.GetMass();
+        constexpr bool can_ionize1 = false;
 
         // Get particles SoA data for species 2
-        auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex());
-        auto& soa2 = ptile2.GetStructOfArrays();
+        //auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex());
+        auto& soa2 = pti.GetStructOfArrays();
         amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux_half_step).data();
         amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy_half_step).data();
         amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi_half_step).data();
@@ -298,8 +298,9 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
         const amrex::Real wp = std::sqrt(static_cast<double>(background_density_SI) *
                                          PhysConstSI::q_e*PhysConstSI::q_e /
                                          (PhysConstSI::ep0*PhysConstSI::m_e));
-        const amrex::Real dt = normalized_units ? geom.CellSize(2)/wp
-                                                : geom.CellSize(2)/PhysConstSI::c;
+        const amrex::Real dt = normalized_units ? Hipace::GetInstance().m_dt/wp
+                                                : Hipace::GetInstance().m_dt;
+
         // 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);
@@ -310,34 +311,34 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
             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]`
-                // PlasmaBins::index_type const cell_start1 = offsets1[i_cell];
-                // PlasmaBins::index_type const cell_stop1  = offsets1[i_cell+1];
-                // // Same for species 2
-                // PlasmaBins::index_type const cell_start2 = offsets2[i_cell];
-                // PlasmaBins::index_type const cell_stop2  = offsets2[i_cell+1];
-                //
-                // // 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_stop1 - cell_start1 < 1 ||
-                //      cell_stop2 - cell_start2 < 1 ) return;
-                // // shuffle
-                // ShuffleFisherYates(indices1, cell_start1, cell_stop1, engine);
-                // ShuffleFisherYates(indices2, cell_start2, cell_stop2, engine);
-                //
-                // // TODO: FIX DT.
-                // // Call the function in order to perform collisions
-                // ElasticCollisionPerez(
-                //     cell_start1, cell_stop1, cell_start2, cell_stop2,
-                //     indices1, indices2,
-                //     ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev1, ion_lev2,
-                //     q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2,
-                //     dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI,
-                //     normalized_units, background_density_SI, is_same_species, engine );
+                // The particles from species1 that are in the cell `i_cell` are
+                // given by the `indices_1[cell_start_1:cell_stop_1]`
+                BeamBins::index_type const cell_start1 = offsets1[i_cell];
+                BeamBins::index_type const cell_stop1  = offsets1[i_cell+1];
+                // Same for species 2
+                PlasmaBins::index_type const cell_start2 = offsets2[i_cell];
+                PlasmaBins::index_type const cell_stop2  = offsets2[i_cell+1];
+
+                // 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_stop1 - cell_start1 < 1 ||
+                     cell_stop2 - cell_start2 < 1 ) return;
+                // shuffle
+                ShuffleFisherYates(indices1, cell_start1, cell_stop1, engine);
+                ShuffleFisherYates(indices2, cell_start2, cell_stop2, engine);
+
+                // TODO: FIX DT.
+                // Call the function in order to perform collisions
+                ElasticCollisionPerez(
+                    cell_start1, cell_stop1, cell_start2, cell_stop2,
+                    indices1, indices2,
+                    ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev2, ion_lev2, // passing ion_lev2 for beam particles, will never be used
+                    q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2,
+                    dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI,
+                    normalized_units, background_density_SI, false, true, engine );
             }
             );
         count++;
diff --git a/src/particles/collisions/ElasticCollisionPerez.H b/src/particles/collisions/ElasticCollisionPerez.H
index 0dfa68c581..36cb1d1a54 100644
--- a/src/particles/collisions/ElasticCollisionPerez.H
+++ b/src/particles/collisions/ElasticCollisionPerez.H
@@ -69,7 +69,7 @@ void ElasticCollisionPerez (
     const bool can_ionize1, const bool can_ionize2,
     T_R const  dt, T_R const   L, T_R const inv_dV, T_R const clight, T_R const inv_c,
     T_R const inv_c_SI, T_R const inv_c2, T_R const inv_c2_SI, const bool normalized_units,
-    const amrex::Real background_density_SI, const bool is_same_species,
+    const amrex::Real background_density_SI, const bool is_same_species, bool is_beam_coll,
     amrex::RandomEngine const& engine)
 {
     using namespace amrex::literals;
@@ -151,24 +151,22 @@ void ElasticCollisionPerez (
             if (can_ionize1) q1 *= ion_lev1[I1[i1]];
             if (can_ionize2) q2 *= ion_lev2[I2[i2]];
             // particle's Lorentz factor
-            amrex::Real g1 = (
-                1.0_rt +
-                u1x[I1[i1]]*u1x[I1[i1]]*inv_c2 + u1y[I1[i1]]*u1y[I1[i1]]*inv_c2 +
-                psi1[I1[i1]]*psi1[I1[i1]]) / (2.0_rt * psi1[I1[i1]] );
+            amrex::Real g1 = is_beam_coll ? std::sqrt( 1._rt
+                + (u1x[I1[i1]]*u1x[I1[i1]] + u1y[I1[i1]]*u1y[I1[i1]] + psi1[I1[i1]]*psi1[I1[i1]])*inv_c2 )
+                : (1.0_rt + u1x[I1[i1]]*u1x[I1[i1]]*inv_c2 + u1y[I1[i1]]*u1y[I1[i1]]*inv_c2 +
+                   psi1[I1[i1]]*psi1[I1[i1]]) / (2.0_rt * psi1[I1[i1]] );
             // particle's Lorentz factor
-            amrex::Real g2 = (
-                1.0_rt +
-                u2x[I2[i2]]*u2x[I2[i2]]*inv_c2 + u2y[I2[i2]]*u2y[I2[i2]]*inv_c2 +
-                psi2[I2[i2]]*psi2[I2[i2]]) / (2.0_rt * psi2[I2[i2]] );
+            amrex::Real g2 = (1.0_rt + u2x[I2[i2]]*u2x[I2[i2]]*inv_c2 + u2y[I2[i2]]*u2y[I2[i2]]*inv_c2 +
+                   psi2[I2[i2]]*psi2[I2[i2]]) / (2.0_rt * psi2[I2[i2]] );
 
             // Convert from pseudo-potential to momentum
-            amrex::Real u1z = clight * (g1 - psi1[I1[i1]]);
+            amrex::Real u1z = is_beam_coll ? psi1[I1[i1]] : clight * (g1 - psi1[I1[i1]]);
             amrex::Real u2z = clight * (g2 - psi2[I2[i2]]);
 
             // In the longitudinal push of plasma particles, the dt is different for each particle.
             // The dt applied for collision probability is the average (in the lab frame) of these
             // dts. This is NOT clean. TODO FIXME.
-            const amrex::Real dt_fac = 0.5_rt * (g1/psi1[I1[i1]] + g2/psi2[I2[i2]]);
+            const amrex::Real dt_fac = is_beam_coll ? 1.0_rt : 0.5_rt * (g1/psi1[I1[i1]] + g2/psi2[I2[i2]]);
             UpdateMomentumPerezElastic(
                 u1x[ I1[i1] ], u1y[ I1[i1] ], u1z, g1,
                 u2x[ I2[i2] ], u2y[ I2[i2] ], u2z, g2,
@@ -177,7 +175,7 @@ void ElasticCollisionPerez (
 
             g1 = std::sqrt( T_R(1.0) +
                             (u1x[I1[i1]]*u1x[I1[i1]]+u1y[I1[i1]]*u1y[I1[i1]]+u1z*u1z)*inv_c2 );
-            psi1[I1[i1]] = g1 - u1z*inv_c;
+            psi1[I1[i1]] = is_beam_coll ? u1z : g1 - u1z*inv_c;
             g2 = std::sqrt( T_R(1.0) +
                             (u2x[I2[i2]]*u2x[I2[i2]]+u2y[I2[i2]]*u2y[I2[i2]]+u2z*u2z)*inv_c2 );
             psi2[I2[i2]] = g2 - u2z*inv_c;

From 5965fbed2711461ec1d969a232595bad0ba1f9e0 Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Thu, 15 Jun 2023 22:23:26 +0200
Subject: [PATCH 05/14] fix doc

---
 src/particles/collisions/ElasticCollisionPerez.H | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/particles/collisions/ElasticCollisionPerez.H b/src/particles/collisions/ElasticCollisionPerez.H
index 36cb1d1a54..3086175175 100644
--- a/src/particles/collisions/ElasticCollisionPerez.H
+++ b/src/particles/collisions/ElasticCollisionPerez.H
@@ -50,6 +50,7 @@
  * @param[in] normalized_units whether normalized units are used
  * @param[in] background_density_SI background plasma density (only needed for normalized units)
  * @param[in] is_same_species whether the collisions happen within the same species
+ * @param[in] is_beam_coll whether species1 is a beam
  * @param[in] engine AMReX engine for the random number generator.
 */
 

From b2931706b2703396bf038f08a0544ab258fd8254 Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Thu, 15 Jun 2023 22:33:33 +0200
Subject: [PATCH 06/14] fix CI

---
 tests/collisions.SI.1Rank.sh | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/collisions.SI.1Rank.sh b/tests/collisions.SI.1Rank.sh
index 6aaa80dedf..6ed72e0111 100755
--- a/tests/collisions.SI.1Rank.sh
+++ b/tests/collisions.SI.1Rank.sh
@@ -29,7 +29,7 @@ TEST_NAME="${FILE_NAME%.*}"
 OMP_NUM_THREADS=1 mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \
         hipace.do_tiling = 0 \
         hipace.file_prefix=$TEST_NAME \
-        plasmas.collisions = collision1 \
+        hipace.collisions = collision1 \
         collision1.species = plasma plasma
 
 # Compare the results with checksum benchmark

From 749042a1acbd89350fe063dae8e93e1eecba411e Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Mon, 19 Jun 2023 14:26:40 +0200
Subject: [PATCH 07/14] fix automatic temperature calculation for beams

---
 src/particles/collisions/ComputeTemperature.H    | 10 ++++++----
 src/particles/collisions/ElasticCollisionPerez.H |  5 +++--
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/particles/collisions/ComputeTemperature.H b/src/particles/collisions/ComputeTemperature.H
index 107ff2e027..413a82d475 100644
--- a/src/particles/collisions/ComputeTemperature.H
+++ b/src/particles/collisions/ComputeTemperature.H
@@ -13,7 +13,7 @@ AMREX_GPU_HOST_DEVICE
 T_R ComputeTemperature (
     T_index const Is, T_index const Ie, T_index const * AMREX_RESTRICT I,
     T_R const * AMREX_RESTRICT ux, T_R const * AMREX_RESTRICT uy, T_R const * AMREX_RESTRICT psi,
-    T_R const m, T_R const clight, T_R const inv_c2 )
+    T_R const m, T_R const clight, T_R const inv_c2, bool is_beam_coll )
 {
     using namespace amrex::literals;
 
@@ -27,9 +27,11 @@ T_R ComputeTemperature (
     for (int i = Is; i < (int) Ie; ++i)
     {
         // particle's Lorentz factor
-        gm = (1.0_rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]])*inv_c2
-              + psi[I[i]]*psi[I[i]]) / (2.0_rt * psi[I[i]] );
-        const amrex::Real uz = clight * (gm - psi[I[i]]);
+        gm = is_beam_coll ? std::sqrt( 1._rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]]
+                                                + psi[I[i]]*psi[I[i]])*inv_c2 )
+                          : (1.0_rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]])*inv_c2
+                                       + psi[I[i]]*psi[I[i]]) / (2.0_rt * psi[I[i]] );
+        const amrex::Real uz = is_beam_coll ? psi[I[i]] : clight * (gm - psi[I[i]]);
         us = ( ux[ I[i] ] * ux[ I[i] ] +
                uy[ I[i] ] * uy[ I[i] ] +
                uz         * uz);
diff --git a/src/particles/collisions/ElasticCollisionPerez.H b/src/particles/collisions/ElasticCollisionPerez.H
index 3086175175..cb33f0df54 100644
--- a/src/particles/collisions/ElasticCollisionPerez.H
+++ b/src/particles/collisions/ElasticCollisionPerez.H
@@ -82,12 +82,13 @@ void ElasticCollisionPerez (
     T_R T1t; T_R T2t;
     if ( T1 <= T_R(0.0) && L <= T_R(0.0) )
     {
-        T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,psi1,m1,clight,inv_c2);
+        T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,psi1,m1,clight,inv_c2,is_beam_coll);
     }
     else { T1t = T1; }
     if ( T2 <= T_R(0.0) && L <= T_R(0.0) )
     {
-        T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,psi2,m2,clight,inv_c2);
+        // second species is always plasma thus is_beam_coll is always false
+        T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,psi2,m2,clight,inv_c2,false);
     }
     else { T2t = T2; }
 

From 3314585705fff5a2e9bc716205e47b9930e78ca4 Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Mon, 19 Jun 2023 15:25:35 +0200
Subject: [PATCH 08/14] add benchmark CI test for collisions with beam

---
 CMakeLists.txt                                |  6 +++
 .../collisions_beam.SI.1Rank.json             | 32 +++++++++++++++
 tests/checksum/reset_all_benchmarks.sh        | 12 ++++++
 tests/collisions_beam.SI.1Rank.sh             | 40 +++++++++++++++++++
 4 files changed, 90 insertions(+)
 create mode 100644 tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json
 create mode 100755 tests/collisions_beam.SI.1Rank.sh

diff --git a/CMakeLists.txt b/CMakeLists.txt
index b7091e3622..a67e2db460 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -345,6 +345,12 @@ if(BUILD_TESTING)
                     WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
             )
 
+            add_test(NAME collisions_beam.SI.1Rank
+                    COMMAND bash ${HiPACE_SOURCE_DIR}/tests/collisions_beam.SI.1Rank.sh
+                            $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}
+                    WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+            )
+
             add_test(NAME ion_motion.SI.1Rank
                     COMMAND bash ${HiPACE_SOURCE_DIR}/tests/ion_motion.SI.1Rank.sh
                             $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}
diff --git a/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json b/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json
new file mode 100644
index 0000000000..088cb26ee5
--- /dev/null
+++ b/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json
@@ -0,0 +1,32 @@
+{
+  "lev=0": {
+    "Bx": 185897.4521776,
+    "By": 185897.45217552,
+    "Bz": 5395.5503546257,
+    "ExmBy": 84161261851902.0,
+    "EypBx": 84161261852290.0,
+    "Ez": 71705050369318.0,
+    "Psi": 1089585554.3813,
+    "Sx": 5219225582596600.0,
+    "Sy": 5219225582743900.0,
+    "chi": 4056826388137200.0,
+    "jx": 1.4609497758209e+16,
+    "jx_beam": 2.1120335304489e-06,
+    "jy": 1.4609497758179e+16,
+    "jy_beam": 4.2519723034246e-06,
+    "jz_beam": 1.0840992316526e+16,
+    "rhomjz": 178176970.1344
+  },
+  "beam": {
+    "charge": 1.1933011570032e-15,
+    "id": 27740076,
+    "mass": 6.7846689808772e-27,
+    "x": 0.03871,
+    "y": 0.03871,
+    "z": 0.2189712,
+    "ux": 0.0,
+    "uy": 0.0,
+    "uz": 14896000.0,
+    "w": 1692775082.5317
+  }
+}
diff --git a/tests/checksum/reset_all_benchmarks.sh b/tests/checksum/reset_all_benchmarks.sh
index ccf8647931..9a4ca5812e 100755
--- a/tests/checksum/reset_all_benchmarks.sh
+++ b/tests/checksum/reset_all_benchmarks.sh
@@ -367,6 +367,18 @@ then
                      --test-name collisions.SI.1Rank
 fi
 
+# collisions_beam.SI.1Rank
+if [[ $all_tests = true ]] || [[ $one_test_name = "collisions_beam.SI.1Rank" ]]
+then
+    cd $build_dir
+    ctest --output-on-failure -R collisions_beam.SI.1Rank \
+        || echo "ctest command failed, maybe just because checksums are different. Keep going"
+    cd $checksum_dir
+    ./checksumAPI.py --reset-benchmark \
+                     --file_name ${build_dir}/bin/collisions_beam.SI.1Rank \
+                     --test-name collisions_beam.SI.1Rank
+fi
+
 # beam_in_vacuum_open_boundary.normalized.1Rank
 if [[ $all_tests = true ]] || [[ $one_test_name = "beam_in_vacuum_open_boundary.normalized.1Rank" ]]
 then
diff --git a/tests/collisions_beam.SI.1Rank.sh b/tests/collisions_beam.SI.1Rank.sh
new file mode 100755
index 0000000000..39c43d9cf6
--- /dev/null
+++ b/tests/collisions_beam.SI.1Rank.sh
@@ -0,0 +1,40 @@
+#! /usr/bin/env bash
+
+# Copyright 2022
+#
+# This file is part of HiPACE++.
+#
+# Authors: MaxThevenet
+# License: BSD-3-Clause-LBNL
+
+
+# This file is part of the HiPACE++ test suite.
+# It runs a Hipace simulation with plasma collisions and compare with benchmark
+# in normalized units
+
+# abort on first encounted error
+set -eu -o pipefail
+
+# Read input parameters
+HIPACE_EXECUTABLE=$1
+HIPACE_SOURCE_DIR=$2
+
+HIPACE_EXAMPLE_DIR=${HIPACE_SOURCE_DIR}/examples/blowout_wake
+HIPACE_TEST_DIR=${HIPACE_SOURCE_DIR}/tests
+
+FILE_NAME=`basename "$0"`
+TEST_NAME="${FILE_NAME%.*}"
+
+# Run the simulation
+OMP_NUM_THREADS=1 mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \
+        hipace.do_tiling = 0 \
+        hipace.file_prefix=$TEST_NAME \
+        hipace.collisions = collision1 \
+        collision1.species = beam plasma
+
+# Compare the results with checksum benchmark
+$HIPACE_TEST_DIR/checksum/checksumAPI.py \
+    --evaluate \
+    --file_name $TEST_NAME \
+    --test-name $TEST_NAME \
+    --skip "{'beam': 'id'}"

From fd649f1f2039eb32eed0d12dc7f1d9a93ee46032 Mon Sep 17 00:00:00 2001
From: Tools <hipace.tools@desy.de>
Date: Mon, 19 Jun 2023 15:59:57 +0200
Subject: [PATCH 09/14] reset benchmark for new CI test

---
 .../collisions_beam.SI.1Rank.json              | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json b/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json
index 088cb26ee5..95d033bfb2 100644
--- a/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json
+++ b/tests/checksum/benchmarks_json/collisions_beam.SI.1Rank.json
@@ -1,19 +1,19 @@
 {
   "lev=0": {
     "Bx": 185897.4521776,
-    "By": 185897.45217552,
-    "Bz": 5395.5503546257,
-    "ExmBy": 84161261851902.0,
-    "EypBx": 84161261852290.0,
-    "Ez": 71705050369318.0,
+    "By": 185897.45217553,
+    "Bz": 5395.5503546258,
+    "ExmBy": 84161261851900.0,
+    "EypBx": 84161261852289.0,
+    "Ez": 71705050369317.0,
     "Psi": 1089585554.3813,
-    "Sx": 5219225582596600.0,
-    "Sy": 5219225582743900.0,
+    "Sx": 5219225582596500.0,
+    "Sy": 5219225582743600.0,
     "chi": 4056826388137200.0,
     "jx": 1.4609497758209e+16,
-    "jx_beam": 2.1120335304489e-06,
+    "jx_beam": 2.1124728260007e-06,
     "jy": 1.4609497758179e+16,
-    "jy_beam": 4.2519723034246e-06,
+    "jy_beam": 4.2513859542859e-06,
     "jz_beam": 1.0840992316526e+16,
     "rhomjz": 178176970.1344
   },

From 497cfe3434f81e462306fe06f72ed904bd0fd37c Mon Sep 17 00:00:00 2001
From: Severin <severin.diederichs@desy.de>
Date: Mon, 19 Jun 2023 16:20:39 +0200
Subject: [PATCH 10/14] cleaning and doc

---
 docs/source/run/parameters.rst       | 126 +++++++++++++++++++++++++++
 src/particles/plasma/MultiPlasma.H   |  10 ---
 src/particles/plasma/MultiPlasma.cpp |   7 --
 3 files changed, 126 insertions(+), 17 deletions(-)

diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst
index 66842a089f..19b77689c4 100644
--- a/docs/source/run/parameters.rst
+++ b/docs/source/run/parameters.rst
@@ -207,6 +207,132 @@ Time step
     Only used when using adaptive time step (see ``hipace.dt`` above).
     Threshold beam momentum, below which the time step is not decreased (to avoid arbitrarily small time steps).
 
+* ``hipace.normalized_units`` (`bool`) optional (default `0`)
+    Using normalized units in the simulation.
+
+* ``hipace.verbose`` (`int`) optional (default `0`)
+    Level of verbosity.
+
+      * ``hipace.verbose = 1``, prints only the time steps, which are computed.
+
+      * ``hipace.verbose = 2`` additionally prints the number of iterations in the
+        predictor-corrector loop, as well as the B-Field error at each slice.
+
+      * ``hipace.verbose = 3`` also prints the number of particles, which violate the quasi-static
+        approximation and were neglected at each slice. It prints the number of ionized particles,
+        if ionization occurred. It also adds additional information if beams
+        are read in from file.
+
+* ``hipace.do_device_synchronize`` (`int`) optional (default `0`)
+    Level of synchronization on GPU.
+
+      * ``hipace.do_device_synchronize = 0``, synchronization happens only when necessary.
+
+      * ``hipace.do_device_synchronize = 1``, synchronizes most functions (all that are profiled
+        via ``HIPACE_PROFILE``)
+
+      * ``hipace.do_device_synchronize = 2`` additionally synchronizes low-level functions (all that
+        are profiled via ``HIPACE_DETAIL_PROFILE``)
+
+* ``hipace.depos_order_xy`` (`int`) optional (default `2`)
+    Transverse particle shape order. Currently, `0,1,2,3` are implemented.
+
+* ``hipace.depos_order_z`` (`int`) optional (default `0`)
+    Longitudinal particle shape order. Currently, only `0` is implemented.
+
+* ``hipace.depos_derivative_type`` (`int`) optional (default `2`)
+    Type of derivative used in explicit deposition. `0`: analytic, `1`: nodal, `2`: centered
+
+* ``hipace.outer_depos_loop`` (`bool`) optional (default `0`)
+    If the loop over depos_order is included in the loop over particles.
+
+* ``hipace.beam_injection_cr`` (`integer`) optional (default `1`)
+    Using a temporary coarsed grid for beam particle injection for a fixed particle-per-cell beam.
+    For very high-resolution simulations, where the number of grid points (`nx*ny*nz`)
+    exceeds the maximum `int (~2e9)`, it enables beam particle injection, which would
+    fail otherwise. As an example, a simulation with `2048 x 2048 x 2048` grid points
+    requires ``hipace.beam_injection_cr = 8``.
+
+* ``hipace.do_beam_jx_jy_deposition`` (`bool`) optional (default `1`)
+    Using the default, the beam deposits all currents ``Jx``, ``Jy``, ``Jz``. Using
+    ``hipace.do_beam_jx_jy_deposition = 0`` disables the transverse current deposition of the beams.
+
+* ``hipace.boxes_in_z`` (`int`) optional (default `1`)
+    Number of boxes along the z-axis. In serial runs, the arrays for 3D IO can easily exceed the
+    memory of a GPU. Using multiple boxes reduces the memory requirements by the same factor.
+    This option is only available in serial runs, in parallel runs, please use more GPU to achieve
+    the same effect.
+
+* ``hipace.openpmd_backend`` (`string`) optional (default `h5`)
+    OpenPMD backend. This can either be ``h5``, ``bp``, or ``json``. The default is chosen by what is
+    available. If both Adios2 and HDF5 are available, ``h5`` is used. Note that ``json`` is extremely
+    slow and is not recommended for production runs.
+
+* ``hipace.file_prefix`` (`string`) optional (default `diags/hdf5/`)
+    Path of the output.
+
+* ``hipace.do_tiling`` (`bool`) optional (default `true`)
+    Whether to use tiling, when running on CPU.
+    Currently, this option only affects plasma operations (gather, push and deposition).
+    The tile size can be set with ``plasmas.sort_bin_size``.
+
+* ``hipace.comms_buffer_on_gpu`` (`bool`) optional (default `0`)
+    Whether the buffers used for MPI communication should be allocated on the GPU (device memory).
+    By default they will be allocated on the CPU (pinned memory).
+    Setting this option to `1` is necessary to take advantage of GPU-Enabled MPI, however for this
+    additional enviroment variables need to be set depending on the system.
+
+* ``hipace.do_beam_jz_minus_rho`` (`bool`) optional (default `0`)
+    Whether the beam contribution to :math:`j_z-c\rho` is calculated and used when solving for Psi (used to caculate the transverse fields Ex-By and Ey+Bx).
+    if 0, this term is assumed to be 0 (a good approximation for an ultra-relativistic beam in the z direction with small transverse momentum).
+
+* ``hipace.deposit_rho`` (`bool`) optional (default `0`)
+    If the charge density ``rho`` of the plasma should be deposited so that it is available as a diagnostic.
+    Otherwise only ``rhomjz`` equal to :math:`\rho-j_z/c` will be available.
+    If ``rho`` is explicitly mentioned in ``diagnostic.field_data``, then the default will become `1`.
+
+* ``hipace.salame_n_iter`` (`int`) optional (default `3`)
+    Number of iterations the SALAME algorithm should do when it is used.
+
+* ``hipace.salame_do_advance`` (`bool`) optional (default `1`)
+    Whether the SALAME algorithm should calculate the SALAME-beam-only Ez field
+    by advancing plasma (if `1`) particles or by approximating it using the chi field (if `0`).
+
+* ``hipace.salame_Ez_target(zeta,zeta_initial,Ez_initial)`` (`string`) optional (default `Ez_initial`)
+    Parser function to specify the target Ez field at the witness beam for SALAME.
+    ``zeta``: position of the Ez field to set.
+    ``zeta_initial``: position where the SALAME algorithm first started.
+    ``Ez_initial``: field value at `zeta_initial`.
+    For `zeta` equal to `zeta_initial`, the function should return `Ez_initial`.
+    The default value of this function corresponds to a flat Ez field at the position of the SALAME beam.
+    Note: `zeta` is always less than or equal to `zeta_initial` and `Ez_initial` is typically below zero for electron beams.
+
+* ``hipace.background_density_SI`` (`float`) optional
+    Background plasma density in SI units. Certain physical modules (collisions, ionization, radiation reactions) depend on the actual background density.
+    Hence, in normalized units, they can only be included, if a background plasma density in SI units is provided using this input parameter.
+
+Binary collisions
+-----------------
+
+WARNING: this module is in development.
+
+HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX,
+for collisions between plasma-plasma and beam-plasma.
+As collisions depend on the physical density, in normalized units `hipace.background_density_SI` must be specified.
+
+* ``plasmas.collisions`` (list of `strings`) optional
+    List of names of binary Coulomb collisions.
+    Each will represent collisions between 2 species.
+
+* ``<collision name>.species`` (two `strings`) optional
+    The name of the two species for which collisions should be included.
+    This can either be plasma-plasma or beam-plasma collisions. For plasma-plasma collisions, the species can be the same to model collisions within a species.
+    The names must be in `plasmas.names` or `beams.names` (for beam-plasma collisions).
+
+* ``<collision name>.CoulombLog`` (`float`) optional (default `-1.`)
+    Coulomb logarithm used for this collision.
+    If not specified, the Coulomb logarithm is determined from the temperature in each cell.
+
 Field solver parameters
 -----------------------
 
diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H
index f22c8ede07..1616f8c2f9 100644
--- a/src/particles/plasma/MultiPlasma.H
+++ b/src/particles/plasma/MultiPlasma.H
@@ -122,16 +122,6 @@ public:
     /** returns number of plasma species */
     int GetNPlasmas() const {return m_nplasmas;}
 
-    // /** Perform binary elastic Coulomb collision, inter- and/or intra-species.
-    //  *
-    //  * The algorithm implemented is that of [Perez et al., Phys. Plasmas 19, 083104 (2012)]
-    //  *
-    //  * \param[in] lev MR level
-    //  * \param[in] bx box on which plasma particles are sorted per-cell and collide together
-    //  * \param[in] geom Corresponding gemetry
-    //  */
-    // void doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom);
-
     /** Reorder particles to speed-up current deposition
      * \param[in] islice zeta slice index
      */
diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp
index 2984a13bb6..a6954be746 100644
--- a/src/particles/plasma/MultiPlasma.cpp
+++ b/src/particles/plasma/MultiPlasma.cpp
@@ -159,13 +159,6 @@ MultiPlasma::TileSort (amrex::Box bx, amrex::Geometry geom)
     }
 }
 
-// void
-// MultiPlasma::doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom)
-// {
-//     HIPACE_PROFILE("MultiPlasma::doCoulombCollision()");
-//
-// }
-
 void
 MultiPlasma::ReorderParticles (const int islice)
 {

From e061f592e19ddaf0e421bd29c7df5b1230d538a7 Mon Sep 17 00:00:00 2001
From: Severin Diederichs <severin.diederichs@desy.de>
Date: Wed, 5 Jul 2023 14:03:57 +0200
Subject: [PATCH 11/14] fix doc

---
 docs/source/run/parameters.rst | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst
index 19b77689c4..0198c8066e 100644
--- a/docs/source/run/parameters.rst
+++ b/docs/source/run/parameters.rst
@@ -936,15 +936,18 @@ Binary collisions
 
 WARNING: this module is in development.
 
-HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX, between plasma species.
+HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX,
+for collisions between plasma-plasma and beam-plasma.
 As collisions depend on the physical density, in normalized units `hipace.background_density_SI` must be specified.
 
 * ``plasmas.collisions`` (list of `strings`) optional
-    List of names of types binary Coulomb collisions.
-    Each will represent collisions between 2 plasma species (potentially the same).
+    List of names of binary Coulomb collisions.
+    Each will represent collisions between 2 species.
 
 * ``<collision name>.species`` (two `strings`) optional
-    The name of the two plasma species for which collisions should be included.
+    The name of the two species for which collisions should be included.
+    This can either be plasma-plasma or beam-plasma collisions. For plasma-plasma collisions, the species can be the same to model collisions within a species.
+    The names must be in `plasmas.names` or `beams.names` (for beam-plasma collisions).
 
 * ``<collision name>.CoulombLog`` (`float`) optional (default `-1.`)
     Coulomb logarithm used for this collision.

From 877491531354510393f6fb77c912ba4a05756696 Mon Sep 17 00:00:00 2001
From: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com>
Date: Wed, 5 Jul 2023 15:22:29 +0200
Subject: [PATCH 12/14] Apply suggestions from code review

---
 docs/source/run/parameters.rst | 112 ---------------------------------
 1 file changed, 112 deletions(-)

diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst
index 0198c8066e..50db07cc2d 100644
--- a/docs/source/run/parameters.rst
+++ b/docs/source/run/parameters.rst
@@ -207,118 +207,6 @@ Time step
     Only used when using adaptive time step (see ``hipace.dt`` above).
     Threshold beam momentum, below which the time step is not decreased (to avoid arbitrarily small time steps).
 
-* ``hipace.normalized_units`` (`bool`) optional (default `0`)
-    Using normalized units in the simulation.
-
-* ``hipace.verbose`` (`int`) optional (default `0`)
-    Level of verbosity.
-
-      * ``hipace.verbose = 1``, prints only the time steps, which are computed.
-
-      * ``hipace.verbose = 2`` additionally prints the number of iterations in the
-        predictor-corrector loop, as well as the B-Field error at each slice.
-
-      * ``hipace.verbose = 3`` also prints the number of particles, which violate the quasi-static
-        approximation and were neglected at each slice. It prints the number of ionized particles,
-        if ionization occurred. It also adds additional information if beams
-        are read in from file.
-
-* ``hipace.do_device_synchronize`` (`int`) optional (default `0`)
-    Level of synchronization on GPU.
-
-      * ``hipace.do_device_synchronize = 0``, synchronization happens only when necessary.
-
-      * ``hipace.do_device_synchronize = 1``, synchronizes most functions (all that are profiled
-        via ``HIPACE_PROFILE``)
-
-      * ``hipace.do_device_synchronize = 2`` additionally synchronizes low-level functions (all that
-        are profiled via ``HIPACE_DETAIL_PROFILE``)
-
-* ``hipace.depos_order_xy`` (`int`) optional (default `2`)
-    Transverse particle shape order. Currently, `0,1,2,3` are implemented.
-
-* ``hipace.depos_order_z`` (`int`) optional (default `0`)
-    Longitudinal particle shape order. Currently, only `0` is implemented.
-
-* ``hipace.depos_derivative_type`` (`int`) optional (default `2`)
-    Type of derivative used in explicit deposition. `0`: analytic, `1`: nodal, `2`: centered
-
-* ``hipace.outer_depos_loop`` (`bool`) optional (default `0`)
-    If the loop over depos_order is included in the loop over particles.
-
-* ``hipace.beam_injection_cr`` (`integer`) optional (default `1`)
-    Using a temporary coarsed grid for beam particle injection for a fixed particle-per-cell beam.
-    For very high-resolution simulations, where the number of grid points (`nx*ny*nz`)
-    exceeds the maximum `int (~2e9)`, it enables beam particle injection, which would
-    fail otherwise. As an example, a simulation with `2048 x 2048 x 2048` grid points
-    requires ``hipace.beam_injection_cr = 8``.
-
-* ``hipace.do_beam_jx_jy_deposition`` (`bool`) optional (default `1`)
-    Using the default, the beam deposits all currents ``Jx``, ``Jy``, ``Jz``. Using
-    ``hipace.do_beam_jx_jy_deposition = 0`` disables the transverse current deposition of the beams.
-
-* ``hipace.boxes_in_z`` (`int`) optional (default `1`)
-    Number of boxes along the z-axis. In serial runs, the arrays for 3D IO can easily exceed the
-    memory of a GPU. Using multiple boxes reduces the memory requirements by the same factor.
-    This option is only available in serial runs, in parallel runs, please use more GPU to achieve
-    the same effect.
-
-* ``hipace.openpmd_backend`` (`string`) optional (default `h5`)
-    OpenPMD backend. This can either be ``h5``, ``bp``, or ``json``. The default is chosen by what is
-    available. If both Adios2 and HDF5 are available, ``h5`` is used. Note that ``json`` is extremely
-    slow and is not recommended for production runs.
-
-* ``hipace.file_prefix`` (`string`) optional (default `diags/hdf5/`)
-    Path of the output.
-
-* ``hipace.do_tiling`` (`bool`) optional (default `true`)
-    Whether to use tiling, when running on CPU.
-    Currently, this option only affects plasma operations (gather, push and deposition).
-    The tile size can be set with ``plasmas.sort_bin_size``.
-
-* ``hipace.comms_buffer_on_gpu`` (`bool`) optional (default `0`)
-    Whether the buffers used for MPI communication should be allocated on the GPU (device memory).
-    By default they will be allocated on the CPU (pinned memory).
-    Setting this option to `1` is necessary to take advantage of GPU-Enabled MPI, however for this
-    additional enviroment variables need to be set depending on the system.
-
-* ``hipace.do_beam_jz_minus_rho`` (`bool`) optional (default `0`)
-    Whether the beam contribution to :math:`j_z-c\rho` is calculated and used when solving for Psi (used to caculate the transverse fields Ex-By and Ey+Bx).
-    if 0, this term is assumed to be 0 (a good approximation for an ultra-relativistic beam in the z direction with small transverse momentum).
-
-* ``hipace.deposit_rho`` (`bool`) optional (default `0`)
-    If the charge density ``rho`` of the plasma should be deposited so that it is available as a diagnostic.
-    Otherwise only ``rhomjz`` equal to :math:`\rho-j_z/c` will be available.
-    If ``rho`` is explicitly mentioned in ``diagnostic.field_data``, then the default will become `1`.
-
-* ``hipace.salame_n_iter`` (`int`) optional (default `3`)
-    Number of iterations the SALAME algorithm should do when it is used.
-
-* ``hipace.salame_do_advance`` (`bool`) optional (default `1`)
-    Whether the SALAME algorithm should calculate the SALAME-beam-only Ez field
-    by advancing plasma (if `1`) particles or by approximating it using the chi field (if `0`).
-
-* ``hipace.salame_Ez_target(zeta,zeta_initial,Ez_initial)`` (`string`) optional (default `Ez_initial`)
-    Parser function to specify the target Ez field at the witness beam for SALAME.
-    ``zeta``: position of the Ez field to set.
-    ``zeta_initial``: position where the SALAME algorithm first started.
-    ``Ez_initial``: field value at `zeta_initial`.
-    For `zeta` equal to `zeta_initial`, the function should return `Ez_initial`.
-    The default value of this function corresponds to a flat Ez field at the position of the SALAME beam.
-    Note: `zeta` is always less than or equal to `zeta_initial` and `Ez_initial` is typically below zero for electron beams.
-
-* ``hipace.background_density_SI`` (`float`) optional
-    Background plasma density in SI units. Certain physical modules (collisions, ionization, radiation reactions) depend on the actual background density.
-    Hence, in normalized units, they can only be included, if a background plasma density in SI units is provided using this input parameter.
-
-Binary collisions
------------------
-
-WARNING: this module is in development.
-
-HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX,
-for collisions between plasma-plasma and beam-plasma.
-As collisions depend on the physical density, in normalized units `hipace.background_density_SI` must be specified.
 
 * ``plasmas.collisions`` (list of `strings`) optional
     List of names of binary Coulomb collisions.

From 42721d5cc9047d73c19d73221af08058fef36d9f Mon Sep 17 00:00:00 2001
From: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com>
Date: Wed, 5 Jul 2023 15:22:43 +0200
Subject: [PATCH 13/14] Update docs/source/run/parameters.rst

---
 docs/source/run/parameters.rst | 14 --------------
 1 file changed, 14 deletions(-)

diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst
index 50db07cc2d..bb752f9019 100644
--- a/docs/source/run/parameters.rst
+++ b/docs/source/run/parameters.rst
@@ -207,20 +207,6 @@ Time step
     Only used when using adaptive time step (see ``hipace.dt`` above).
     Threshold beam momentum, below which the time step is not decreased (to avoid arbitrarily small time steps).
 
-
-* ``plasmas.collisions`` (list of `strings`) optional
-    List of names of binary Coulomb collisions.
-    Each will represent collisions between 2 species.
-
-* ``<collision name>.species`` (two `strings`) optional
-    The name of the two species for which collisions should be included.
-    This can either be plasma-plasma or beam-plasma collisions. For plasma-plasma collisions, the species can be the same to model collisions within a species.
-    The names must be in `plasmas.names` or `beams.names` (for beam-plasma collisions).
-
-* ``<collision name>.CoulombLog`` (`float`) optional (default `-1.`)
-    Coulomb logarithm used for this collision.
-    If not specified, the Coulomb logarithm is determined from the temperature in each cell.
-
 Field solver parameters
 -----------------------
 

From 775840dccee64dece167b5d928bc7acf3442c455 Mon Sep 17 00:00:00 2001
From: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com>
Date: Wed, 5 Jul 2023 15:23:39 +0200
Subject: [PATCH 14/14] Update docs/source/run/parameters.rst

---
 docs/source/run/parameters.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst
index bb752f9019..dbb863e586 100644
--- a/docs/source/run/parameters.rst
+++ b/docs/source/run/parameters.rst
@@ -814,7 +814,7 @@ HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (
 for collisions between plasma-plasma and beam-plasma.
 As collisions depend on the physical density, in normalized units `hipace.background_density_SI` must be specified.
 
-* ``plasmas.collisions`` (list of `strings`) optional
+* ``hipace.collisions`` (list of `strings`) optional
     List of names of binary Coulomb collisions.
     Each will represent collisions between 2 species.