finiteT calculation of 2d lattice

+1 vote

Dear Itensor community

I have tried to extend the spin-1/2 chain finiteT example in the tutorial folder to a finiteT temperature calculation of 2D lattice (e.g: triangular lattice). Comparing with small 4*4 tiangular lattice with ED results, what I am getting is totally wrong.

Bellow you can find the code:

===========

include "itensor/util/print_macro.h"

using namespace std;
using namespace itensor;

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

//Get parameter file
if(argc != 2)
{
printfln("Usage: %s inputfile.",argv[0]);
return 0;
}
auto input = InputGroup(argv[1],"input");

//auto yperiodic = false;
auto Nx = input.getInt("Nx",4);
auto Ny = input.getInt("Ny",4);
auto N = Nx*Ny;

auto beta = input.getReal("beta",1);
auto tau = input.getReal("tau",0.01);

auto maxdim = input.getInt("maxdim",1000);
auto cutoff = input.getReal("cutoff",1E-11);

auto verbose = input.getYesNo("verbose",false);

Args args;

//
// Initialize the site degrees of freedom.
//
auto sites = SpinHalf(2*N,{"ConserveQNs=",false});

//
// Use the AutoMPO feature to create the
// next-neighbor Heisenberg model.
//
auto ampo = AutoMPO(sites);

// 4*4 triangular lattice  : bound between nearest neighbor sites


// sites with even number reserved for "ancilla", so lattice sites numbered in odd
auto site1 = {1,1,1,3,3,3,5,5,5,7,9,9,9,11,11,11,13,13,13,15,17,17,17,19,19,19,21,21,21,23};
auto site2 = {9,11,3,11,13,5,13,15,7,15,17,19,11,19,21,13,21,23,15,23,25,27,19,27,29,21,29,31,23,31};

 for(auto bnd : zip(site1,site2)){
auto [s1,s2] = bnd;
//printfln(" s1 : %i , s2 : %i ",bnd.s1,bnd.s2);
ampo += 0.5,"S+",s1,"S-",s2;
ampo += 0.5,"S-",s1,"S+",s2;
ampo +=       ,"Sz",s1,"Sz",s2;
}

auto H    = toMPO(ampo);
auto expH = toExpH(ampo,tau);
auto S2   = makeS2(sites,{"SkipAncilla=",true});
//auto Sz2 = makeTotSz2(sites,{"SkipAncilla=",true});

//
// Make initial 'wavefunction' which is a product
// of perfect singlets between neighboring sites
//
auto psi = MPS(sites);
for(int n = 1; n <= 2*N; n += 2)
{
auto s1 = sites(n);
auto s2 = sites(n+1);
auto wf = ITensor(s1,s2);
wf.set(s1(1),s2(2), ISqrt2);
wf.set(s1(2),s2(1), -ISqrt2);

ITensor D;
psi.ref(n) = ITensor(s1);
psi.ref(n+1) = ITensor(s2);
svd(wf,psi.ref(n),D,psi.ref(n+1));
psi.ref(n) *= D;
}

auto obs = TStateObserver(psi);
auto ttotal = beta/2.;
const int nt = int(ttotal/tau+(1e-9*(ttotal/tau)));
if(fabs(nt*tau-ttotal) > 1E-9)
{
Error("Timestep not commensurate with total time");
}
printfln("Doing %d steps of tau=%f",nt,tau);

auto targs = args;

auto En     = Vector(nt);
auto Sus    = Vector(nt);
auto SzSz   = Vector(nt);
auto Betas  = Vector(nt);
Real tsofar = 0;

std::ofstream enf("en_cutoffE-10_dmrg.txt");
std::ofstream susf("sus_cutoffE-10_dmrg.txt");
for(int tt = 1; tt <= nt; ++tt)
{
psi = applyMPO(expH,psi,args);

//Normalize wavefunction
psi.ref(1) /= norm(psi(1));
//psi.noPrime().normalize();

tsofar += tau;
obs.measure(targs);

//Record beta value
auto bb = (2*tsofar);
Betas(tt-1) = bb;

//
// Measure Energy
//
auto en = inner(psi,H,psi);
En(tt-1) = real(en)/N;

SzSz(tt-1) = 0.0;//(SziSzj*bb/3.)/N;
printfln("\n Temp : %.4f,  E/N :%.20f, Sus/N :%.20f",bb,en/N,Sus(tt-1));
enf  << format("%.14f %.14f %.14f\n",Betas(tt-1), 1.0/Betas(tt-1), En(tt-1));
}

enf.close();
susf.close();

return 0;
}


===========

Any comments or suggestion are welcome.

Regards

commented by (70.1k points)
Hi Javad, it could be any number of things. Is the lattice exactly the same, including whether the boundaries are periodic or not? What happens when you adjust the cutoff and maxdim settings? What are you setting those to?

What properties are you comparing to ED and how are you obtaining them?

Miles
commented by (420 points)
Dear Miles

boundary is "Open"
I compared "E/N"
I doubled checked lattice configuration in both my "ED" and "dmrg" scripts, it seems to be identical.
I considered two sets of input:

{
Nx      = 4
Ny      = 4
beta   = 10.0
tau    = 0.01
cutoff = 1E-10
maxdim = 1000

}

and
{
Nx      = 4
Ny      = 4
beta   = 10.0
tau    = 0.01
cutoff = 1E-12
maxdim = 2000

}

but E/N did not change at all buy changing input parameters!!!!

moreover, energy is not changing with temperature (just a  fix number)!!!!

so, something is wrong but I was not able to figure it out!!

commented by (70.1k points)
I copied your code and ran it. One line didn't compile for me so I'm surprised it compiled for you. It was this line:

ampo +=      , "Sz",s1,"Sz",s2;

where I needed to remove the leading comma before the first "Sz" to get it to compile.

But other than that, I do see the energy (E/N) changing with every time step that is applied. So when you say it does not change, do you mean the energy that is reached at the target imaginary time of beta? And when you say does not change, do you mean as a function of beta, or of cutoff and maxdim?

Thanks,
Miles
commented by (420 points)
Dear Miles

The problem is solved and I have great tresults.!!!

Regards
commented by (70.1k points)
Really glad to hear it! If you don't mind mentioning what the fix was I would be curious. Happy New Year to you.

Miles
commented by (420 points)
Dear Miles

Indeed, I had a stupid mistake in the "for loop" of the code were ruuning on cluster.

Insted of having loop over

for(auto bnd : zip(site1,site2)){
auto [s1,s2] = bnd;
ampo += 0.5,"S+",s1,"S-",s2;
ampo += 0.5,"S-",s1,"S+",s2;
ampo +=        "Sz",s1,"Sz",s2;
}

I was iterating over the follwoing one!!!:

for(auto bnd : zip(site1,site1)){
auto [s1,s2] = bnd;
ampo += 0.5,"S+",s1,"S-",s2;
ampo += 0.5,"S-",s1,"S+",s2;
ampo +=       "Sz",s1,"Sz",s2;
}

both where refereing to site1. ;-(

Many thanks for your kind help and consideration.
Happy New Year to you.