Reduced density matrix of nonadjacent sites

[Hank]

Hello everyone.

With the hope you are doing well, I would like to address a concern regarding some calculations I am currently doing.

In fact, I have a 14-site alternating mixed-spin chain of spin-1/2 and spin-1 and I would like to get the reduced density matrix of say 4 sites (not necessarily adjacent). Let A, B, C, and D be the sites of interest.

Following some ideas proposed here before as well as suggestions, I wrote some lines but I currently have some issue due to the non-adjacent character of the sites of interest.

Here, I provide with an example of the code I wrote:

let
    c=1
  N = 14
  sites = siteinds("S=1/2",N)

  os = OpSum()
  for j=1:N-1
    os += 0.5,"S+",j,"S-",j+1
    os += 0.5,"S-",j,"S+",j+1
    os += "Sz",j,"Sz",j+1
  end
  H = MPO(os,sites)

  nsweeps = 5 # number of sweeps is 5
  maxdim = [10, 100, 1000, 10000, 25580] # gradually increase states kept
  cutoff = [1E-10] # desired truncation error

  psi0 = random_mps(sites;linkdims=10)

  energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff)
# Reduced density matrix between sites 1, 2, 3, 4 in the order 1234
    
psidag=prime(dag(psi),linkinds(psi))

rho1=prime(psi[1],sites[1])*dag(psi[1])
for i=2:4
    rho1 *= psi[i]
    rho1 *= psidag[i]
end

rho1_PT = swapinds(rho1,inds(rho1)[1],inds(rho1)[3])
    
    
############################################################################################
    
# Reduced density matrix between sites 1, 2, 3, 4 in the order 4123
orthogonalize!(psi,c+3)
ket2=psi[c+3]
ket2 *= psi[c]*psi[c+1]*psi[c+2]   
bra2 = dag(ket2)
rho2 = prime(ket2,"Site")*bra2 
rho2_PT = swapinds(rho2,inds(rho2)[1],inds(rho2)[3])
    
    
    
############################################################################################

# Reduced density matrix between sites 1, 2, 3, 4 in the order 3124
orthogonalize!(psi,c+4)
ket3= psi[c+4]
ket3 *= psi[c]*psi[c+1]*psi[c+3]   
bra3 = dag(ket3)
rho3 = prime(ket3,"Site")*bra3 
rho3_PT = swapinds(rho3,inds(rho3)[1],inds(rho3)[3])  
    

The last lines of each subpart represents a part transpose in each matrix, with the intention to check whether they are well-defined (diagonal summing up to 1).

The lines I wrote, following the first part of the code is the following:

a1=0
a2=0
a3=0

d1,u1 = eigen(rho1_PT)
d2,u2 = eigen(rho2_PT)
d3,u3 = eigen(rho3_PT)

for j=1:length(diag(real(d1)))
   a1+=diag(real(d1))[j]
end
for j=1:length(diag(real(d2)))
   a2+=diag(real(d2))[j]
end
for j=1:length(diag(real(d3)))
   a3+=diag(real(d3))[j]
end

As outputs, I obtain the following:

a1 = 0.9999999999999999
a2 = 0.8323756817707685
a3 = 31.99999999999999

Meaning that something wrong is happening.

Therefore, I would like to ask if there is a way you might know to deal with this situation because my main goal is to compute the log-negativity based on these 4-sites RDMs.

Thanks in advance for your time and consideration,

Sincerely,
Idriss

[ryanlevy]

Hi Idriss,
Could you try and get your post formatting working please? For code blocks, you’ll want 3 ` ticks, instead of 2

[Hank]

let
    c=1
  N = 14
  sites = siteinds("S=1/2",N)

  os = OpSum()
  for j=1:N-1
    os += 0.5,"S+",j,"S-",j+1
    os += 0.5,"S-",j,"S+",j+1
    os += "Sz",j,"Sz",j+1
  end
  H = MPO(os,sites)

  nsweeps = 5 # number of sweeps is 5
  maxdim = [10, 100, 1000, 10000, 25580] # gradually increase states kept
  cutoff = [1E-10] # desired truncation error

  psi0 = random_mps(sites;linkdims=10)

  energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff)
        
orthogonalize!(psi,2)

psidag=prime(dag(psi),linkinds(psi))

rho1=prime(psi[1],sites[1])*dag(psi[1])
for i=2:4
    rho1 *= psi[i]
    rho1 *= psidag[i]
end
rho1_PT = swapinds(rho1,inds(rho1)[1],inds(rho1)[3])

orthogonalize!(psi,c+3)
ket2=psi[c+3]
ket2 *= psi[c]*psi[c+1]*psi[c+2]   
bra2 = dag(ket2)
rho2 = prime(ket2,"Site")*bra2 
rho2_PT = swapinds(rho2,inds(rho2)[1],inds(rho2)[3])

orthogonalize!(psi,c+4)
ket3=psi[c+4]
ket3 = psi[c]*psi[c+1]*psi[c+3]   
bra3 = dag(ket3)
rho3 = prime(ket3,"Site")*bra3 
rho3_PT = swapinds(rho3,inds(rho3)[1],inds(rho3)[3])  


a1=0
a2=0
a3=0

d1,u1 = eigen(rho1_PT)
d2,u2 = eigen(rho2_PT)
d3,u3 = eigen(rho3_PT)

for j=1:length(diag(real(d1)))
   a1+=diag(real(d1))[j]
end
for j=1:length(diag(real(d2)))
   a2+=diag(real(d2))[j]
end
for j=1:length(diag(real(d3)))
   a3+=diag(real(d3))[j]
end

[ryanlevy]

Thanks, its much more helpful to read and copy/paste now.

The first thing I would as ask, is before the partial trace is it a normalized density matrix? So I remove your swapinds and did

d1,u1 = eigen(rho1,ishermitian=true)
d2,u2 = eigen(rho2,ishermitian=true)
d3,u3 = eigen(rho3,ishermitian=true)

d1 and d3 did not add up to 1 as I would expect.
For rho1, when you write

rho1=prime(psi[1],sites[1])*dag(psi[1])

you’ve traced over the link indices, forming a 1 site RDM. The rest of the multiplications of psi[i] and psidag[i] aren’t doing anything. I would also recommend doing the orthogonalize at the left most site, here site 1.

Another small comment, when you do swapinds(rho1,inds(rho1)[1],inds(rho1)[3]) you assume something of the order of the indices. This is dangerous as the order is not guaranteed.

rho2 seems to give a normalized density matrix.

For rho3, you’ll want to orthogonalize to the left most (or middle) site in your RDM, and you can’t skip over psi[c+2]. You’ll need to include that tensor and make sure it is traced out when applying primed “bra”.

In general I would recommend drawing a picture of of the tensor operations you are trying to implement and print out what is going on each step to make sure the tensors have the indices you expect.

[Hank]

Hi sir,

Thank you for the help as well as the suggestions. In fact, when I run the code, only a1 returns 1 whereas a2 and a3 do not.

Anyway, I would follow your suggestions and come back with a feedback asap.

Sincerely,
Idriss

[adityab999]

Not sure if it’s relevant anymore (perhaps you already figured it out), but in case it is (and also for record as this thread might come up in search results for many of us looking for similar questions) –

for nonadjacent sites such as for example a 3-site RDM with config. AB-C (where C is separated from AB by a distance of one site), just do the following :

orthogonalize!(psi1,A)
ket1 = psi1[A]
for i=A+1:C
    ket1 *= psi1[i]
end
bra1 = dag(ket1)
rho = prime(ket1,"Site")*bra1
#@show inds(rho)  <-- this shows that the indices of rho corresponding to the site between B and C are the 3rd and 7th element of the inds(rho) vector. 
rho_ABC = rho*delta(inds(rho)[3],inds(rho)[7])     # delta(inds,inds) traces out the given inds (see the documentation)

This will yield the correct sum of eigenvalues = 1 and the correct number of eigenvalues 2^i, where i=3 in this example.
This is now readily generalizable to setups like AB–CD, AB–C–D and so on.

[vyason]

Hi @adityab999 I am still confused. Lets say I have 4 sites AB – CD. A and B are consecutive and C and D are consecutive. but they are separted by J number of sites. Could you please specify a clean code for the reduced density matrix?

A

[adityab999]

Hi, the code example I posted above is readily generalizable to your question as well… first create a reduced density matrix for sites A to D, then use show(inds) to see which indices correspond to the J number of sites between B and C and then use the delta(inds,inds) function to trace over them so that what remains is the reduced density matrix of AB - CD.

Just note that the J or |D-A| lengths can not be large, as constructing reduced density matrices is exponentially costly in the length.

[vyason]

Thank you @adityab999

This is the only solution that finally worked. It is true that the distance between the pair of sites must be small, but it is working very well for my problem

Just to mention, there is a way to get the RDM for sites that are arbitrarily far away, but it only works for two-site RDM rho_AB. I am attaching the code here:

orthogonalize!(ket,A)
bra = prime(dag(ket),linkinds(ket))
  
rho = prime(ket[A],linkinds(ket,A-1)) * prime(bra[A],sites[A])
for for j in A+1:B-1
    rho *= ket[j]*bra[j]
end
rho *= prime(ket[B],linkinds(ket,B)) * prime(bra[B],sites[B])  

Can you see a solution for extending it for 3-site or 4-site RDMs? If someone can figure it out, it will be wonderful!

Ankit

[adityab999]

Hi,

Sorry for the late reply. Yes I am familiar with that approach, it was posted originally by someone long ago in this forum… For quite some time I’d tried to get around to making it work for RDMs for 3 or more spins with arbitrary separation but it never worked completely correctly… I’d managed to modify it for 3 or more spins to give out results in numerical agreement with numbers obtained by other methods, but there always remained two unresolved issues with it whenever dealing with non-contiguous series of spins –

1. the number of eigenvalues was always wrong, yet their total sum was always = 1 as it should be… this was because, in addition to the correct number of correct eigenvalues, a number of “zeros” (1e-10 or lower) were getting added to the eigenvalue vector for some reason that I never fully understood or could rectify…

2. the indices corresponding to the tensors at each site were a mess and showed many spurious-looking objects at the intermediate sites that were supposed to be traced over (I believe the appended zeros mentioned above were due to these things)… As a result, while you could get agreeable results for, for example, von Neumann entropy or mutual info and such objects, you couldn’t do any operations such as the partial transpose, etc with it.

It would be indeed great if this method can be made to work for 3 or more spins with arbitrary separation.

[vyason]

Hi @adityab999

Yes. the original is posted here:

Is this a correct way to calculate one-site and two-site entanglement spectrum and entropy?

Hello everyone! Here I want to calculate the one-site and two-site entropy to obtain the mutual information. With the guides of these disccusion 1-RDM from Julia version - ITensor Julia Questions - ITensor Discourse, Mutual Information matrix - DMRG and Numerical Methods - ITensor Discourse and Two-site reduced density matrix - ITensor Julia Questions - ITensor Discourse, I wrote some test codes to ahieve my goal. But nearly all of the returned eigenvalues are zero, except for one value that …

I also tried a lot for extending it to 3-site RDM. However:

1. The final tensor only has the indices the indices of the sites on the very left and the very right, but not for the intermediate one.
2. There are alwyas leftover links in the resulting tensor object.

Hope somebody will figure it out.

A

[miles]

@ryanlevy recently introduced a new ITensor Entropy Tools package. One of its main features is efficiently computing rdms. Does it help with the case you wanted?

See post below:

Introducing ITensorEntropyTools.jl

Introducing a new 3rd party package to compute entanglement entropy ITensorEntropyTools.jl: Features The idea behind ITensorEntropyTools.jl is to automate some common calculations/tools: Entanglement Entropy calculations von Neumann and Rényi Bipartite regions Arbitrary (potentially-disconnected) regions Automatic calculation of (potentially) minimal rank RDMs Mutual information Tripartite mutual information Experimental Support ITensorNetworks.jl (see network branch) Entanglement negat…