Yes this is possible to do, but it can get rather technical to do it efficiently. The easiest thing you can do with ITensor at the moment is the use the AutoMPO system and just input all N^2 terms connecting a site i to a site j, weighted by an 1/|i-j| coefficient. This would work for small to medium size systems, and with our new AutoMPO engine (post version 2.0.11) the resulting Hamiltonian ought to be compressed as much as possible.
For larger systems there are approaches in the literature you can implement using ITensor to help. One approach is to fit the 1/r interaction by a sum of exponentials. An exponentially decaying interaction can be represented exactly within an MPO. Two good papers discussing this approach are:
1. Crosswhite, Doherty, and Vidal, Phys. Rev. B 78, 035116 (2008)
2. Pirvu, Murg, Cirac, and Verstraete, New J. Phys. 12, 025012 (2010)
We (myself and Steve White) are currently writing up a paper with a different approach to long-range interactions that is too technical to describe here. But please look out for our new paper and I can share with you some pretty simple and nice code you could adapt for your own work if you need to study long range interactions.