anchor.diff:       end if
   
        call varia%tic()
  -     call compute_rho_xy
  +     if (store_rho_xy) then
  +        call compute_rho_xy
  +        call rho_xy_el%append(rho_xy)
  +     end if
        call varia%tac()
  -     call rho_xy_el%append(rho_xy)
   
        temperature = compute_temperature(solvent, solvent_cells)
        if (dimer) then
  
positive.code:   call solvent_cells%init(L, 1.d0,has_walls = .true.)
    call vx% init(0.d0, solvent_cells% edges(3), L(3))
  
    allocate(rho_xy(N_species, L(2), L(1)))
    call h5gcreate_f(hfile%id, 'fields', fields_group, error)
    call rho_xy_el%create_time(fields_group, 'rho_xy', rho_xy, ior(H5MD_LINEAR,H5MD_STORE_TIME), &
         step=N_MD_steps, time=N_MD_steps*dt)
    call elem_vx% create_time(fields_group, 'vx', vx% data, ior(H5MD_TIME, H5MD_STORE_TIME))
    call elem_vx_count% create_time(fields_group, 'vx_count', vx% count, ior(H5MD_TIME, H5MD_STORE_TIME))
  
negative.code: 
    call h5gclose_f(timers_group, error)
  
    if (store_rho_xy) call rho_xy_el%close()
    call elem_vx% close()
    call elem_vx_count% close()
    call dimer_io%close()
  