Energy of Bose Hubbard System does not converge to the exact value.

[alfagalileo]

Hi, actually I am trying to reproduce the Saito's paper with NeuralQuantum.jl. However, I am getting energy values that do not match with the values reported in this paper.
Here is the code

`
nsites = 11
nbosons = 9
U = 2.0
J = 1.0
V = (collect(0:10).-5).^2

hilb = HomogeneousFock( nsites, nbosons )
H = LocalOperator( hilb )

for i in 1:nsites
    ni = number( hilb, i)
    ai = destroy( hilb, mod1(i,nsites) )
    adip1 = create( hilb, mod1(i+1,nsites))
    adim1 = create( hilb, mod1(i-1,nsites))
    H +=  V[i]*ni + U/2*(ni*ni - ni)
    H -= J*ai*adip1
    H -= J*ai*adim1
end

ai = LocalOperator( hilb )
for i in 1:nsites
    ai += destroy( hilb, mod1(i,nsites) )
end

net = RBM( Float64, nsites, 1, af_logcosh )

init_random_pars!( net, sigma=0.1 )
sampler = MetropolisSampler( LocalRule(), 1000, nsites, burn=100 )
algo = SR( ϵ = (0.01), algorithm = sr_cg, precision = 1e-3 )

is = BatchedSampler( net, sampler, H, algo; batch_sz=16 )
add_observable!( is, "ai", ai )
optimizer = Optimisers.Descent(0.01)

Evalues = Float64[];
Eerr = Float64[];

for i=1:epochs
    ldata, prec = sample!(is)
    global ob = compute_observables(is)

    push!( Evalues, real(ldata.mean) )
    push!( Eerr, ldata.error )
    grad = precondition!( prec, algo, i )
    Optimisers.update!( optimizer, net, grad )

    @show i, real( ldata.mean )
end

Evalues

`

I would like to know if I am making a mistake.

[PhilipVinc]

Hey! Thanks for your interest.

So... there are a bunch of issues if you want to reproduce the results by Saito...
1 - Are you sure that you are writing the exact-same Hamiltonian? It might be possible that there are a bunch of factors that don't match...
If you then do (not tested...)

using QuantumOptics
H_qo = SparseOperator(H)
???.groundstate_energy(H_qo)

do you get the correct (Saito) result? Note that you'll probably need to reduce the system size for this to work.
If you indeed get the same...

2 - The network: you are using a RBM/trivial NeuralQuantumState. Saito uses a more complicated, effectively a 2-layer network. You can imitate it (but it won't be exactly the same) by using a chain and a weighted-sum layer.

alpha = 1 # from your example...
ch = Chain(Dense(ComplexF64, N, alpha*N, af_logcosh), WSum(ComplexF64, alpha*N))
net = PureStateAnsatz(ch, N)

3 - Sampling: If I recall correctly Saito performs the computation by fixing the total number of bosons/excitations. (Can you confirm? or am I wrong?). This is particularly important when dealing with bosonic systems... To do that, you need

• to use an Hilbert space with a fixed number of photons. I just made a quick PR to implement this (check the doctoring of HomogeneousFock. you'll need to pass a kwarg excitations=Nbosons).
• use a different markov-chain transition rule. One that preserves the total number of excitations. You should use ExchangeRule(hamiltonian) in lieu of LocalRule.

[PhilipVinc]

Ah, I see! By the way, your way of encoding Nbosons=9 as the local cutoff is... well, not wrong, but... it's completely different from his point of view that there should only be 9 particles in the system.

If you update to the last version of NeuralQuantum you should really just do

cutoff = 3
hilb = HomogeneousFock( nsites, cutoff, excitations=nbosons)

maybe even cutoff = 2 to start with...