anchor.diff:  
        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
  +     if (sampling) then
  +        n_plus = 0
  +        n_minus = 0
           call reset_enzyme_region_bit
           call select_substrate
        end if
  
positive.code:            if (threefry_double(state(1)) < rate_release_s / total_rate) then
                ! release substrate
                call unbind_molecule(enzyme_i, 1)
             else
                ! release product
                call unbind_molecule(enzyme_i, 2)
             end if
             next_reaction_time(enzyme_i) = huge(next_reaction_time(enzyme_i))
          end if
  
negative.code:            n_solvent(j) = n_solvent(j) + 1
          end do
          call n_solvent_el%append(n_solvent)
          call n_plus_el%append(n_plus)
          call n_minus_el%append(n_minus)
          call link_angle_el%append(link_angle)
          if (modulo(i, colloid_sampling)==0) then
             call dimer_io%position%append(colloids%pos)
  