Dear ITensor team,

Recently, I have construct Hubbard model in honeycomb lattice, but I execute the code, I find such error :

libc++abi.dylib: terminating with uncaught exception of type std::out*of*range: vector

zsh: abort ./hubbard_2d

I guess this error is due to : auto H = toMPO(ampo);

when I convert other lattice (for example square lattice or triangular lattice), the code is works! I just imitate square.h formalism to construct honeycomb lattice. Is there any subtle in the constructing ? the following is my naive code :

## include "itensor/all.h"

## include "itensor/mps/lattice/honeycomb.h"

using namespace itensor;

int main(int argc, char *argv[])

{

int Nx = 4;

int Ny = 8;

double U = 8.0;

if(argc > 3)

U = std::stof(argv[3]);

if(argc > 2)

Ny = std::stoi(argv[2]);

if(argc > 1)

Nx = std::stoi(argv[1]);

auto N = Nx * Ny;

auto sites = Electron(N);

auto t = 1.0;

auto ampo = AutoMPO(sites);

auto lattice = honeycombLattice(Nx, Ny, {"YPeriodic = ", true});

//println(lattice);

for(auto j : lattice)

{

ampo += -t, "Cdagup", j.s1, "Cup", j.s2;

ampo += -t, "Cdagup", j.s2, "Cup", j.s1;

ampo += -t, "Cdagdn", j.s1, "Cdn", j.s2;

ampo += -t, "Cdagdn", j.s2, "Cdn", j.s1;

}

for(auto j : range1(N))

{

ampo += U, "Nupdn", j;

}

println(ampo);

auto H = toMPO(ampo);

auto state = InitState(sites);

for(auto j : range1(N))

{

state.set(j, (j % 2 == 1 ? "Up" : "Dn"));

}

auto sweeps = Sweeps(15);

sweeps.maxdim() = 20, 60, 100, 100, 200, 400, 800, 2000, 3000;

sweeps.noise() = 1E-7, 1E-8, 1E-10, 0;

sweeps.cutoff() = 1E-6;

PrintData(sweeps);

auto psi0 = randomMPS(state);

auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet=",true});

PrintData(Nx);

PrintData(Ny);

PrintData(U);

PrintData(t);

PrintData(totalQN(psi));

PrintData(maxLinkDim(psi));

PrintData(energy);

return 0;

}

and this is my naive constructing honeycomb lattice code:

## ifndef __ITENSOR*LATTICE*HONEYCOMB*H*

## define __ITENSOR*LATTICE*HONEYCOMB*H*

## include "itensor/mps/lattice/latticebond.h"

namespace itensor {

LatticeGraph inline

honeycombLattice(int Nx,

int Ny,

Args const& args = Args::global())

{

auto yperiodic = args.getBool("YPeriodic",true);

// Periodicity on y is meaningless for one dimensional chain or a ladder

yperiodic = yperiodic && (Ny > 2);

auto N = Nx*Ny;
//auto N = 2*Nx

*Ny; // transfer to x direction 2*Nx; y direction Ny

auto Nbond = 3*N-Ny + (yperiodic ? 0 : -Nx);

LatticeGraph latt;

latt.reserve(Nbond);

```
for(int I = 1; I <= N; ++I)
{
int xu = (I-1)/Ny+1;
int yu = (I-1)%Ny+1;
int ni = 2*(I-1);
latt.emplace_back(ni+2,ni+1,2*xu-1,yu,2*xu,yu); // inter-unitcell
//X-direction bond
if(xu < Nx) latt.emplace_back(ni+1,ni+2+2*Ny,2*xu,yu,2*(xu+1)-1,yu);
if(Ny > 1)
{
//diagonal bond
if(yu < Ny) latt.emplace_back(ni+2,ni+1*2+1,2*xu-1,yu,2*xu,yu+1);
//Periodic bond
if(yperiodic && yu == 1) latt.emplace_back(ni+1,ni+2*(Ny-1)+2,2*xu,yu,2*xu-1,yu+Ny);
}
}
//println(latt.size());
//println(Nbond);
if(int(latt.size()) != Nbond) Error("Honeycomb latt wrong number of bonds");
return latt;
}
```

} //namespace itensor

## endif

Many thanks!