Contracting ITenor with MPO

ITensor Julia Questions
1

Hi all,

I am trying to do the contraction

Tr(H^\dagger c_m^\dagger H c_n)

where H is an MPO and C and C dagger are fermionic annihilation / creation operators. My problem is I have to the contraction in this order. I cannot make an MPO from single C or C dagger operators because that would not preserve particle number and I am using symmetry preserving MPOs.

Whenever I tr to diligently do the above contraction using the apply function, I have run into a myriad of different issues. For instnace, I am pretty sure the contraction H c_n requires c_n to act on the right side or unprimed side of H. However, whenever I try to do this I run into so many issues with indices and primes. I have tried many things but it just seems to become a complicated mess. Any help is appreciated. Here is my current code.

            O1 = op("C",siteinds(H,i),1)
            O2 = swapprime(op("Cdag",siteinds(H,j),2), 0 => 2)

            println(O2)
            println(siteinds(H,j))
            tmp1 = apply(O2,H)
            tmp1 = swapprime(tmp1, 2=>0)
            tmp2 = apply(O1,tmp1)
            value = tr(apply(Hdag, tmp2))

Right now I get the error: ERROR: In product, must have common indices with prime level 0.
This is frustrating because I dont know how else to multiply the operator on the right of my MPO, H.
It seems most of my trouble is figuring out how to right multiply an MPO

2

Hi Ben, try this

using ITensors
using ITensorMPS

function create_or_anihilate(op_name::String, i::Int, sites::Vector{<:Index})::MPO
  H = MPO(sites)
  
  for j in 1:(i - 1)
    H[j] = op("F", sites[j])
  end

  H[i] = op(op_name, sites[i])

  for j in (i + 1):length(sites)
    H[j] = op("I", sites[j])
  end

  # This puts in the dummy link indices.
  orthogonalize!(H, 1)

  return H
end

function my_trace(H::MPO, m::Int, n::Int)::Number
  sites = noprime(siteinds(first, H))
  Cdag_m = create_or_anihilate("Cdag", m, sites)
  C_n = create_or_anihilate("C", n, sites)

  HC_n = apply(H, C_n; alg="naive", truncate=false)
  Cdag_mHC_n = apply(Cdag_m, HC_n; alg="naive", truncate=false)

  return inner(H, Cdag_mHC_n)
end

let 
  sites = siteinds("Fermion", 10; conserve_qns=true)

  H = MPO(sites, "I")
  
  @show my_trace(H, 4, 8)
  @show my_trace(H, 4, 4)
end

nothing

Output:

my_trace(H, 4, 8) = 0.0
my_trace(H, 4, 4) = 512.0000000000011

Which seems correct to me since \operatorname{tr}(I^{\otimes 9} \otimes n) = 2^9 = 512.

No promises that this is research ready, but I think the general idea is correct. Since the creation and annihilation operators are non-local due to the Jordan-Wigner string you need to represent them as bond-dimension one MPOs and not just local operators.

If H is a Hamiltonian, as a sanity check I would construct c^\dagger_m H c_n directly as an MPO via the OpSum machinery and compare the results with my_trace above. If H is not a Hamiltonian, I would cook up some example Hamiltonians to make sure the implementation is correct.

3

Brilliant! Yes this seems to work. I did not think of this since opsum seems to usually autoapply the JW strings. Thank you very much!

4

Great. It agree it can be confusing when the JW strings are handled for you and when you need to be aware of them.