anchor.diff:    type(particle_system_t) :: solvent
   
     type(profile_t) :: tz
  +  type(histogram_t) :: rhoz
   
     type(h5md_file_t) :: datafile
     type(h5md_element_t) :: elem
     type(h5md_element_t) :: e_solvent, e_solvent_v
     type(h5md_element_t) :: elem_tz, elem_tz_count
  +  type(h5md_element_t) :: elem_rhoz
  +  type(h5md_element_t) :: elem_T
     integer(HID_T) :: box_group, solvent_group
   
     integer, parameter :: N = 10800
  
positive.code:      call solvent% sort(solvent_cells)
       call solvent_cells%count_particles(solvent% pos)
  
       write(13,*) compute_temperature(solvent, solvent_cells, tz), sum(solvent% vel**2)/(3*solvent% Nmax), v_com
  
       if (modulo(i, 50) == 0) then
          call tz% norm()
          call elem_tz% append(tz% data, i, i*1.d0)
          call elem_tz_count% append(tz% count, i, i*1.d0)
          call tz% reset()
       end if
  
    end do
  
    call e_solvent% append(solvent% pos, i, i*1.d0)
    call e_solvent_v% append(solvent% vel, i, i*1.d0)
  
    clock = 0
    do i = 1 , solvent_cells% N
  
negative.code: 
    call datafile% create('data_setup_simple_fluid.h5', 'RMPCDMD:setup_simple_fluid', '0.0 dev', 'Pierre de Buyl')
  
    call tz% init(0.d0, solvent_cells% edges(3), L(3))
    call rhoz% init(0.d0, solvent_cells% edges(3), L(3))
    call h5gcreate_f(datafile% id, 'observables', datafile% observables, error)
  
    call h5gcreate_f(datafile% particles, 'solvent', solvent_group, error)
  