anchor.diff:       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
  +     T = compute_temperature(solvent, solvent_cells, tz)
  +     write(13,*) T, sum(solvent% vel**2)/(3*solvent% Nmax), v_com
  +     call elem_T% append(T, i, i*tau)
  +
  +     call compute_rho(solvent, rhoz)
   
        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 elem_tz% append(tz% data, i, i*tau)
  +        call elem_tz_count% append(tz% count, i, i*tau)
           call tz% reset()
  +        rhoz% data = rhoz% data / (50.d0 * rhoz% dx)
  +        call elem_rhoz% append(rhoz% data, i, i*tau)
  +        rhoz% data = 0
        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)
  +  call e_solvent% append(solvent% pos, i, i*tau)
  +  call e_solvent_v% append(solvent% vel, i, i*tau)
   
     clock = 0
     do i = 1 , solvent_cells% N
  
positive.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 solvent% sort(solvent_cells)
    call solvent_cells%count_particles(solvent% pos)
  
  
negative.code:   integer, allocatable :: seed(:)
  
    double precision :: v_com(3), wall_v(3,2), wall_t(2)
    double precision :: T
  
    double precision, parameter :: tau = 0.1d0
  
  