Index Contraction Problem (Urgent!!)

DMRG and Numerical Methods
, ,
1

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

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

2

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

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.

3

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

4

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:

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

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.

1 Like
5

Thanks a lot Miles. You just saved my day!! :slight_smile:

1 Like