anchor.diff:      integer :: list_s(list_size), list_p(list_size)
       logical :: reaction_happens
   
  -    ! Do nothing if enzyme is already bound
  -    if (enzyme_bound) return
  -
  -    n_s = 0
  -    n_p = 0
  -
  -    x_enzyme = modulo(colloids%pos(:,2), solvent_cells%edges)
  -
  -    !$omp parallel private(thread_id)
  -    thread_id = omp_get_thread_num() + 1
  -    !$omp do private(i, m, s_sp, dist)
  -    do i = 1, solvent_cells%N
  -       select_loop: do m = solvent_cells%cell_start(i), solvent_cells%cell_start(i) + solvent_cells%cell_count(i) - 1
  +    integer :: enzyme_i, enz_2
  +
  +
  +    select_substrate_loop: do enzyme_i = 1, N_enzymes
  +
  +       if (enzyme_bound(enzyme_i)) cycle select_substrate_loop
  +
  +       enz_2 = (enzyme_i-1)*3 + 2
  +
  +       n_s = 0
  +       n_p = 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 == 3) cycle select_loop
  +
             dist = norm2(rel_pos(x_enzyme, solvent%pos(:,m), solvent_cells%edges))
  +
             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(2)) ) ) then
  +             if ( (dist <= enzyme_capture_radius) .and. (dist > solvent_colloid_lj%cut(s_sp, colloids%species(enz_2)) ) ) then
                   ! select for current round
                   solvent%flags(m) = ibset(solvent%flags(m), ENZYME_REGION_BIT)
                   if (s_sp==1) then
  
positive.code: 
    call h5fflush_f(hfile%id, H5F_SCOPE_GLOBAL_F, error)
    write(*,*) 'Running for', equilibration_loops, '+', N_loop, 'loops'
    write(*,*) 'mass', mass 
    call main%tic()
    sampling = .false.
    enzyme_bound = .false.
    next_reaction_time = huge(next_reaction_time)
  
    do i = 0, N_loop+equilibration_loops
       if (i==equilibration_loops) sampling = .true.
  
negative.code:         call md_pos(solvent, dt)
          do k=1, colloids% Nmax
             colloids% pos(:,k) = colloids% pos(:,k) + dt * colloids% vel(:,k) + &
                  dt**2 * colloids% force(:,k) / (2 * colloids%mass(k))
          end do
          
          so_max = solvent% maximum_displacement()
  