Extract/Project IQTensor on a QN subspace

Question

Hi all,

I am fairly new to ITensor and C++.
Thus if you have any advice that can be easier implemented in Julia I could take it into account, too.

First let me describe the system I am dealing with:
I consider a Bosonic chain of N sites: Boson(N,{"ConserveQNs=",true,"MaxOcc",1})).
I would like to construct the reduced density matrix of the system and project it on a given subspace.
In formulas, given the reduced density matrix @@\rho_A=Tr(|\psi><\psi|)@@
I would like
$$
\rho_A(q)=\frac{P_q \rho_A P_q}{Tr(\rho_A P_q)}
$$
Where @@P_q@@ is the projector on the subspace with a given number of particles
$$
\sum_{i=1}^{l}n_i=q.
$$
Here @@n_i@@ is the number of Bosons at site i and l is the length of the subsystem.

The question is the following: is there a way to do obtain @@\rho_A(q)@@, using IQTensor, without computing explicitly the reduced density matrix and the projector, since it is very inefficient?

I saw in the source code of qutil.h the function getBlock, but I'm not used to C++ and I could not deduce its functionalities.

Thank you for your time.

Best,
Vittorio

Comments

Hi Vittorio, yes I think there should be a good way to do this, with mostly tools we have available in the code already. (Yes it would be easier to do some of this with the Julia version since it is more flexible and offers more low-level tools to manipulate QN blocks.)

I assume your state |psi> is an MPS correct?

Best,
Miles

---

Yes, It is an MPS.

Vittorio

Answer 1

Hi all,

I don't know if it is allowed to answer in the same post but I reconsidered the problem after a while an came with the following idea.

Considering a chain of L sites, with two species of bosons, I wanted to construct the reduced density matrix of the first 4 sites with a given flux QN(("Na",2),("Nb",2)), Na and Nb are the number of bosons of type a and type b. The values of Na and Nb are taken arbitrarly in this example.

I tried this:

ψ=orthogonalize(starting_MPS,1)
ψdag = dag(ψ)
r1=linkind(ψ,1)
site1=1 #first site
site2=4 #last site
ρA = prime(ψ[site1],r1)*prime(ψdag[site1],"Site")
for j in site1+1:site2-1
lj = linkind(ψ,j-1)
rj = linkind(ψ,j)
ρA *= prime(prime(ψ[j],lj),rj)
ρA *= prime(ψdag[j],"Site")
end
l2 = linkind(ψ,site2-1)
l2p = prime(linkind(ψdag,site2-1))
r2 = prime(linkind(ψ,site2))
r2p = prime(linkind(ψdag,site2))
ρA *= prime(ψ[site2],l2)
ρA *= delta(QN(("Na",2),("Nb",2)),l2,l2p)
ρA *= delta(QN(("Na",2),("Nb",2)),r2p,r2)
ρA *= prime(ψdag[site2])

and realized that delta(QN(("Na",2),("Nb",2)),l2,l2p) and delta(QN(("Na",2),("Nb",2)),r2,r2p) were empty.

At first I thought I was selecting an incorrect values of the flux but it is compatible with l2, which is the following

(dim=64|id=983|"Link,fact")
1: QN(("Na",4),("Nb",4)) => 1
2: QN(("Na",4),("Nb",3)) => 3
3: QN(("Na",4),("Nb",2)) => 3
4: QN(("Na",4),("Nb",1)) => 1
5: QN(("Na",3),("Nb",4)) => 3
6: QN(("Na",3),("Nb",3)) => 9
7: QN(("Na",3),("Nb",2)) => 9
8: QN(("Na",3),("Nb",1)) => 3
9: QN(("Na",2),("Nb",4)) => 3
10: QN(("Na",2),("Nb",3)) => 9
11: QN(("Na",2),("Nb",2)) => 9
12: QN(("Na",2),("Nb",1)) => 3
13: QN(("Na",1),("Nb",4)) => 1
14: QN(("Na",1),("Nb",3)) => 3
15: QN(("Na",1),("Nb",2)) => 3
16: QN(("Na",1),("Nb",1)) => 1

for this reason I assumed the problem was in l2p in the sense that dag(index) inverts the sign of the flux as I understood correctly looking at the documentation.

However delta(QN(("Na",2),("Nb",-2)),l2,l2p) does not work as well and I'm a bit confused about the construction of this delta tensor.

For example delta(QN(("Na",2),("Nb",-2)),l2,l2) works and selects the right sector but the indexes are both outwards and the contraction with the other tensors becomes cumbersome.

Thank you for your time and I apologize if this answer goes against the policy of this forum. It is not written in the guidelines.

Best,
Vittorio

EDIT: The idea of multiplying by the delta tensor is to have a tensor which is different from zero only in the subspace I want to construct the reduced density matrix.

EDIT2: I managed to do what I wanted using setelt in some fashion like this one:

temp=setelt(r2p=>96,r2=>96)
for n in 97:124
temp+=setelt(r2p=>n,r2=>n)
end

It is not pretty nor automatic but it works. I put it here in case anyone is interested.

Comments

Hi Vittorio, I would say that it is great that you wrote an answer, even an especially on your own post. Good to have this information here!

Answer 2

To add on to Vittorio's post: as he deduced, @@delta@@ isn't quite the right tool for this case, since it is really meant for making a diagonal tensor with a uniform value along the diagonal (and right now I think only works if all of the diagonal blocks exist, though in principle it could be generalized). Here is a more general approach using some internal and custom made tools:

using ITensors

# Get the range of values of a block
blockrange(i::Index, b::Block) = _blockrange(i, Int(b))
function _blockrange(i::Index, b)
  start = 1
  for j in 1:b-1
    start += ITensors.blockdim(i, j)
  end
  return start:start+ITensors.blockdim(i, b)-1
end

function blockrange(i::Index, find_qn::QN)
  # The block the QN is in
  b = ITensors.findfirstblock(x -> qn(x) == find_qn, i)
  return blockrange(i, b)
end

space = [QN(("Na",4),("Nb",4)) => 1, QN(("Na",4),("Nb",3)) => 3, QN(("Na",4),("Nb",2)) => 3, QN(("Na",4),("Nb",1)) => 1, QN(("Na",3),("Nb",4)) => 3]

i = Index(space)

# The QN we are projecting into
proj_qn = QN(("Na",4),("Nb",2))

r = blockrange(i, proj_qn)

T = emptyITensor(i', dag(i))
for n in r
  T[n, n] = 1.0
end

# Alternative, remove the blocks completely
i2 = Index([proj_qn => length(r)])
T2 = emptyITensor(i2', dag(i))
for n in 1:length(r)
  T2[n, n+first(r)-1] = 1.0
end

It is essentially the same as the @@setelt@@ version but uses indexing instead which should be faster. In addition, this handles an arbitrary QN, where the range of the QN is determined by the @@blockrange@@ function. This would be a good function to have in ITensors.jl to make operations like this easier.

Comments

Thank you Matt!

I finally wrote this to construct the tensor which selects only the QNs present in QNlist:

temp=ITensor(r2,r2p)
for n in 1:dim(r2)
       if qn(splitblocks(r2),n) in QNlist
            temp+=setelt(r2p=>n,r2=>n)
       end
end

I guess yours is faster so I'll try as soon as I can.

Best,
Vittorio