anchor.diff:      integer :: m
       double precision :: x(3)
   
  -    !$omp parallel do private(x, dist_to_C_sq, dist_to_N_sq)
  +    catalytic_change = 0
  +    !$omp parallel do private(x, dist_to_C_sq, dist_to_N_sq) reduction(+:catalytic_change)
       do m = 1, solvent% Nmax
          if (solvent% flag(m) == 1) then
             x = rel_pos(colloids% pos(:,1), solvent% pos(:,m), solvent_cells% edges)
  
positive.code:   allocate(rho_xy(N_species, L(2), L(1)))
    call h5gcreate_f(hfile%id, 'fields', fields_group, error)
    call rho_xy_el%create_time(fields_group, 'rho_xy', rho_xy, ior(H5MD_LINEAR,H5MD_STORE_TIME), &
         step=N_loop*N_MD_steps, time=N_loop*N_MD_steps*dt)
    call h5gclose_f(fields_group, error)
  
    colloids% pos(3,:) = solvent_cells% edges(3)/2.d0
    if (order) then
  
negative.code:      double precision, intent(in) :: edges(3)
       integer, intent(in) :: bufferlength
    
       integer :: k, s
  
       bulk_change = 0
       do k = 1, particles% Nmax
          s = particles% species(k)
          if (s <= 0) continue
          if (particles% pos(1,k) < bufferlength) then
             if (particles% pos(2,k) < edges(2)/2.d0) then
                bulk_change(s) = bulk_change(s) - 1
                particles% species(k) = 1
                s = 1
                bulk_change(s) = bulk_change(s) + 1
             else
                bulk_change(s) = bulk_change(s) - 1
                particles% species(k) = 3
                s = 3
                bulk_change(s) = bulk_change(s) + 1
             end if  
          end if 
       end do
  