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)

psidag*j = 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 entropy*von*neumann(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 = entropy*von*neumann(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