about orthogonalize! function

ITensor Julia Questions
1

Dear there,
I have a quick question about calculating expectation of a two-site operator like simply c^{\dag}_jc_{j+1} with a MPS got from DMRG. I tried two ways as follows, and my questions is I find they give totally different results:

# Method 1 -- use MPO
value=0
for j in 1:1:(N-1)
    hop = OpSum()
    hop += 1.0, "Cdagup", j, "Cup", j + 1
    hop_mpo = MPO(hop,sites) 
    value = inner(psi0',hop_mpo,psi0)
    println("j=$j, C+upCup=$value")
end
       
# Method 2 -- use orthogonalize  
value2 = 0
for j in 1:1:(N-1)
	orthogonalize!(psi0, j)
	psi0_j = psi0[j]*psi0[j+1]
	psidag0_j = dag(prime(psi0[j]*psi0[j+1], "Site"))
	ampo = op(sites, "Cdagup", j) * op(sites, "Cup", j + 1)
	value2 = scalar(psidag0_j * ampo * psi0_j) 
	println("j=$j, C+upCup=$value2")
end

So in the 1st way, I made c^{\dag}_jc_{j+1} as a MPO, and in the 2nd way I orthogonalize at site j and use op. But I’m worried about if the 2nd way is valid because for each \langle c^{\dag}_jc_{j+1}\rangle I use both site j and j+1, but I only orthogonalize j not j+1.

So I’m just wondering could you help to see what’s wrong in 1 or 2 so that the results are different, and if the 2nd way is indeed wrong due to my orthoginalize, then should I orthogonalize two sites simultaneously or maybe this way can only be used for single-site operator expectation?

Thank you so much!!

2

Hi zz, The MPO method should give you correct answer, but it is very inefficient. I would also recommend that you add a third method and that is using the correlation_matrix function. It does what you are trying to do in method 2. It orthogonalizes on site j and the calculates what is effectively a left envorinment for the operator at site j, and sweeps that environment through to site j+n.
You can view the code here to see how it works. This should be the most efficient way to calculate multi site expectations.
Let me know if you have more questions.
I hope this helps.
JR

1 Like
3

To answer the original question as to why method 2 is wrong, when applying ops to the MPS like that you have to keep track of the fermion strings yourself, which the helper function Jan posted does. To fix your method 2 code, just change "Cdagup" → "Cdagup * F".

As for the orthogonalize, this works because your operator is acting on sites (j,j+1) and the rest of the MPS contracts into an identity. It can be helpful to draw this contraction out in order to see why you only need either j or j+1 called in orthogonalize!

2 Likes