anchor.diff:  
          n_s = 0
          n_p = 0
  +       s_found = 0
  +       p_found = 0
   
          x_enzyme = modulo(colloids%pos(:,enz_2), solvent_cells%edges)
   
          select_loop: do i = 1, neigh%n(enz_2)
             m = neigh%list(i, enz_2)
   
  -          ! skip substrates undergoing MD
  -          if (btest(solvent%flags(m), MD_BIT)) cycle select_loop
  -
             s_sp = solvent%species(m)
   
             ! skip neutral fluid particles
             if ((s_sp == 0) .or. (s_sp == 3)) cycle select_loop
   
             dist = norm2(rel_pos(x_enzyme, solvent%pos(:,m), solvent_cells%edges))
  +
  +          if ( (dist <= enzyme_capture_radius) .and. (dist > solvent_colloid_lj%cut(s_sp, colloids%species(enz_2)) ) ) then
  +             if (s_sp==1) s_found = s_found + 1
  +             if (s_sp==2) p_found = p_found + 1
  +          end if
   
             if (.not. btest(solvent%flags(m), ENZYME_REGION_BIT)) then
                if ( (dist <= enzyme_capture_radius) .and. (dist > solvent_colloid_lj%cut(s_sp, colloids%species(enz_2)) ) ) then
  
positive.code:                    end if
                     list_p(n_p) = m
                  end if
               end if
            else
               if (dist > enzyme_capture_radius) then
                  solvent%flags(m) = ibclr(solvent%flags(m), ENZYME_REGION_BIT)
               end if
            end if
         end do select_loop
  
negative.code:   sigma_cut = sigma*2**(1.d0/6.d0)
    max_cut = maxval(sigma_cut)
  
    if (max_cut > enzyme_capture_radius) then
       write(*,*) 'enzyme_capture_radius must be greater than max_cut'
       stop
    end if
  
    call solvent_colloid_lj% init(epsilon, sigma, sigma_cut)
  
    epsilon(:,:) = PTread_d(config, 'epsilon_colloid', loc=params_group)
  