From 456444957ad5099f7cbb917bdcf5eb51e7ffa7af Mon Sep 17 00:00:00 2001
From: Frederic Perez <frederic.perez@polytechnique.edu>
Date: Wed, 24 Apr 2024 15:08:55 +0200
Subject: [PATCH] fix collisions: wrong check when 0 particles

---
 src/Collisions/BinaryProcesses.cpp | 26 ++++++++++++++------------
 1 file changed, 14 insertions(+), 12 deletions(-)

diff --git a/src/Collisions/BinaryProcesses.cpp b/src/Collisions/BinaryProcesses.cpp
index 0b93d0284..107d12c8f 100644
--- a/src/Collisions/BinaryProcesses.cpp
+++ b/src/Collisions/BinaryProcesses.cpp
@@ -238,22 +238,26 @@ void BinaryProcesses::apply( Params &params, Patch *patch, int itime, vector<Dia
         if( shuffle1 ) {
             swap( npartmin, npartmax );
         }
-        // skip to next bin if no particles
-        if( npartmax == 0 ) {
+        
+        // skip to next bin if not enough pairs
+        if( npartmin == 0 || ( intra_ && npartmin < 2 ) ) {
             continue;
         }
         
-        size_t npairs, N2max;
+        size_t npairs, npairs_not_repeated;
+        double weight_correction_1, weight_correction_2;
         if( intra_ ) { // In the case of pairing within one species
-            if( npartmax < 2 ) {
-                continue;
-            }
             npairs = ( npartmax + 1 ) / 2; // half as many pairs as macro-particles
-            N2max = npartmax - npairs; // number of not-repeated particles (in group with less particles)
+            npairs_not_repeated = npartmax - npairs;
+            weight_correction_1 = 1.;
+            weight_correction_2 = 0.5;
         } else { // In the case of pairing between two species
-            npairs = npartmax; // as many pairs as macro-particles (in group with more particles)
-            N2max = npartmin; // number of not-repeated particles (in group with less particles)
+            npairs = npartmax; // as many pairs as macro-particles in group with more particles
+            npairs_not_repeated = npartmin;
+            weight_correction_1 = 1. / (double)( npairs / npairs_not_repeated );
+            weight_correction_2 = 1. / (double)( npairs / npairs_not_repeated + 1 );
         }
+        
         RandomShuffle shuffler( *patch->rand_, npartmax );
         
         // Calculate the densities
@@ -286,8 +290,6 @@ void BinaryProcesses::apply( Params &params, Patch *patch, int itime, vector<Dia
         // Pre-calculate some numbers before the big loop
         unsigned int ncorr = intra_ ? 2*npairs-1 : npairs;
         double dt_corr = every_ * params.timestep * ((double)ncorr) * inv_cell_volume;
-        double weight_correction_1 = 1. / (double)( (npairs-1) / N2max );
-        double weight_correction_2 = 1. / (double)( (npairs-1) / N2max + 1 );
         n1  *= inv_cell_volume;
         n2  *= inv_cell_volume;
         D.n123 = pow( n1, 2./3. );
@@ -343,7 +345,7 @@ void BinaryProcesses::apply( Params &params, Patch *patch, int itime, vector<Dia
             
             // Calculate the timestep correction
             D.dt_correction = D.maxW * dt_corr;
-            if( i % N2max <= (npairs-1) % N2max ) {
+            if( i % npairs_not_repeated < npairs % npairs_not_repeated ) {
                 D.dt_correction *= weight_correction_2 ;
             } else {
                 D.dt_correction *= weight_correction_1;