|
12 | 12 | "using QuadGK: quadgk\n", |
13 | 13 | "using SpecialFunctions: ellipe\n", |
14 | 14 | "using Plots\n", |
15 | | - "using LinearAlgebra" |
| 15 | + "using LinearAlgebra\n", |
| 16 | + "using BenchmarkFreeFermions" |
16 | 17 | ], |
17 | 18 | "metadata": {}, |
18 | 19 | "execution_count": null |
|
25 | 26 | "This example shows how to simulate the finite temperature behavior of the XY model in 1D.\n", |
26 | 27 | "Importantly, the Hamiltonian can be diagonalized in terms of fermionic creation and annihilation operators.\n", |
27 | 28 | "As a result, many properties have analytical expressions that can be used to verify our results.\n", |
| 29 | + "Here, we use [BenchmarkFreeFermions.jl](https://github.com/Qiaoyi-Li/BenchmarkFreeFermions.jl/) to compare our results.\n", |
28 | 30 | "\n", |
29 | 31 | "$$\n", |
30 | 32 | " H = J \\sum_{i=1}^{N} \\left( \\sigma^x_i \\sigma^x_{i+1} + \\sigma^y_i \\sigma^y_{i+1} \\right)\n", |
|
46 | 48 | "cell_type": "code", |
47 | 49 | "source": [ |
48 | 50 | "J = 1 / 2\n", |
49 | | - "N = 30\n", |
| 51 | + "T = ComplexF64\n", |
| 52 | + "symmetry = U1Irrep\n", |
50 | 53 | "\n", |
51 | 54 | "function XY_hamiltonian(::Type{T}=ComplexF64, ::Type{S}=Trivial; J=1 / 2,\n", |
52 | 55 | " N) where {T<:Number,S<:Sector}\n", |
|
85 | 88 | "outputs": [], |
86 | 89 | "cell_type": "code", |
87 | 90 | "source": [ |
88 | | - "single_particle_energy(k, J, N) = J * cos(k * 2π / (N + 0))\n", |
89 | 91 | "function groundstate_energy(J, N)\n", |
90 | | - " return isfinite(N) ?\n", |
91 | | - " -sum(n -> abs(single_particle_energy(n, J, N)), 1:N) / 2N : -J / π\n", |
| 92 | + " isfinite(N) || return -J / π\n", |
| 93 | + " T = diagm(1 => J / 2 * ones(N - 1), -1 => J / 2 * ones(N - 1))\n", |
| 94 | + " ϵ = SingleParticleSpectrum(T)\n", |
| 95 | + " return Energy(ϵ, Inf, 0) / N\n", |
92 | 96 | "end" |
93 | 97 | ], |
94 | 98 | "metadata": {}, |
|
107 | 111 | "outputs": [], |
108 | 112 | "cell_type": "code", |
109 | 113 | "source": [ |
110 | | - "H = periodic_boundary_conditions(XY_hamiltonian(; J, N=Inf), 6)\n", |
| 114 | + "N_exact = 6\n", |
| 115 | + "H = open_boundary_conditions(XY_hamiltonian(T, symmetry; J, N=Inf), N_exact)\n", |
111 | 116 | "H_dense = convert(TensorMap, H);\n", |
112 | | - "vals = eigvals(H_dense)[Trivial()] ./ 6\n", |
113 | | - "groundstate_energy(J, 6)\n", |
| 117 | + "vals = eigvals(H_dense)[one(symmetry)] ./ N_exact\n", |
| 118 | + "groundstate_energy(J, N_exact)\n", |
114 | 119 | "\n", |
115 | | - "println(\"Exact (N=6):\\t\", groundstate_energy(J, 6))\n", |
116 | | - "println(\"Exact (N=Inf):\\t\", groundstate_energy(J, Inf))\n", |
117 | | - "println(\"Numerical:\\t\", minimum(real(vals)))" |
| 120 | + "println(\"Numerical:\\t\", minimum(real(vals)))\n", |
| 121 | + "println(\"Exact (N=$(N_exact)):\\t\", groundstate_energy(J, N_exact))\n", |
| 122 | + "println(\"Exact (N=Inf):\\t\", groundstate_energy(J, Inf))" |
118 | 123 | ], |
119 | 124 | "metadata": {}, |
120 | 125 | "execution_count": null |
|
132 | 137 | "outputs": [], |
133 | 138 | "cell_type": "code", |
134 | 139 | "source": [ |
135 | | - "H = XY_hamiltonian(; J, N)\n", |
| 140 | + "N = 32\n", |
| 141 | + "H = XY_hamiltonian(T, symmetry; J, N)\n", |
136 | 142 | "D = 64\n", |
137 | | - "psi_init = FiniteMPS(N, physicalspace(H, 1), ℂ^D)\n", |
138 | | - "psi, envs, = find_groundstate(psi_init, H, DMRG(; maxiter=10));\n", |
| 143 | + "V_init = symmetry === Trivial ? ℂ^32 : U1Space(i => 10 for i in -1:(1 // 2):1)\n", |
| 144 | + "psi_init = FiniteMPS(N, physicalspace(H, 1), V_init)\n", |
| 145 | + "trscheme = truncdim(D)\n", |
| 146 | + "psi, envs, = find_groundstate(psi_init, H, DMRG2(; trscheme, maxiter=5));\n", |
139 | 147 | "E_0 = expectation_value(psi, H, envs) / N\n", |
140 | 148 | "\n", |
| 149 | + "println(\"Numerical:\\t\", real(E_0))\n", |
141 | 150 | "println(\"Exact (N=$N):\\t\", groundstate_energy(J, N))\n", |
142 | | - "println(\"Exact (N=Inf):\\t\", groundstate_energy(J, Inf))\n", |
143 | | - "println(\"Numerical:\\t\", real(E_0))" |
| 151 | + "println(\"Exact (N=Inf):\\t\", groundstate_energy(J, Inf))" |
144 | 152 | ], |
145 | 153 | "metadata": {}, |
146 | 154 | "execution_count": null |
|
197 | 205 | "cell_type": "code", |
198 | 206 | "source": [ |
199 | 207 | "function partition_function(β::Number, J::Number, N::Number)\n", |
200 | | - " return prod(k -> (1 + exp(-β * single_particle_energy(k, J, N))), 1:N)^(1 / N)\n", |
| 208 | + " T = diagm(1 => J / 2 * ones(N - 1), -1 => J / 2 * ones(N - 1))\n", |
| 209 | + " ϵ = SingleParticleSpectrum(T)\n", |
| 210 | + " return LogPartition(ϵ, β, 0) / N\n", |
201 | 211 | "end\n", |
202 | 212 | "function free_energy(β, J, N)\n", |
203 | | - " return -1 / β * log(partition_function(β, J, N))\n", |
204 | | - "end" |
| 213 | + " T = diagm(1 => J / 2 * ones(N - 1), -1 => J / 2 * ones(N - 1))\n", |
| 214 | + " ϵ = SingleParticleSpectrum(T)\n", |
| 215 | + " return FreeEnergy(ϵ, β, 0) / N\n", |
| 216 | + "end\n", |
| 217 | + "\n", |
| 218 | + "βs = 0.0:0.2:8.0\n", |
| 219 | + "\n", |
| 220 | + "Z_analytic = partition_function.(βs, J, N);\n", |
| 221 | + "F_analytic = free_energy.(βs, J, N);" |
205 | 222 | ], |
206 | 223 | "metadata": {}, |
207 | 224 | "execution_count": null |
|
224 | 241 | "outputs": [], |
225 | 242 | "cell_type": "code", |
226 | 243 | "source": [ |
227 | | - "βs = 0.0:0.2:8.0\n", |
228 | 244 | "expansion_orders = 1:3\n", |
229 | 245 | "\n", |
230 | | - "function partition_function_taylor(β, H; expansion_order)\n", |
| 246 | + "function logpartition_taylor(β, H; expansion_order)\n", |
231 | 247 | " dτ = im * β\n", |
232 | 248 | " expH = make_time_mpo(H, dτ,\n", |
233 | 249 | " TaylorCluster(; N=expansion_order))\n", |
234 | | - " return real(tr(expH))^(1 / N)\n", |
| 250 | + " return log(real(tr(expH))) / length(H)\n", |
235 | 251 | "end\n", |
236 | 252 | "\n", |
237 | 253 | "Z_taylor = map(Iterators.product(βs, expansion_orders)) do (β, expansion_order)\n", |
238 | 254 | " @info \"Computing β = $β at order $expansion_order\"\n", |
239 | | - " return partition_function_taylor(β, H; expansion_order)\n", |
| 255 | + " return logpartition_taylor(β, H; expansion_order)\n", |
240 | 256 | "end\n", |
241 | | - "F_taylor = -(1 ./ βs) .* log.(Z_taylor)\n", |
| 257 | + "F_taylor = -(1 ./ βs) .* Z_taylor\n", |
242 | 258 | "\n", |
243 | 259 | "p_taylor = let\n", |
244 | 260 | " labels = reshape(map(expansion_orders) do N\n", |
245 | 261 | " return \"Taylor N=$N\"\n", |
246 | 262 | " end, 1, :)\n", |
247 | | - " p1 = plot(βs, partition_function.(βs, J, N); label=\"analytic\",\n", |
| 263 | + " p1 = plot(βs, Z_analytic; label=\"analytic\",\n", |
248 | 264 | " title=\"Partition function\",\n", |
249 | 265 | " xlabel=\"β\", ylabel=\"Z(β)\")\n", |
250 | | - " plot!(p1, βs, real.(Z_taylor); label=labels)\n", |
251 | | - " p2 = plot(βs, free_energy.(βs, J, N); label=\"analytic\", title=\"Free energy\",\n", |
| 266 | + " plot!(p1, βs, Z_taylor; label=labels)\n", |
| 267 | + " p2 = plot(βs, F_analytic; label=\"analytic\", title=\"Free energy\",\n", |
252 | 268 | " xlabel=\"β\", ylabel=\"F(β)\")\n", |
253 | | - " plot!(p2, βs, real.(F_taylor); label=labels)\n", |
| 269 | + " plot!(p2, βs, F_taylor; label=labels)\n", |
254 | 270 | " plot(p1, p2)\n", |
255 | 271 | "end" |
256 | 272 | ], |
|
298 | 314 | "outputs": [], |
299 | 315 | "cell_type": "code", |
300 | 316 | "source": [ |
| 317 | + "expansion_orders = 2:3\n", |
| 318 | + "Z_taylor = Z_taylor[:, 2:end]\n", |
| 319 | + "F_taylor = F_taylor[:, 2:end]\n", |
| 320 | + "\n", |
301 | 321 | "p_taylor_diff = let\n", |
302 | | - " labels = reshape(map(expansion_orders[2:end]) do N\n", |
| 322 | + " labels = reshape(map(expansion_orders) do N\n", |
303 | 323 | " return \"Taylor N=$N\"\n", |
304 | 324 | " end, 1, :)\n", |
305 | | - " p1 = plot(βs, abs.(real.(Z_taylor[:, 2:end]) .- partition_function.(βs, J, N));\n", |
| 325 | + " p1 = plot(βs, abs.(Z_taylor .- Z_analytic);\n", |
306 | 326 | " label=labels, title=\"Partition function error\",\n", |
307 | 327 | " xlabel=\"β\", ylabel=\"ΔZ(β)\", legend=:topleft)\n", |
308 | | - " p2 = plot(βs, abs.(real.(F_taylor[:, 2:end]) .- free_energy.(βs, J, N)); label=labels,\n", |
| 328 | + " p2 = plot(βs, abs.(F_taylor .- F_analytic); label=labels,\n", |
309 | 329 | " xlabel=\"β\", ylabel=\"ΔF(β)\", title=\"Free energy error\", legend=:topleft)\n", |
310 | 330 | " plot(p1, p2)\n", |
311 | 331 | "end" |
|
341 | 361 | "outputs": [], |
342 | 362 | "cell_type": "code", |
343 | 363 | "source": [ |
344 | | - "function partition_function_taylor2(β, H; expansion_order)\n", |
| 364 | + "double_logpartition(ρ₁, ρ₂=ρ₁) = log(real(dot(ρ₁, ρ₂))) / length(ρ₁)\n", |
| 365 | + "\n", |
| 366 | + "function logpartition_taylor2(β, H; expansion_order)\n", |
345 | 367 | " dτ = im * β / 2\n", |
346 | 368 | " expH = make_time_mpo(H, dτ, TaylorCluster(; N=expansion_order))\n", |
347 | | - " return real(dot(expH, expH))^(1 / N)\n", |
| 369 | + " return double_logpartition(expH)\n", |
348 | 370 | "end\n", |
349 | 371 | "\n", |
350 | | - "Z_taylor2 = map(Iterators.product(βs, expansion_orders[2:end])) do (β, expansion_order)\n", |
| 372 | + "Z_taylor2 = map(Iterators.product(βs, expansion_orders)) do (β, expansion_order)\n", |
351 | 373 | " @info \"Computing β = $β at order $expansion_order\"\n", |
352 | | - " return partition_function_taylor2(β, H; expansion_order)\n", |
| 374 | + " return logpartition_taylor2(β, H; expansion_order)\n", |
353 | 375 | "end\n", |
354 | | - "F_taylor2 = -(1 ./ βs) .* log.(Z_taylor2)\n", |
| 376 | + "F_taylor2 = -(1 ./ βs) .* Z_taylor2\n", |
355 | 377 | "\n", |
356 | 378 | "p_taylor2_diff = let\n", |
357 | 379 | " labels = reshape(map(expansion_orders[2:end]) do N\n", |
358 | 380 | " return \"Taylor N=$N\"\n", |
359 | 381 | " end, 1, :)\n", |
360 | | - " p1 = plot(βs, abs.(real.(Z_taylor2) .- partition_function.(βs, J, N));\n", |
| 382 | + " p1 = plot(βs, abs.(Z_taylor2 .- Z_analytic);\n", |
361 | 383 | " label=labels, title=\"Partition function error\",\n", |
362 | 384 | " xlabel=\"β\", ylabel=\"ΔZ(β)\", legend=:topleft)\n", |
363 | | - " p2 = plot(βs, abs.(real.(F_taylor2) .- free_energy.(βs, J, N)); label=labels,\n", |
| 385 | + " p2 = plot(βs, abs.(F_taylor2 .- F_analytic); label=labels,\n", |
364 | 386 | " xlabel=\"β\", ylabel=\"ΔF(β)\", title=\"Free energy error\", legend=:topleft)\n", |
365 | 387 | " plot(p1, p2)\n", |
366 | 388 | "end" |
|
389 | 411 | "To achieve this, we can reinterpret the density matrix as an MPS with two physical indices.\n", |
390 | 412 | "Then, we have some control over the approximations we make by tuning the maximal bond dimension.\n", |
391 | 413 | "\n", |
| 414 | + "Here, we swich to a logarithmic scale for the errors to better illustrate the results.\n", |
| 415 | + "\n", |
392 | 416 | "> **Warning**\n", |
393 | 417 | ">\n", |
394 | 418 | "> Using MPS techniques to approximate the multiplication of density matrices does not necessarily inherit all of the nice properties of approximating MPS.\n", |
|
408 | 432 | "# first iteration: start from high order Taylor expansion\n", |
409 | 433 | "ρ₀ = make_time_mpo(H, im * βs[2] / 2, TaylorCluster(; N=3))\n", |
410 | 434 | "Z_mpo_mul[1] = Z_taylor[1]\n", |
411 | | - "Z_mpo_mul[2] = real(dot(ρ₀, ρ₀))^(1 / N)\n", |
| 435 | + "Z_mpo_mul[2] = double_logpartition(ρ₀)\n", |
412 | 436 | "\n", |
413 | 437 | "# subsequent iterations: multiply by ρ₀\n", |
414 | 438 | "ρ_mps = convert(FiniteMPS, ρ₀)\n", |
|
417 | 441 | " @info \"Computing β = $(βs[i])\"\n", |
418 | 442 | " ρ_mps, = approximate(ρ_mps, (ρ₀, ρ_mps),\n", |
419 | 443 | " DMRG2(; trscheme=truncdim(D_max), maxiter=10))\n", |
420 | | - " Z_mpo_mul[i] = real(dot(ρ_mps, ρ_mps))^(1 / N)\n", |
| 444 | + " Z_mpo_mul[i] = double_logpartition(ρ_mps)\n", |
421 | 445 | "end\n", |
422 | | - "F_mpo_mul = -(1 ./ βs) .* log.(Z_mpo_mul)\n", |
| 446 | + "F_mpo_mul = -(1 ./ βs) .* Z_mpo_mul\n", |
423 | 447 | "\n", |
424 | 448 | "p_mpo_mul_diff = let\n", |
425 | | - " labels = reshape(map(expansion_orders[2:end]) do N\n", |
| 449 | + " labels = reshape(map(expansion_orders) do N\n", |
426 | 450 | " return \"Taylor N=$N\"\n", |
427 | 451 | " end, 1, :)\n", |
428 | | - " p1 = plot(βs, abs.(real.(Z_taylor2) .- partition_function.(βs, J, N));\n", |
| 452 | + " p1 = plot(βs, abs.(Z_taylor2 .- Z_analytic);\n", |
429 | 453 | " label=labels, title=\"Partition function error\",\n", |
430 | | - " xlabel=\"β\", ylabel=\"ΔZ(β)\", legend=:topleft)\n", |
431 | | - " plot!(p1, βs, abs.(real.(Z_mpo_mul) .- partition_function.(βs, J, N));\n", |
| 454 | + " xlabel=\"β\", ylabel=\"ΔZ(β)\", legend=:bottomright, yscale=:log10)\n", |
| 455 | + " plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic);\n", |
432 | 456 | " label=\"MPO multiplication\")\n", |
433 | | - " p2 = plot(βs, abs.(real.(F_taylor2) .- free_energy.(βs, J, N)); label=labels,\n", |
434 | | - " xlabel=\"β\", ylabel=\"ΔF(β)\", title=\"Free energy error\", legend=:topleft)\n", |
435 | | - " plot!(p2, βs, abs.(real.(F_mpo_mul) .- free_energy.(βs, J, N));\n", |
| 457 | + " p2 = plot(βs, abs.(F_taylor2 .- F_analytic); label=labels,\n", |
| 458 | + " xlabel=\"β\", ylabel=\"ΔF(β)\", title=\"Free energy error\", legend=nothing,\n", |
| 459 | + " yscale=:log10)\n", |
| 460 | + " plot!(p2, βs, abs.(F_mpo_mul .- F_analytic);\n", |
436 | 461 | " label=\"MPO multiplication\")\n", |
437 | 462 | " plot(p1, p2)\n", |
438 | 463 | "end" |
|
484 | 509 | "cell_type": "code", |
485 | 510 | "source": [ |
486 | 511 | "βs_exp = 2.0 .^ (-3:3)\n", |
| 512 | + "Z_analytic_exp = partition_function.(βs_exp, J, N)\n", |
| 513 | + "F_analytic_exp = free_energy.(βs_exp, J, N)\n", |
| 514 | + "\n", |
487 | 515 | "Z_mpo_mul_exp = zeros(length(βs_exp))\n", |
488 | 516 | "\n", |
489 | 517 | "# first iteration: start from high order Taylor expansion\n", |
490 | 518 | "ρ₀ = make_time_mpo(H, im * first(βs_exp) / 2, TaylorCluster(; N=3))\n", |
491 | | - "Z_mpo_mul_exp[1] = real(dot(ρ₀, ρ₀))^(1 / N)\n", |
| 519 | + "Z_mpo_mul_exp[1] = double_logpartition(ρ₀)\n", |
492 | 520 | "\n", |
493 | 521 | "# subsequent iterations: square\n", |
494 | 522 | "ρ = ρ₀\n", |
|
498 | 526 | " @info \"Computing β = $(βs_exp[i])\"\n", |
499 | 527 | " ρ_mps, = approximate(ρ_mps, (ρ, ρ_mps),\n", |
500 | 528 | " DMRG2(; trscheme=truncdim(D_max), maxiter=10))\n", |
501 | | - " Z_mpo_mul_exp[i] = real(dot(ρ_mps, ρ_mps))^(1 / N)\n", |
| 529 | + " Z_mpo_mul_exp[i] = double_logpartition(ρ_mps)\n", |
502 | 530 | " ρ = convert(FiniteMPO, ρ_mps)\n", |
503 | 531 | "end\n", |
504 | | - "F_mpo_mul_exp = -(1 ./ βs_exp) .* log.(Z_mpo_mul_exp)\n", |
| 532 | + "F_mpo_mul_exp = -(1 ./ βs_exp) .* Z_mpo_mul_exp\n", |
505 | 533 | "\n", |
506 | 534 | "p_mpo_mul_exp_diff = let\n", |
507 | 535 | " labels = reshape(map(expansion_orders[2:end]) do N\n", |
508 | 536 | " return \"Taylor N=$N\"\n", |
509 | 537 | " end, 1, :)\n", |
510 | | - " p1 = plot(βs, abs.(real.(Z_taylor2) .- partition_function.(βs, J, N));\n", |
511 | | - " label=labels, title=\"Partition function error\",\n", |
512 | | - " xlabel=\"β\", ylabel=\"ΔZ(β)\", legend=:topleft)\n", |
513 | | - " plot!(p1, βs, abs.(real.(Z_mpo_mul) .- partition_function.(βs, J, N));\n", |
514 | | - " label=\"MPO multiplication\")\n", |
515 | | - " plot!(p1, βs_exp, abs.(real.(Z_mpo_mul_exp) .- partition_function.(βs_exp, J, N));\n", |
516 | | - " label=\"MPO multiplication exp\")\n", |
517 | | - "\n", |
518 | | - " p2 = plot(βs, abs.(real.(F_taylor2) .- free_energy.(βs, J, N)); label=labels,\n", |
519 | | - " xlabel=\"β\", ylabel=\"ΔF(β)\", title=\"Free energy error\", legend=:topleft)\n", |
520 | | - " plot!(p2, βs, abs.(real.(F_mpo_mul) .- free_energy.(βs, J, N));\n", |
521 | | - " label=\"MPO multiplication\")\n", |
522 | | - " plot!(p2, βs_exp, abs.(real.(F_mpo_mul_exp) .- free_energy.(βs_exp, J, N));\n", |
523 | | - " label=\"MPO multiplication exp\")\n", |
| 538 | + " p1 = plot(βs, abs.(Z_taylor2 .- Z_analytic);\n", |
| 539 | + " label=labels, title=\"Partition function error\", xlabel=\"β\", ylabel=\"ΔZ(β)\",\n", |
| 540 | + " legend=:bottomright, yscale=:log10)\n", |
| 541 | + " plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic); label=\"MPO multiplication\")\n", |
| 542 | + " plot!(p1, βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp); label=\"MPO multiplication exp\")\n", |
| 543 | + "\n", |
| 544 | + " p2 = plot(βs, abs.(F_taylor2 .- F_analytic); label=labels, xlabel=\"β\", ylabel=\"ΔF(β)\",\n", |
| 545 | + " title=\"Free energy error\", legend=nothing, yscale=:log10)\n", |
| 546 | + " plot!(p2, βs, abs.(F_mpo_mul .- F_analytic); label=\"MPO multiplication\")\n", |
| 547 | + " plot!(p2, βs_exp, abs.(F_mpo_mul_exp .- F_analytic_exp); label=\"MPO multiplication exp\")\n", |
524 | 548 | " plot(p1, p2)\n", |
525 | 549 | "end" |
526 | 550 | ], |
|
567 | 591 | "\n", |
568 | 592 | "# first iteration: start from infinite temperature state\n", |
569 | 593 | "ρ₀ = infinite_temperature_density_matrix(H)\n", |
570 | | - "Z_tdvp[1] = real(dot(ρ₀, ρ₀))^(1 / N)\n", |
| 594 | + "Z_tdvp[1] = double_logpartition(ρ₀)\n", |
571 | 595 | "\n", |
572 | 596 | "# subsequent iterations: evolve by H\n", |
573 | 597 | "ρ_mps = convert(FiniteMPS, ρ₀)\n", |
574 | 598 | "for i in 2:length(βs)\n", |
575 | 599 | " global ρ_mps\n", |
576 | 600 | " @info \"Computing β = $(βs[i])\"\n", |
577 | 601 | " ρ_mps, = timestep(ρ_mps, H, βs[i - 1] / 2, -im * (βs[i] - βs[i - 1]) / 2,\n", |
578 | | - " TDVP2(; trscheme=truncdim(D_max)))\n", |
579 | | - " Z_tdvp[i] = real(dot(ρ_mps, ρ_mps))^(1 / N)\n", |
| 602 | + " TDVP2(; trscheme=truncdim(64)))\n", |
| 603 | + " Z_tdvp[i] = double_logpartition(ρ_mps)\n", |
580 | 604 | "end\n", |
581 | | - "F_tdvp = -(1 ./ βs) .* log.(Z_tdvp)\n", |
| 605 | + "F_tdvp = -(1 ./ βs) .* Z_tdvp\n", |
582 | 606 | "\n", |
583 | 607 | "p_mpo_mul_diff = let\n", |
584 | | - " labels = reshape(map(expansion_orders[2:end]) do N\n", |
| 608 | + " labels = reshape(map(expansion_orders) do N\n", |
585 | 609 | " return \"Taylor N=$N\"\n", |
586 | 610 | " end, 1, :)\n", |
587 | | - " p1 = plot(βs, abs.(real.(Z_taylor2) .- partition_function.(βs, J, N));\n", |
588 | | - " label=labels, title=\"Partition function error\",\n", |
589 | | - " xlabel=\"β\", ylabel=\"ΔZ(β)\", legend=:topleft)\n", |
590 | | - " plot!(p1, βs, abs.(real.(Z_mpo_mul) .- partition_function.(βs, J, N));\n", |
591 | | - " label=\"MPO multiplication\")\n", |
592 | | - " plot!(p1, βs, abs.(real.(Z_tdvp) .- partition_function.(βs, J, N));\n", |
593 | | - " label=\"TDVP\")\n", |
594 | | - "\n", |
595 | | - " p2 = plot(βs, abs.(real.(F_taylor2) .- free_energy.(βs, J, N)); label=labels,\n", |
596 | | - " xlabel=\"β\", ylabel=\"ΔF(β)\", title=\"Free energy error\", legend=:topleft)\n", |
597 | | - " plot!(p2, βs, abs.(real.(F_mpo_mul) .- free_energy.(βs, J, N));\n", |
598 | | - " label=\"MPO multiplication\")\n", |
599 | | - " plot!(p2, βs, abs.(real.(F_tdvp) .- free_energy.(βs, J, N));\n", |
600 | | - " label=\"TDVP\")\n", |
| 611 | + " p1 = plot(βs, abs.(Z_taylor2 .- Z_analytic); label=labels,\n", |
| 612 | + " title=\"Partition function error\", xlabel=\"β\", ylabel=\"ΔZ(β)\",\n", |
| 613 | + " legend=:bottomright, yscale=:log10)\n", |
| 614 | + " plot!(p1, βs, abs.(Z_mpo_mul .- Z_analytic); label=\"MPO multiplication\")\n", |
| 615 | + " plot!(p1, βs_exp, abs.(Z_mpo_mul_exp .- Z_analytic_exp); label=\"MPO multiplication exp\")\n", |
| 616 | + " plot!(p1, βs, abs.(Z_tdvp .- Z_analytic); label=\"TDVP\")\n", |
| 617 | + "\n", |
| 618 | + " p2 = plot(βs, abs.(F_taylor2 .- F_analytic); label=labels, xlabel=\"β\", ylabel=\"ΔF(β)\",\n", |
| 619 | + " title=\"Free energy error\", legend=nothing, yscale=:log10)\n", |
| 620 | + " plot!(p2, βs, abs.(F_mpo_mul .- F_analytic); label=\"MPO multiplication\")\n", |
| 621 | + " plot!(p2, βs_exp, abs.(F_mpo_mul_exp .- F_analytic_exp); label=\"MPO multiplication exp\")\n", |
| 622 | + " plot!(p2, βs, abs.(F_tdvp .- F_analytic); label=\"TDVP\")\n", |
601 | 623 | "\n", |
602 | 624 | " plot(p1, p2)\n", |
603 | 625 | "end" |
|
0 commit comments