anchor.diff:  
     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
  
positive.code: 
    end function compute_temperature
  
  end module mpcd
  
negative.code:   private
  
    public :: compute_temperature, simple_mpcd_step
    public :: wall_mpcd_step
    public :: mpcd_stream
  
  contains
  
  