anchor.diff:      solvent%species(idx) = 0
   
       ! compute conformational energy difference when changing the angle
  -    en1 = compute_bead_energy(enzyme_idx)
  +    en1 = compute_bead_energy(enzyme_idx, factor=-1)
       link_angle(enzyme_idx) = link_angles(2)
  -    en2 = compute_bead_energy(enzyme_idx)
  +    en2 = compute_bead_energy(enzyme_idx, factor=1)
   
       call add_energy_to_cell(cell_idx, excess_kinetic_energy + en1 - en2)
   
  
positive.code: 
    end function compute_bead_force
  
    function compute_bead_energy(enzyme_idx) result(en)
      integer, intent(in) :: enzyme_idx
      double precision :: en
  
      double precision :: r12(3), d12, r32(3), d32, costheta
  
      r12 = rel_pos(colloids%pos(:,(enzyme_idx-1)*3+1), colloids%pos(:,(enzyme_idx-1)*3+2), solvent_cells%edges)
      d12 = norm2(r12)
      r32 = rel_pos(colloids%pos(:,(enzyme_idx-1)*3+3), colloids%pos(:,(enzyme_idx-1)*3+2), solvent_cells%edges)
      d32 = norm2(r32)
      costheta = dot_product(r12, r32)/(d12*d32)
      en = angle_k * (acos(costheta)-link_angle(enzyme_idx))**2 / 2
  
    end function compute_bead_energy
  
negative.code:     solvent%vel(:,i) = colloids%vel(:,enz_2)
  
      ! compute conformational energy difference when changing the angle
      en1 = compute_bead_energy(enzyme_idx, factor=-1)
      link_angle(enzyme_idx) = link_angles(1)
      en2 = compute_bead_energy(enzyme_idx, factor=1)
  
      p = solvent_cells%cartesian_indices(solvent%pos(:, i))
  
  