anchor.diff:    call e_solvent% append(solvent% pos, 0, 0.d0)
     call e_solvent_v% append(solvent% vel, 0, 0.d0)
   
  +  wall_v = 0
     do i = 1, 20
  -     call simple_mpcd_step(solvent, solvent_cells)
  +     call wall_mpcd_step(solvent, solvent_cells, mt, [ 0.8d0, 1.2d0 ], wall_v, [10,10])
        v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
   
  -     solvent% pos = solvent% pos + 0.1 * solvent% vel
  -     do j= 1, 3
  -        solvent% pos(:, j) = modulo( solvent% pos(:, j), L(j)*1.d0 )
  -     end do
  +     call mpcd_stream(solvent, solvent_cells, 0.1d0)
   
        call solvent% sort(solvent_cells)
        call e_solvent% append(solvent% pos, i, i*1.d0)
  
positive.code:   !!
    !! Use the method from Marsaglia & Bray (1964) to generate numbers from a
    !! normal distribution.
    function localnormal() result(r)
      double precision :: r
  
      double precision :: u1, u2, radius
  
negative.code:     else
         found = .false.
         do while (.not. found)
            u1 = 2*genrand_real1(state) - 1
            u2 = 2*genrand_real1(state) - 1
            radius = (u1**2+u2**2)
            if ( ( radius < 1 ) .and. (radius > 0) ) found = .true.
         end do
  