ITensor can calculate imaginary time correlation function?

Question

Dear ITensor,
I follow tutorial in order to calculate imaginary time correlation function such as <S^{z}{1}(L/2)S^{z}{1}>. here is my naive code. it can run but I get wrong answer (compare my QMC datas). so please help me check my naive code, Thanks!!!

int main()
{
ofstream outfile;
outfile.open("cor_t.dat",ios::app);

int length[7]={12,16,24,32,48,64,80};

for(int &N : length)
{
auto sites = SpinHalf(N);
auto state = InitState(sites);
for(int i = 1; i <= N; ++i)
{
//Neel state
state.set(i,i%2==1 ? "Up" : "Dn");
}
auto psi0 = MPS(state);
auto psi = MPS(state);

auto ampo = AutoMPO(sites);
//Make the Heisenberg Hamiltonian
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);

//
// Begin the DMRG calculation
//
auto energy = dmrg(psi,H,sweeps,"Quiet");
//

auto Sz_1=sites.op("Sz",1);

auto tau =0.1;
auto expHL0 = toExpH(ampo,tau);
auto expHR0 = toExpH(ampo,-tau);

auto args = Args("Cutoff=",1E-9,"Maxm=",3000);
auto ttotal = N/4;
auto nt = int(ttotal/tau+(1e-9*(ttotal/tau)));

// left part <tau(t)|S^{z}
for(int n = 1; n <= nt; ++n)
{
psi = exactApplyMPO(expHL0,psi,args);
normalize(psi);
}

auto newpsi= Sz_1*psi.A(1);
newpsi.noprime();
psi.setA(1,newpsi);

// Right part exp(-H\tau)*S^{z}|g.s>

auto newpsi0 = Sz_1*psi0.A(1);
newpsi0.noprime();

psi0.setA(1,newpsi0);

for(int n = 1; n <= nt; ++n)
{
psi0 = exactApplyMPO(expHR0,psi0,args);
normalize(psi0);
}

//Print(psi0);
auto result=overlapC(psi,psi0).real();
//Print(result);
outfile <<N<<" "<<abs(result)<< endl;
}
outfile.close();
return 0;
}

Comments

Hi, thanks for the question. I could imagine there are a number of different things about the code that could be leading to the incorrect results. One of them is that you consider a time step of tau=0.1 which is kind of large. You could also try tau=0.05, tau=0.02, or tau=0.0.1 to see if those work better.

But one potentially very serious issue is that on this line:

    auto newpsi0 = Sz_1*psi0.A*(1);

you are making the state Sz_1 * |0> where |0> is the Neel product state, not the ground state. Is this correct for what you want to do? In your comments, you write “Right part exp(-H\tau) Sz|g.s>” so it seems you meant to act Sz on the ground state, not on the Neel state.

Finally, given that you are studying a nearest-neighbor, one-dimensional system, I’d actually recommend you use the Trotter gate method for time evolution instead of the expH / AutoMPO approach to making the time evolution operator. The Trotter approach is generally more accurate.

Please let me know if that’s helpful.

Miles

---

Hi Miles,
     Thank you very much! your answer is great! but I have some other questions:
     1) According to your answer, I find L=12, imaginary time correlation function is good, but L=16,24...etc, it very close to zero (such as 10^(-34) ). So it has some optimize parameters I should choose??
     2)  if my one-dimensional hamiltonian have i and i +2 interaction term(for example \sigma_{I}\sigma_{I+2}). should I use Trotter approach? or I still use MPO time evolution?
                                                                                                                                                                         Best,
                                                                                                                                                                         Sugar

---

Hi,
So first of all, did you see my comment about psi0 not being the ground state? Is this what you want in the code or do you think it is an error?

(1) I’m not sure why you would see such a different behavior for different L. So I’d have to look at the code and results myself probably to understand this, but please do answer the question about psi0 first before I get to that.

(2) What is the Hamiltonian you wish to study eventually? Is it still basically one dimensional? If so, I would recommend using Trotter gates but with swap gates to handle further-neighbor interactions. The Trotter approach is very reliable for one-dimensional systems and will probably get you the best results. It’s also ok to use the MPO approach, but it generally has worse time-step errors and is less reliable for getting highly accurate results.

Miles

---

Hi Miles,
     1) I follow you comment (and correct my error ) but our L=16,24,....result close to zero. so maybe Trotter gate is better.
     2) yes, I want to study one-dimensional chain. for example cluster Ising model have next-nearest neighbor interaction(I and I+2), in this case we can use Trotter gates right?

Sugar

---

Hi,Miles,
    This is my correct Trotter gate code(follow tutorial and I want to calculate imaginary time evolution: <g.s|exp(H\tau)S^{z}exp(-H\tau)S^{z}|g.s>):

int main()
    {
ofstream outfile;
outfile.open("cor_t.dat",ios::app);

int length[7]={12,16,24,32,48,64,80};

for(int &N : length)
    {        
auto sites = SpinHalf(N);
auto state = InitState(sites);
for(int i = 1; i <= N; ++i)
    {
    //Neel state
    state.set(i,i%2==1 ? "Up" : "Dn");
    }
auto psi0 = MPS(state);
auto psi =  MPS(state);
Complex tstep=-Cplx_i*0.02;
Complex ttotal=-Cplx_i*N/4;
Real cutoff=1E-10;


auto gates = vector<Gate>();
auto gates0 = vector<Gate>();

auto ampo = AutoMPO(sites);
//Make the Heisenberg Hamiltonian
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 g = Gate(sites,b,b+1,Gate::tCReal,tstep/2,ampo);
    gates.push_back(g);

    auto g0 = Gate(sites,b,b+1,Gate::tCReal,-tstep/2,ampo);
    gates0.push_back(g0);

    }

for(int b = N-1; b >= 1; --b)
    {
    ampo += 0.5,"S+",b,"S-",b+1;
    ampo += 0.5,"S-",b,"S+",b+1;
    ampo +=     "Sz",b,"Sz",b+1;

    auto g = Gate(sites,b,b+1,Gate::tReal,tstep/2,ampo);
    gates.push_back(g);

    auto g0 = Gate(sites,b,b+1,Gate::Real,-tstep/2,ampo);
    gates0.push_back(g0);

    }

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);

//
// Begin the DMRG calculation
//
auto energy = dmrg(psi,H,sweeps,"Quiet");
//
auto energy0=  dmrg(psi0,H,sweeps,"Quiet");

auto Sz_1=sites.op("Sz",1);


// Time evolve

gateTEvol(gates,ttotal,tstep,psi,{"Cutoff=",cutoff,"Verbose=",true});

auto newpsi= Sz_1*psi.A(1);
newpsi.noprime();
psi.setA(1,newpsi);


// Right part exp(-H\tau)*S^{z}|g.s>

auto newpsi0 = Sz_1*psi0.A(1);
newpsi0.noprime();

psi0.setA(1,newpsi0);

gateTEvol(gates0,ttotal,tstep,psi0,{"Cutoff=",cutoff,"Verbose=",true});


//Print(psi0);
auto result=overlapC(psi,psi0).real();
//Print(result);
outfile <<N<<"  "<<abs(result)<< endl;
}
outfile.close();
return 0;
}

---

Hi, so yes, you can study any one-dimensional chain with Trotter gates if you also use "swap gates" to move the sites around so that they temporarily become nearest-neighbor. If you want to see a somewhat detailed discussion of this, you can see this paper: https://arxiv.org/abs/1002.1305

Thanks for posting your code, and please let me know if you have additional questions.

---

Hi Miles,
       Thanks for your reply. my posted naive code can work but results is not correct, is there any problem?? please help me, Thanks!

---

Hi, so I see a number of possible issues with your code. But it's a little hard to be of help right away, because you haven't stated very clearly what you are trying to calculate.

The main issue I can see is that the gates you are creating to do the time evolution are of type Gate::tReal whereas you said you are interested in doing imaginary time evolution. So would you not wish to use imaginary time evolution gates? To do that, set the type to Gate::tImag. For more information on BondGate objects please see this documentation page:
http://itensor.org/docs.cgi?vers=cppv3&page=classes/bondgate

One other thing I noticed is that you are computing the ground state MPS twice: once for psi and once for psi0. But these should be the same ground state, so you could just make a copy:
auto psi0 = psi;
and that way you don't need to run DMRG twice.

Finally, you should run your code on a very small system and compare in detail to an exact-diagonalization code (which doesn't have to be sophisticated for a small enough system) to debug it. You could find a lot of bugs just testing on a 4-site system, which is what I would do even with my experience. Be sure to print out all of your intermediate objects and as much information about them as you can, such as the energy of your wavefunctions, and their norms, at each step.

Hope that helps -

Miles

---

Hi Miles,
       your answer is great. I have last question(sorry about that). my ITensor is v2, if I use MPO imaginary time evolution (see http://itensor.org/docs.cgi?vers=cppv2&page=formulas/tevol_mps_mpo)
auto expH = toExpH<ITensor>(ampo,tau); (this term mean imaginary time evolution operator exp(-H\tau)) right, but I want to calculate imaginary time correlation function <g.s| exp(H\atu)S^{z}exp(-H\tau)S^{z} |g.s>. so I want to build exp(H\tau). I guess it is :
auto expH1=toExpH<ITensor>(amp,-tau); right?  Thank you very much !!!

---

Hi, yes if your question is whether inputting a tau of the opposite sign means the exponent will also have the flipped sign then the answer is yes. The code does not take the absolute value of your input. It does put a negative sign on whatever input tau you give it though, consistent with what you wrote.

Also please do test your code thoroughly against some exact results for small systems. As you know this is a good idea not only when using ITensor but when setting up a new numerical code more generally.

Best,
Miles

Answer 1

(Marking this question as answered: please see the discussion above.)