anchor.diff:         z_pos = dot_product(z, solvent_pos(:,o))
          solvent_pos(:,o) = (/x_pos,y_pos,z_pos/)
       end do
  +    !$omp parallel do private(r, theta, check)
       do o = 1, solvent% Nmax
          r = sqrt(solvent_pos(1,o)**2 + solvent_pos(2,o)**2)
          theta = atan(solvent_pos(2,o)/solvent_pos(1,o))
  
positive.code:     integer :: m
      double precision :: x(3)
  
      do m = 1, solvent% Nmax
         if (solvent% flag(m) == 1) then
            x = rel_pos(colloids% pos(:,1), solvent% pos(:,m), solvent_cells% edges)
  
negative.code: 
      far = (L(1)*0.45)**2
  
      !$omp parallel do private(x, dist_to_C_sq, dist_to_N_sq)
      do n = 1,solvent% Nmax
         if (solvent% species(n) == 2) then
            x = rel_pos(colloids% pos(:,1), solvent% pos(:,n), solvent_cells% edges)
  