anchor.diff:             com_pos = sum((colloids%pos(:,:)+ &
                   colloids%image(:,:)*spread(solvent_cells%edges, dim=2, ncopies=colloids%Nmax) ) &
                   , dim=2)/ colloids%Nmax
  -           head_pos = colloids%pos(:,head)+colloids%image(:,head)*solvent_cells%edges
  -           unit_r = rel_pos(head_pos, com_pos, solvent_cells%edges)
  -           unit_r = unit_r / norm2(unit_r)
  +           if (do_rattle .and. (.not. do_quaternion)) then
  +              rigid_janus%pos = com_pos
  +              rigid_janus%vel = v_com
  +           end if
  +           unit_r = get_unit_r()
              call axial_cf%add_fast((i-equilibration_loops)*N_MD_steps+j-1, v_com, unit_r)
              call janus_pos_el%append(rigid_janus%pos)
              call janus_vel_el%append(rigid_janus%vel)
  
positive.code: !! \param epsilon_N             interaction parameter of N sphere with both solvent species (2 elements)
  !! \param data_filename         filename for input Janus coordinates
  !! \param epsilon_colloid       interaction parameter for colloid-colloid interactions
  !! \param link_treshold         distance criterion for finding rigid-body links
  !! \param do_read_links         read link information from a file
  !! \param links_file            filename for the links data
  
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
  