TeNPy Forum

How can we add terms like \(O_k=\Pi_{<i,j>} (1-n_i n_j)B_k \quad i,j \neq k\) to Hamitonian ?

At a glimplse, this looks horrible for DMRG, because each term is non-local, giving you an immense MPO bond dimension if you implement it naively.

Do you have only those terms, i.e. is \(H=\left(\prod_{\langle i, j \rangle} (1-n_i n_j) \right) H^0 \) with a "usual" \(H^0\)?
It looks like you're trying to project out nearest-neighbour excitations in a constrained Rydberg-type model, similar to the PXP model.
Can you actually re-write it such that you only have local few-body terms, \(O_k = \prod_{\langle i, k\rangle} (1-n_i n_k ) B_k\) where the product is only over nearest neighbors i of k? Possibly applied to an initial state fulfilling the constraints? Then you have a chance to implement it as an MPO using the multi_coupling_terms....

Yes, true. That is why I first do the projection and have few-body terms where I use multi_coupling_terms to implement the Hamiltonian. Just curiosity, how would we implement it naively? How would we implement this sort of product terms given in my original post?

very naively, you can multiply out all the products such that you're left with a sum of terms, many of which are long-range.

Less naively, you'd try to rewrite it in the nearest-neighbor version \(O_k = \prod_{\langle i, k\rangle} (1-n_i n_k ) B_k\) mentioned above, and multiply that out to get sums of terms - if you have a reasonable number of nearest neighbors, that should be somewhat okay.

If I had to implment it, I would try to be more clever, decompose the \( H = P H^0\), and try to somehow implement the P by modifying the matvec of the (effective) Hamiltonian, projecting out the "invalid" states. I'm absolutely not sure about the details, but it might be possible to implement something along the lines of this paper.