|
1 | 1 | # [Hubbard Model](@id Hubbard) |
2 | 2 |
|
3 | | -TODO |
| 3 | +## Ground state |
| 4 | +In this section, we use Hubbard model as an example to demonstrate how to simulate fermionic systems and compute multi-site correlations. We consider the square-lattice Hubbard model on a $W\times L$ cylinder. Note we use a relatively small bond dimension here to avoid too much online build time, which is not sufficient to obtain converged results. |
| 5 | + |
| 6 | +The Hamiltonian of Hubbard model reads\ |
| 7 | +$H = -t\sum_{\langle i, j\rangle\sigma} c_{i\sigma}^\dagger c_{j\sigma} + U\sum_i n_{i\uparrow}n_{i\downarrow}$. Here we set $t = 1$ as energy unit and use charge U(1) and spin SU(2) symmetries. |
| 8 | +```@example Hubbard |
| 9 | +using FiniteMPS |
| 10 | +using FiniteLattices |
| 11 | +using CairoMakie, Statistics |
| 12 | +
|
| 13 | +
|
| 14 | +mkpath("figs_Hubbard") # save figures |
| 15 | +
|
| 16 | +# parameters |
| 17 | +L = 24 |
| 18 | +W = 4 |
| 19 | +U = 12.0 |
| 20 | +Ntot = 88 # total particle number |
| 21 | +Dmin = 128 # min bond dimension |
| 22 | +Dmax = 1024 # max bond dimension |
| 23 | +
|
| 24 | +# generate lattices with FiniteLattices.jl |
| 25 | +Latt = YCSqua(L, W) |> Snake! |
| 26 | +
|
| 27 | +# generate the Hamiltonian MPO |
| 28 | +Tree = InteractionTree(size(Latt)) |
| 29 | +# hopping |
| 30 | +for (i, j) in neighbor(Latt; ordered = true) |
| 31 | + addIntr!(Tree, U1SU2Fermion.FdagF, (i, j), (true, true), -1.0; name = (:Fdag, :F), Z = U1SU2Fermion.Z) |
| 32 | +end |
| 33 | +
|
| 34 | +# U terms |
| 35 | +for i in 1:size(Latt) |
| 36 | + addIntr!(Tree, U1SU2Fermion.nd, i, U; name = :nd) |
| 37 | +end |
| 38 | +H = AutomataMPO(Tree) |
| 39 | +``` |
| 40 | +The above code generates the Hamiltonian MPO for the Hubbard model, the only new usage is setting `(true, true)` to indicate both operators are fermionic and passing the parity operator `Z` when adding hopping terms. Next we generate a random density distribution and obtain the corresponding MPS. |
| 41 | +```@example Hubbard |
| 42 | +# random density distribution |
| 43 | +lsn = zeros(Int64, size(Latt)) |
| 44 | +for _ in 1:Ntot |
| 45 | + i = rand(findall(==(minimum(lsn)), lsn)) |
| 46 | + lsn[i] += 1 |
| 47 | +end |
| 48 | +
|
| 49 | +# charge quantum numbers of horizontal bonds |
| 50 | +lsqc = [0] |
| 51 | +for n in reverse(lsn) # from right to left |
| 52 | + if n == 0 |
| 53 | + push!(lsqc, lsqc[end] + 1) |
| 54 | + elseif n == 1 |
| 55 | + push!(lsqc, lsqc[end]) |
| 56 | + else |
| 57 | + push!(lsqc, lsqc[end] - 1) |
| 58 | + end |
| 59 | +end |
| 60 | +lsqc = reverse(lsqc[2:end]) |
| 61 | +
|
| 62 | +# spaces of horizontal bonds |
| 63 | +lsspaces = map(1:size(Latt)) do i |
| 64 | + # the left boundary bond |
| 65 | + i == 1 && return Rep[U₁×SU₂]((lsqc[1], iseven(Ntot) ? 0 : 1/2) => 1) |
| 66 | + return Rep[U₁×SU₂]((lsqc[i], s) => 1 for s in 0:1/2:1/2) |
| 67 | +end |
| 68 | +Ψ = randMPS(fill(U1SU2Fermion.pspace, size(Latt)), lsspaces) |
| 69 | +``` |
| 70 | +Then we prepare the `ObservableTree` used to calculate observables. |
| 71 | +```@example Hubbard |
| 72 | +ObsTree = ObservableTree(size(Latt)) |
| 73 | +
|
| 74 | +# on-site terms |
| 75 | +for i in 1:size(Latt) |
| 76 | + addObs!(ObsTree, U1SU2Fermion.n, i; name = :n) # density |
| 77 | + addObs!(ObsTree, U1SU2Fermion.nd, i; name = :nd) # double occupancy |
| 78 | +end |
| 79 | +
|
| 80 | +# all-to all 2-site correlations |
| 81 | +for i in 1:size(Latt), j in i:size(Latt) |
| 82 | + # spin correlation |
| 83 | + addObs!(ObsTree, U1SU2Fermion.SS, (i, j), (false, false); name = (:S, :S)) |
| 84 | + # density correlation |
| 85 | + addObs!(ObsTree, (U1SU2Fermion.n, U1SU2Fermion.n), (i, j), (false, false); name = (:n, :n)) |
| 86 | + # single particle correlation |
| 87 | + addObs!(ObsTree, U1SU2Fermion.FdagF, (i, j), (true, true); name = (:Fdag, :F), Z = U1SU2Fermion.Z) |
| 88 | +end |
| 89 | +
|
| 90 | +# singlet pairing correlations of y-bonds |
| 91 | +Pairs = [(Latt[x, y], Latt[x, y % W + 1]) for x in 1:L for y in 1:W] |
| 92 | +for idx1 in 1:length(Pairs), idx2 in idx1:length(Pairs) |
| 93 | + addObs!(ObsTree, U1SU2Fermion.ΔₛdagΔₛ, (Pairs[idx1]..., Pairs[idx2]...), (true, true, true, true); name = (:Fdag, :Fdag, :F, :F), Z = U1SU2Fermion.Z) |
| 94 | +end |
| 95 | +
|
| 96 | +# observables of the initial random state |
| 97 | +calObs!(ObsTree, Ψ) |
| 98 | +Obs = convert(Dict, ObsTree) |
| 99 | +``` |
| 100 | +We can verify that the initial density distribution matches the random configuration via |
| 101 | +```@example Hubbard |
| 102 | +# check initial density distribution |
| 103 | +@assert all(i -> Obs["n"][(i,)] ≈ lsn[i], 1:length(lsn)) |
| 104 | +``` |
| 105 | +Then, we perform the CBE-DMRG to obtain the ground state. |
| 106 | +```@example Hubbard |
| 107 | +Env = Environment(Ψ', H, Ψ) |
| 108 | +lsD = [Dmin] |
| 109 | +lsE = Float64[] |
| 110 | +lsObs = Dict[] |
| 111 | +etol = 1e-5 # energy tolerance |
| 112 | +MaxSweeps = 100 # avoid dead loops |
| 113 | +
|
| 114 | +for i in 1:MaxSweeps |
| 115 | + D = lsD[end] |
| 116 | + info, _ = DMRGSweep1!(Env; K = 16, trunc = truncdim(D), |
| 117 | + CBEAlg = NaiveCBE(D + div(D, 4), 1e-8; rsvd = true), |
| 118 | + ) |
| 119 | + |
| 120 | + push!(lsE, info[2].dmrg[1].Eg) |
| 121 | + # @show D, lsE[end] |
| 122 | +
|
| 123 | + if i > 1 && (lsE[end-1] - lsE[end])/size(Latt) < etol |
| 124 | + calObs!(ObsTree, Ψ) |
| 125 | + push!(lsObs, convert(Dict, ObsTree)) |
| 126 | + |
| 127 | + D *= 2 # increase bond dimension |
| 128 | + D > Dmax && break # converged |
| 129 | + end |
| 130 | +
|
| 131 | + push!(lsD, D) |
| 132 | +end |
| 133 | +
|
| 134 | +# plot the energy vs nsweep |
| 135 | +fig = Figure(size = (480, 240)) |
| 136 | +ax = Axis(fig[1, 1]; |
| 137 | + xlabel = "nsweep", |
| 138 | + ylabel = "energy per site") |
| 139 | +scatterlines!(ax, lsE / size(Latt)) # per site |
| 140 | +save("figs_Hubbard/GS_energy.png", fig) |
| 141 | +``` |
| 142 | + |
| 143 | +The parameters are chosen as $U = 12$ and $\delta = 1/12$, where the ground state belongs to [LE1](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.2.033073) phase that exhibits half-filled charge density wave with $\lambda_\textrm{CDW}= 1/(2\delta) = 6$ and power-law pairing correlation. |
| 144 | + |
| 145 | +```@example Hubbard |
| 146 | +# illustrate CDW |
| 147 | +fig = Figure(size = (480, 500)) |
| 148 | +ax = Axis(fig[1, 1]; |
| 149 | + xlabel = L"x", |
| 150 | + ylabel = L"n(x)") |
| 151 | +map(lsObs, unique(lsD), range(0.2, 1.0;length = length(lsObs))) do Obs, D, α |
| 152 | + lsnx = map(1:L) do x |
| 153 | + return mean(1:W) do y |
| 154 | + Obs["n"][(Latt[x, y],)] |
| 155 | + end |
| 156 | + end |
| 157 | + scatterlines!(ax, 1:L, lsnx; |
| 158 | + color = (:blue, α), |
| 159 | + label = L"D = %$(D)" |
| 160 | + ) |
| 161 | +end |
| 162 | +
|
| 163 | +# pairing correlation in 1/4 to 3/4 bulk region |
| 164 | +ax2 = Axis(fig[2, 1]; |
| 165 | + xlabel = L"r", |
| 166 | + ylabel = L"\Phi_{yy}(r)", |
| 167 | + xscale = log10, |
| 168 | + yscale = log10, |
| 169 | +) |
| 170 | +lsr = 1:div(L, 2) - 1 |
| 171 | +map(lsObs, unique(lsD), range(0.2, 1.0;length = length(lsObs))) do Obs, D, α |
| 172 | + lsΦyy = map(lsr) do r |
| 173 | + mean([(x, y) for x in div(L, 4)+1 : div(3*L, 4)-r for y in 1:W]) do (x, y) |
| 174 | + Obs["FdagFdagFF"][(Latt[x, y], Latt[x, y % W + 1], Latt[x + r, y], Latt[x + r, y % W + 1])] |
| 175 | + end |
| 176 | + end |
| 177 | + scatterlines!(ax2, lsr, lsΦyy; |
| 178 | + color = (:blue, α), |
| 179 | + label = L"D = %$(D)" |
| 180 | + ) |
| 181 | +end |
| 182 | +axislegend(ax2; position = (0, 0), rowgap = 0) |
| 183 | +
|
| 184 | +save("figs_Hubbard/GS_Obs.png", fig) |
| 185 | +``` |
| 186 | + |
| 187 | +Finally, we illustrate the CDW and the pairing correlation. Emphasize again this is only a demonstrative example, a much larger bond dimension is required to reproduce the power-law pairing correlation. |
| 188 | + |
| 189 | +## Finite temperature |
| 190 | +TODO |
0 commit comments