Fermionic PEPS with flat Z2 storage

This example shows how to set up a fermionic PEPS with Z2 symmetry and the ‘flat’ backend using symmray and quimb, then optimize using imaginary time simple update, then evaluate the energy using boundary contraction. Because we are using Z2 symmetry we rely on setting a chemical potential to reach half-filling.

%config InlineBackend.figure_formats = ['svg']
import quimb as qu
import quimb.tensor as qtn

import symmray as sr

symmetry = "Z2"

Lx = 4
Ly = 4
D = 6
seed = 1

# note fully random initialization will not be a very
# good initial state, used as a demonstration only here
peps = sr.PEPS_fermionic_rand(
    symmetry=symmetry,
    Lx=Lx,
    Ly=Ly,
    bond_dim=D,
    seed=seed,
    phys_dim=4,
    flat=True,
    subsizes="equal",
)

# verify the spin structure using a basic cluster expectation
Ospin = sr.fermi_spin_operator_local_array(symmetry, flat=True)

peps.compute_local_expectation_cluster(
    {(site,): Ospin for site in peps.gen_site_coos()},
    max_distance=1,
    return_all=True,
)

{((0, 0),): np.float64(0.0327419300969182),
 ((0, 1),): np.float64(-0.01598635686807416),
 ((0, 2),): np.float64(-0.056781828870607676),
 ((0, 3),): np.float64(-0.04050117272143959),
 ((1, 0),): np.float64(0.032587690022989765),
 ((1, 1),): np.float64(0.00016372219003649158),
 ((1, 2),): np.float64(-0.006227043206973079),
 ((1, 3),): np.float64(-0.033253735745883584),
 ((2, 0),): np.float64(0.056101826385629536),
 ((2, 1),): np.float64(0.012242773893809566),
 ((2, 2),): np.float64(0.0028582463963720273),
 ((2, 3),): np.float64(0.0128546755014783),
 ((3, 0),): np.float64(0.08064024594552714),
 ((3, 1),): np.float64(-0.022580207467178273),
 ((3, 2),): np.float64(0.0018963528877691538),
 ((3, 3),): np.float64(-0.05837888912565602)}

terms = sr.ham_fermi_hubbard_from_edges(
    symmetry=symmetry,
    edges=tuple(peps.gen_bond_coos()),
    t=1.0,
    U=8.0,
    mu=4.0,
    flat=True,
)
ham = qtn.LocalHamGen(terms)

su = qtn.SimpleUpdateGen(
    peps,
    ham,
    # flat only supports zero cutoff
    cutoff=0.0,
    second_order_reflect=True,
    # SimpleUpdateGen computes cluster energies by default
    # which might not be accurate
    compute_energy_every=10,
    compute_energy_opts=dict(max_distance=1),
    compute_energy_per_site=True,
    # use a fixed trotterization order
    ordering="sort",
    # if the gauge difference drops below this, we consider the PEPS converged
    tol=1e-9,
)

# run the evolution, these are reasonable defaults
tau = 0.5 * D ** (-3 / 2)
steps = round(50 / tau)

su.evolve(steps, tau=tau)

n=1470, D=6, tau=0.034, max|dS|=2.43e-07, energy≈-4.4002: 100%|###########################################################| 1470/1470 [03:49<00:00,  6.40it/s]

su.plot();

gs = su.get_state()

Check energy properly with boundary contraction and increasing \(\chi\).

terms_without_mu = sr.ham_fermi_hubbard_from_edges(
    symmetry=symmetry,
    edges=tuple(gs.gen_bond_coos()),
    t=1.0,
    U=8.0,
    mu=0.0,
)

terms_without_mu = {k: v.to_flat() for k, v in terms_without_mu.items()}

full_energies = {}

for chi in [16, 24, 32, 48, 64]:
    en = (
        gs.compute_local_expectation(
            terms_without_mu,
            normalized=True,
            max_bond=chi,
            cutoff=0.0,
        )
        / gs.nsites
    )
    print(f"chi={chi}: E/N = {en}")
    full_energies[chi] = en

chi=16: E/N = -0.4129848610534672
chi=24: E/N = -0.4128461676264311
chi=32: E/N = -0.4128453545547058
chi=48: E/N = -0.4128488960852359
chi=64: E/N = -0.41284968352636936

qu.plot(full_energies.keys(), full_energies.values(), marker=".");

# verify the spin structure using a basic cluster expectation
Ospin = sr.fermi_spin_operator_local_array(symmetry).to_flat()
# Nelec = sr.fermi_number_operator_spinful_local_array(symmetry).to_flat()

gs.compute_local_expectation_cluster(
    {(site,): Ospin for site in peps.gen_site_coos()},
    max_distance=2,
    return_all=True,
)

{((0, 0),): np.float64(-0.2953250799887245),
 ((0, 1),): np.float64(0.2941976047163778),
 ((0, 2),): np.float64(-0.29425867264051664),
 ((0, 3),): np.float64(0.29549669709151427),
 ((1, 0),): np.float64(0.2939438711703014),
 ((1, 1),): np.float64(-0.3136221430478892),
 ((1, 2),): np.float64(0.3134570999248664),
 ((1, 3),): np.float64(-0.29392501511436775),
 ((2, 0),): np.float64(-0.2941208664940671),
 ((2, 1),): np.float64(0.3137381559209075),
 ((2, 2),): np.float64(-0.3136845444304519),
 ((2, 3),): np.float64(0.2941651876360549),
 ((3, 0),): np.float64(0.2955559103339454),
 ((3, 1),): np.float64(-0.29411440423158275),
 ((3, 2),): np.float64(0.29414247661918214),
 ((3, 3),): np.float64(-0.2955821994756589)}