Calculatation code for Von Neumann entropy

DMRG and Numerical Methods
, , ,
1

Hi all,

I want to calculate the entanglement entropy of a fermion chain, and I find the following code.
(Entanglement entropy (Julia) - ITensor Support Q&A).
The code is

function entropy_von_neumann(psi::MPS, b::Int)
  s = siteinds(psi)  
  orthogonalize!(psi, b)
  _,S = svd(psi[b], (linkind(psi, b-1), s[b]))
  SvN = 0.0
  for n in 1:dim(S, 1)
    p = S[n,n]^2
    SvN -= p * log(p)
  end
  return SvN
end

I would like to ask if the b mentioned in this program refers to the entanglement entropy of the two subsystems before and after the truncation at b position

But then I see a code that adds a “for loop”. So what does b for i in 2:N-1 mean in this code?

function entropyvonneumann(psi::MPS, N::Int)
    SvN = fill(0.0, N)
    for b in 2:N-1
    orthogonalize!(psi, b)
    _,S = svd(psi[b], (linkind(psi, b-1), s[b]))
    for n in dim(S, 1)
    p = S[n,n]^2
    SvN[b] -= p * log(p)
    end
    println("$(SvN[b])")
    end
    return SvN
end

Please tell me if I need to add a “for loop” if I want to calculate the entanglement entropy of a subsystem that cuts off the first part and the second part at parameter b.

2

I’m not sure what that for loop is doing either, and I did not see any such for loop at the link you provided.

I would recommend you use the code example from our most current documentation here:
<MPS and MPO Examples · ITensors.jl >
In that code, “b” refers to a bond number. So like bond number 1,2,3, … of the MPS. The entanglement computed is between sites 1,2,…b and sites b+1,b+2,…,N. It is the “bipartite von Neumann” entanglement entropy between those two groups of sites.

3

Thank you very much for your answer. I think I understand.

4

Hi Miles,

Sorry if I should have started a new post, but this is just a quick technical question.

I noticed that the svd function returns spec, where I believe spec.eigs are the abs2 of the svdvals. This allows for a simplified code formula:

function entanglement(ψ, b)
    orthogonalize!(ψ, b)
    _, _, _, spec = svd(ψ[b], linkind(ψ, b))    # cut the bond between b & b+1
    SvN = sum(p -> p > 0.0 ? -log(p)*p : 0.0 , spec.eigs)
    SvN = SvN/log(2)
    return SvN
end

Please let me know if I’ve missed any subtleties.

Thanks!

5

I agree that code should work correctly, and it is nice and short. It looks correct to me, but please do verify it on some test cases of your own to be 100% sure.