Add entanglement entropy to examples (#80)

Author: ryanlevy

examples/vumps/src/entropy.jl

@@ -0,0 +1,33 @@
+using ITensors
+using ITensorInfiniteMPS
+
+# calculate -\sum_i rho_i log(rho_i)
+function entropy(S::ITensor)
+  @assert order(S) == 2
+  SvN, norm = 0.0, 0.0
+  for n in 1:dim(S, 1)
+    p = S[n, n]^2
+    SvN -= p * log(p)
+    norm += p
+  end
+  return SvN, norm
+end
+
+# calculate von Neumann entropy of a MPS cut between sites b-1 and b
+function entropy(ψ::MPS, b)
+  ψ = orthogonalize(ψ, b)
+  U, S, V = svd(ψ[b], (linkind(ψ, b - 1), siteind(ψ, b)))
+  SvN, norm = entropy(S)
+  @assert norm ≈ 1.0
+  return SvN
+end
+
+# calculate von Neumann entropy of an infinite MPS at cut b of the unit cell
+function entropy(ψ::InfiniteCanonicalMPS, b)
+  #calculate entropy
+  C = ψ.C[b]
+  Ũ, S, Ṽ = svd(C, inds(C)[1])
+  SvN, norm = entropy(S)
+  @assert norm ≈ 1.0
+  return SvN
+end

examples/vumps/vumps_hubbard_extended.jl

@@ -1,11 +1,11 @@
 using ITensors
 using ITensorInfiniteMPS
 
-include(
-  joinpath(
-    pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl"
-  ),
-)
+base_path = joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src")
+src_files = ["vumps_subspace_expansion.jl", "entropy.jl"]
+for f in src_files
+  include(joinpath(base_path, f))
+end
 
 ##############################################################################
 # VUMPS parameters
@@ -105,7 +105,7 @@ maxdims =
 ## setcutoff!(sweeps, cutoff)
 
 println("\nRun DMRG on $Nfinite sites")
-energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; nsweeps, maxdims, cutoff)
+energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite; nsweeps, maxdim=maxdims, cutoff)
 println("\nEnergy density")
 @show energy_finite_total / Nfinite
 
@@ -125,6 +125,9 @@ corr_finite = correlation_matrix(
   ψfinite, "Cdagup", "Cup"; sites=Int(Nfinite / 2):Int(Nfinite / 2 + 9)
 )
 
+S_finite = [entropy(ψfinite, b) for b in (Nfinite ÷ 2):(Nfinite ÷ 2 + N - 1)]
+S_infinite = [entropy(ψ, b) for b in 1:N]
+
 println("\nResults from VUMPS")
 @show energy_infinite
 @show energy_exact
@@ -133,6 +136,7 @@ println("\nResults from VUMPS")
 @show Nup .+ Ndn
 @show Sz
 @show corr_infinite
+@show S_infinite
 
 println("\nResults from DMRG")
 @show energy_finite
@@ -141,5 +145,6 @@ println("\nResults from DMRG")
 @show Nup_finite .+ Ndn_finite
 @show Sz_finite
 @show corr_finite
+@show S_finite
 
 nothing

examples/vumps/vumps_ising.jl

@@ -1,11 +1,11 @@
 using ITensors
 using ITensorInfiniteMPS
 
-include(
-  joinpath(
-    pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl"
-  ),
-)
+base_path = joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src")
+src_files = ["vumps_subspace_expansion.jl", "entropy.jl"]
+for f in src_files
+  include(joinpath(base_path, f))
+end
 
 ##############################################################################
 # VUMPS/TDVP parameters
@@ -102,3 +102,10 @@ Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz")
 
 @show Sz1_finite, Sz2_finite
 @show Sz1_infinite, Sz2_infinite
+
+S_finite = [entropy(ψ_finite, b) for b in (nsite_finite ÷ 2):(nsite_finite ÷ 2 + nsite - 1)]
+S_infinite = [entropy(ψ, b) for b in 1:nsite]
+@show S_finite
+@show S_infinite
+
+nothing

test/test_vumps.jl

@@ -36,7 +36,7 @@ end
   ψfinite = randomMPS(sfinite, initstate)
   nsweeps = 20
   energy_finite_total, ψfinite = dmrg(
-    Hfinite, ψfinite; nsweeps, maxdims=10, cutoff=1e-10, outputlevel=0
+    Hfinite, ψfinite; nsweeps, maxdim=10, cutoff=1e-10, outputlevel=0
   )
 
   @testset "VUMPS/TDVP with: multisite_update_alg = $multisite_update_alg, conserve_qns = $conserve_qns, nsites = $nsites, time_step = $time_step, localham_type = $localham_type" for multisite_update_alg in
@@ -142,7 +142,7 @@ end
   ψfinite = randomMPS(sfinite, initstate)
   nsweeps = 20
   energy_finite_total, ψfinite = dmrg(
-    Hfinite, ψfinite; nsweeps, maxdims=10, cutoff=1e-10, outputlevel=0
+    Hfinite, ψfinite; nsweeps, maxdim=10, cutoff=1e-10, outputlevel=0
   )
 
   multisite_update_alg = "sequential"

test/test_vumpsmpo.jl

@@ -105,7 +105,7 @@ end
   ψfinite = randomMPS(sfinite, initstate)
   nsweeps = 20
   energy_finite_total, ψfinite = dmrg(
-    Hfinite, ψfinite; nsweeps, maxdims=10, cutoff=1e-10, outputlevel=0
+    Hfinite, ψfinite; nsweeps, maxdim=10, cutoff=1e-10, outputlevel=0
   )
   Szs_finite = expect(ψfinite, "Sz")