Site index and labels while doing TEBD

Hi,

I a using a standard TEBD code from ITensor and I have a small question in relation to that.

Essentially the code I am using is given in MPS Time Evolution · ITensors.jl

Now, when I start with a simple product state (as in the example given in the previous link I pasted), and I print the initial mps , I get

MPS
[1] ((dim=2|id=624|"S=1/2,Site,n=1"), (dim=1|id=408|"Link,l=1"))
[2] ((dim=1|id=408|"Link,l=1"), (dim=2|id=106|"S=1/2,Site,n=2"), (dim=1|id=538|"Link,l=2"))
[3] ((dim=1|id=538|"Link,l=2"), (dim=2|id=808|"S=1/2,Site,n=3"), (dim=1|id=878|"Link,l=3"))
[4] ((dim=1|id=878|"Link,l=3"), (dim=2|id=339|"S=1/2,Site,n=4"), (dim=1|id=157|"Link,l=4"))
[5] ((dim=1|id=157|"Link,l=4"), (dim=2|id=366|"S=1/2,Site,n=5"), (dim=1|id=197|"Link,l=5"))
[6] ((dim=1|id=197|"Link,l=5"), (dim=2|id=289|"S=1/2,Site,n=6"), (dim=1|id=434|"Link,l=6"))
[7] ((dim=1|id=434|"Link,l=6"), (dim=2|id=167|"S=1/2,Site,n=7"), (dim=1|id=147|"Link,l=7"))
[8] ((dim=1|id=147|"Link,l=7"), (dim=2|id=386|"S=1/2,Site,n=8"), (dim=1|id=116|"Link,l=8"))
[9] ((dim=1|id=116|"Link,l=8"), (dim=2|id=134|"S=1/2,Site,n=9"), (dim=1|id=79|"Link,l=9"))
[10] ((dim=1|id=79|"Link,l=9"), (dim=2|id=801|"S=1/2,Site,n=10"))

which is all fine in terms of index structure that gives me link and site informations. Nonetheless after I have run the code (i.e., finite number of TEBD ‘iterations’), I get the mps structure as

[1] ((dim=2|id=624|"S=1/2,Site,n=1"), (dim=2|id=213|"Link,n=1"))
[2] ((dim=2|id=213|"Link,n=1"), (dim=2|id=106|"S=1/2,Site,n=2"), (dim=4|id=897|"Link,n=1"))
[3] ((dim=2|id=808|"S=1/2,Site,n=3"), (dim=8|id=926|"Link,n=1"), (dim=4|id=897|"Link,n=1"))
[4] ((dim=2|id=339|"S=1/2,Site,n=4"), (dim=16|id=13|"Link,n=1"), (dim=8|id=926|"Link,n=1"))
[5] ((dim=2|id=366|"S=1/2,Site,n=5"), (dim=32|id=701|"Link,n=1"), (dim=16|id=13|"Link,n=1"))
[6] ((dim=2|id=289|"S=1/2,Site,n=6"), (dim=16|id=99|"Link,n=1"), (dim=32|id=701|"Link,n=1"))
[7] ((dim=2|id=167|"S=1/2,Site,n=7"), (dim=8|id=531|"Link,n=1"), (dim=16|id=99|"Link,n=1"))
[8] ((dim=2|id=386|"S=1/2,Site,n=8"), (dim=4|id=292|"Link,n=1"), (dim=8|id=531|"Link,n=1"))
[9] ((dim=2|id=134|"S=1/2,Site,n=9"), (dim=2|id=47|"Link,n=1"), (dim=4|id=292|"Link,n=1"))
[10] ((dim=2|id=801|"S=1/2,Site,n=10"), (dim=2|id=47|"Link,n=1"))

as you can clearly see all the “Link,n=1”. I am not sure why it is now behaving like this. Would be nice if someone can point me out how to get rid of this.

The reason I am asking this question is as follows. I have a vectorized density matrix (which is an MPS) and then I want to convert it to MPO structure and compute the log negativity (see Figure 2 in https://journals.aps.org/prb/pdf/10.1103/PhysRevB.102.064304) and for that I need some contractions (as you can see in fig. 2 of the paper I just referred to) conducted correctly. Therefore, it is important for me that all these link and site index infos are structured in a correct and comprehensive way. In case, I am missing something please point that out.

Thanks in advance for help.

Best,
Sourav.

Hi Sourav,
Thanks for letting us know about this. It is indeed a bug in the apply function behavior (though probably due to a lower-level routine that apply calls down to). So we will check if there is an open Issue for it our on Github tracker, and if not start one.

Importantly, though, this is very likely just a “cosmetic” issue for the applications you are doing. Note that the Index ID numbers are unique and different for the links you see in the printout. So the indices are actually different. It is only the tags which are ending up the same. But unless you are using the tags explicitly to ‘grab’ or modify certain link indices, this should not affect your code’s correctness at all. Are you needing to do that? If so, there might be a better pattern we can recommend to you.

Hi Miles,

Thanks for your reply. Indeed, it is a cosmetic issue and for most of the things I have been doing I never needed these details. But in one of my projects, I am trying to compute log negativity for which I need to compute from my density matrix MPO something as shown in the attached figure (original reference is https://arxiv.org/pdf/2110.07384).

Now for this computation, I guess I need to do some very specific contraction and in this context I was asking this question. Although I have (frankly speaking) not thought in detail, I had the impression that for this computation I might need links and everything specified precisely–maybe I am missing something. In any case, is there any routine/example code in ITensor for computing such quantity from MPOs? I was trying to search some old chats in this discussion forum but did not find something very relevant/useful.

Thanks again!

Best,
Sourav

There is not a pre-coded routine to compute negativity. (It would be a nice thing to have as an “overlay package”, something like ITensorEntanglement.jl but we don’t have that at this time, so you must compute entanglement measures manually.)

I would encourage you to start coding it up and see what you end up needing to do in terms of manipulating indices. One pattern that we like a lot is to try to use “graph” information to obtain indices (when you need to obtain them at all), such as calling commonind(A,B) on a pair of ITensors A and B. If these have a single common Index between them, this function will return it. Also there is a function uniqueind(A,B) that returns the unique Index of A that does not connect to B. There are also commoninds and uniqueinds when you expect that there may be more than one common or unique index.

Hi Miles,
Thanks a lot. Let me try it, if there is some issue/problem I will possibly write here–

best,
sourav.

Great – I think it would be easier to answer any questions that may arise after you got some way along trying it. Glad you will give it a go yourself first.

Hi Miles,

I have managed to write a code for computing log negativity. I am not going into details but there is one problem I have got stuck. To simplify the scenario let me give me an example. Let us consider a simple TEBD as follows:

N, cf, tau, ttotal= 4, 1e-10, 0.1, 5.0
s = siteinds("S=1/2", N)
chi = 60

gates = ITensor[]
for j in 1:(N - 1)
    s1 = s[j]
    s2 = s[j + 1]
    hj =
    op("Sz", s1) * op("Sz", s2) +
    1 / 2 * op("S+", s1) * op("S-", s2) +
    1 / 2 * op("S-", s1) * op("S+", s2)
    Gj = exp(-im * tau / 2 * hj)
    push!(gates, Gj)
end

append!(gates, reverse(gates))
psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
c = div(N, 2) 
for t in 0.0:tau:ttotal
    t≈ttotal && break
    psi = apply(gates, psi; cutoff = cf, maxdim = chi )
    normalize!(psi)
end

Now, this spits out an MPS psi. In the next part of my code I use something like

ψ1 = deepcopy(MPO(psi))
ψ2 = deepcopy(MPO(psi))
ψ3 = deepcopy(MPO(psi))

contraction_left_side_1 = contract(prime(ψ2[i], "Site"), ψ1[i], cutoff=cf)

again it is just a part of a routine and I am not going to any details. Problem is as follows:

From the line " contraction_left_side_1 = contract(prime(ψ2[i], “Site”), ψ1[i], cutoff=cf)" if I remove the cutoff=cf, everything works smoothly. But when I put the cutoff I get the following error message

ERROR: StackOverflowError:
Stacktrace:
     [1] contract(::ITensor, ::Vararg{ITensor}; kwargs::@Kwargs{cutoff::Float64})
       @ ITensors ~/.julia/packages/ITensors/23eRw/src/tensor_operations/tensor_algebra.jl:168
     [2] (::ITensors.var"#213#215"{@Kwargs{cutoff::Float64}})(A::ITensor, B::ITensor)
       @ ITensors ~/.julia/packages/ITensors/23eRw/src/tensor_operations/tensor_algebra.jl:158
     [3] (::Base.BottomRF{ITensors.var"#213#215"{@Kwargs{cutoff::Float64}}})(acc::ITensor, x::ITensor)
       @ Base ./reduce.jl:86
     [4] afoldl(::Base.BottomRF{ITensors.var"#213#215"{@Kwargs{cutoff::Float64}}}, ::Base._InitialValue, ::ITensor, ::ITensor)
       @ Base ./operators.jl:545
     [5] _foldl_impl(op::Base.BottomRF{ITensors.var"#213#215"{…}}, init::Base._InitialValue, itr::Tuple{ITensor, ITensor})
       @ Base ./reduce.jl:68
     [6] foldl_impl(op::Base.BottomRF{ITensors.var"#213#215"{@Kwargs{…}}}, nt::Base._InitialValue, itr::Tuple{ITensor, ITensor})
       @ Base ./reduce.jl:48
     [7] mapfoldl_impl(f::typeof(identity), op::ITensors.var"#213#215"{…}, nt::Base._InitialValue, itr::Tuple{…})
       @ Base ./reduce.jl:44
     [8] mapfoldl(f::Function, op::Function, itr::Tuple{ITensor, ITensor}; init::Base._InitialValue)
       @ Base ./reduce.jl:175
     [9] mapfoldl
       @ ./reduce.jl:175 [inlined]
    [10] foldl(op::Function, itr::Tuple{ITensor, ITensor}; kw::@Kwargs{})
       @ Base ./reduce.jl:198
    [11] contract(As::Tuple{ITensor, ITensor}; sequence::String, kwargs::@Kwargs{cutoff::Float64})
       @ ITensors ~/.julia/packages/ITensors/23eRw/src/tensor_operations/tensor_algebra.jl:158
--- the last 11 lines are repeated 3084 more times ---
 [33936] contract(::ITensor, ::Vararg{ITensor}; kwargs::@Kwargs{cutoff::Float64})
       @ ITensors ~/.julia/packages/ITensors/23eRw/src/tensor_operations/tensor_algebra.jl:168
Some type information was truncated. Use `show(err)` to see complete types.

Not sure what’s going. A help will be appreciated :slight_smile:

Best,
Sourav

It’s hard to say what’s going wrong from that error message. I would suggesting printing out the indices (@show inds(T)) for the tensors that you are contracting, then making sure that the resulting tensor will be of a reasonable size, and also making sure that the indices which match (and thus will be automatically summed over) are the correct ones you are trying to contract.

Hi @sourav - just wanted to say that if your main focus or aim is to calculate the negativity (or log-negativity) and its cousin quantities obtained via partial transposition, there are other more simple ways to do so within itensor (and they have been discussed elsewhere in this forum - just wished to point this out in case you missed it)… of course if your aim is not to just calculate negativity anyhow but with only this specific approach you’ve described above, my comment is mute : )

Hi @adityab999 , I indeed require to compute log negativity. Can you point me out what other easier method(s) one can use for that in ITensor?

Best,
Sourav.

see for instance the discussions in the following links - Calculating the Log Negativity - #2 by adityab999 and Reduced density matrix of nonadjacent sites - #7 by adityab999, and Retrieve bipartite reduced density matrix - #6 by miles