anchor.diff:      if (thermostat) error stop 'thermostatting not implemented'
   
       call particles%time_step%tic()
  +    sin_alpha = sin(alpha)
  +    cos_alpha = cos(alpha)
       !$omp parallel
       thread_id = omp_get_thread_num() + 1
       !$omp do private(start, n, local_v, i, vec, omega)
  
positive.code:           vec = rand_sphere(state(thread_id))
            omega = &
                 reshape( (/ &
                 cos(alpha)+vec(1)**2*(1-cos(alpha)), vec(1)*vec(2)*(1-cos(alpha)) - vec(3)*sin(alpha), &
                 vec(1)*vec(3)*(1-cos(alpha)) + vec(2)*sin(alpha) ,&
                 vec(2)*vec(1)*(1-cos(alpha)) + vec(3)*sin(alpha) , cos(alpha)+vec(2)**2*(1-cos(alpha)) ,&
                 vec(2)*vec(3)*(1-cos(alpha)) - vec(1)*sin(alpha),&
                 vec(3)*vec(1)*(1-cos(alpha)) - vec(2)*sin(alpha), vec(3)*vec(2)*(1-cos(alpha)) + vec(1)*sin(alpha),&
                 cos(alpha) + vec(3)**2*(1-cos(alpha)) &
                 /), (/3, 3/))
  
            do i = start, start + n - 1
  
negative.code:     integer :: i, start, n
      integer :: cell_idx
      double precision :: local_v(3), omega(3,3), vec(3)
      double precision :: sin_alpha, cos_alpha
      logical :: thermostat
      integer :: thread_id
  
  