anchor.diff:  
     end subroutine read_links
   
  +  function get_unit_r() result(u_r)
  +    double precision :: u_r(3)
  +    double precision :: r1(3), r2(3), r3(3)
  +
  +    r1 = rel_pos(colloids%pos(:,1), com_pos, solvent_cells%edges)
  +    r2 = rel_pos(colloids%pos(:,2), com_pos, solvent_cells%edges)
  +    r3 = rel_pos(colloids%pos(:,3), com_pos, solvent_cells%edges)
  +
  +    u_r = (r1 + alpha*r2 + beta*r3) / z0
  +
  +  end function get_unit_r
  +
   end program single_janus_pbc
  
positive.code:   call radial_hist_el%close()
    call n_solvent_el%close()
    call q_el%close()
    call janus_pos_el%close()
    call janus_vel_el%close()
    call hfile%close()
  
negative.code:     end do
    end subroutine refuel
  
    subroutine read_links(filename, n_links, links, links_d)
      character(len=*), intent(in) :: filename
      integer, intent(out) :: n_links
  