I want to code the following operator (1) in Julia to calculate correlations (2):

where, \alpha =a,b,c represent different nearest neighboring bonds as illustrated.

I wrote a function for this as:

```
function Delta_beta(r,beta)
x=0
if beta=='a'
x=Ly
end
if beta=='b'
x=1
end
if beta=='c'
x=1-Ly
end
op1=op("Cup",sites,r)
op2=op("Cdn",sites,r+x)
op3=op("Cdn",sites,r)
op4=op("Cup",sites,r+x)
return ((op1*op2)-(op3*op4))/sqrt(2)
end
#Defining Pairing Operators
function Delta_alpha_dagger(r,alpha)
x=0
if alpha=='a'
x=Ly
end
if alpha=='b'
x=1
end
if alpha=='c'
x=1-Ly
end
op1=op("Cdagup",sites,r)
op2=op("Cdagdn",sites,r+x)
op3=op("Cdagdn",sites,r)
op4=op("Cdagup",sites,r+x)
return ((op2*op1)-(op4*op3))/sqrt(2)
end
function P_alpha_beta(r0,r,alpha,beta)
return Delta_alpha_dagger(r0,alpha)*Delta_beta(r0+r,beta)
end
```

But when I run the code: `P_bb=[expect(psi,P_alpha_beta(1,r,'b','b')) for r in 1:Ly:L]`

I get

```
Attempting to contract IndexSet:
((dim=3|id=141|"Site,n=2,tJ")' <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1, (dim=3|id=141|"Site,n=2,tJ") <In>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1, (dim=3|id=523|"Site,n=1,tJ")' <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1, (dim=3|id=523|"Site,n=1,tJ") <In>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1)
with IndexSet:
((dim=3|id=141|"Site,n=2,tJ")' <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1, (dim=3|id=141|"Site,n=2,tJ") <In>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1, (dim=3|id=582|"Site,n=3,tJ")' <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1, (dim=3|id=582|"Site,n=3,tJ") <In>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1)
QN indices must have opposite direction to contract, but indices:
(dim=3|id=141|"Site,n=2,tJ")' <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1
and:
(dim=3|id=141|"Site,n=2,tJ")' <Out>
1: QN(("Nf",0,-1),("Sz",0)) => 1
2: QN(("Nf",1,-1),("Sz",1)) => 1
3: QN(("Nf",1,-1),("Sz",-1)) => 1
do not have opposite directions.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] compute_contraction_labels(Ais::NTuple{4, Index{Vector{Pair{QN, Int64}}}}, Bis::NTuple{4, Index{Vector{Pair{QN, Int64}}}})
@ ITensors C:\Users\HP\.julia\packages\ITensors\5dcHw\src\indexset.jl:648
[3] _contract(A::NDTensors.BlockSparseTensor{Float64, 4, NTuple{4, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 4}}, B::NDTensors.BlockSparseTensor{Float64, 4, NTuple{4, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 4}})
@ ITensors C:\Users\HP\.julia\packages\ITensors\5dcHw\src\itensor.jl:1853
[4] _contract(A::ITensor, B::ITensor)
@ ITensors C:\Users\HP\.julia\packages\ITensors\5dcHw\src\itensor.jl:1860
[5] contract(A::ITensor, B::ITensor)
@ ITensors C:\Users\HP\.julia\packages\ITensors\5dcHw\src\itensor.jl:1964
[6] *(A::ITensor, B::ITensor)
@ ITensors C:\Users\HP\.julia\packages\ITensors\5dcHw\src\itensor.jl:1951
[7] P_alpha_beta(r0::Int64, r::Int64, alpha::Char, beta::Char)
@ Main .\In[40]:43
[8] (::var"#13#14")(r::Int64)
@ Main .\none:0
[9] iterate
@ .\generator.jl:47 [inlined]
[10] collect(itr::Base.Generator{StepRange{Int64, Int64}, var"#13#14"})
@ Base .\array.jl:724
[11] top-level scope
@ In[41]:1
[12] eval
@ .\boot.jl:373 [inlined]
[13] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1196
```

Please help me out. It’s Urgent!!

P.S. I’ve used `sites = siteinds("tJ",L,conserve_qns=true)`

Best

Kartikeya