string correlation function of several sites

Hi,
I am new to ITensor, so please forgive me if my question is naive. I want to get the results of string correlation function like <psi | N_i N_i+1 N_i+2 … N_L|psi>. From searching the categories, I finded it that I can achieve my goal as far as using the OpSum interface:

os = OpSum()
os *= "N",i,"N",2,"N",3,...,"N",L

But, at the end, I encountered some errors. I could not find the problems from the tensors. Here are my codes:

using ITensors
L=10;cutoff=1E-8;δτ=0.1;beta_max=2.0;
function ITensors.op(::OpName"expτnn", ::SiteType"Boson", adag::Index,amis::Index; τ)
    h =
      1 / 2 * op("a†", adag) * op("A", amis) +
      1 / 2 * op("A", amis) * op("a†", adag) +
      op("N", adag) * op("N",amis)
    return exp(τ * h)
end
s = siteinds("Boson", L; conserve_number=false);
gates = ops([("expτnn", (e, e + 1), (τ=-δτ / 2,)) for e in 1:(L - 1)], s);
append!(gates, reverse(gates));
rho = MPO(s, "Id") ./ √2;
terms = OpSum();
U=4;V=4;
for j in 1:(L - 1)
    terms += "a†", j, "A", j + 1
    terms += "A", j, "a†", j + 1
    terms += U / 2, "N", j, "N", j
    terms += V, "N", j, "N", j + 1
end
H = MPO(terms, s)
for w in 0:δτ:beta_max
    energy = inner(rho, H)
    @printf("β = %.2f energy = %.8f\n", w, energy)
    rho = apply(gates, rho; cutoff)
    rho = rho / tr(rho)
end
τ_range = δτ:δτ:(beta_max / 2)
psi = randomMPS(s);
for τ in τ_range
    psi = apply(gates, psi; cutoff)
    normalize!(psi)
end
n_2=op("N",s[2]);

n_3=op("N",s[3]);

n_4=op("N",s[4]);

n_5=op("N",s[5]);

orthogonalize!(psi,2)

c=psi[2];

c *= n_2;

ir2=commonindex(psi[2],psi[3],"Link");

c *=dag(prime(prime(psi[2],"Site"),ir2));

c *=psi[3];

c *=n_3;

c *=dag(prime(prime(psi[3],"Site"),"Link"));

c *=psi[4];

c *=n_4;

c *=dag(prime(prime(psi[4],"Site"),"Link"));

c *=psi[5];

c *=n_5;

il5=commonindex(psi[5],psi[4],"Link");

c *=dag(prime(prime(psi[5],"Site"),il5));

scalar(c)

The error tolds me that “In scalar(T) or T, ITensor T is not a scalar.” Please tell me the reasons and give me a hand.
Thanks a lot!

Hi, thanks for the question. Your code looks really good, so I think you are close to getting it working.

If you could print out inds(c) at the end, right before you call scalar on it, what indices does c have? Based on the error, it has one or more indices and is not actually a scalar. But your code seems correct i.e. I would also have thought it would give a scalar. So printing out inds(c) would help to see which indices were not contracted.

1 Like

HI,miles.
Thanks a lot for your quick reply. I worked hard to do the contraction, but I still could not. I think there might be some minor errors. Here are my codes:

orthogonalize!(psi,2)

n_2=op("N",s[2]) /ITensor ord=2 (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=862|"Boson,Site,n=2")
n_3=op("N",s[3])
n_4=op("N",s[4])
n_5=op("N",s[5])

c=psi[2] /ITensor ord=3 (dim=2|id=574|"Link,n=1") (dim=4|id=952|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")'
c *= n_2 /ITensor ord=3 (dim=2|id=574|"Link,n=1") (dim=4|id=952|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")'
ir=commonind(psi[2],psi[3],"Link")
c *=dag(prime(prime(psi[2],"Site"),ir)) /ITensor ord=6 (dim=2|id=574|"Link,n=1") (dim=4|id=952|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=4|id=952|"Link,n=1")'

c *=psi[3] /ITensor ord=6 (dim=2|id=574|"Link,n=1") (dim=4|id=952|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=4|id=952|"Link,n=1")'
c *=n_3 /ITensor ord=7 (dim=2|id=574|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=4|id=952|"Link,n=1")' (dim=5|id=720|"Link,n=1") (dim=2|id=253|"Boson,Site,n=3")'
c *=dag(prime(prime(psi[3],"Site")))  /ITensor ord=8 (dim=2|id=574|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=5|id=720|"Link,n=1") (dim=2|id=253|"Boson,Site,n=3")' (dim=2|id=253|"Boson,Site,n=3")'' (dim=5|id=720|"Link,n=1")'

c *=psi[4] /ITensor ord=9 (dim=2|id=574|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=2|id=253|"Boson,Site,n=3")' (dim=2|id=253|"Boson,Site,n=3")'' (dim=5|id=720|"Link,n=1")' (dim=2|id=546|"Boson,Site,n=4") (dim=4|id=508|"Link,n=1")
c *=n_4 /ITensor ord=9 (dim=2|id=574|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=2|id=253|"Boson,Site,n=3")' (dim=2|id=253|"Boson,Site,n=3")'' (dim=5|id=720|"Link,n=1")' (dim=4|id=508|"Link,n=1") (dim=2|id=546|"Boson,Site,n=4")'
c *=dag(prime(prime(psi[4],"Site"))) /ITensor ord=10 (dim=2|id=574|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=2|id=253|"Boson,Site,n=3")' (dim=2|id=253|"Boson,Site,n=3")'' (dim=4|id=508|"Link,n=1") (dim=2|id=546|"Boson,Site,n=4")' (dim=2|id=546|"Boson,Site,n=4")'' (dim=4|id=508|"Link,n=1")'

c *=psi[5] /ITensor ord=11 (dim=2|id=574|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=2|id=253|"Boson,Site,n=3")' (dim=2|id=253|"Boson,Site,n=3")'' (dim=2|id=546|"Boson,Site,n=4")' (dim=2|id=546|"Boson,Site,n=4")'' (dim=4|id=508|"Link,n=1")' (dim=2|id=223|"Boson,Site,n=5") (dim=5|id=977|"Link,n=1")
c *=n_5  /ITensor ord=11 (dim=2|id=574|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=2|id=253|"Boson,Site,n=3")' (dim=2|id=253|"Boson,Site,n=3")'' (dim=2|id=546|"Boson,Site,n=4")' (dim=2|id=546|"Boson,Site,n=4")'' (dim=4|id=508|"Link,n=1")' (dim=5|id=977|"Link,n=1") (dim=2|id=223|"Boson,Site,n=5")'
il5=commonindex(psi[5],psi[4],"Link")
c *=dag(prime(prime(psi[5],"Site"),il5)) /ITensor ord=12 (dim=2|id=574|"Link,n=1") (dim=2|id=862|"Boson,Site,n=2")' (dim=2|id=574|"Link,n=1")' (dim=2|id=862|"Boson,Site,n=2")'' (dim=2|id=253|"Boson,Site,n=3")' (dim=2|id=253|"Boson,Site,n=3")'' (dim=2|id=546|"Boson,Site,n=4")' (dim=2|id=546|"Boson,Site,n=4")'' (dim=5|id=977|"Link,n=1") (dim=2|id=223|"Boson,Site,n=5")' (dim=2|id=223|"Boson,Site,n=5")'' (dim=5|id=977|"Link,n=1")'

With Regards,
Yuanyu

Thanks for providing that output – it’s helpful. I found the problem, which is a bit subtle.

If you print out the variable ir returned from commonindex you will see that it is equal to nothing rather than an Index. This is because you should actually call commonind this way:

ir = commonind(psi[2],psi[3];tags="Link")

where I have put a semicolon then passed "Link" as the value of the keyword argument named tags. This way commonind knows that "Link" is referring to a specific tag value.

Please try that and let me know if it works! We will improve the documentation for commonind to make it clearer what keyword arguments you can pass to it. It’s a bit confusing right now that it lets you pass the string "Link" without giving an error.

1 Like

Is it necessary to use the tag filtering at all, or is ir = commonind(psi[2], psi[3]) sufficient? I usually prefer to avoid using tags to filter when possible, usually it is avoidable by using set/network functions like siteinds, linkinds, commoninds, uniqueinds`, etc.

1 Like

Agreed, it would be totally sufficient here to use

ir = commonind(psi[2],psi[3])

also, since for an MPS by definition there is only one common index between neighboring tensors.

I guess that isn’t really constrained anywhere and in principle someone could make an MPS with more than one index shared between neighboring tensors, but usually there is only one.

The safest thing would be to use:

ir = only(commoninds(psi[2], psi[3]))

which errors if there are either zero or more than one shared Index.

Hi
The results are a bit amusing,

commonind(psi[2],psi[3],tags="Link") /(dim=4|id=507|"Link,n=1")
commonind(psi[2],psi[3]) /(dim=4|id=507|"Link,n=1")
commonind(psi[2],psi[3],"Link") /nothing

The tags are avoidable or need to be correct.

In the end, I got:

ir=commonind(psi[2],psi[3],tags="Link") /(dim=4|id=507|"Link,n=1")
...
il5=commonindex(psi[5],psi[4],tags="Link") / il5=commonindex(psi[5],psi[4],tags="Link")
...
c *=dag(prime(prime(psi[5],tags="Site"),il5)) /ITensor ord=4 (dim=4|id=507|"Link,n=1")' (dim=4|id=507|"Link,n=1") (dim=4|id=263|"Link,n=1") (dim=4|id=263|"Link,n=1")'

Is it possible to change the id number?

At this point I think you should have all the pieces you need to have a working code. In terms of how to combine these pieces correctly, I would strongly recommend drawing all of the tensors as diagrams on paper or on a board, then marking which indices are primed and by how much, in order to reason about how the next contraction will go and what priming to do next. Then I think you will find a pattern that works such that c always has the correct number of indices that you want.

1 Like

Hi miles,

Thanks for your patiently reply, you have helped me save a lot time. I think I will try to code in C.

With regards
Yuanyu

Note that you can try out this (slightly) experimental package to compute general correlators: GitHub - ITensor/ITensorCorrelators.jl.

EDIT: There is a slightly more extensive example here: https://github.com/ITensor/ITensorCorrelators.jl/blob/main/examples/4_point_correlator.jl which compares against using a more naive approach with OpSum/MPO, which you can run by including the example file and running main(; N=20) where N is the system size.

2 Likes

Hi,
Thanks for your reply. A few days ago, I have tried to do this update but got no response from system. Luckily, I success this time.
With regards,
Yuanyu