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/all.h"
include "TStateObserver.h"
include "S2.h"
include <math.h> /* sqrt */
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;
args.add("MaxDim",maxdim);
args.add("Cutoff",cutoff);
args.add("Verbose",verbose);
args.add("Method","DensityMatrix");
//args.add("Method","Fit");
//
// 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;
targs.add("TimeStepNum",tt);
targs.add("Time",tsofar);
targs.add("TotalTime",ttotal);
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
Javad