Bose-Hubbard model fails in extended Hilbert space for newer versions

Created by: RomanArmenta

Hello, we are working on an extended Hilbert space in order to calculate the ground state energy of the Bose-Hubbard model. We added spin 1/2 sites between every boson site in the one-dimensional chain, but we left the Hamiltonian unchanged. We expected that the ground state energy should be the same in this extended Hilbert space, but in the newer versions of Netket (post 3.3.3), this energy differs from the actual solution. Even more, every time we run the exact diagonalization of the Hamiltonian matrix (either with the full_ed or lanczos_ed method), its eigenvalues change; we also diagonalized the same Hamiltonian matrix with the linalg package of numpy and the result also differs every time we do it. This doesn't happen in the older versions; the problem begins in the 3.3.4 version of Netket.

We claim that the problem resides in the way that the newer versions of Netket interpret the Hamiltonian matrices with mixed Hilbert spaces. When we do the same Bose-Hubbard model in a Fock-Boson space, the calculation gives the correct value of the ground state energy.

Here is the code that we run in both versions of Netket, and a video of the discrepancy.

Notebook ran in Nektet 3.3.3: https://www.youtube.com/watch?v=ENhRJfYg7dg Notebook ran in Netket 3.4.1: https://www.youtube.com/watch?v=Q3XfWrnR7LU

You can download the Jupyter Script here: https://www.dropbox.com/s/dza4kbyem2ycg6v/BoseHubbardNetket.ipynb?dl=0

# Extended Hilbert space
n_max = 3
L = 3 # Number of spin sites
hil_gen = nk.hilbert.Fock(n_max = n_max)*nk.hilbert.Spin(1/2)
hi_ext = hil_gen
for i in range(L-1):
    hi_ext = hi_ext * hil_gen 
hi_ext = hi_ext * nk.hilbert.Fock(n_max = n_max)
# Boson-Fock Hilbert space
N = 4 # Number of bosons in the chain
# Chain graph
g = nk.graph.Chain(length=N, pbc=False)
hi_fock = nk.hilbert.Fock(n_max=n_max, N=N)

# Bose-Hubbard Hamiltonian for extended Hilbert space
J = 0.5
U = 0.5
h = create(hi_ext,0)*destroy(hi_ext,2) + create(hi_ext,2)*destroy(hi_ext,4) + create(hi_ext,4)*destroy(hi_ext,6)
h_hc = h.H.collect()
h_u = (number(hi_ext,0)*(number(hi_ext,0)-1) + number(hi_ext,2)*(number(hi_ext,2)-1) + number(hi_ext,4)*(number(hi_ext,4)-1) + 
      number(hi_ext,6)*(number(hi_ext,6)-1))
H_bh = -J*(h + h_hc) + (U/2)*h_u

# Bose-Hubbard Hamiltonian for the Fock-Bose Hilbert space
H_bh2 = nk.operator.BoseHubbard(hilbert=hi_fock,U=0.5,J=0.5,graph=g)

# lanczos computation of the eigenvalues of the Hamiltonian with the extended Hilbert space
E_bh, ket_bh = nk.exact.lanczos_ed(H_bh, compute_eigenvectors=True)
print("Exact ground state energy = {0:.3f}".format(E_bh[0]))

# full computation of the eigenvalues of the Hamiltonian with the extended Hilbert space
E_bhfull = nk.exact.full_ed(H_bh)
print("Exact ground state energy = {0:.3f}".format(E_bhfull[0]))

# numpy.linalg computation of the eigenvalues of the Hamiltonian with the extended Hilbert space
M_bh = H_bh.to_dense()
E_bhla = min(np.linalg.eig(M_bh)[0])
print("Ground state energy:", E_bhla.real)

# scipy sparse computation of the eigenvalues of the Hamiltonian with the extended Hilbert space
H_bhsparse = H_bh.to_sparse()
eig_vals, eig_vecs = eigsh(H_bhsparse, k=1, which="SA")
print("eigenvalues with scipy sparse:", eig_vals)

# lanczos computation of the eigenvalues of the Hamiltonian with the Fock Hilbert space
E_bh2, ket_bh2 = nk.exact.lanczos_ed(H_bh2, compute_eigenvectors=True)
print("Exact ground state energy = {0:.3f}".format(E_bh2[0]))

# full computation of the eigenvalues of the Hamiltonian with the Fock Hilbert space
E_bhfull2= nk.exact.full_ed(H_bh2)
print("Exact ground state energy = {0:.3f}".format(E_bhfull2[0]))

# numpy.linalg computation of the eigenvalues of the Hamiltonian with the Fock Hilbert space
M_bh2 = H_bh2.to_dense()
E_bhla2 = min(np.linalg.eig(M_bh2)[0])
print("Ground state energy:", E_bhla2.real)

# scipy sparse computation of the eigenvalues of the Hamiltonian with the extended Hilbert space
H_bhsparse2 = H_bh2.to_sparse()
eig_vals2, eig_vecs2 = eigsh(H_bhsparse2, k=1, which="SA")
print("eigenvalues with scipy sparse:", eig_vals2)