anchor.diff:    private
   
     public :: compute_temperature, simple_mpcd_step
  +  public :: wall_mpcd_step
  +  public :: mpcd_stream
   
   contains
   
  
positive.code: 
    end function compute_temperature
  
  end module mpcd
  
negative.code: 
    end subroutine simple_mpcd_step
  
    !! Lamura, Gompper, Ihle and Kroll, EPL 56, 319-325 (2001)
    !! http://dx.doi.org/10.1209/epl/i2001-00522-9
    subroutine wall_mpcd_step(particles, cells, state, wall_temperature, wall_v, wall_n)
      use hilbert
      class(particle_system_t), intent(in) :: particles
      class(cell_system_t), intent(in) :: cells
      type(mt19937ar_t), intent(inout) :: state
      double precision, intent(in) :: wall_temperature(2)
      double precision, intent(in) :: wall_v(3,2)
      integer, intent(in) :: wall_n(2)
  
      integer :: i, start, n
      integer :: cell_idx, count
      double precision :: local_v(3), omega(3,3), vec(3)
      integer :: n_virtual
      integer :: cell(3)
      integer :: wall_idx
      double precision :: virtual_v(3)
  
      do cell_idx = 1, cells% N
         if (cells% cell_count(cell_idx) <= 1) cycle
  
         start = cells% cell_start(cell_idx)
         n = cells% cell_count(cell_idx)
         count = count + 1
         n_virtual = 0
  
         ! Find whether we are in a wall cell
         cell = compact_h_to_p(cell_idx, cells% M)
         if (cell(3) == 1) then
            wall_idx = 1
            vec = wall_v(:,wall_idx)
         else if (cell(3) == cells% L(3)) then
            wall_idx = 2
            vec = wall_v(:,wall_idx)
         else
            wall_idx = -1
            vec = local_v
         end if
  
         if ( (wall_idx > 0) .and. (n < wall_n(wall_idx)) ) then
            n_virtual = wall_n(wall_idx) - n
            call mt_normal_data(virtual_v, state)
            virtual_v = virtual_v * sqrt(n_virtual*wall_temperature(wall_idx))
         end if
  
         local_v = 0
         do i = start, start + n - 1
            local_v = local_v + particles% vel(:, i)
         end do
         if (n_virtual > 0) then
            local_v = (local_v + virtual_v) / wall_n(wall_idx)
         else
            local_v = local_v / n
         end if
  
         vec = rand_sphere(state)
         omega = &
              reshape( (/ &
              vec(1)**2, vec(1)*vec(2) + vec(3), vec(1)*vec(3) - vec(2) ,&
              vec(2)*vec(1) - vec(3) , vec(2)**2 , vec(2)*vec(3) + vec(1),&
              vec(3)*vec(1) + vec(2), vec(3)*vec(2) - vec(1), vec(3)**2 &
              /), (/3, 3/))
  
         do i = start, start + n - 1
            particles% vel(:, i) = local_v + matmul(omega, (particles% vel(:, i)-local_v))
         end do
  
      end do
  
    end subroutine wall_mpcd_step
  
    function compute_temperature(particles, cells) result(te)
      class(particle_system_t), intent(in) :: particles
      class(cell_system_t), intent(in) :: cells
  