Reduced density matrix of nonadjacent sites

DMRG and Numerical Methods
1

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 somehting 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

2

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

4 
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
5

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.

6

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

7

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.