+1 vote
asked by (270 points)
edited by

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(psidag
j * 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

1 Answer

+1 vote
answered by (70.1k points)
selected by
 
Best answer

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
commented by (270 points)
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?
commented by (70.1k points)
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.
commented by (270 points)
I see! Yes I also didn't get error message from this "s" so a little unsure about the difference. Thanks for explaining this! :)
commented by (70.1k points)
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.
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.

Categories

...