anchor.diff:             if (threefry_double(state(1)) < rate_release_s / total_rate) then
                 ! release substrate
                 call unbind_molecule(enzyme_i, 1)
  +              n_plus(1) = n_plus(1) + 1
              else
                 ! release product
                 call unbind_molecule(enzyme_i, 2)
  +              n_plus(2) = n_plus(2) + 1
              end if
              next_reaction_time(enzyme_i) = huge(next_reaction_time(enzyme_i))
           end if
  
positive.code: 
       call neigh% make_stencil(solvent_cells, enzyme_capture_radius)
       call neigh% update_list(colloids, solvent, enzyme_capture_radius, solvent_cells)
       if (i>equilibration_loops) then
          call reset_enzyme_region_bit
          call select_substrate
       end if
  
negative.code:              ! pick in substrates
               idx = floor(xi * n_s / total_p) + 1
               m = list_s(idx)
               n_minus(1) = n_minus(1) + 1
               ! use list_s(idx)
            else
               ! pick in products
               idx = floor(xi * n_p / total_p) + 1
               m = list_p(idx)
               n_minus(2) = n_minus(2) + 1
            end if
            call bind_molecule(enzyme_i, m)
         end if
  