anchor.diff:    call solvent_io%image%append(solvent%image)
     call solvent_io%species%append(solvent%species)
   
  +  call rho_xy_el%close()
     call dimer_io%close()
     call hfile%close()
     call h5close_f(error)
  
positive.code: 
    epsilon = 1.d0
    sigma(1,:) = [sigma_C, sigma_N]
    sigma(2,:) = [sigma_C, sigma_N]
    sigma_cut = sigma*3**(1.d0/6.d0)
    shift = maxval(colloid_lj% cut) +0.25
    call walls_colloid_lj% init(epsilon(1:2,:), sigma(1:2,:), sigma_cut(1:2,:), shift)
  
    mass(1) = rho * sigma_C**3 * 4 * 3.14159265/3
    mass(2) = rho * sigma_N**3 * 4 * 3.14159265/3
  
negative.code:   end do
  
    call solvent_cells%init(L, 1.d0,has_walls = .true.)
  
    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_loop*N_MD_steps, time=N_loop*N_MD_steps*dt)
    call h5gclose_f(fields_group, error)
  
    colloids% pos(3,:) = solvent_cells% edges(3)/2.d0
    if (order) then
  