anchor.diff:         n = cells% cell_count(cell_idx)
          count = count + 1
   
  +       n_effective = 0
          local_v = 0
          local_kin = 0
          do i = start, start + n - 1
  -          local_v = local_v + particles% vel(:, i)
  +          if (particles%species(i) > 0) then
  +             local_v = local_v + particles% vel(:, i)
  +             n_effective = n_effective + 1
  +          end if
          end do
          local_v = local_v / n
          do i = start, start + n - 1
  -          local_kin = local_kin + sum((particles% vel(:, i)-local_v)**2)
  -       end do
  -       local_kin = local_kin / dble(3*(n-1))
  +          if (particles%species(i) > 0) then
  +             local_kin = local_kin + sum((particles% vel(:, i)-local_v)**2)
  +          end if
  +       end do
  +       local_kin = local_kin / dble(3*(n_effective-1))
          te = te + local_kin
   
          if (do_tz) then
  
positive.code:     call particles%time_step%tic()
      !$omp parallel private(thread_id)
      thread_id = omp_get_thread_num() + 1
      !$omp do private(start, n, local_v, i, vec, omega)
      do cell_idx = 1, cells% N
         if (cells% cell_count(cell_idx) <= 1) cycle
  
         start = cells% cell_start(cell_idx)
         n = cells% cell_count(cell_idx)
  
         local_v = 0
         do i = start, start + n - 1
            local_v = local_v + particles% vel(:, i)
         end do
         local_v = local_v / n
  
         vec = rand_sphere(state(thread_id))
         omega = &
  
negative.code: 
      double precision :: te
  
      integer :: i, start, n, n_effective
      integer :: cell_idx, count
      integer :: cell(3)
      double precision :: local_v(3), local_kin
  