anchor.diff:             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)
  
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:   type(h5md_element_t) :: n_solvent_el
    type(h5md_element_t) :: link_angle_el
  
    integer, dimension(N_species) :: n_plus, n_minus
    type(h5md_element_t) :: n_plus_el, n_minus_el
  
    !! Histogram variable
    type(histogram_t), allocatable :: radial_hist(:)
    double precision, allocatable :: dummy_hist(:,:,:)
  