Reduced density matrix of nonadjacent sites

Hank · July 1, 2024, 11:37am

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 · July 1, 2024, 11:49am

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

Hank · July 1, 2024, 12:07pm

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 · July 1, 2024, 1:06pm

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 · July 1, 2024, 1:47pm

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 · August 3, 2024, 12:47am

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 · December 25, 2024, 1:54am

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 · December 27, 2024, 9:13am

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 · December 30, 2024, 3:03pm

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

Topic		Replies	Views
Calculating the Log Negativity DMRG and Numerical Methods	3	139	July 1, 2024
Multiorbital electron mixed Hamiltonian DMRG and Numerical Methods julia	3	161	February 13, 2024
Larger linkdim in initial state leads to unexpected behavior in excited state DMRG DMRG and Numerical Methods	1	135	February 23, 2024
Mixed-site (Boson,spin) model Expected value of <ni> DMRG and Numerical Methods julia , mpo	5	82	October 5, 2024
Errors in the Ground-State Energy in a Next NN 1D Spin Chain DMRG and Numerical Methods dmrg , mpo	2	220	July 28, 2023

Reduced density matrix of nonadjacent sites

Related topics