anchor.diff:       call e_solvent_v% append(solvent% vel, i, i*1.d0)
   
        call solvent_cells%count_particles(solvent% pos)
  -     write(*,*) compute_temperature(solvent, solvent_cells), sum(solvent% vel**2)/(3*solvent% Nmax), v_com
  +     write(13,*) compute_temperature(solvent, solvent_cells, tz), sum(solvent% vel**2)/(3*solvent% Nmax), v_com
  +
     end do
   
     call e_solvent% close()
  
positive.code:   solvent% species = 1
    call solvent% random_placement(L*1.d0)
  
    call solvent_cells%init(L, 1.d0, .true.)
  
    call solvent_cells%count_particles(solvent% pos)
  
    call datafile% create('data_setup_simple_fluid.h5', 'RMPCDMD:setup_simple_fluid', '0.0 dev', 'Pierre de Buyl')
  
    call h5gcreate_f(datafile% particles, 'solvent', solvent_group, error)
    call h5gcreate_f(solvent_group, 'box', box_group, error)
  
negative.code: 
    call solvent% sort(solvent_cells)
  
    wall_v = 0
    do i = 1, 200
       call simple_mpcd_step(solvent, solvent_cells, mt)
       call mpcd_stream(solvent, solvent_cells, 0.1d0)
       call solvent% sort(solvent_cells)
       call solvent_cells%count_particles(solvent% pos)
    end do
  
    do i = 1, 200
       call simple_mpcd_step(solvent, solvent_cells, mt)
       v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
  
       call mpcd_stream(solvent, solvent_cells, 0.1d0)
  