From exact diagonalization to MPS

[sebastianvromero]

Hi, hope you are doing well.

Given an op in the form of OpSum, I would like to obtain its whole set of eigenstates and turn them into MPSs to compute separately the overlaps between each eigenstate of op with a given MPS psi, thus \langle\text{each eigst of op}|\psi\rangle. I have been following different posts regarding exact diagonalization with ITensor

• Converting the eigenstates obtained from the exact diagonalization to a MPS
• ED with Itensor

but, unfortunately, they apparently do not answer my question.

Here there is a snippet of code where the op to diagonalize is prepared and exact diagonalized using eigen, returning an 2^N\times 2^N matrix V for the eigenstates, as expected.

using ITensors

let
  N = 4 # number of sites

  s = siteinds("Qubit", N)
  op = OpSum()
  for i in 1:N
    op += "σy", i
  end

  cb=combiner(dag(s)...)
  cbp=combiner(s'...)
  H = prod(MPO(op, s))*cb*cbp
  E, V = eigen(H,combinedind(cb),combinedind(cbp), ishermitian=true)

  psi = randomMPS(s;linkdims=2)
  return
end

The problem comes when I try to access each individual eigenstate of V, since I do not know how to manipulate the NDTensors.Dense object. As I said, I would also like to convert them into MPSs to extract each inner(each eigst of op', psi).

Thanks a lot for your time devoted!

[ryanlevy]

First note that we have separated ITensor and ITensorMPS so you’ll need using ITensorMPS as well

Unfortunately that snippet isn’t the best way of doing things (and we should really only have one combiner), but with some modification here is how you can do the overlaps

using ITensors, ITensorMPS

let
  N = 4 # number of sites

  s = siteinds("Qubit", N)
  os = OpSum()
  for i in 1:N
    os += "σy", i
  end

  cb=combiner(dag(s)...)
  H = prod(MPO(os, s))*cb*dag(cb)'
  E, V = eigen(H,ishermitian=true)

  psi = randomMPS(s;linkdims=2)

  psic = prod(conj(psi))*cb # conj psi in case complex
  l = noncommonind(V,psic)
  @assert hastags(l,"eigen")

  # overlap of single vector
  i = 7
  vec = V*onehot(l=>i)
  overlap = (psic*vec)[]
  @show overlap

  # all at once
  overlaps = V*psic
  overlaps = [(overlaps*onehot(l=>i))[] for i=1:2^N]
  @show overlaps
  @show overlap ≈ overlaps[i]
  return
end

[sebastianvromero]

It works! Thanks a lot