anchor.diff:  
     call datafile% create('data_setup_simple_fluid.h5', 'RMPCDMD:setup_simple_fluid', '0.0 dev', 'Pierre de Buyl')
   
  -  call tz% init(&
  -       solvent_cells% origin(3), &
  -       solvent_cells% origin(3) + solvent_cells% L(3) * solvent_cells% a, &
  -       solvent_cells% L(3))
  +  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)
  
positive.code:   integer, allocatable :: seed(:)
  
    double precision :: v_com(3), wall_v(3,2), wall_t(2)
  
    double precision, parameter :: tau = 0.1d0
  
  
negative.code:   call e_solvent_v% create_time(solvent_group, 'velocity', solvent% vel, ior(H5MD_TIME, H5MD_STORE_TIME))
    call elem_tz% create_time(datafile% observables, 'tz', tz% data, ior(H5MD_TIME, H5MD_STORE_TIME))
    call elem_tz_count% create_time(datafile% observables, 'tz_count', tz% count, ior(H5MD_TIME, H5MD_STORE_TIME))
    call elem_rhoz% create_time(datafile% observables, 'rhoz', rhoz% data, ior(H5MD_TIME, H5MD_STORE_TIME))
    call elem_T% create_time(datafile% observables, 'temperature', T, ior(H5MD_TIME, H5MD_STORE_TIME))
  
    call solvent% sort(solvent_cells)
    call solvent_cells%count_particles(solvent% pos)
  
  