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 lli = leftLinkIndex(psi,ii);
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);
psidag.prime("Link");
auto lLI = leftLinkIndex(psi,ii);
// 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){
auto rLI = rightLinkIndex(psi,ii);
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