diff --git a/src/modules/common/moduleRefParam.f90 b/src/modules/common/moduleRefParam.f90
index e8819ce..051190b 100644
--- a/src/modules/common/moduleRefParam.f90
+++ b/src/modules/common/moduleRefParam.f90
@@ -3,7 +3,7 @@ MODULE moduleRefParam
   !Parameters that define the problem (inputs)
   REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigmaVrel_ref
   !Reference parameters for non-dimensional problem
-  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref, e_ref
+  REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref
 
 END MODULE moduleRefParam
 
diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index 7af67db..c8beb0f 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -144,7 +144,6 @@ MODULE moduleInput
       EF_ref    = qe*n_ref*L_ref/eps_0           !reference electric field
       Volt_ref  = EF_ref*L_ref                   !reference voltage
       B_ref     = m_ref / (ti_ref * qe)          !reference magnetic field
-      e_ref     = qe/(m_ref*v_ref**2)            !reference charge
 
       !If a reference cross section is given, it is used
       CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found)
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index c80eec3..a62cfb2 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -975,10 +975,11 @@ MODULE moduleMesh
       REAL(8):: delta_par, delta_par_square, delta_per, delta_per_square
       REAL(8):: W(3), dW(2), normW !Relative velocity between particle and species and its increment
       REAL(8):: l2, l, lW, AW
+      REAL(8):: deltaV(1:3), totalDeltaV_ij, normDeltaV
       REAL(8):: rnd
 
 
-      !$OMP SINGLE
+      !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
       DO e = 1, self%numCells
         cell => self%cells(e)%obj
         cellNodes = cell%getNodes(cell%nNodes)
@@ -1003,6 +1004,7 @@ MODULE moduleMesh
           END DO
 
           !Loop over particles of species_i
+          totalDeltaV_ij = 0.D0
           partTemp => cell%listPart_in(i)%head
           DO WHILE(ASSOCIATED(partTemp))
             density     = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
@@ -1036,16 +1038,18 @@ MODULE moduleMesh
             e2(3) =  0.D0
             e2 = normalize(e2)
             !Third one is perpendicular to the other two
-            e3 = crossProduct(e2, e1)
+            e3 = crossProduct(e1, e2)
             e3 = normalize(e3)
 
             !Random number for direction
             rnd = PI2*random()
 
-            !Change particle velocity
-            partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + &
-                                                                  SIN(rnd)*e3)
+            deltaV =  dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3)
 
+            totalDeltaV_ij = totalDeltaV_ij + NORM2(deltaV)
+            !Change particle velocity
+            partTemp%part%v = partTemp%part%v + deltaV
+                                               
             partTemp => partTemp%next
 
           END DO
@@ -1062,6 +1066,8 @@ MODULE moduleMesh
   
             END DO
 
+            !Divide total momentum exchanged among all the particles of species j
+            normDeltaV = totalDeltaV_ij / REAL(cell%listPart_in(j)%amount) * (self%sp_i%weight*self%sp_i%m)/(self%sp_j%weight*self%sp_j%m)
             !Loop over particles of species_j
             partTemp => cell%listPart_in(j)%head
             DO WHILE(ASSOCIATED(partTemp))
@@ -1086,6 +1092,10 @@ MODULE moduleMesh
               dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin)
               dW(2) =                ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
 
+              !Normalize with average exchange per particle
+              !TODO: This is a dirty trick, it should not be neccessary but I don't lnow why it is.
+              dW = normDeltaV*dW/NORM2(dW)
+
               !System of reference for the velocity change
               !First one is parallel to the relative velocity
               e1 = normalize(W)
@@ -1095,15 +1105,16 @@ MODULE moduleMesh
               e2(3) =  0.D0
               e2 = normalize(e2)
               !Third one is perpendicular to the other two
-              e3 = crossProduct(e2, e1)
+              e3 = crossProduct(e1, e2)
               e3 = normalize(e3)
   
               !Random number for direction
               rnd = PI2*random()
   
+              deltaV =  dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3)
+
               !Change particle velocity
-              partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + &
-                                                                    SIN(rnd)*e3)
+              partTemp%part%v = partTemp%part%v + deltaV
   
               partTemp => partTemp%next
   
@@ -1114,7 +1125,7 @@ MODULE moduleMesh
         END DO
 
       END DO
-      !$OMP END SINGLE
+      !$OMP END DO
 
     END SUBROUTINE doCoulomb
 
diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90
index 85f4f49..ab931e0 100644
--- a/src/modules/moduleCoulomb.f90
+++ b/src/modules/moduleCoulomb.f90
@@ -54,6 +54,7 @@ MODULE moduleCoulomb
       CLASS(interactionsCoulomb), INTENT(out):: self
       INTEGER, INTENT(in):: i, j
       REAL(8):: Z_i, Z_j
+      REAL(8):: scaleFactor
 
       self%sp_i => species(i)%obj
       self%sp_j => species(j)%obj
@@ -81,8 +82,10 @@ MODULE moduleCoulomb
 
       self%lnCoulomb = 10.0 !Make this function dependent
 
-      self%A_i = 8.D0*PI*Z_i**2*Z_j**2*self%lnCoulomb / self%sp_i%m**2
-      self%A_j = 8.D0*PI*Z_j**2*Z_i**2*self%lnCoulomb / self%sp_j%m**2
+      scaleFactor = (n_ref * qe**4) / (eps_0**2 * m_ref**2 * v_ref**3) * ti_ref
+
+      self%A_i = 2.D0*Z_i**2*Z_j**2*self%lnCoulomb / self%sp_i%m**2 * scaleFactor
+      self%A_j = 2.D0*Z_j**2*Z_i**2*self%lnCoulomb / self%sp_j%m**2 * scaleFactor
 
       self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent
       self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent