Define an initial not product MPS with SiteType

[Giuseppe]

Hi I am performing a time evolution for a 13 level qudit model where I define the following sitetype:

ITensors.space(::SiteType"SU2") = 13
sites = siteinds(“SU2”,N)

I want to use as initial state a MPS with bond dimension 2 where I know each element.
The problem is that sitetype seems to work only for the definition of product states not for MPS with higher bond dimension, is that true?
I have tried to define the MPS by end setting on each element the right tensor.
e.g.
psi=MPS(sites;linkdims)
psi[n]=T
where T is an ITensor(i,j,k)

The problem is that such state gives error when used as initial state for the time evolution defined in the sitetype environment.
For instance while a random mps (for three sites) with bond dimension 2 reads:

[1] ((dim=13|id=913|“SU2,Site,n=1”), (dim=2|id=69|“Link,l=1”))
[2] ((dim=2|id=69|“Link,l=1”), (dim=13|id=67|“SU2,Site,n=2”), (dim=2|id=518|“Link,l=2”))
[3] ((dim=2|id=518|“Link,l=2”), (dim=13|id=459|“SU2,Site,n=3”))

The MPS constructed by scratch reads:
1] ((dim=13|id=321|“index_j”), (dim=2|id=406|“index_i”))
[2] ((dim=2|id=406|“index_i”), (dim=13|id=321|“index_j”), (dim=2|id=861|“index_k”))
[3] ((dim=2|id=861|“index_k”), (dim=13|id=321|“index_j”))

There is some way to convert the second to the same sitetype style?
Many thanks in advance!

[kmp5]

Hi Giuseppe,

I am sorry you’re running into an issue. Unfortunately, I think I am confused. Can you please write out a full code block to run. What is linkdims and what is i, j, k? How are you constructing an MPS from scratch?
It would be helpful if you could please use the preformatted text option in the editor to make your code more readable.
Thanks!
Karl

[Giuseppe]

Hi Karl thank you for your time.
Here the part of the code under discussion:

#start--------#
N=3 #number of sites
linkdims=2#bond dimension

#define site type
ITensors.space(::SiteType"SU2") = 13
sites = siteinds(“SU2”,N)

#tensor definition

i = Index(2,“index_i”)
j = Index(13,“index_j”)
k = Index(2,“index_k”)

T1 = ITensor(j,i)#left edge
T1[i=>1,j=>5] = 1
T1[i=>2,j=>12] = 1

T2 = ITensor(i,j,k)#bulk
T2[i=>1,j=>5,k=>1] = 1
T2[i=>1,j=>12,k=>2] = 1
T2[i=>2,j=>13,k=>1] = 1
T2[i=>2,j=>6,k=>2] = 1

T3=ITensor(k,j)#right edge
T3[j=>1,k=>1] = 1
T3[j=>8,k=>2] = 1

#define state
psi=MPS(sites;linkdims)
psi[1]=T1
psi[2]=T2
psi[3]=T3

#end--------#

If I print psi this is the output I get
1] ((dim=13|id=321|“index_j”), (dim=2|id=406|“index_i”))
[2] ((dim=2|id=406|“index_i”), (dim=13|id=321|“index_j”), (dim=2|id=861|“index_k”))
[3] ((dim=2|id=861|“index_k”), (dim=13|id=321|“index_j”))

The problem is that I get error if I do operation or time evolution on this state with operators defined with sitetype.
The problem does not occur If I just define for instance a random mps with bond dimension 2
that is printed as

[1] ((dim=13|id=913|“SU2,Site,n=1”), (dim=2|id=69|“Link,l=1”))
[2] ((dim=2|id=69|“Link,l=1”), (dim=13|id=67|“SU2,Site,n=2”), (dim=2|id=518|“Link,l=2”))
[3] ((dim=2|id=518|“Link,l=2”), (dim=13|id=459|“SU2,Site,n=3”))

or if I start with whatever initial product state.

The question is how I can convert the first state to something that can be worked out with site type define operators. Or if there is a way to define not product states MPS with site type operators. A solution I guess is building it applying gates from an initial product state but once I know how to write it down the state I wonder if there is a more straightforward approach.
Thanks again,
Giuseppe

[miles]

Hi Giuseppe,
I’m not totally sure what your question is. Maybe I can partially answer it, though, by shedding some light on what the sitetype feature in ITensor does.

Mainly all that the sitetype system does is define mappings from indices to operators. So if a tensor or MPS has a certain Index with a special tag, such as “SU2” in your case, then our sitetype system lets ITensor know the answer to questions such as “what is the definition of the S^z operator for this Index?” because recall of course than an Index mathematically really represents a vector space, and there can be linear operators mapping that vector space to itself.

To do time evolution by applying gates on an MPS, technically you don’t need to use the sitetype system at all. You could instead just make your gates as ITensors, and apply them to the appropriate sites using tensor contraction and other helpful operations like using “prime levels” on indices to make it easy to make maps from a space back into itself.

In fact, I’m pretty sure that our apply function does not use the sitetype system at all, though it may look for indices to have the “Site” tag (I’d have to check though… it may even be agnostic to that).

So finally, when you say things still aren’t working with the way you defined an MPS in your example code, could you give a very specific example of precisely one thing you want to do that isn’t working? Then we could help you debug that one and then from that example maybe it would clarify your other questions.

[Giuseppe]

Hi Miles thanks a lot for your time. Here a simplified version of my code where I changed the local basis from 13 to 3 and defined all the tensors with zeros entries.

#----start----
using LinearAlgebra
using ITensors

let
cutoff = 1E-8

N=3 #number of sites
linkdims=2#bond dimension

#tensor definition
i = Index(2,“index_i”)
j = Index(3,“index_j”)
k = Index(2,“index_k”)

T1 = ITensor(j,i)#left edge
T2 = ITensor(i,j,k)
T3=ITensor(k,j)#right edge

ITensors.space(::SiteType"3level") = 3
sites = siteinds(“3level”,N)

U=complex(zeros(3,3))
ITensors.op(::OpName"U",::SiteType"3level") = U

psi=MPS(sites;linkdims)

psi[1]=T1
psi[2]=T2
psi[3]=T3

#two sites operator
hj = op(“U”, sites[1])*op(“U”, sites[2])

psie = apply(hj, psi; cutoff, maxdim=10)

end
#----end

This code gives the following error:
ERROR: BoundsError: attempt to access 0-element Vector{Int64} at index [1]
Stacktrace:
[1] getindex
@ ./essentials.jl:13 [inlined]
[2] product(o::ITensor, ψ::MPS, ns::Vector{Int64}; move_sites_back::Bool, apply_dag::Bool, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:cutoff, :maxdim), Tuple{Float64, Int64}}})
@ ITensors ~/.julia/packages/ITensors/WMeVS/src/mps/abstractmps.jl:2108
[3] product
@ ~/.julia/packages/ITensors/WMeVS/src/mps/abstractmps.jl:2093 [inlined]
[4] top-level scope
@ ~/Documents/Julia_codes/Su2_1D_julia/domanda_itensor.jl:42

The same does not happen if I apply the operator hj on the state psi=randomMPS(sites,2) that has the same local and bond dimension.

My question was if I could still define gates and operator using sitetype with the MPS defined as ITensor. If I understood correctly from your reply seems that the answer is just to do not use sitetype and define gates and operator as ITensors. Is that right?
Thanks again for the reply and for your great work with the library!

[kmp5]

Giuseppe,

One issue I see with your most recent code is that the default constructor to ITensor creates a tensor with an EmptyStorage type.

julia> T1 = ITensor(j,i)
ITensor ord=2 (dim=3|id=706|"index_j") (dim=2|id=321|"index_i")
EmptyStorage{NDTensors.EmptyNumber, Dense{NDTensors.EmptyNumber, Vector{NDTensors.EmptyNumber}}}

You need to have data in your tensors and we suggest adding an element type to your constructor

julia> T1 = ITensor(Float64, j,i)
julia> fill!(T1, zero(Float64))

Another issue is that your hj tensor has no matching indices with your wavefunction when you change psi to [T1, T2, T3].

[Giuseppe]

Hi Karl,
thanks for your message. It seems this was the issue:

“Another issue is that your hj tensor has no matching indices with your wavefunction when you change psi to [T1, T2, T3]”

I was misinterpreting the index labelling, thanks to point it out.
By the way, in parallel I wrote down the state just applying gates on an initial product state and it turned out to be the most convenient choice. Maybe it can be useful to someone.
Many thanks again,
Giuseppe