Hi, thanks for the question. Are you interested in doing this for the C++ or Julia version? I’ll assume C++ and for Julia it’s actually quite similar.
We have an example code of doing a 2D system you can use here:
http://itensor.org/docs.cgi?vers=cppv3&page=formulas/2d_dmrg
The only difference from a 1D code in that example is that the AutoMPO part of the code runs over site numberings which are obtained from an array of “LatticeBond” objects that are returned from the triangularLattice function. You can see the code for this function here:
https://github.com/ITensor/ITensor/blob/v3/itensor/mps/lattice/triangular.h
Despite the apparent complexity of that function (triangularLattice) it is actually doing something really simple: it is making an array of pairs of integers, such as (n,n+1) and (n,n+Ny-1), etc. which define the sites on which the Hamiltonian terms should act, using a 1D ordering of the sites. This is because an MPS always has a 1D ordering of its sites, so to implement a 2D Hamiltonian you have to define some mapping of the 2D bonds (pairs of sites) to pairs of sites in the 1D, MPS ordering. The AutoMPO system takes care of the rest. Usually the 2D to 1D mapping used is just a “zig zag” pattern that goes up each column of the 2D system.
To implement the honeycomb lattice, then, all you have to do is come up with such a 2D to 1D mapping. I would recommend reshaping the honeycomb lattice first to a square lattice but with some bonds missing and perhaps some bonds added if necessary. Then on paper enumerate the sites from 1,2,...,N where N=Nx*Ny. Now make a function which generates the integer pairs of all the 2D Hamiltonian bonds following this enumeration and returns an array of structs (data structures) of these integer pairs.
A helpful practice step is to modify the 2D DMRG code linked above to just print out the values bnd.s1 and bnd.s2. Compare these values to the numberings of sites of the triangular lattice following the zig zag convention I mentioned above.
Best regards,
Miles