Minor changes

mbtrack2/pic/grid2d.f90

 subroutine deposit_cic(nx, ny, x_min, y_min, dx, dy, &
     x, y, q, n_particles, rho)
+    use omp_lib
     implicit none
 
     integer, intent(in) :: nx, ny       ! Number of grid points in x and y directions
@@ -13,6 +14,7 @@ subroutine deposit_cic(nx, ny, x_min, y_min, dx, dy, &
 
     integer :: i, i_grid, j_grid
     real(8) :: wx, wy
+    !$omp parallel do private(i, i_grid, j_grid, wx, wy)
     do i = 1, n_particles
 
         i_grid = int((x(i) - x_min) / dx)
@@ -22,18 +24,24 @@ subroutine deposit_cic(nx, ny, x_min, y_min, dx, dy, &
         wy = ((y(i) - y_min) / dy) - dble(j_grid)
 
         if (i_grid >= 0 .and. i_grid < nx-1 .and. j_grid >= 0 .and. j_grid < ny-1) then
+            !$omp atomic
             rho(i_grid+1, j_grid+1) = rho(i_grid+1, j_grid+1) + q * (1.0d0 - wx) * (1.0d0 - wy)                
+            !$omp atomic
             rho(i_grid+2, j_grid+1) = rho(i_grid+2, j_grid+1) + q * wx * (1.0d0 - wy)
+            !$omp atomic
             rho(i_grid+1, j_grid+2) = rho(i_grid+1, j_grid+2) + q * (1.0d0 - wx) * wy
+            !$omp atomic
             rho(i_grid+2, j_grid+2) = rho(i_grid+2, j_grid+2) + q * wx * wy
         end if
     end do
+    !$omp end parallel do 
 
 end subroutine deposit_cic
 
 subroutine interpolate_electric_field(nx, ny, n_particles, x_min, x_max, y_min, y_max, &
     x_particles, y_particles, Ex_grid, Ey_grid, &
     Ex_interp, Ey_interp)
+    use omp_lib
     implicit none
 
     ! Input parameters
@@ -54,6 +62,7 @@ subroutine interpolate_electric_field(nx, ny, n_particles, x_min, x_max, y_min,
     dy = (y_max - y_min) / dble(ny - 1)
 
     ! Loop over each particle to interpolate electric fields
+    !$omp parallel do private(i, j, px, x_rel, y_rel, wx, wy, ex, ey)
     do px = 1, n_particles
         ! Find the indices of the lower-left grid point
         i = int((x_particles(px) - x_min) / dx)
@@ -86,5 +95,6 @@ subroutine interpolate_electric_field(nx, ny, n_particles, x_min, x_max, y_min,
             Ey_interp(px) = 0.0d0
         end if
     end do
+    !$omp end parallel do
 
 end subroutine interpolate_electric_field

mbtrack2/pic/pic2d.py

@@ -23,7 +23,9 @@ class PICSolver2D:
         '''Implemented in fortran pic2d.f90'''
         grid2d.deposit_cic(self.x_min, self.y_min, self.dx, self.dy, x, y, q,
                            self.rho)
-        self.rho /= (epsilon_0 * self.dx * self.dy * self.rho.sum())
+        total_charge = self.rho.sum()
+        if total_charge != 0:
+            self.rho /= (epsilon_0 * self.dx * self.dy * total_charge)
 
 
 # Original Python implementation (only to quickly check that fortran version is much faster before the merge request)