I'm not completely sure if I understand what you want to do and what you mean with "Gutzwiller projection".
As far as I understand, you have found an MPS of SpinHalfFermionSite -based model with both Sz and N conservation.
Now you want to "project" into the subspace of fixed particle number N=1 on each site, right?
The projection operator \(P_1 = \mathbb{1} - |vac\rangle \langle vac| - | \uparrow \downarrow \rangle \langle \uparrow \downarrow| = |\uparrow\rangle \langle \uparrow | + | \downarrow \rangle \langle \downarrow | \) is
diagonal, hence it commutes with both the Sz and N operator! Hence, you actually
can write it down while conserving the charges!
So the simple approach would be to just add a local operator doing the projection and applying it to psi.
Doing this will make it trivial to do measurements afterwards.
Code: Select all
psi = .... # from DMRG, assuming that it has
# add a new operator "P_1" projection operator to each site of the MPS
# (you can also do this already when you initialize the site for a custom model)
for site in psi.sites:
if "P_1" in site.opnames:
# we already added the operator, as the Site instances are re-used
continue
# add the P_1 operator
P_1 = np.zeros([site.dim, site.dim])
for i in site.state_indices(['up', 'down']):
P_1[i, i] = 1.
site.add_op("P_1", P_1, hc="P_1")
# apply the product of single-site projectors P = \product_x (P_1)_x
for x in range(psi.L - 1):
# hack: unitary=True just avoids calling `psi.canonical_form`
psi.apply_local_op(x, 'P_1', unitary=True)
psi.apply_local_op(psi.L - 1, unitary=False, renormalize=True)
# the last call calls psi.canonical_form() in the end
# measure psi expecaton values etc as usual
Yesterday, I implemented MPS.
apply_product_op, a slight generalization of MPS.
apply_local_op.
Using that function, the last loop becomes a bit easier. If you have a TenPy version after
ac42eed45, you can use:
Code: Select all
# apply the product of single-site projectors P = \product_x (P_1)_x
psi.apply_product_op(["P_1"] * psi.L, renormalize=True)
Alternatively, if you want to be really pedantic (or need to optimize this for a looong-running subsequent algorithm) you can in principle also project out the non-used local states completely, basically using Array.
iproject, but that's a bit more complicated; you need to take care to do this consistently on all the B tensors and sites at the same time. In that case, you could even drop the N charge completely. Let me know if you would like to see code for that.