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