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:         call vx% reset()
       end if
  
       call solvent% sort(solvent_cells)
       call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells)
  
  
negative.code:         end if 
       end do
    end subroutine buffer_particles
  
    subroutine compute_rho_xy
      integer :: i, s, ix, iy
  
      rho_xy = 0
      do i = 1, solvent%Nmax
         s = solvent%species(i)
         if (s <= 0) continue
         ix = modulo(floor(solvent%pos(1,i)/solvent_cells%a), L(1)) + 1
         iy = modulo(floor(solvent%pos(2,i)/solvent_cells%a), L(2)) + 1
         rho_xy(s, iy, ix) = rho_xy(s, iy, ix) + 1
      end do
  
    end subroutine compute_rho_xy
  
  end program setup_single_dimer
  