[Julia] Regarding index contraction issue in correlation_matrix function for conserved QN case

Question

Hi,

Recently I am playing with the Julia version of ITensor. And I think the Julia version is just amazing! It is much easier to use once one gets used to the language. And there are a lot of implemented functions, which can save lots of time.

However, I recently meet with some problems when using the correlation_matrix() function which is implemented

I would like to calculate the correlation functions on a cylinder with PBC and conserved QNs, say

Nx = 30
Ny = 2
sites = siteinds("S=1", N; conserve_qns=true)
state = [isodd(n) ? "Up":"Dn" for n=1:N]
psi = randomMPS(state,sites,20)

And calculate the correlation functions starting from the middle

Czz = correlation_matrix(psi, "Sz", "Sz"; site_range=28:45)

But it gives the following errors.

    Attempting to contract IndexSet:
((dim=1|id=875|"Link,l=28") <Out>
 1: QN("Sz",0) => 1, (dim=1|id=875|"Link,l=28")' <Out>
 1: QN("Sz",0) => 1)with IndexSet:
((dim=3|id=838|"S=1,Site,n=29") <Out>
 1: QN("Sz",2) => 1
 2: QN("Sz",0) => 1
 3: QN("Sz",-2) => 1, (dim=1|id=875|"Link,l=28") <Out>
 1: QN("Sz",0) => 1, (dim=3|id=839|"Link,l=29") <Out>
 1: QN("Sz",2) => 1
 2: QN("Sz",0) => 1
 3: QN("Sz",-2) => 1)QN indices must have opposite direction to contract.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] compute_contraction_labels(::Tuple{Index{Array{Pair{QN,Int64},1}},Index{Array{Pair{QN,Int64},1}}}, ::Tuple{Index{Array{Pair{QN,Int64},1}},Index{Array{Pair{QN,Int64},1}},Index{Array{Pair{QN,Int64},1}}}) at /home/slowlight/.julia/packages/ITensors/OLHU9/src/indexset.jl:585
 [3] _contract(::ITensors.NDTensors.Tensor{Float64,2,ITensors.NDTensors.DiagBlockSparse{Float64,Float64,2},Tuple{Index{Array{Pair{QN,Int64},1}},Index{Array{Pair{QN,Int64},1}}}}, ::ITensors.NDTensors.Tensor{Float64,3,ITensors.NDTensors.BlockSparse{Float64,Array{Float64,1},3},Tuple{Index{Array{Pair{QN,Int64},1}},Index{Array{Pair{QN,Int64},1}},Index{Array{Pair{QN,Int64},1}}}}) at /home/slowlight/.julia/packages/ITensors/OLHU9/src/itensor.jl:1627
 [4] _contract(::ITensor, ::ITensor) at /home/slowlight/.julia/packages/ITensors/OLHU9/src/itensor.jl:1634
 [5] contract(::ITensor, ::ITensor) at /home/slowlight/.julia/packages/ITensors/OLHU9/src/itensor.jl:1732
 [6] *(::ITensor, ::ITensor) at /home/slowlight/.julia/packages/ITensors/OLHU9/src/itensor.jl:1720
 [7] correlation_matrix(::MPS, ::String, ::String; kwargs::Base.Iterators.Pairs{Symbol,UnitRange{Int64},Tuple{Symbol},NamedTuple{(:site_range,),Tuple{UnitRange{Int64}}}}) at /home/slowlight/.julia/packages/ITensors/OLHU9/src/mps/mps.jl:598
 [8] top-level scope at /home/slowlight/xuyi/julia_projects/eval_obs.jl:167
 [9] include(::Function, ::Module, ::String) at ./Base.jl:380
 [10] include(::Module, ::String) at ./Base.jl:368
 [11] exec_options(::Base.JLOptions) at ./client.jl:296
 [12] _start() at ./client.jl:506
in expression starting at /home/slowlight/xuyi/julia_projects/eval_obs.jl:6

So, I follow the Stacktrace and find that it is the following line in the correlation_matrix that has some problems. (line 598 in mps.jl)

598   Li = L * psi[i]

But the same function works well for site_range=1:N. This leads me to think that the problem has something to do with this if-else loop.

  if start_site == 1
    L = ITensor(1.0)
  else
    lind = commonind(psi[start_site], psi[start_site - 1])
    L = delta(lind, lind')
  end

And I did a test with the same state and siteset.

lind1 = commonind(psi[i_start-1], psi[i_start])
@show lind1
lind2 = commonind(psi[i_start], psi[i_start-1])
@show lind2

L1 = delta(lind1, lind1')
Li1 = L1 * psi[i_start]
println("1 works!")

L2 = delta(lind2, lind2')
Li2 =  psi[i_start] * L2
println("2 works!")

It turns out that the previous contraction is alright, while the second has the same problem. I guess it has something to do with the diretion of QN flux when doing contraction.

    lind1 = (dim=1|id=154|"Link,l=28") <In>
 1: QN("Sz",0) => 1
lind2 = (dim=1|id=154|"Link,l=28") <Out>
 1: QN("Sz",0) => 1
1 works!
ERROR: LoadError: Attempting to contract IndexSet:
((dim=3|id=610|"Link,l=29") <Out>
 1: QN("Sz",-2) => 1
 2: QN("Sz",0) => 1
 3: QN("Sz",2) => 1, (dim=3|id=838|"S=1,Site,n=29") <Out>
 1: QN("Sz",2) => 1
 2: QN("Sz",0) => 1
 3: QN("Sz",-2) => 1, (dim=1|id=154|"Link,l=28") <Out>
 1: QN("Sz",0) => 1)with IndexSet:
((dim=1|id=154|"Link,l=28") <Out>
 1: QN("Sz",0) => 1, (dim=1|id=154|"Link,l=28")' <Out>
 1: QN("Sz",0) => 1)QN indices must have opposite direction to contract.

Maybe I miss something.. Could anyone help me recognize the problem?

Best,
Yi

Answer 1

Hi Yi,

Thanks for looking into this, probably the line L = delta(lind, lind') should be changed to L = delta(dag(lind), lind').

Could you try that out and see if it helps?

A simple way to modify the file is to check out the ITensors.jl package in "development mode" by doing:

julia> using Pkg

julia> Pkg.dev("ITensors")

and then modify the file ~/.julia/dev/ITensors/src/mps/mps.jl. To have the changes you make reflected the next time you run the correlation_matrix function, either start a new Julia session and rerun your code or install the Revise package:

https://timholy.github.io/Revise.jl/stable/

which will reflect changes you make automatically without having to start a new Julia session.

If that code fixes it, it would also be great if you could make a pull request for the fix on Github (for example by making a fork that includes the change and making a PR using the online interface: https://github.com/ITensor/ITensors.jl/compare )!

-Matt

Comments

The same instructions for making local changes to the ITensors.jl library can be found here:

https://itensor.github.io/ITensors.jl/stable/AdvancedUsageGuide.html#Developing-ITensors.jl

by the way.

---

Hi Matt,

Thank you ver much for your help! Now it works perfectly! And I have pulled the request about this change.

By the way, I was wondering if there is a substitue for the nmultMPO function in Julia version of ITensor, which can multiply two MPOs that have some overlap sites.

Best,
Yi

---

Hi Yi, since this is a different question, could you post it as a new question? It will help us to be able to see what questions have been answered or not.