anchor.diff:  
           call rattle_body_vel(colloids, links, links_d, dt, solvent_cells% edges, rattle_pos_tolerance)
   
  +        if (sampling) then
  +           v_com = sum(colloids%vel, dim=2)/colloids%Nmax
  +           com_pos = sum((colloids%pos(:,:)+ &
  +                colloids%image(:,:)*spread(solvent_cells%edges, dim=2, ncopies=colloids%Nmax) ) &
  +                , dim=2)/ colloids%Nmax
  +           unit_r = rel_pos(colloids%pos(:,head), com_pos, solvent_cells%edges)
  +           unit_r = unit_r / norm2(unit_r)
  +           call axial_cf%add_fast((i-equilibration_loops)*N_MD_steps+j-1, v_com, unit_r)
  +        end if
  +
  +        if ((sampling) .and. (modulo(j, colloid_sampling)==0)) then
  +           call janus_io%position%append(colloids%pos)
  +           call janus_io%velocity%append(colloids%vel)
  +           call janus_io%image%append(colloids%image)
  +        end if
  +
           call time_flag%tic()
           call flag_particles
           call time_flag%tac()
  
positive.code:   solvent% force_old = solvent% force
    colloids% force_old = colloids% force
  
    write(*,*) 'Running for', N_loop, 'loops'
    call main%tic()
    do i = 1, N_loop
       if (modulo(i,5) == 0) write(*,'(i05)',advance='no') i
       md_loop: do j = 1, N_MD_steps
          call md_pos(solvent, dt)
  
  
negative.code:   double precision :: link_treshold
    integer :: i_link
    double precision :: dist, rattle_pos_tolerance, rattle_vel_tolerance
    double precision :: unit_r(3)
    double precision :: com_pos(3)
    integer :: head
  
    integer, parameter :: block_length = 8
    type(axial_correlator_t) :: axial_cf
    integer(HID_T) :: correlator_group
  
    integer :: equilibration_loops
    integer :: colloid_sampling
    logical :: sampling
  
    type(args_t) :: args
    character(len=144) :: data_filename
  