# Index Contraction Problem (Urgent!!)

**URL:** https://itensor.discourse.group/t/index-contraction-problem-urgent/220
**Category:** DMRG and Numerical Methods
**Tags:** julia, dmrg, mpo
**Created:** [June 20, 2022, 11:17am UTC](https://itensor.discourse.group/t/index-contraction-problem-urgent/220 "2022-06-20T11:17:37Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![Kartikeya](https://avatars.discourse-cdn.com/v4/letter/k/e9c0ed/32.png) [@Kartikeya](https://itensor.discourse.group/u/Kartikeya)
#### Post date: [June 20, 2022, 11:17am UTC](https://itensor.discourse.group/t/index-contraction-problem-urgent/220/1 "2022-06-20T11:17:37Z")

</div>

I want to code the following operator (1) in Julia to calculate correlations (2):  
 ![4](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/391d9dcde2b1183c92fcec48f8d1fb5df6414734.png)

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

I wrote a function for this as:

```auto
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

```auto
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

---

<div class="post-metadata">

### Author: ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)
#### Post date: [June 20, 2022, 8:24pm UTC](https://itensor.discourse.group/t/index-contraction-problem-urgent/220/2 "2022-06-20T20:24:06Z")

</div>

Unfortunately this won’t be possible to do very efficiently with the existing tools that we offer. We have a function `correlation_matrix` that computes correlations involving operators on pairs of sites but here you have operators spread out across four different sites in general.

A very convenient, but possibly slow, option is to use `OpSum` to make each operator that you want to measure, turn it into an MPO, then compute the expectation value of that MPO. Please let me know if you need help about how to use `OpSum` to do this but it is pretty similar to how you use it to make a Hamiltonian.

The other option, and most ideal way (computational efficiency wise) to do it is to write a custom loop that computes the correlation function by contracting MPS tensors along the MPS path and sticking in the operators on the four specific sites that you need. Here what one wants to do additionally is to save tensors that can be reused.

Two resources that could help you understand what I mean about this second strategy are:

1. code for the correlation\_matrix function, which would need to be generalized to four-site correlation functions for your case:  
[https://github.com/ITensor/ITensors.jl/blob/6d4e77cdd7694ef525584dfbc5466e8e85bafacb/src/mps/mps.jl#L602](https://github.com/ITensor/ITensors.jl/blob/6d4e77cdd7694ef525584dfbc5466e8e85bafacb/src/mps/mps.jl#L602)
2. an older tutorial (based on the C++ version of ITensor) that is about the strategy for computing correlation functions efficiently. The code won’t be directly applicable, but the general code pattern discussed there could be useful to you.  
[http://itensor.org/docs.cgi?vers=cppv3&page=tutorials/correlations](http://itensor.org/docs.cgi?vers=cppv3&page=tutorials/correlations)

Finally, the strategy you were trying of using `expect` won’t work because `expect` only accepts an operator that acts on a single site then computes its expectation value on each site of the MPS.

---

<div class="post-metadata">

### Author: ![Kartikeya](https://avatars.discourse-cdn.com/v4/letter/k/e9c0ed/32.png) [@Kartikeya](https://itensor.discourse.group/u/Kartikeya)
#### Post date: [June 21, 2022, 1:18am UTC](https://itensor.discourse.group/t/index-contraction-problem-urgent/220/3 "2022-06-21T01:18:59Z")

</div>

Hi Miles,

Thanks a lot for the quick reply. Could you please explain the `OpSum` approach a bit?  
It’d be helpful if you could explain with some rough code that I could expand on (because to calculate the required expectation value (1), four operatos will be needed to be multiplied together and I don’t know how to do that).

Best  
Kartikeya

---

<div class="post-metadata">

### Author: ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)
#### Post date: [June 22, 2022, 1:58pm UTC](https://itensor.discourse.group/t/index-contraction-problem-urgent/220/4 "2022-06-22T13:58:52Z")

</div>

Sure thing. Say you want to make the operator O = \frac{1}{\sqrt{2}} (c\_{i \uparrow} c\_{j \downarrow} - c\_{i \downarrow} c\_{j \uparrow}) \times (c\_{k \uparrow} c\_{l \downarrow} - c\_{k \downarrow} c\_{l \uparrow}) . Expanding the terms out, we have the terms:

O = \frac{1}{\sqrt{2}} c\_{i \uparrow} c\_{j \downarrow} c\_{k \uparrow} c\_{l \downarrow} - \frac{1}{\sqrt{2}} c\_{i \uparrow} c\_{j \downarrow} c\_{k \downarrow} c\_{l \uparrow} - \frac{1}{\sqrt{2}} c\_{i \downarrow} c\_{j \uparrow} c\_{k \uparrow} c\_{l \downarrow} + \frac{1}{\sqrt{2}} c\_{i \downarrow} c\_{j \uparrow} c\_{k \downarrow} c\_{l \uparrow} 

Then you can make O like this:

```auto
os = OpSum()
os += (1/sqrt(2)),"Cup",i,"Cdn",j,"Cup",k,"Cdn",l
os += (-1/sqrt(2)),"Cup",i,"Cdn",j,"Cdn",k,"Cup",l
os += (-1/sqrt(2)),"Cdn",i,"Cup",j,"Cup",k,"Cdn",l
os += (1/sqrt(2)),"Cdn",i,"Cup",j,"Cdn",k,"Cup",l
O = MPO(os,sites)

```

where `sites` is your array of site indices.

Then if you have an MPS wavefunction `psi`, you can compute

```auto
inner(psi,O,psi)

```

to get \langle \psi | O | \psi \rangle.

Please check my work above! I might have messed up the algebra somewhere like missing a minus sign.

---

<div class="post-metadata">

### Author: ![Kartikeya](https://avatars.discourse-cdn.com/v4/letter/k/e9c0ed/32.png) [@Kartikeya](https://itensor.discourse.group/u/Kartikeya)
#### Post date: [June 22, 2022, 4:17pm UTC](https://itensor.discourse.group/t/index-contraction-problem-urgent/220/5 "2022-06-22T16:17:37Z")

</div>

Thanks a lot Miles. You just saved my day!! 🙂