Tracing out sites inbetween two four site subsets of a 1D-quantum-spin system

+1 vote
retagged

Hi,

I have calculated the ground state and its energy using the implemented DMRG algorithm.
Is there a possibility to trace out the sites inbetween subsets A and B of an open 1D quantum spin chain. In my specific case A and B consist of each 4 spins on the edges of the chain.
I would like to calculate the entropy of these two subsets combined (S_{AB}).

Thanks, Mathias

(I am using Julia with ITensors v0.2.3)

commented by (480 points)
edited
The documentation I've found for this is for the C++ code, so you'd need to translate it to Julia: http://www.itensor.org/docs.cgi?page=formulas/mps_two_rdm&vers=cppv3
For your purpose, you'd need to replace the i and j with code for the portion of the reduced density matrix for blocks A and B, but summing over the intermediate sites and external sites works the same as in that (you should gauge the MPS to a point in the space of interest, say the leftmost site of region A).

The questioner here https://www.itensor.org/support/229/evaluate-block-entanglement-blocks-that-extend-the-lattice?show=230#a230 has possibly not the best motivations, but the content from Miles is still really useful.

Another question on this that's in my bookmarks: https://www.itensor.org/support/1907/entanglement-entropy-of-a-finite-block-of-a-spin-chain?show=1907#q1907

If your reduced density matrix can be well-approximated with a product state of the separate reduced density matrices of A and B, then do that. If you're looking for something like the mutual information of A and B, this can be a very long calculation in my naive attempts to code it. For the Kitaev Wire model (a very simple model), the code stalled out (took too dang long) at L=23 sites because of the joint density matrix needed for mutual information calculations on my Linux laptop. Something that I did not test with the Kitaev model is reducing the number of bonds to make such calculations more feasible - my last edit of that code has the final maxBondDimension = 250 and a cutoff = 1E-10, it could certainly be tweaked and perhaps find better results (or perhaps make a copy to post-process and reduced in size to capture the entropy but maybe a less-accurate energy... Untested by me).

Miles has mentioned some sampling method to approximate the entanglement entropy for the reduced joint density matrix of AB, but to my knowledge he hasn't shared that code. I'm not sure what other tricks can be done with MPSs for things on opposite ends of the lattice.

Best of luck,
Jared

Edit: assuming you have qubits/two-state local spaces, that gives a 2^8=256 for the local space size for the reduced density matrix of AB - that should be within the realm of calculation.
commented by (170 points)
Hi Jared,

There's one point that still worries me: So far as I have understood it, I need to diagonalize rho = psi^2 of the sites I would like to trace out. This will scale exponentionally with the length of the spin chain, since the amount of sites I would like to keep is constant(?).

I guess you've encountered that in the Kitaev Wire model and Miles mentioned it in one of his answers as well.

I will implement it anyways and see how it performs.

Thanks and kind regards,
Mathias
commented by (480 points)
edited
Hello Mathias, it *can* be a little better, but it's not something that I've found a clean answer to. I'll post an official answer below shortly to get indentation correct (to make it more readable), but it'll be for the C++ version... Not ideal, but I'm not currently able to take the extra time for the Julia side.

+1 vote

Here is some mildly adjusted code from some C++ code that would need to be translated to Julia. It's posted as an answer to get the code portion to be nicely formatted.

The main trick is that assuming you have a globally pure state for a system divided into three regions [A|B|C], and you'd like to compute the entanglement entropy of AC, then you can also look at just B, and take whichever problem is smaller since the Schmidt values are the same inside/outside. This does produce a few if-thens, but it speeds up calculations.

• La, Lb, Lc are the lengths of regions A,B,C (La=Lc for convenience)
• Bl, Br are the leftmost and rightmost site indices of the interior region, B.
• Lp1by1m1 is just a way of calculating how many non-trivial ways we can divide the lattice into [A|B|C] such that La=Lc.

While, it does seem to indicate the topological phase transition in the Kitaev wire model for initial testing, I've not tested this thoroughly, so there could certainly be errors.

Here's the snippet of code:

// Define Lattice Length
int L = 20 ;
int Lp1by2m1 = (L+1)/2 - 1 ; // need integer division here. Use a floor function if needed.

println("L is \n");
std::cout << L << " \n ";
println("Lp1by2m1 is \n");
std::cout << Lp1by2m1 << " \n ";

int NumberOfStates = 5;

float SvN1array[NumberOfStates][L-1];
float SvN1fracfoundarray[NumberOfStates][L-1];
float SvN2array[NumberOfStates][Lp1by2m1];
float SvN2fracfoundarray[NumberOfStates][Lp1by2m1];
float MIarray[NumberOfStates][Lp1by2m1];

//
//
//    Your code to find the states of interest here
//
//

//
// Entanglement spectra & related quantities
//

// Single-cut entanglement entropies (bi-partite)

float SvN1vec[L-1];
float SvN1fracfound[L-1];

// Single-cuts are very fast. May as well compute the entanglement entropy for all single-interior cuts across the lattice
for( int ii : range1(L-1))
{
// Gauge the MPS to site ii
psi.position(ii);

// Define the two-site wavefunction for sites (ii,ii+1)
//// with the bond between them being the bond for which you wish to compute entanglement entropy
auto si  = siteIndex(psi,ii) ;
auto [U1,S1,V1] = svd(psi(ii),{lli,si}) ;
auto u1 = commonIndex(U1, S1);
double SvN1 = 0.0 ;

double fractionfound = 0.0 ;
for( auto n : range1(dim(u1)) )
{
auto Sn = elt(S1,n,n);
auto p1 = sqr(Sn) ;
if( p1 > 1E-12) {
SvN1 += -p1*log(p1);
fractionfound += p1;
}
}
SvN1vec[ii-1] = SvN1;
SvN1fracfound[ii-1] = fractionfound;
SvN1array[nn][ii-1] = SvN1;
SvN1fracfoundarray[nn][ii-1] = fractionfound;

}

// Fraction found
//// sum(spectrum.eigs())
// Linear Entanglement Measure
//// maxval(spectrum.eigs())
// Separability
//// - spectrum.eigs *dot* log(spectrum.eigs() )
// Second Renyi
//// - log( spectrum.eigs ^ 2 )
// Minimum entropy -- (infinite order) Renyi entropy
//// - log (separability)

// Two-cut entanglement entropies for systems divided into regions [A|B|C]
///// Want to calculate S(AC) such that the lengths of A and C are the same
//// (still bi-partite, but with a less intuitive basis, not as efficient with MPS methods)
//// Able to be done w/ usual Schmidt decomp. if starting with a pure state across the lattice
////// Entanglement spectrum of a pure state is the same for AC and complement of AC (that is, region B)
////// Take the space of AC or B, and solve whichever problem is smaller

float SvN2vec[Lp1by2m1];
float SvN2fracfound[Lp1by2m1];
double evalcutoff = 1E-2;
// // Uncomment the for loop if wanting to loop over a range of lengths of region A
//  for( int La : range1(Lp1by2m1))
//  {
int Lc = La;
int Bl = La + 1; // Bleft - leftmost site of region B
int Br = L  - Lc;
int Cl = Br + 1;
int Cr = L ;

double SvN2 = 0.0 ;
double fractionfound = 0.0 ;
// Remove any primes previously applied
psi.noPrime();
if ( Bl != Br ) {

if ( Br-Bl < L/2 ) {
// construct Psi for region B, then make density matrix

// Gauge the MPS to site Bl
psi.position(Bl);
auto psiB = psi(Bl) ;
// Loop over sites between Bl and Br, exclusive
for(int ii = Bl+1; ii < Br; ++ii)
{
psiB *= psi(ii);
}
psiB *= psi(Br) ;

auto rhoB = prime(psiB,"Site") * dag(psiB);

valcutoff = 1E-4;
auto [Q,D] = diagHermitian(rhoB, {"Cutoff=",evalcutoff,"ShowEigs=",false, "Truncate=", true});
// This is of rho, not just psi! No need to square!!!
auto u2 = commonIndex(Q, D);
// println("dim(u2) is \n") ;
// std::cout << dim(u2) << " \n ";
fractionfound = 0.0;

for( auto n : range1(dim(u2)) )
{
auto p2 = elt(D,n,n);
auto p2 = sqr(Sn) ;
if( p2 > evalcutoff) {
SvN2          += -p2*log(p2);
fractionfound += p2;
// println("fraction found so far is \n");
// std::cout << fractionfound << " \n ";
}
}

SvN2vec[La-1] = SvN2;
SvN2fracfound[La-1] = fractionfound;
SvN2array[nn][La-1] = SvN2;
SvN2fracfoundarray[nn][La-1] = fractionfound;

} else {   // Bl != Br   and   (Br-Bl) > L/2

// Gauge the MPS to site 1
psi.noPrime();
psi.position(1);
// Remove any primes previously applied
auto psidag = dag(psi);
// Initiallize reduced density matrix to begin construction over subsystem A
auto rhoAC = prime(psi(ii),lLI) * prime(psidag(ii),"Site");

for(int ii : range1(L) ){
// Construct rho over region A
if( ii < Bl){
if (ii > 1) {
rhoAC *= psi(ii) * prime(psidag(ii),"Site");
} // Region A execpt site 1
} else if (ii < Cl) { // construct such that region B is summed over
rhoAC *= psi(ii) * psidag(ii);
} else { // ii >= Cl Write over region C
if(ii < Cr){
rhoAC *= psi(ii) * prime(psidag(ii),"Site");
} // C region Sites except for final site region C
if(ii == Cr){
rhoAC *= prime(psi(ii),rLI) * prime(psidag(ii),"Site");
} // Final site region C
} // A, B, C region if
} // ii dummy loop

valcutoff = 1E-4;
auto [Q,D] = diagHermitian(rhoAC, {"Cutoff=",evalcutoff,"ShowEigs=",false, "Truncate=", true});
// This is of rho, not just psi! No need to square!!!
auto u2 = commonIndex(Q, D);
// println("dim(u2) is \n") ;
// std::cout << dim(u2) << " \n ";
fractionfound = 0.0;

for( auto n : range1(dim(u2)) )
{
auto p2 = elt(D,n,n);
auto p2 = sqr(Sn) ;
if( p2 > evalcutoff) {
SvN2          += -p2*log(p2);
fractionfound += p2;
// println("fraction found so far is \n");
// std::cout << fractionfound << " \n ";
}
}

SvN2vec[La-1] = SvN2;
SvN2fracfound[La-1] = fractionfound;
SvN2array [nn][La-1] = SvN2;
SvN2fracfoundarray [nn][La-1] = fractionfound;

} //  Bl != Br   and   (Br-Bl) > L/2

} else { // if Bl == Br  // then just regular old single cut entanglement entropy

psi.noPrime();
// Gauge the MPS to site Bl
psi.position(Bl);
auto rhoB = prime(psi(Bl),"Site") * dag(psi(Bl));

evalcutoff = 1E-4;
auto [Q,D] = diagHermitian(rhoB, {"Cutoff=",evalcutoff,"ShowEigs=",false, "Truncate=", true});

auto u2 = commonIndex(Q, D);
fractionfound = 0.0 ;
for( auto n : range1(dim(u2)) )
{
auto p2 = elt(D,n,n);
auto p2 = sqr(Sn) ;
if( p2 > evalcutoff) {
SvN2          += -p2*log(p2);
fractionfound += p2;
}
}
SvN2vec[La-1] = SvN2;
SvN2fracfound[La-1] = fractionfound;

SvN2array[nn][La-1] = SvN2;
SvN2fracfoundarray[nn][La-1] = fractionfound;

} // Bl = Br

//  }   // END FOR over La

// Put this in a loop if varying La
MIarray[nn][La-1] = SvN1array[nn][La-1] + SvN1array[nn][L-(La-1)-2] - SvN2array[nn][La-1] ;

} // nn state of interest loop

+1 vote