anchor.diff:  
        ! Check if unbinding should occur
        do enzyme_i = 1, N_enzymes
  -        if (i*N_MD_steps*dt > next_reaction_time(enzyme_i)) then
  +        if (current_time > next_reaction_time(enzyme_i)) then
              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
  +              call kinetics_data(enzyme_i)%release_substrate%append(current_time)
              else
                 ! release product
                 call unbind_molecule(enzyme_i, 2)
                 n_plus(2) = n_plus(2) + 1
  +              call kinetics_data(enzyme_i)%release_product%append(current_time)
              end if
              next_reaction_time(enzyme_i) = huge(next_reaction_time(enzyme_i))
           end if
  
positive.code:     colloids%mass(enz_2) = colloids%mass(enz_2) + 1
  
      bound_molecule_id(enzyme_idx) = solvent%id(idx)
      solvent%species(idx) = 0
  
      ! compute conformational energy difference when changing the angle
  
negative.code:   type(histogram_t), allocatable :: radial_hist(:)
    double precision, allocatable :: dummy_hist(:,:,:)
  
    !! Kinetics data
    type(enzyme_kinetics_t), allocatable :: kinetics_data(:)
  
    integer :: equilibration_loops
    integer :: colloid_sampling
    logical :: sampling
  