|
| 1 | +using Test |
| 2 | +using Random |
| 3 | +using PEPSKit |
| 4 | +using TensorKit |
| 5 | +using KrylovKit |
| 6 | +using OptimKit |
| 7 | + |
| 8 | +using MPSKit: add_physical_charge |
| 9 | + |
| 10 | +# This example demonstrates the simulation of the two-dimensional Bose-Hubbard model using |
| 11 | +# PEPSKit.jl. In particular, it showcases the use of internal symmetries and finite particle |
| 12 | +# densities in PEPS simulations. |
| 13 | + |
| 14 | +## Defining the model |
| 15 | + |
| 16 | +# We will construct the Bose-Hubbard model Hamiltonian through the |
| 17 | +# [`bose_hubbard_model` function from MPSKitModels.jl](https://quantumkithub.github.io/MPSKitModels.jl/dev/man/models/#MPSKitModels.bose_hubbard_model), |
| 18 | +# as reexported by PEPSKit.jl. We'll simulate the model in its Mott-insulating phase |
| 19 | +# where the ratio U/t is large, since in this phase we expect the ground state to be well |
| 20 | +# approximated by a PEPS with a manifest global U(1) symmetry. Furthermore, we'll impose |
| 21 | +# a cutoff at 2 bosons per site, set the chemical potential to zero and use a simple 1x1 |
| 22 | +# unit cell. |
| 23 | + |
| 24 | +t = 1.0 |
| 25 | +U = 30.0 |
| 26 | +cutoff = 2 |
| 27 | +mu = 0.0 |
| 28 | +lattice = InfiniteSquare(1, 1) |
| 29 | + |
| 30 | +# We'll impose an explicit global U(1) symmetry as well as a fixed particle number density |
| 31 | +# in our simulations. We can do this by setting the `symmetry` keyword argument to `U1Irrep` |
| 32 | +# and passing one as the particle number density keyword argument `n`. |
| 33 | + |
| 34 | +symmetry = U1Irrep |
| 35 | +n = 1 |
| 36 | + |
| 37 | +# We can then construct the Hamiltonian, and inspect the corresponding lattice of physical |
| 38 | +# spaces. |
| 39 | + |
| 40 | +H = bose_hubbard_model(ComplexF64, symmetry, lattice; cutoff, t, U, n) |
| 41 | +Pspaces = H.lattice |
| 42 | + |
| 43 | +# Note that the physical space contains U(1) charges -1, 0 and +1. Indeed, imposing a |
| 44 | +# particle number density of +1 corresponds to shifting the physical charges by -1 to |
| 45 | +# 're-center' the physical charges around the desired density. When we do this with a cutoff |
| 46 | +# of two bosons per site, i.e. starting from U(1) charges 0, 1 and 2 on the physical level, |
| 47 | +# we indeed get the observed charges. |
| 48 | + |
| 49 | +## Characterizing the virtual spaces |
| 50 | + |
| 51 | +# When running PEPS simulations with explicit internal symmetries, specifying the structure |
| 52 | +# of the virtual spaces of the PEPS and its environment becomes a bit more involved. For the |
| 53 | +# environment, one could in principle allow the virtual space to be chosen dynamically |
| 54 | +# during the boundary contraction using CTMRG by using a truncation scheme that allows for |
| 55 | +# this (e.g. using alg=:truncdim or alg=:truncbelow to truncate to a fixed total bond |
| 56 | +# dimension or singular value cutoff respectively). For the PEPS virtual space however, the |
| 57 | +# structure has to be specified before the optimization. |
| 58 | + |
| 59 | +# While there are a host of techniques to do this in an informed way (e.g. starting from |
| 60 | +# a simple update result), here we just specify the virtual space manually. Since we're |
| 61 | +# dealing with a model at unit filling our physical space only contains integer U(1) irreps. |
| 62 | +# Therefore, we'll build our PEPS and environment spaces using integer U(1) irreps centered |
| 63 | +# around the zero charge. |
| 64 | + |
| 65 | +Vpeps = U1Space(0 => 2, 1 => 1, -1 => 1) |
| 66 | +Venv = U1Space(0 => 6, 1 => 4, -1 => 4, 2 => 2, -2 => 2) |
| 67 | + |
| 68 | +## Finding the ground state |
| 69 | + |
| 70 | +# Having defined our Hamiltonian and spaces, it is just a matter of pluggin this into the |
| 71 | +# optimization framework in the usual way to find the ground state. |
| 72 | + |
| 73 | +# specify algorithms and tolerances |
| 74 | +boundary_alg = (; tol=1e-8, alg=:simultaneous, verbosity=2, trscheme=(; alg=:fixedspace)) |
| 75 | +gradient_alg = (; tol=1e-6, maxiter=10, alg=:eigsolver, iterscheme=:diffgauge) |
| 76 | +optimizer_alg = (; tol=1e-4, alg=:lbfgs, verbosity=3, maxiter=200, ls_maxiter=2, ls_maxfg=2) |
| 77 | +reuse_env = true |
| 78 | + |
| 79 | +# initialize state |
| 80 | +Nspaces = fill(Vpeps, size(lattice)...) |
| 81 | +Espaces = fill(Vpeps, size(lattice)...) |
| 82 | +Random.seed!(2928528935) # for reproducibility |
| 83 | +ψ₀ = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) |
| 84 | +env₀ = CTMRGEnv(ψ₀, Venv) |
| 85 | +env₀, = leading_boundary(env₀, ψ₀; boundary_alg...) |
| 86 | + |
| 87 | +# optimize |
| 88 | +ψ, env, E, info = fixedpoint( |
| 89 | + H, ψ₀, env₀; boundary_alg, gradient_alg, optimizer_alg, reuse_env |
| 90 | +) |
| 91 | + |
| 92 | +## Check the result |
| 93 | + |
| 94 | +# We can compare our PEPS result to the energy obtained using a cylinder-MPS calculation |
| 95 | +# using a cylinder circumference of Ly = 7 and a bond dimension of 446, which yields |
| 96 | +# E = -0.273284888 |
| 97 | + |
| 98 | +E_ref = -0.273284888 |
| 99 | +@test E ≈ E_ref rtol = 1e-3 |
0 commit comments