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::outofrange: 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 __ITENSORLATTICEHONEYCOMBH
define __ITENSORLATTICEHONEYCOMBH
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 = NxNy;
//auto N = 2NxNy; // transfer to x direction 2Nx; 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!