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

0 votes
asked

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

## 1 Answer

0 votes
answered by (14.1k points)
selected by

Best answer

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

commented by (14.1k points)
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.
commented by (400 points)
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
commented by (70.1k points)
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.