von Neumann entanglement entropy of 1D hubbard model (Julia)

Question

Hi all,
I got a question when calculating von Neumann entanglement entropy of a simple 1D hubbard model... I use the sample code (https://github.com/ITensor/ITensors.jl/blob/main/examples/dmrg/1d_hubbard_extended.jl), keep only t1 and U and change sweep setting a little bit (attached below in code) because I have a longer chain (100sites). And I combined it with another sample code for the entropy calculation which I find here (http://itensor.org/support/2423/entanglement-entropy-julia).

I use a loop for that index b which labels a cut at bond b, so that I can plot this entropy as a function of the site number. I expect the order of magnitude of S(b) to be O(1), but the results I got are all very samll numbers ~10^(-6)... I'm confused...Did I do anything wrong in my code?

Thank you so much! !
(I attched my complete code and printed entropy at each b below.)

Code:
using ITensors

let
N = 100
Npart = 100
t1 = 1.0
t2 = 0.0
U = -5.0
V1 = 0.0

s = siteinds("Electron", N; conserve_qns=true)

ampo = OpSum()
for b in 1:(N - 1)
ampo += -t1, "Cdagup", b, "Cup", b + 1
ampo += -t1, "Cdagup", b + 1, "Cup", b
ampo += -t1, "Cdagdn", b, "Cdn", b + 1
ampo += -t1, "Cdagdn", b + 1, "Cdn", b
ampo += V1, "Ntot", b, "Ntot", b + 1
end
for b in 1:(N - 2)
ampo += -t2, "Cdagup", b, "Cup", b + 2
ampo += -t2, "Cdagup", b + 2, "Cup", b
ampo += -t2, "Cdagdn", b, "Cdn", b + 2
ampo += -t2, "Cdagdn", b + 2, "Cdn", b
end
for i in 1:N
ampo += U, "Nupdn", i
end
H = MPO(ampo, s)

sweeps = Sweeps(30)
setmaxdim!(sweeps, 500)
setcutoff!(sweeps, 1E-7)
@show sweeps

state = ["Emp" for n in 1:N]
p = Npart
for i in N:-1:1
if p > i
println("Doubly occupying site $i")
state[i] = "UpDn"
p -= 2
elseif p > 0
println("Singly occupying site $i")
state[i] = (isodd(i) ? "Up" : "Dn")
p -= 1
end
end
# Initialize wavefunction to be bond
# dimension 10 random MPS with number
# of particles the same as state
psi0 = randomMPS(s, state, 10)

# Check total number of particles:
@show flux(psi0)

# Start DMRG calculation:
energy, psi = dmrg(H, psi0, sweeps)

upd = fill(0.0, N)
dnd = fill(0.0, N)
for j in 1:N
orthogonalize!(psi, j)
psidagj = dag(prime(psi[j], "Site"))
upd[j] = scalar(psidagj * op(s, "Nup", j) * psi[j])
dnd[j] = scalar(psidag_j * op(s, "Ndn", j) * psi[j])
end

println("Up Density:")
for j in 1:N
println("$j $(upd[j])")
end
println()

println("Dn Density:")
for j in 1:N
println("$j $(dnd[j])")
end
println()

println("Total Density:")
for j in 1:N
println("$j $(upd[j]+dnd[j])")
end
println()

println("\nGround State Energy = $energy")

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

SvN = entropyvonneumann(psi, N)

end

---

The entropy S(b):
0.09070691669706957
0.001332649961200171
1.3686167960001409e-6
4.359976386344567e-6
6.733373654214966e-6
1.3372604646791153e-7
1.5001040623576317e-5
5.058686612647704e-7
2.5543410516658077e-5
1.0639655773856194e-6
3.829868192360102e-5
1.4027735830778272e-7
2.0115153047422772e-7
2.2972119655530066e-7
3.4283980977130424e-7
3.115574369307899e-7
4.871244534235853e-7
3.9846497809687774e-7
6.440428978804487e-7
4.885839791828351e-7
8.217896181088025e-7
5.822938107843419e-7
1.060520215982716e-6
6.831068268753088e-7
1.3465796889154753e-6
1.4128831887793594e-7
1.6578331240586298e-6
1.6817179418251726e-7
1.9086191014162828e-6
1.9218002525787778e-7
2.1334820606118912e-6
2.1459105376943014e-7
2.3313186179033887e-6
2.3442546255436683e-7
2.5091554302506395e-6
2.5148073806805365e-7
2.6629455105823636e-6
2.675461252705469e-7
2.7969288395425305e-6
2.820721836141702e-7
2.908469131486962e-6
2.940051897395952e-7
2.9878629561903115e-6
3.0171249585454165e-7
3.0612731469971742e-6
3.083795037091508e-7
3.100147849310539e-6
3.123873318518695e-7
3.1099061931736388e-6
3.123841279545293e-7
3.0987763981597106e-6
3.0976940496368236e-7
3.063501520406781e-6
3.0449799942869857e-7
3.001787667600418e-6
2.957587377413038e-7
2.918827460496549e-6
2.847481294234705e-7
2.808596119674664e-6
2.7068863487085386e-7
2.6786544449685913e-6
2.5486386260203346e-7
2.524835813245228e-6
2.3709999310016173e-7
2.3417164843292646e-6
2.154108041515205e-7
2.1471156146958798e-6
1.9392643033737067e-7
1.9277674229420978e-6
1.6959635573688854e-7
1.6764798541578795e-6
1.4214995738633714e-7
1.3460961878017196e-6
6.8199231820002e-7
1.0490130295538297e-6
5.814456742233713e-7
8.18106307245907e-7
4.878329051121316e-7
6.425305910731801e-7
3.9783346478151117e-7
4.842645529017947e-7
3.105588577603732e-7
3.2992276174983606e-7
2.099190336426272e-7
1.9297394199410213e-7
1.3973463812131638e-7
3.825571380443112e-5
1.059329165404299e-6
2.554198355682492e-5
4.920291374342437e-7
1.4994953866517857e-5
1.3062726493860295e-7
6.736631969987577e-6
4.323008567401327e-6
1.3668060657135116e-6
0.0013324848707008954
0.09070625418213862
0.3597747329567827

Answer 1

Hi, so I found the issue and it's due to a bug in the entropy measuring function you were using. The bug is that the for loop should be written like this:

    for n in 1:dim(S, 1)

instead of writing for n in dim(S,1) as in the original code. As I think you can see, in the original buggy version, dim(S,1) is just a number, whereas a for loop should go over a range.

The function in the other question you linked to had that bug plus a second bug, so I went and edited that question reply to fix both bugs. Here is a corrected version of your implementation of it below:

function entropyvonneumann(psi::MPS, N::Int)
  s = siteinds(psi)
  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 1:dim(S, 1)
      p = S[n,n]^2
      SvN[b] -= p * log(p)
    end
    println("$(SvN[b])")
  end
  return SvN
end

Comments

I see! Thank you so much! !

And just a quick check,  is "s" not defined in function the second bug you mentioned? That's why you add this line, "s = siteinds(psi)", in function?

---

Yes, exactly: that was the other bug. The function linked in the other question actually did run without erroring, but it only worked because the code below it was in global scope, so the "s" variable was a global variable. It goes to show that working in global scope is risky and can lead to a lot of bugs. So I usually prefer to work in a local scope (inside a function or a let..end block) and also define most of my functions outside of that scope. I had an issue recently where I didn't and had a similar bug in one of my own codes.

---

I see! Yes I also didn't get error message from this "s" so a little unsure about the difference. Thanks for explaining this! :)

---

Right, I think you got it but in your code above, the function "captures" the variable s from the outer scope which is why it works. But capturing variables can be a dangerous feature to use because it can lead to unexpected behaviors if variables have the same name by accident. So I usually define my functions outside of the scope where I define my variables, unless they are small "lambda" functions where the capturing behavior is intentional and very useful.