anchor.diff:      end do
       kin = kin / 2
   
  +    write(31,*) 'com_v', com_v, 'kin', kin
  +
  +    if ((kin + energy) <= 0) then
  +       stop 'negative energy in add_energy_to_cell'
  +    end if
       factor = sqrt((kin + energy)/kin)
       do i = start, start + n - 1
          solvent%vel(:,i) = com_v + factor*(solvent%vel(:,i) - com_v)
  
positive.code:   double precision :: enzyme_capture_radius
    double precision :: proba_p, proba_s
    double precision :: rate_release_p, rate_release_s, total_rate, next_reaction_time
    double precision :: substrate_fraction
  
    integer, dimension(N_species) :: n_solvent
    type(h5md_element_t) :: n_solvent_el
  
negative.code:     en = en + angle_k * (acos(costheta)-link_angle)**2 / 2
  
    end function compute_bead_force
  
    function compute_bead_energy() result(en)
      double precision :: en
  
      double precision :: r12(3), d12, r32(3), d32, costheta
  
      r12 = rel_pos(colloids%pos(:,1), colloids%pos(:,2), solvent_cells%edges)
      d12 = norm2(r12)
      r32 = rel_pos(colloids%pos(:,3), colloids%pos(:,2), solvent_cells%edges)
      d32 = norm2(r32)
      costheta = dot_product(r12, r32)/(d12*d32)
      en = angle_k * (acos(costheta)-link_angle)**2 / 2
    end function compute_bead_energy
  
    !> derivative of the angular harmonic arccos term
    function fprime(x) result(r)
  