AutoMPO using operators with more arguments

Question

[Using Julia]

Hi! I am hoping to construct an ampo using operators that contain extra arguments. For example,

function ITensors.op(::OpName"Nup", st::SiteType"HubHolst"; max_phonons::Int=1)

Initially I just had a global MAX_PHONONS variable defined in the HubHolst site type file, but for the sake of making the code more hygienic I was hoping to explicitly pass in max_phonons as an argument to specify the dimension of the phonon Hilbert space whenever I construct the operators. Alternatively (and preferably), the st variable could also carry this information. However, I am not sure whether it is possible (that is to say, internally compatible) to modify the SiteType struct to contain more fields, such as the dimension of the phonon hilbert space on each site.

Thanks so much for the help! I really appreciate it

Answer 1

Hi thanks for the question. So if the main parameter you are wanting to pass is the dimension of one of the Hilbert space types (the bosons representing the phonons) then the best way to do that is the idea you suggest at the end of your question: just by making your sites in a different way.

The "Boson" site type included in ITensors.jl recognizes a dim keyword argument whose default is 2, but can be increased to any number.

So then instead of using the siteinds function to make your sites, you can do the following for your case

dim = 3
s = [isodd(n) ? siteind("Fermion", n) : siteind("Boson", n; dim) for n in 1:N]

which will make an array of site indices and pass the dim keyword argument in the case that the site is a "Boson". (There is a Julia trick here where instead of writing dim=dim you can just write dim.)

Would that address the issue you're facing? You can also print out the array s of site indices afterward to verify that your "Boson" sites have the correct dimension.

Finally, when you then make operators for these sites, such as by using the OpSum/AutoMPO system, the correct operators will be made to match the dimension of the "Boson" sites that you make.

Happy to discuss more -

Miles

Comments

Thank you so much! For anybody following the question, here is some example code for what ended up working for me. I am doing on-site phonon interactions, so the dimension of the Hilbert space is actually 4 * (number of phonon degrees of freedom). This information is captured in how I define the ITensors.space(...) function, and can be read out directly from the index I pass into ITensors.(...)

Working code:

function ITensors.op(::OpName"I", ::SiteType"HubHolst", s::Index)
  M = Matrix(I, Int(dim(s)), Int(dim(s)))
  return ITensor(M, prime(s), dag(s))
end

Small follow-up question, but by no means urgent: as you can see in the last line of the function above, I have to (for some unknown reason) specifically convert M into an ITensor with the appropriate indices before returning. However, this was not necessary when I had, say,

function ITensors.op(::OpName"I", ::SiteType"HubHolst")
  return Matrix(I, Int(dim(s)), Int(dim(s)))
end

Why is this the case? I assume it has something to do with the extra s::Index argument, but I am afraid I don't know enough about the inner mysteries of Julia to answer this question for myself.

Thanks again for all the help!

---

Glad it helped and thanks for posting the code!

Your follow up question is a good one. The answer first of all is that the difference in behavior there isn't a Julia language thing, it's just about how we designed the op function overloads in ITensors.jl.

What we did was first to define the versions of op that take an Index at the end, expecting or requring them to return an ITensor.

Later, we added the version of op that doesn't take an Index, and only takes an OpName and a SiteType and set things up so that version of op is expected to just return a matrix. So it was just a deliberate design choice because we realized it would be nice to have the matrix version to use, to allow things to be as simple as possible.

So you can think of the version taking Index as kind of an older interface that we are keeping for backwards compatibility, or as a more advanced interface for cases where you need to do special things to the ITensor before returning it (say, if you wanted to make it a kind of sparse ITensor or something like that).