Peierl substitution/custom complex hopping parameters in Bose-Hubbard (or Fermi-Hubbard) models

[jeohjee]

I am interested in studying Bose-Hubbard model with complex hoppings. To this end I wrote the following code for my Hamiltonian (definition of the lattice geometry is not included):

positions = graph.positions # lattice positions
q_vec = np.array([0,0]) #vector potential 

def a(site):
    return nk.operator.boson.destroy(hi, site)

def adag(site):
    return nk.operator.boson.create(hi, site)

def nc(site):
    return nk.operator.boson.number(hi, site)

U = 4.0
t = 1.0
ha = LocalOperator(hi, dtype=complex)
for i in range(pos_shape[0]):
    for j in range(pos_shape[0]):
        if i==j:
            continue;
        r_tmp_ij = [positions[i,0] - positions[j,0],positions[i,1] - positions[j,1] ]
        exp_term = np.exp(1j*(q_vec[0]*r_tmp_ij[0] + q_vec[1]*r_tmp_ij[1]))
        ha += t *exp_term* adag(i) * a(j)

    ha += U * nc(i) * (nc(i) -1.0)/2

With this Hamiltonian my code crashes due to high memory cost. If I just use the usual ha = nk.operator.BoseHubbard(hilbert=hi, graph=graph, U=4.0, J= -1.0), the program runs just fine. But with the default BoseHubbard constructor I cannot input finite vector potentials to my lattice model. I tried to build my own customized Bose-Hubbard class based on the existing BoseHubbard class but this seems bit involved and I got errors.

Any ideas how to proceed?

[gcarleo]

Hi, how many sites do you have? I am surprised this crashes

[gcarleo]

it doesn't crash on my computer

[gcarleo]

what version of netket are you using?

[jeohjee]

Answering to your question, I'm running the code in Kaggle and Colab via command %pip install --quiet netket.
But,I think I was bit unclear in my first message (sorry about that). So the code does not crash after declaring the Hamiltonian but when I try to do an actual run. Here is my full code so you can understand better what I'm doing

import netket as nk
import json
import numpy as np
import matplotlib.pyplot as plt

# sqrt(3)x sqrt(3) unit cell arrangement for kagome:
basis_vectors = [[3/2,np.sqrt(3)/2],[3/2,-1*np.sqrt(3)/2]]
atom_positions = [[0.5,0],[0,0],[0.25,np.sqrt(3)/4],[1.5,0],[1,0],[1.25,np.sqrt(3)/4],[2.5,0],[2,0],[2.25,np.sqrt(3)/4]]
dimensions = [2,2]
M = len(atom_positions)

#Define the graph
graph = nk.graph.Lattice(basis_vectors=basis_vectors,
                         atoms_coord = atom_positions,
                         extent = dimensions,
                         max_neighbor_order = 1,
                         pbc = True
                        )

positions = graph.positions # lattice positions
symmetries = graph.automorphisms()
pos_shape =np.shape(positions)

from jax import numpy as jnp

# Boson Hilbert Space
hi = nk.hilbert.Fock(n_particles=12, N=graph.n_nodes, n_max = 3)


#custom Hamiltonian with finite Peierl substitution

q_vec = np.array([0,0])

def a(site):
    return nk.operator.boson.destroy(hi, site)

def adag(site):
    return nk.operator.boson.create(hi, site)

def nc(site):
    return nk.operator.boson.number(hi, site)

U = 4.0
t = 1.0
ha = nk.operator.LocalOperator(hi, dtype=complex)

for i in range(pos_shape[0]):
    for j in range(pos_shape[0]):
        if i==j:
            continue;
        r_tmp_ij = [positions[i,0] - positions[j,0],positions[i,1] - positions[j,1] ]
        exp_term = np.exp(1j*(q_vec[0]*r_tmp_ij[0] + q_vec[1]*r_tmp_ij[1]))
        ha += t *exp_term* adag(i) * a(j)

    ha += U * nc(i) * (nc(i) -1.0)/2


ma = nk.models.RBMSymm(symmetries=symmetries, alpha=10, param_dtype=complex)

sa = nk.sampler.MetropolisExchange(hilbert = hi, graph=graph, d_max=2)
op = nk.optimizer.Sgd(learning_rate=3e-3)
sr = nk.optimizer.SR(diag_shift=0.01,holomorphic=True)


vs = nk.vqs.MCState(sampler = sa, model = ma, n_discard_per_chain=100, n_samples=3520)

gs = nk.VMC(hamiltonian=ha, optimizer=op, variational_state=vs, preconditioner=sr)
gs.run(n_iter=600, out="test")

The code crashes (due to lack of RAM) in the last line.

[gcarleo]

ok but then this likely crashes because your model does not fit in the memory when computing the local energy. the main issue I see, compared to standard Bose Hubbard, is that here you have all-to-all hoppings, so the off-diagonal elements on which you need to compute the value of the wave functions in the local energy are O(N^2) instead of O(N) as in the case of short range hopping. you can try to set chunk_size to a smaller value

[jeohjee]

Thanks, yeah I also just realized my custom Hamiltonian has a bug (all sites connected). Sorry! I'll try to fix that first =)

[jeohjee]

I fixed it the Hamiltonian to include only NN hoppings and now it works without crashing. Thanks Giuseppe for your time and sorry for such silly issue.