anchor.diff:         end do
       end do placement_loop
   
  +    ! Retrieve location in particle array
       i = solvent%id_to_idx(bound_molecule_id(enzyme_idx))
   
       solvent%species(i) = to_species
  
positive.code:         do k = 1, colloids%Nmax
             call compute_radial_histogram(radial_hist(k), colloids%pos(:,k), solvent_cells%edges, solvent)
          end do
  
       end if
  
       call varia%tac()
  
       if (sampling) then
          n_solvent = 0
          do k = 1, solvent%Nmax
             j = solvent%species(k)
             if (j <= 0) cycle
  
negative.code: 
      enzyme_bound(enzyme_idx) = .true.
  
      ! sample from Poisson process
      next_reaction_time(enzyme_idx) = current_time - log(threefry_double(state(1)))/total_rate
  
    end subroutine bind_molecule
  
    !> Bind solvent molecule idx to enzyme
    subroutine unbind_molecule(enzyme_idx, to_species)
      integer, intent(in) :: enzyme_idx
      integer, intent(in) :: to_species
  