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.(...)
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))
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)))
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!