anchor.diff:    call random_seed(put = seed)
     deallocate(seed)
   
  +  call system_clock(count=clock)
  +  call init_genrand(mt, int(clock, c_long))
  +
     call h5open_f(error)
   
     L = [8, 5, 5]
   
     call solvent% init(N)
   
  -  do j=1, solvent% Nmax
  -     solvent% vel(1, j) = localnormal()
  -     solvent% vel(2, j) = localnormal()
  -     solvent% vel(3, j) = localnormal()
  -  end do
  +  call mt_normal_data(solvent% vel, mt)
     v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
     solvent% vel = solvent% vel - spread(v_com, dim=2, ncopies=size(solvent% vel, dim=2))
   
  
positive.code:   use hdf5
    use h5md_module
    use mpcd
    implicit none
  
    type(cell_system_t) :: solvent_cells
    type(particle_system_t) :: solvent
  
negative.code:   !!
    !! Use the method from Marsaglia & Bray (1964) to generate numbers from a
    !! normal distribution.
    function localnormal(state) result(r)
      type(mt19937ar_t), intent(inout) :: state
      double precision :: r
  
      double precision :: u1, u2, radius
  