anchor.diff:      do i = this%i_start, this%i_stop
          f = ps%force(:,i)
          this%force = this%force + f
  -       this%torque = this%torque + cross(ps%pos(:,i)+ps%image(:,i)*edges-this%pos, f)
  +       pos = qrot(this%q, this%pos_body(:,i))
  +       this%torque = this%torque + cross(pos, f)
       end do
   
     end subroutine rigid_body_compute_force_torque
  
positive.code:     this%vel = this%vel + this%force*dt/(2*this%mass)
  
      do i = this%i_start, this%i_stop
         ps%vel(:,i) = this%vel + cross(omega, ps%pos(:,i)+ps%image(:,i)*edges-this%pos)
      end do
  
    end subroutine rigid_body_vv2
  
negative.code:   end subroutine rigid_body_init
  
    !! Compute torque and force on the rigid body in the lab frame
    subroutine rigid_body_compute_force_torque(this, ps)
      class(rigid_body_t), intent(inout) :: this
      type(particle_system_t), intent(in) :: ps
  
      double precision :: pos(3), f(3)
      integer :: i
  
      this%force = 0
  