Spinful fermionic amplitudes¶
import numpy as np
import quimb.operator as qop
import quimb.tensor as qtn
import symmray as sr
Lx = 2
Ly = 3
D = 4
seed = 42
# only the flat backend is compatible with jax.jit
flat = True
peps = sr.networks.PEPS_fermionic_rand(
"Z2",
Lx,
Ly,
D,
phys_dim=4,
# put an odd number of odd sites in, for testing
site_charge=lambda site: int(site in [(0, 0), (0, 1), (1, 0)]),
subsizes="equal",
flat=flat,
seed=seed,
)
edges = qtn.edges_2d_square(Lx, Ly)
sites = [(i, j) for i in range(Lx) for j in range(Ly)]
# exact energy via local expectation contraction
terms = sr.hamiltonians.ham_fermi_hubbard_from_edges(
"Z2",
edges=edges,
U=3.71,
mu=0.119,
)
if flat:
terms = {k: v.to_flat() for k, v in terms.items()}
eref = peps.compute_local_expectation_exact(terms, normalized=True)
eref
np.float64(6.439453406305376)
# compute the spin expectation at site (1,1)
G_s11 = sr.fermi_spin_operator_local_array("Z2")
if flat:
G_s11 = G_s11.to_flat()
s11ref = peps.local_expectation_exact(G_s11, where=[(1, 1)], normalized=True)
s11ref
np.float64(0.03469485193879984)
# get a symbolic representation of the Hamiltonian
H = qop.fermi_hubbard_from_edges(
edges,
U=3.71,
mu=0.119,
# this ordering pairs spins together, as with the fermionic TN
order=lambda site: (site[1], site[0]),
sector=int(sum(ary.charge for ary in peps.arrays) % 2),
symmetry="Z2",
)
hs = H.hilbert_space
hs.sites
(('↑', (0, 0)),
('↓', (0, 0)),
('↑', (0, 1)),
('↓', (0, 1)),
('↑', (0, 2)),
('↓', (0, 2)),
('↑', (1, 0)),
('↓', (1, 0)),
('↑', (1, 1)),
('↓', (1, 1)),
('↑', (1, 2)),
('↓', (1, 2)))
if flat:
import jax
import jax.numpy as jnp
peps.apply_to_arrays(jnp.array)
def flat_amplitude(fx):
# convert neighboring pairs (up, down) to single index 0..3
# these should match up with the phys_dim ordering above
fx = 2 * (fx[::2] != fx[1::2]) + fx[1::2]
selector = {peps.site_ind(site): val for site, val in zip(peps.sites, fx)}
tnb = peps.isel(selector)
return tnb.contract()
if flat:
flat_amplitude = jax.jit(flat_amplitude)
# compute the full exact energy at the amplitude level
O = 0.0
p = 0.0
fcs = []
for i in range(hs.size):
fx = hs.rank_to_flatconfig(i)
xpsi = flat_amplitude(fx)
if not xpsi:
continue
pi = abs(xpsi) ** 2
p += pi
Oloc = 0.0
for fy, hxy in zip(*H.flatconfig_coupling(fx)):
ypsi = flat_amplitude(fy)
Oloc = Oloc + hxy * ypsi / xpsi
O += Oloc * pi
O / p
Array(6.43945, dtype=float32)
eref
np.float64(6.439453406305376)
# define symbolic spin operator at site (1,1)
S = qop.SparseOperatorBuilder(
terms=[
(+1 / 2, ("n", ("↑", (1, 1)))),
(-1 / 2, ("n", ("↓", (1, 1)))),
],
hilbert_space=H.hilbert_space,
)
# compute the full exact energy at the amplitude level
s = 0.0
p = 0.0
fcs = []
for i in range(hs.size):
fx = hs.rank_to_flatconfig(i)
xpsi = flat_amplitude(fx)
if not xpsi:
continue
pi = abs(xpsi) ** 2
p += pi
sloc = 0.0
for fy, hxy in zip(*S.flatconfig_coupling(fx)):
ypsi = flat_amplitude(fy)
sloc = sloc + hxy * ypsi / xpsi
s += sloc * pi
s / p
Array(0.03469491, dtype=float32)
s11ref
np.float64(0.03469485193879984)
If one only has access to the hamiltonian in the another ordering, e.g. grouped by spin species, there will be extra signs to account for:
H = qop.fermi_hubbard_from_edges(
edges,
U=3.71,
mu=0.119,
# this ordering
# order=lambda site: (site[1], site[0]),
sector=1,
symmetry="Z2",
)
hs = H.hilbert_space
Ha = qop.fermi_hubbard_from_edges(
edges,
U=3.71,
mu=0.119,
# this ordering
order=lambda site: (site[1], site[0]),
sector=1,
symmetry="Z2",
)
hsa = Ha.hilbert_space
print(
" spin first ordering Hij : site-first ordering Hij"
)
for i in [633, 1846]:
fx = hs.rank_to_flatconfig(i)
fys, hxys = H.flatconfig_coupling(fx)
fxa = np.empty_like(fx)
fxa[::2] = fx[: hs.nsites // 2]
fxa[1::2] = fx[hs.nsites // 2 :]
fyas, hxyas = Ha.flatconfig_coupling(fxa)
print(
f"{i:>4} {fx[: hs.nsites // 2]} {fx[hs.nsites // 2 :]} {fxa}"
)
print("---- couples to ... with ...")
for i in range(len(fys)):
fy = fys[i]
fya = fyas[i]
hxy = hxys[i]
hxya = hxyas[i]
print(
" ->",
fy[: hs.nsites // 2],
fy[hs.nsites // 2 :],
# f"{hxy:>6.4g}", " : ", fya[::2], fya[1::2], f"{hxya:>6.4g}"
f"{hxy:>6.4g}",
" : ",
fya,
f"{hxya:>6.4g}",
)
print("----")
print()
spin first ordering Hij : site-first ordering Hij
633 [0 1 0 0 1 1] [1 1 0 0 1 1] [0 1 1 1 0 0 0 0 1 1 1 1]
---- couples to ... with ...
-> [1 0 0 0 1 1] [1 1 0 0 1 1] -1 : [1 1 0 1 0 0 0 0 1 1 1 1] 1
-> [0 1 0 0 1 1] [0 1 0 1 1 1] 1 : [0 0 1 1 0 0 0 1 1 1 1 1] -1
-> [0 0 1 0 1 1] [1 1 0 0 1 1] -1 : [0 1 0 1 1 0 0 0 1 1 1 1] 1
-> [0 1 0 0 1 1] [1 0 1 0 1 1] -1 : [0 1 1 0 0 1 0 0 1 1 1 1] -1
-> [0 1 1 0 1 0] [1 1 0 0 1 1] 1 : [0 1 1 1 1 0 0 0 1 1 0 1] -1
-> [0 1 0 0 1 1] [1 1 1 0 1 0] 1 : [0 1 1 1 0 1 0 0 1 1 1 0] 1
-> [0 1 0 1 0 1] [1 1 0 0 1 1] -1 : [0 1 1 1 0 0 1 0 0 1 1 1] -1
-> [0 1 0 0 1 1] [1 1 0 1 0 1] -1 : [0 1 1 1 0 0 0 1 1 0 1 1] 1
-> [0 1 0 0 1 1] [1 1 0 0 1 1] 10.3 : [0 1 1 1 0 0 0 0 1 1 1 1] 10.3
----
1846 [1 1 1 0 0 1] [1 0 1 1 0 0] [1 1 1 0 1 1 0 1 0 0 1 0]
---- couples to ... with ...
-> [1 1 1 0 0 1] [0 1 1 1 0 0] -1 : [1 0 1 1 1 1 0 1 0 0 1 0] 1
-> [0 1 1 1 0 1] [1 0 1 1 0 0] -1 : [0 1 1 0 1 1 1 1 0 0 1 0] -1
-> [1 1 1 0 0 1] [1 1 0 1 0 0] -1 : [1 1 1 1 1 0 0 1 0 0 1 0] 1
-> [1 0 1 0 1 1] [1 0 1 1 0 0] 1 : [1 1 0 0 1 1 0 1 1 0 1 0] 1
-> [1 1 1 0 0 1] [1 0 0 1 0 1] 1 : [1 1 1 0 1 0 0 1 0 0 1 1] -1
-> [1 1 1 0 0 1] [1 0 1 0 1 0] -1 : [1 1 1 0 1 1 0 0 0 1 1 0] -1
-> [1 1 1 0 1 0] [1 0 1 1 0 0] -1 : [1 1 1 0 1 1 0 1 1 0 0 0] -1
-> [1 1 1 0 0 1] [1 0 1 1 0 0] 6.587 : [1 1 1 0 1 1 0 1 0 0 1 0] 6.587
----