anchor.diff:    type(polar_fields_t) :: polar
     integer(HID_T) :: polar_id
   
  -  logical :: do_rattle, do_read_links, do_lennard_jones, do_elastic
  +  logical :: do_rattle, do_read_links, do_lennard_jones, do_elastic, do_quaternion
     integer, allocatable :: links(:,:)
     double precision, allocatable :: links_d(:)
     double precision :: link_treshold
  
positive.code:   double precision :: unit_r(3)
    double precision :: com_pos(3), head_pos(3)
    integer :: head
  
    integer, parameter :: block_length = 8
    type(axial_correlator_t) :: axial_cf
  
negative.code: 
          call md_vel(solvent, dt)
  
          if (do_quaternion) then
             call rigid_janus%compute_force_torque(colloids, solvent_cells%edges)
             call rigid_janus%vv2(colloids, dt, solvent_cells%edges)
          else
             call varia%tic()
             do k = 1, colloids%Nmax
                colloids% vel(:,k) = colloids% vel(:,k) + &
                     dt * ( colloids% force(:,k) + colloids% force_old(:,k) ) / (2 * colloids%mass(colloids%species(k)))
             end do
             call varia%tac()
          end if
  
          if (do_rattle) then
             call rattle_body_vel(colloids, links, links_d, dt, solvent_cells% edges, rattle_pos_tolerance)
          end if
  
          if (sampling) then
             v_com = sum(colloids%vel, dim=2)/colloids%Nmax
  