Dear ITensor team,
Thanks for your wonderful tool about Tensor Network. I have some silly questions about spectrum function. I want to repeat PHYSICAL REVIEW B 77, 134437 2008 Fig.4(t=10 (vertical green dot line ) results S^{zz}(0,10) about 0.02, spin 1 Heisenberg chain)
According to definition: C^{z}(x,t)=<g.s|exp{iHt}S^{z}(x)exp(-iHt)S^{z}(0)|g.s> right?
my naive strategy as follow:
1) S^{z}(0) act on ground state |g.s>(single site operator act on MPS get new MPS psi1)
2) exp(-iHt) act new state psi1(Trotter evolution) get new MPS state psi2
3) S^{z}(x) act new psi2 to get new MPS psi3(same as step 1 )
4) exp(iHt) act psi3 get new MPS psi4(Note that in the code I want to construct exp(Itstephterm) gate. Naively, I replaced tstep to -tstep because if BondGate argument is tstep we get exp(-Itstephterm), is that correct? (please see my naive code))
5) Finally, C^{z}(x,t)=innerC(psi,psi4) right?
According to this strategy, I write a naive code, but my answer is not consistent with that paper. So I don't know is there any wrong with my naive code? Thank you very much !
int main()
{
ofstream outfile;
outfile.open("cor_t.dat",ios::app);
int N = 50; //number of sites
Real tt = 10.0;
Real tstep = 0.05; //time step (smaller is generally more accurate)
Real ttotal = tt; //total time to evolve
Real cutoff = 1E-8; //truncation error cutoff when restoring MPS form
auto sites = SpinOne(N);
auto state = InitState(sites);
auto ampo = AutoMPO(sites);
for(auto j : range1(N))
{
state.set(j,j%2==1?"Up":"Dn");
}
auto psi = MPS(state);
for(int b = 1; b < N; ++b)
{
ampo += 0.5,"S+",b,"S-",b+1;
ampo += 0.5,"S-",b,"S+",b+1;
ampo += "Sz",b,"Sz",b+1;
}
auto H = MPO(ampo);
auto sweeps = Sweeps(30);
sweeps.maxm() = 100,200,800,800,1000;
sweeps.cutoff() = 1E-10;
sweeps.niter() = 2;
sweeps.noise() = 1E-7,1E-8,1E-9;
println(sweeps);
auto energy = dmrg(psi,H,sweeps,"Quiet");
println(energy);
auto psi0 = psi;
// S^{z}{1} act on RHS
auto Sz1 = op(sites,"Sz",1);
println(Sz1);
auto newpsi = Sz1*psi(1);
newpsi.noPrime();
psi.set(1,newpsi);
//Create a std::vector (dynamically sizeable array)
//to hold the Trotter gates
auto gates = vector();
//Create the gates exp(-itstep/2hterm)
//and add them to gates
for(int b = 1; b <= N-1; ++b)
{
auto hterm = op(sites,"Sz",b)op(sites,"Sz",b+1);
hterm += 0.5op(sites,"S+",b)op(sites,"S-",b+1);
hterm += 0.5op(sites,"S-",b)*op(sites,"S+",b+1);
auto g = BondGate(sites,b,b+1,BondGate::tReal,tstep/2.,hterm);
gates.push_back(g);
}
//Create the gates exp(-itstep/2hterm) in reverse
//order (to get a second order Trotter breakup which
//does a time step of "tstep") and add them to gates
for(int b = N-1; b >= 1; --b)
{
auto hterm = op(sites,"Sz",b)op(sites,"Sz",b+1);
hterm += 0.5op(sites,"S+",b)op(sites,"S-",b+1);
hterm += 0.5op(sites,"S-",b)*op(sites,"S+",b+1);
auto g = BondGate(sites,b,b+1,BondGate::tReal,tstep/2.,hterm);
gates.push_back(g);
}
//Time evolve, overwriting psi when done
gateTEvol(gates,ttotal,tstep,psi,{"Cutoff=",cutoff,"Verbose=",true});
// S^{z}{1} act RHS
auto newpsi1 = Sz1*psi(1);
newpsi1.noPrime();
psi.set(1,newpsi1);
// exp(itstepH) act RHS
//Create the gates exp(itstep/2hterm)
//and add them to gates
for(int b = 1; b <= N-1; ++b)
{
auto hterm = op(sites,"Sz",b)op(sites,"Sz",b+1);
hterm += 0.5op(sites,"S+",b)op(sites,"S-",b+1);
hterm += 0.5op(sites,"S-",b)*op(sites,"S+",b+1);
auto g = BondGate(sites,b,b+1,BondGate::tReal,-tstep/2.,hterm);
gates.push_back(g);
}
//Create the gates exp(itstep/2*hterm) in reverse
//order (to get a second order Trotter breakup which
//does a time step of "tstep") and add them to gates
for(int b = N-1; b >= 1; --b)
{
auto hterm = op(sites,"Sz",b)op(sites,"Sz",b+1);
hterm += 0.5op(sites,"S+",b)op(sites,"S-",b+1);
hterm += 0.5op(sites,"S-",b)*op(sites,"S+",b+1);
auto g = BondGate(sites,b,b+1,BondGate::tReal,-tstep/2.,hterm);
gates.push_back(g);
}
//Time evolve, overwriting psi when done
gateTEvol(gates,ttotal,tstep,psi,{"Cutoff=",cutoff,"Verbose=",true});
printfln("Maximum MPS bond dimension after time evolution is %d",maxLinkDim(psi));
outfile <<tt<<" "<<real(innerC(psi0,psi))<< endl;
outfile.close();
return 0;
}