Entanglment entropy calculation for SSH model

DMRG and Numerical Methods
1

Hi,
I have been trying to calculate the entanglement entropy for the SSH model with open boundary conditions. I’m very new to the Julia language but to the best of my knowledge, I believe I’ve coded it correctly. Below is the full code that I’ve written following different discussions here as well as in the ITensor documentation.
Im trying to calculate the bipartite entropy for a fixed lattice size, but with one of the hopping parameters varying between a certain range.

using ITensors
using Plots

@time let
  N = 60
  t1 = collect(0.5:0.01:2.5)

  
  function SSH_OBC(N,t1)
    sites = siteinds("Fermion", N; conserve_qns=false)
    t2 = 1.8
    os = OpSum()
    
    for b in 1:(N-1)
      if isodd(b)
        os += -t1, "Cdag", b, "C", b + 1
        os += -t1, "Cdag", b+1, "C", b

      elseif iseven(b)
        os += -t2, "Cdag", b, "C", b + 1
        os += -t2, "Cdag", b+1, "C", b
      end
    end
    H = MPO(os, sites)
    return sites,H
  end

  function ent(psi::MPS; kwargs...)
    N = length(psi)
    b = ÷(N,2)
    orthogonalize!(psi, b)
    s = siteinds(psi)
    U,S,V = 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

  energies = []
  states = MPS[]
  SvN = 0.0
  SvN_loop = []

  
  for t1 in t1
    sites,Hlist = SSH_OBC(N,t1)
    state = []
     
    for n in 1:N
      if isodd(n)
        push!(state,"1")
      elseif iseven(n)
        push!(state,"0")
      end
    end

    psi0 = randomMPS(sites,state,10)
    nsweeps = 15
    maxdim = [50, 100, 200, 400, 800, 800]
    cutoff = [1E-12]
    println("---------------------------------------------------------")
    println("Running DMRG for t1 = $t1")
    energy, psi = dmrg(Hlist, psi0; nsweeps, maxdim, cutoff)
    println("---------------------------------------------------------")

    SvN = ent(psi)


    push!(energies,energy)
    push!(states,psi)
    push!(SvN_loop,SvN)

  end
  println("t1\t G.S Energy\t Entropy")
  for j in 1:length(energies)
    println("$(t1[j])\t $(energies[j])\t $(SvN_loop[j])")
  end

  p = plot(t1,SvN_loop,seriestype=:scatter,show = true)
  savefig(p, "myplot.pdf")

end

The ground state found by the DMRG algorithm matches with my exact diagonalization code. There are a few things I’m getting wrong in this calculation.

  1. When I run this code inside the REPL environment vs. when I run it in cmd with "julia *.jl ", I see slightly different results in my final plot, particularly along the regions where the t1 approaches the value of t2. There is a topological phase transition at t1=t2 so I’m guessing the inclusion of more singular values in the truncation part causes slight variation. Is that what is happening?
    There are also abrupt changes in the entropy values along the transition regime. Might that also be a numerical artifact?

  2. The entanglement entropy behaviour is different from the analytical calculations I’ve done for the same system. I’m unsure if there is some syntax error or if there is some physical argument I’m missing.

TIA!

2

So entanglement is one of the last quantities to converge in a DMRG calculation – it is a very sensitive quantity, especially near a phase transition as you correctly guess. Therefore if you are getting wrong or inconsistent entanglement results, a first thing to try is to do more DMRG sweeps to see if your results change. Also you’re correct that entanglement can be sensitive to the level of truncation, so after you are more sure your state is converged for a given cutoff level, check also that your entanglement results are converged with choice of cutoff.

If you are getting different results each time you run the program (I doubt it’s the command line versus repl; more likely just that each run of the program is giving different results since you are using a random initialization of your MPS), then the thing to do is to converge your results more fully, again mainly by doing more DMRG sweeps.

In general, I’d recommend not to put a DMRG calculation inside a for loop over parameters. It’s better to focus on a single choice of Hamiltonian parameters and check carefully that your results are converged before doing a second or third calculation.

About it not matching your analytical calculations, the two possibilities I can think of are either insufficient convergence of the numerics (see above), or else possibly a bug in your code when computing entanglement.

One thing I would caution against in your code is putting functions inside other functions, because in Julia nested functions automatically “capture” variables from the enclosing scope. So it’s safer to put functions into global scope and pass all of their arguments to them.

3

Hi Miles,
Thanks for the help. You were right that the DMRG wasn’t converging. I increased the sweeps, and it seems to work, although, for certain parameters, the convergence isn’t as great as I want it to be. Even after 500 sweeps, there is still variation in the G.S. energy for certain t1 values that are not near the transition point. I have also defined the helper functions globally.
Now, coming to the entanglement entropy calculation, I followed the discussion in the ITensor documentation. I defined the following function to calculate it:

function ent(psi::MPS; kwargs...)
    N = length(psi)
    b = ÷(N,2)
    orthogonalize!(psi, b)
    s = siteinds(psi)
    U,S,V = 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

If I have understood it correctly, I define my orthogonality center at the middle of the chain. The orthogonalize! Function left normalizes the tensors to the left of b = 30, and right normalizes the other side. Then the SVD retrieves the singular values along the middle cut, which are used to construct the von Neumann entropy. I can’t see if there is a bug in my code for the entanglement calculation.
Thanks a lot for the pointers in your answer!

1 Like
4

Hi, glad to hear things are working better for you. Here’s two responses:

  1. if your DMRG is taking 500 sweeps to converge well, another parameter you can consider adjusting is the eigsolve_krylovdim parameter. To adjust it, just pass it as a keyword argument to the dmrg function. The default is 3, which can be optimal for easy systems with a large gap and short correlation length, but in some cases taking it to be larger can help a lot with convergence accuracy and speed. Try out some different values to see if it helps.

  2. about your ent function, it looks good to me, and seems to be based on the code example from the ITensor documentation, which is great. Did you have a specific question about it? (As a minor detail, no need to take kwargs... in the function definition unless you are going to pass the kwargs... to another function inside or use the ent function in a generic coding context where you might pass some unused keyword arguments.)

5

Hi Miles,
Thanks so much for that! The eigensolve_krylovdim parameter solved the convergence issue!
The entanglement entropy now matches with the correlation matrix calculations for small lattices, but for larger lattice sizes, say from N=32 onward, they do not match. Since edge modes are present in the SSH model, is it possible that for larger lattice sizes, their effect isn’t being considered by the dmrg algorithm?

6

Glad the eigsolve_krylovdim parameter helped!

I don’t think I fully understand the issue you’re seeing though. Could you explain a bit more of the steps you are doing? Is the correlation matrix calculation done for the exact same MPS that you did the entanglement entropy calculation for? Or is the correlation matrix computed separately from, say, an exact diagonalization or free fermion calculation?

About edge modes, they are always considered by DMRG but DMRG can often treat them differently from how other numerical methods treat them. It’s a long thing to explain in detail, but basically edge modes lead to degeneracies (multiple ground states) and DMRG will return one of these ground states while other methods such as ED usually returns specific linear combinations of two or more ground states. Neither result is technically “right” or “wrong” but just a different behavior of the numerical method or a different basis, basically, for representing the ground state manifold.

7

Apologies for the ambiguity, but I calculated the correlation matrix via free fermion calculation. This system not only has edge modes, but they have zero energy as well. So if there are multiple ground states, they could, in principle, have different entanglement entropies, and DMRG is picking out the one with less entanglement. That’s probably what I see in my calculations as well.
Is there any literature regarding DMRG and dealing with degenerate ground states?

8

Good question, but I don’t think there is all that much literature, at least not very pedagogical literature, describing best practices for DMRG on systems with zero-energy edge states. There is a general set of results (e.g. papers by Tarun Grover) about how DMRG tends to choose states with the least entanglement for systems with ground state degeneracies.

Basically, the rule of thumb here is that among all possible linear combinations of the ground states, DMRG will pick the least entangled one or one of the least entangled states. Then you can use “excited state DMRG” as well as techniques involving total quantum numbers to usually obtain the other ground states. In some systems, like the AKLT chain, there are also tricks involving special boundary conditions that can “quench” or remove the edge states.