Fermi-Hubbard model with $f\mathbb{Z}_2 \boxtimes U(1)$ symmetry, at large $U$ and half-filling

In this example, we will demonstrate how to handle fermionic PEPS tensors and how to optimize them. To that end, we consider the two-dimensional Hubbard model

\[H = -t \sum_{\langle i,j \rangle} \sum_{\sigma} \left( c_{i,\sigma}^+ c_{j,\sigma}^- - c_{i,\sigma}^- c_{j,\sigma}^+ \right) + U \sum_i n_{i,\uparrow}n_{i,\downarrow} - \mu \sum_i n_i\]

where $\sigma \in \{\uparrow,\downarrow\}$ and $n_{i,\sigma} = c_{i,\sigma}^+ c_{i,\sigma}^-$ is the fermionic number operator. As in previous examples, using fermionic degrees of freedom is a matter of creating tensors with the right symmetry sectors - the rest of the simulation workflow remains the same.

First though, we make the example deterministic by seeding the RNG, and we make our imports:

using Random
using TensorKit, PEPSKit
using MPSKit: add_physical_charge
Random.seed!(2928528937);

Defining the fermionic Hamiltonian

Let us start by fixing the parameters of the Hubbard model. We're going to use a hopping of $t=1$ and a large $U=8$ on a $2 \times 2$ unit cell:

t = 1.0
U = 8.0
lattice = InfiniteSquare(2, 2);

In order to create fermionic tensors, one needs to define symmetry sectors using TensorKit's FermionParity. Not only do we want use fermion parity but we also want our particles to exploit the global $U(1)$ symmetry. The combined product sector can be obtained using the Deligne product, called through ⊠ which is obtained by typing \boxtimes+TAB. We will not impose any extra spin symmetry, so we have:

fermion = fℤ₂
particle_symmetry = U1Irrep
spin_symmetry = Trivial
S = fermion ⊠ particle_symmetry

TensorKitSectors.ProductSector{Tuple{TensorKitSectors.FermionParity, TensorKitSectors.U1Irrep}}

The next step is defining graded virtual PEPS and environment spaces using S. Here we also use the symmetry sector to impose half-filling. That is all we need to define the Hubbard Hamiltonian:

D, χ = 1, 1
V_peps = Vect[S]((0, 0) => 2 * D, (1, 1) => D, (1, -1) => D)
V_env = Vect[S](
    (0, 0) => 4 * χ, (1, -1) => 2 * χ, (1, 1) => 2 * χ, (0, 2) => χ, (0, -2) => χ
)
S_aux = S((1, 1))
H₀ = hubbard_model(ComplexF64, particle_symmetry, spin_symmetry, lattice; t, U)
H = add_physical_charge(H₀, fill(S_aux, size(H₀.lattice)...));

Finding the ground state

Again, the procedure of ground state optimization is very similar to before. First, we define all algorithmic parameters:

boundary_alg = (; tol = 1.0e-8, alg = :simultaneous, trscheme = (; alg = :fixedspace))
gradient_alg = (; tol = 1.0e-6, alg = :eigsolver, maxiter = 10, iterscheme = :diffgauge)
optimizer_alg = (; tol = 1.0e-4, alg = :lbfgs, maxiter = 80, ls_maxiter = 3, ls_maxfg = 3)

(tol = 0.0001, alg = :lbfgs, maxiter = 80, ls_maxiter = 3, ls_maxfg = 3)

Second, we initialize a PEPS state and environment (which we converge) constructed from symmetric physical and virtual spaces:

physical_spaces = physicalspace(H)
virtual_spaces = fill(V_peps, size(lattice)...)
peps₀ = InfinitePEPS(randn, ComplexF64, physical_spaces, virtual_spaces)
env₀, = leading_boundary(CTMRGEnv(peps₀, V_env), peps₀; boundary_alg...);

[ Info: CTMRG init:	obj = +5.484842275412e+04 +4.469243203539e+04im	err = 1.0000e+00
[ Info: CTMRG conv 26:	obj = +8.371681846538e+04 -3.789700713241e-07im	err = 7.4963851544e-09	time = 14.25 sec

And third, we start the ground state search (this does take quite long):

peps, env, E, info = fixedpoint(
    H, peps₀, env₀; boundary_alg, gradient_alg, optimizer_alg, verbosity = 3
)
@show E;

[ Info: LBFGS: initializing with f = 6.680719803101, ‖∇f‖ = 9.5849e+00
┌ Warning: Linesearch not converged after 1 iterations and 4 function evaluations:
│ α = 2.50e+01, dϕ = -1.49e-01, ϕ - ϕ₀ = -2.88e+00
└ @ OptimKit ~/.julia/packages/OptimKit/G6i79/src/linesearches.jl:148
[ Info: LBFGS: iter    1, time  874.29 s: f = 3.801380487744, ‖∇f‖ = 2.3456e+01, α = 2.50e+01, m = 0, nfg = 4
┌ Warning: Linesearch not converged after 1 iterations and 4 function evaluations:
│ α = 2.50e+01, dϕ = -5.73e-03, ϕ - ϕ₀ = -3.81e+00
└ @ OptimKit ~/.julia/packages/OptimKit/G6i79/src/linesearches.jl:148
[ Info: LBFGS: iter    2, time  963.43 s: f = -0.009727650286, ‖∇f‖ = 3.2049e+00, α = 2.50e+01, m = 0, nfg = 4
[ Info: LBFGS: iter    3, time  977.14 s: f = -0.115210826428, ‖∇f‖ = 2.7846e+00, α = 1.00e+00, m = 1, nfg = 1
[ Info: LBFGS: iter    4, time  990.71 s: f = -0.616412169228, ‖∇f‖ = 2.3680e+00, α = 1.00e+00, m = 2, nfg = 1
[ Info: LBFGS: iter    5, time 1003.13 s: f = -0.817801148604, ‖∇f‖ = 1.9111e+00, α = 1.00e+00, m = 3, nfg = 1
[ Info: LBFGS: iter    6, time 1015.30 s: f = -0.990286615265, ‖∇f‖ = 2.3790e+00, α = 1.00e+00, m = 4, nfg = 1
[ Info: LBFGS: iter    7, time 1026.31 s: f = -1.142787566798, ‖∇f‖ = 1.5680e+00, α = 1.00e+00, m = 5, nfg = 1
[ Info: LBFGS: iter    8, time 1036.16 s: f = -1.238274330219, ‖∇f‖ = 3.5015e+00, α = 1.00e+00, m = 6, nfg = 1
[ Info: LBFGS: iter    9, time 1046.93 s: f = -1.438136282421, ‖∇f‖ = 1.3366e+00, α = 1.00e+00, m = 7, nfg = 1
[ Info: LBFGS: iter   10, time 1056.67 s: f = -1.523107107396, ‖∇f‖ = 1.3496e+00, α = 1.00e+00, m = 8, nfg = 1
[ Info: LBFGS: iter   11, time 1077.81 s: f = -1.619305193101, ‖∇f‖ = 1.1951e+00, α = 1.72e-01, m = 9, nfg = 2
[ Info: LBFGS: iter   12, time 1097.26 s: f = -1.681451834691, ‖∇f‖ = 9.4848e-01, α = 2.37e-01, m = 10, nfg = 2
[ Info: LBFGS: iter   13, time 1107.99 s: f = -1.720734533825, ‖∇f‖ = 1.4216e+00, α = 1.00e+00, m = 11, nfg = 1
[ Info: LBFGS: iter   14, time 1117.66 s: f = -1.770831967062, ‖∇f‖ = 6.2747e-01, α = 1.00e+00, m = 12, nfg = 1
[ Info: LBFGS: iter   15, time 1127.31 s: f = -1.807572162185, ‖∇f‖ = 5.1320e-01, α = 1.00e+00, m = 13, nfg = 1
[ Info: LBFGS: iter   16, time 1138.09 s: f = -1.859768355558, ‖∇f‖ = 7.1320e-01, α = 1.00e+00, m = 14, nfg = 1
[ Info: LBFGS: iter   17, time 1147.68 s: f = -1.893160382125, ‖∇f‖ = 6.7323e-01, α = 1.00e+00, m = 15, nfg = 1
[ Info: LBFGS: iter   18, time 1157.44 s: f = -1.923092489408, ‖∇f‖ = 5.5419e-01, α = 1.00e+00, m = 16, nfg = 1
[ Info: LBFGS: iter   19, time 1168.12 s: f = -1.948142544093, ‖∇f‖ = 4.7661e-01, α = 1.00e+00, m = 17, nfg = 1
[ Info: LBFGS: iter   20, time 1177.91 s: f = -1.969512080077, ‖∇f‖ = 4.1608e-01, α = 1.00e+00, m = 18, nfg = 1
[ Info: LBFGS: iter   21, time 1187.87 s: f = -1.982557838199, ‖∇f‖ = 4.5138e-01, α = 1.00e+00, m = 19, nfg = 1
[ Info: LBFGS: iter   22, time 1198.67 s: f = -1.994007805763, ‖∇f‖ = 3.1538e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   23, time 1208.61 s: f = -2.002836016203, ‖∇f‖ = 3.0511e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   24, time 1218.59 s: f = -2.014062852739, ‖∇f‖ = 3.3491e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   25, time 1229.70 s: f = -2.022023251661, ‖∇f‖ = 4.3758e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   26, time 1239.66 s: f = -2.030112566345, ‖∇f‖ = 2.0509e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   27, time 1249.67 s: f = -2.035073683394, ‖∇f‖ = 1.6307e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   28, time 1261.85 s: f = -2.038663850455, ‖∇f‖ = 1.6880e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   29, time 1271.84 s: f = -2.041323592429, ‖∇f‖ = 2.4114e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   30, time 1281.66 s: f = -2.044997390530, ‖∇f‖ = 1.2115e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   31, time 1292.49 s: f = -2.046747469529, ‖∇f‖ = 9.5108e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   32, time 1302.32 s: f = -2.048741416292, ‖∇f‖ = 1.0509e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   33, time 1312.17 s: f = -2.049793769908, ‖∇f‖ = 1.7378e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   34, time 1323.04 s: f = -2.051022900848, ‖∇f‖ = 6.4055e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   35, time 1332.89 s: f = -2.051499900828, ‖∇f‖ = 4.9307e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   36, time 1342.73 s: f = -2.051918795787, ‖∇f‖ = 6.2013e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   37, time 1353.60 s: f = -2.052357188363, ‖∇f‖ = 9.4494e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   38, time 1363.45 s: f = -2.052855317283, ‖∇f‖ = 4.8219e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   39, time 1373.40 s: f = -2.053138284528, ‖∇f‖ = 3.5599e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   40, time 1384.12 s: f = -2.053404037719, ‖∇f‖ = 4.1844e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   41, time 1393.97 s: f = -2.053605747242, ‖∇f‖ = 5.7514e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   42, time 1403.69 s: f = -2.053822345457, ‖∇f‖ = 3.1996e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   43, time 1414.45 s: f = -2.054015631924, ‖∇f‖ = 3.1314e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   44, time 1424.21 s: f = -2.054206835742, ‖∇f‖ = 4.1588e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   45, time 1434.10 s: f = -2.054349141892, ‖∇f‖ = 6.7905e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   46, time 1444.85 s: f = -2.054531571463, ‖∇f‖ = 2.9227e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   47, time 1454.61 s: f = -2.054628027248, ‖∇f‖ = 2.5100e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   48, time 1464.37 s: f = -2.054735541814, ‖∇f‖ = 3.1538e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   49, time 1475.16 s: f = -2.054896782689, ‖∇f‖ = 3.4823e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   50, time 1494.75 s: f = -2.055018285181, ‖∇f‖ = 5.2680e-02, α = 5.17e-01, m = 20, nfg = 2
[ Info: LBFGS: iter   51, time 1505.51 s: f = -2.055214629205, ‖∇f‖ = 3.0513e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   52, time 1515.33 s: f = -2.055401907932, ‖∇f‖ = 2.8740e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   53, time 1525.32 s: f = -2.055643036846, ‖∇f‖ = 4.1540e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   54, time 1536.28 s: f = -2.055979753449, ‖∇f‖ = 6.0310e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   55, time 1546.19 s: f = -2.056292876566, ‖∇f‖ = 6.4503e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   56, time 1556.12 s: f = -2.056764405334, ‖∇f‖ = 4.5709e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   57, time 1567.05 s: f = -2.057301128966, ‖∇f‖ = 5.8535e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   58, time 1577.08 s: f = -2.057684443651, ‖∇f‖ = 7.0407e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   59, time 1587.05 s: f = -2.058273607978, ‖∇f‖ = 6.4287e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   60, time 1598.01 s: f = -2.058991887288, ‖∇f‖ = 8.8941e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   61, time 1607.94 s: f = -2.059459011130, ‖∇f‖ = 1.1553e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   62, time 1617.89 s: f = -2.060066395726, ‖∇f‖ = 6.9440e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   63, time 1628.87 s: f = -2.060520108822, ‖∇f‖ = 8.4931e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   64, time 1648.98 s: f = -2.060815447626, ‖∇f‖ = 1.2115e-01, α = 5.26e-01, m = 20, nfg = 2
[ Info: LBFGS: iter   65, time 1669.96 s: f = -2.060925751714, ‖∇f‖ = 8.3903e-02, α = 5.47e-01, m = 20, nfg = 2
[ Info: LBFGS: iter   66, time 1679.97 s: f = -2.061209735729, ‖∇f‖ = 5.4010e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   67, time 1690.95 s: f = -2.061580165208, ‖∇f‖ = 5.5975e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   68, time 1701.12 s: f = -2.062036980891, ‖∇f‖ = 7.8898e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   69, time 1711.14 s: f = -2.062251708572, ‖∇f‖ = 1.1537e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   70, time 1722.19 s: f = -2.062519628412, ‖∇f‖ = 1.2953e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   71, time 1732.19 s: f = -2.063059957021, ‖∇f‖ = 7.2868e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   72, time 1742.20 s: f = -2.063313168886, ‖∇f‖ = 5.2530e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   73, time 1753.23 s: f = -2.063715484908, ‖∇f‖ = 5.0261e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   74, time 1763.31 s: f = -2.064332647622, ‖∇f‖ = 7.7210e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   75, time 1773.36 s: f = -2.064773364480, ‖∇f‖ = 1.2451e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   76, time 1784.47 s: f = -2.065371916869, ‖∇f‖ = 6.8015e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   77, time 1794.59 s: f = -2.065945928950, ‖∇f‖ = 7.6664e-02, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   78, time 1804.62 s: f = -2.066640737009, ‖∇f‖ = 1.1191e-01, α = 1.00e+00, m = 20, nfg = 1
[ Info: LBFGS: iter   79, time 1815.78 s: f = -2.067648349055, ‖∇f‖ = 2.3836e-01, α = 1.00e+00, m = 20, nfg = 1
┌ Warning: LBFGS: not converged to requested tol after 80 iterations and time 1826.00 s: f = -2.069253135434, ‖∇f‖ = 2.0413e-01
└ @ OptimKit ~/.julia/packages/OptimKit/G6i79/src/lbfgs.jl:197
E = -2.069253135433506

Finally, let's compare the obtained energy against a reference energy from a QMC study by Qin et al.. With the parameters specified above, they obtain an energy of $E_\text{ref} \approx 4 \times -0.5244140625 = -2.09765625$ (the factor 4 comes from the $2 \times 2$ unit cell that we use here). Thus, we find:

E_ref = -2.09765625
@show (E - E_ref) / E_ref;

(E - E_ref) / E_ref = -0.013540404709539051

This page was generated using Literate.jl.