anchor.diff:      part_count(binx, biny, binz) = part_count(binx, biny, binz) + 1
       binpart(part_count(binx,biny,binz),binx,biny,binz) = i
     end do
  -
  +  !$OMP END SINGLE
  +
  +  !$OMP WORKSHARE
     nbs = 0
     nc = 0
  -
  -  !$OMP PARALLEL DO default(none) shared(nbs,nc,nbinx,nbiny,nbinz,px,py,pz,h2,part_count,binpart) &
  -  !$OMP private(q,i,binx,biny,binz,stx,sty,stz,cbinx,cbiny,cbinz,nb,distx,disty,distz,tdist2) &
  -  !$OMP NUM_THREADS(procs)
  -  do i = 1, n
  +  !$OMP END WORKSHARE
  +
  +  !$OMP DO private(binx,biny,binz,cbinx,cbiny,cbinz,nb,distx,disty,distz,tdist2,q)
  +  do i = 1,n
       binx = NINT((px(i))/(BOX_SIZE)*(real(nbinx,WP)-1.0_WP)) + 1
       biny = NINT((py(i))/(BOX_SIZE)*(real(nbiny,WP)-1.0_WP)) + 1
       binz = NINT((pz(i))/(BOX_SIZE)*(real(nbinz,WP)-1.0_WP)) + 1
  
positive.code:   use state
    use neighbors
    use posvel
    implicit none
    integer :: i, j, q, p, nb, l
    integer :: binx, biny, binz
  
negative.code:   twrite = frame_freq
    call cpu_time(init_stop)
  
  
    !$OMP PARALLEL NUM_THREADS(procs) default(shared)
    iter = 0
    time = 0.0_WP
    do while (time .lt. tfin)
      !$OMP MASTER
      iter = iter + 1
      !$OMP END MASTER
  
      ! Simulate the fluid flow
      call step
      !$OMP MASTER
      time = time + dt
      !$OMP END MASTER
  
      !$OMP MASTER
        if(time .ge. twrite) then
          call cpu_time(write_start)
          twrite = twrite + frame_freq
  