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:     end do
    end subroutine refuel
  
    !> Get index of the head bead in the Janus particle
    function get_head_idx(filename) result(head)
      character(len=*), intent(in) :: filename
      integer :: head
  
      integer(HID_T) :: file_id, obj_id, attr_id, space_id
      integer :: error
      integer(HSIZE_T) :: dims(1)
      logical :: is_simple
  
      call h5fopen_f(filename, H5F_ACC_RDONLY_F, file_id, error)
      call h5oopen_f(file_id, 'particles/janus/position', obj_id, error)
      call h5aopen_f(obj_id, 'head', attr_id, error)
      call h5aget_space_f(attr_id, space_id, error)
      call h5sis_simple_f(space_id, is_simple, error)
      if (.not. is_simple) error stop 'head attribute not simple'
      call h5aread_f(attr_id, H5T_NATIVE_INTEGER, head, dims, error)
      call h5sclose_f(space_id, error)
      call h5aclose_f(attr_id, error)
      call h5oclose_f(obj_id, error)
      call h5fclose_f(file_id, error)
  
      ! The index in the H5MD file is 0-based.
      head = head + 1
  
    end function get_head_idx
  
    subroutine read_links(filename, n_links, links, links_d)
      character(len=*), intent(in) :: filename
      integer, intent(out) :: n_links
  
negative.code: 
             call rattle_body_pos(colloids, links, links_d, dt, solvent_cells% edges, rattle_pos_tolerance)
  
             com_pos = modulo(sum((colloids%pos(:,:)+ &
                  colloids%image(:,:)*spread(solvent_cells%edges, dim=2, ncopies=colloids%Nmax) ) &
                  , dim=2)/ colloids%Nmax, solvent_cells%edges)
  
          else if (do_quaternion) then
             call rigid_janus%vv1(colloids, dt, quaternion_treshold, solvent_cells%edges)
             com_pos = modulo(rigid_janus%pos, solvent_cells%edges)
          end if
  
          so_max = solvent% maximum_displacement()
  