Partial Contraction

I was looking for a partial contraction of MPS, similar to what is described here:

Essentially I’m looking for something like this (more simply we can assume \psi = \phi)

I’m sure this is some very simple oversight, but when I use the following code

function partial_contract(ψ::MPS, leftend::Int, rightend::Int)
  if leftend > rightend 
    error("leftend cannot be greater than rightend")
  end 

  # if leftend >= rightend - 1
  #   return inner(ψ', ψ)
  # end 

  sys = deepcopy(ψ[ leftend + 1: rightend - 1])

  L = R = ITensor(1.)

  for left = 1:leftend

    L *= ψ[left]
    L *= prime(dag(ψ[left]))

  end 

  for right =rightend:length(ψ)

    R *= prime(dag(ψ[right]))
    R *= ψ[right]
  end

  sys[1] *= L
  sys[end] *= R

  .....
  ......

  return sys
end 

For some states that I calculated, I got the following dimensions for the product (for every site step)
size(L) = (2, 2, 2, 2)
size(L) = (2, 2, 4, 2, 4, 2)
size(L) = (2, 2, 2, 2, 8, 2, 8, 2)
size(L) = (2, 2, 2, 2, 2, 2, 16, 2, 16, 2)
size(L) = (2, 2, 2, 2, 2, 2, 2, 2, 23, 2, 23, 2)
So obviously they are not contracting properly. I’m wondering if there’s something I missed.

I think the problem is that you are priming both the link and site indices when you write things like prime(dag(ψ[right])), so your contraction is probably more like an outer product of the state with itself rather than a partial inner product. Try just priming the link indices. A shorthand for priming the link indices of the entire MPS is:

julia> using ITensors

julia> s = siteinds("S=1/2", 4)
4-element Vector{Index{Int64}}:
 (dim=2|id=957|"S=1/2,Site,n=1")
 (dim=2|id=666|"S=1/2,Site,n=2")
 (dim=2|id=76|"S=1/2,Site,n=3")
 (dim=2|id=930|"S=1/2,Site,n=4")

julia> ψ = randomMPS(s)
MPS
[1] ((dim=2|id=957|"S=1/2,Site,n=1"), (dim=1|id=522|"Link,l=1"))
[2] ((dim=1|id=522|"Link,l=1"), (dim=2|id=666|"S=1/2,Site,n=2"), (dim=1|id=552|"Link,l=2"))
[3] ((dim=1|id=552|"Link,l=2"), (dim=2|id=76|"S=1/2,Site,n=3"), (dim=1|id=475|"Link,l=3"))
[4] ((dim=1|id=475|"Link,l=3"), (dim=2|id=930|"S=1/2,Site,n=4"))

julia> ψdag = prime(linkinds, dag(ψ))
MPS
[1] ((dim=2|id=957|"S=1/2,Site,n=1"), (dim=1|id=522|"Link,l=1")')
[2] ((dim=1|id=522|"Link,l=1")', (dim=2|id=666|"S=1/2,Site,n=2"), (dim=1|id=552|"Link,l=2")')
[3] ((dim=1|id=552|"Link,l=2")', (dim=2|id=76|"S=1/2,Site,n=3"), (dim=1|id=475|"Link,l=3")')
[4] ((dim=1|id=475|"Link,l=3")', (dim=2|id=930|"S=1/2,Site,n=4"))