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:     integer :: i, start, n
      integer :: cell_idx
      double precision :: local_v(3), omega(3,3), vec(3)
      integer :: n_virtual
      integer :: cell(3)
      integer :: wall_idx
  
negative.code:     end if
  
      call particles%time_step%tic()
      sin_alpha = sin(alpha)
      cos_alpha = cos(alpha)
      !$omp parallel private(thread_id)
      thread_id = omp_get_thread_num() + 1
      !$omp do private(cell_idx, start, n, n_virtual, virtual_v, cell, wall_idx, local_v, i, vec, omega)
  