Error in fermi-hubbard example
Created by: GiantElephant123
Hello, there.
It seems that the definition of create and annihilation operators are defined reversely in fermi-hubbard model, and results are wrong. link
Activity
Created by: imi-hub
By putting the correct definitions of the creation and annihilation operators, the model seems to converge towards the ground state energy. We can compare the variational state (red curve) with the ED results (black lines).
Could you perhaps explain what type of errors you are getting?
mentioned in merge request !1487 (merged)
-
Created by: imi-hub
PR: !1487 (merged)
-
Created by: GiantElephant123
Thanks for the prompt response and sorry for causing confusion. I thought the wrong hopping term is already crucial so I didn't explain more. I should be more specific:
I actually changed the setting in the original file to 1d hubbard chain model (L=2, D=1, t=2, U=0.01, n_fermion=(1,1)), and find VMC results has an overall minus sign compared to exact results
E_0.I think the reason is the wrong hopping term in orignal file, which reads
H_{\rm wrong}=-t\sum_i a_{i+1}a_i^\dagger + a_ia_{i+1}^\daggerSince
a_ia_{i+1}^\dagger=-a_{i+1}^\dagger a_i, I think there should be a overall minus sign compared to the correct hopping term and this explains the sign problem in VMC and exact values.After I swap the name, the results becomes reasonable:
Strangely, for D=2, 2d hubbard model
H_{\rm wrong}will still gives the right answer and I don't know why. Maybe due to the particle-hole symmetry for half filling system? Or my setting could be wrong?Below are the modified file, I only changed the value of L, D, t, U and n_fermion.
import netket as nk import numpy as np import matplotlib.pyplot as plt import json from netket import experimental as nkx L = 2 # take a 2x2 lattice D = 1 t = 2 # tunneling/hopping U = 0.01 # coulomb # create the graph our fermions can hop on g = nk.graph.Hypercube(length=L, n_dim=D, pbc=True) n_sites = g.n_nodes # create a hilbert space with 2 up and 2 down spins hi = nkx.hilbert.SpinOrbitalFermions(n_sites, s=1 / 2, n_fermions=(1, 1)) # create an operator representing fermi hubbard interactions # -t (i^ j + h.c.) + U (i^ i j^ j) # we will create a helper function to abbreviate the creation, destruction and number operators # each operator has a site and spin projection (sz) in order to find the right position in the hilbert space samples def c(site, sz): return nkx.operator.fermion.create(hi, site, sz=sz) def cdag(site, sz): return nkx.operator.fermion.destroy(hi, site, sz=sz) def nc(site, sz): return nkx.operator.fermion.number(hi, site, sz=sz) up = +1 / 2 down = -1 / 2 ham = 0.0 for sz in (up, down): for u, v in g.edges(): ham += -t * cdag(u, sz) * c(v, sz) - t * cdag(v, sz) * c(u, sz) for u in g.nodes(): ham += U * nc(u, up) * nc(u, down) ed_energies = np.linalg.eigvalsh(ham.to_dense()) print("Hamiltonian =", ham.operator_string()) # metropolis exchange moves fermions around according to a graph # the physical graph has LxL vertices, but the computational basis defined by the # hilbert space contains (2s+1)*L*L occupation numbers # by taking a disjoint copy of the lattice, we can # move the fermions around independently for both spins # and therefore conserve the number of fermions with up and down spin # g.n_nodes == L*L --> disj_graph == 2*L*L disj_graph = nk.graph.disjoint_union(g, g) sa = nk.sampler.MetropolisExchange(hi, graph=disj_graph, n_chains=16) # since the hilbert basis is a set of occupation numbers, we can take a general RBM # we take complex parameters, since it learns sign structures more easily, and for even fermion number, the wave function might be complex ma = nk.models.RBM(alpha=1, param_dtype=complex, use_visible_bias=False) vs = nk.vqs.MCState(sa, ma, n_discard_per_chain=100, n_samples=512) # we will use sgd with Stochastic Reconfiguration opt = nk.optimizer.Sgd(learning_rate=0.01) sr = nk.optimizer.SR(diag_shift=0.1) gs = nk.driver.VMC(ham, opt, variational_state=vs, preconditioner=sr) # now run the optimization # first step will take longer in order to compile exp_name = "fermions_test" gs.run(500, out=exp_name) ############## plot ################# with open(f"{exp_name}.log", "r") as f: data = json.load(f) x = data["Energy"]["iters"] y = data["Energy"]["Mean"]["real"] # plot the energy levels plt.axhline(ed_energies[0], color="red", label="E0") for e in ed_energies[1:]: plt.axhline(e, color="black") plt.plot(x, y, color="blue", label="VMC") plt.xlabel("step") plt.ylabel("E") plt.legend() plt.show() -
Created by: gcarleo
The naming issue was not affecting the number operator, that was defined separately and in the correct way, so since the hopping is invariant under renaming of c and cdagger, there was no actual problem with the energy... On the other hand I'm not sure I understand the context of the minus sign issue you mention
-
Created by: GiantElephant123
I'm afraid the hopping is not invariant under renaming of fermion operator
a_ianda_i^\dagger. According to the original wrong naming and codedef c(site, sz): return nkx.operator.fermion.create(hi, site, sz=sz) def cdag(site, sz): return nkx.operator.fermion.destroy(hi, site, sz=sz) def nc(site, sz): return nkx.operator.fermion.number(hi, site, sz=sz) up = +1 / 2 down = -1 / 2 ham = 0.0 for sz in (up, down): for u, v in g.edges(): ham += -t * cdag(u, sz) * c(v, sz) - t * cdag(v, sz) * c(u, sz) for u in g.nodes(): ham += U * nc(u, up) * nc(u, down)The hamiltonian used is actually given as below:
H_{\rm wrong}=-t\sum_i a_{i}a_{i+1}^\dagger+a_{i+1}a_{i}^\daggerBut the correct hopping term reads
H_{\rm correct}=-t\sum_i a_{i}^\dagger a_{i+1}+a_{i+1}^\dagger a_{i} = t\sum_{i}a_{i+1}a_i^\dagger+a_ia_{i+1}^\dagger={\color{red} -}H_{\rm wrong}Sorry if I have made some stupid mistakes. Maybe the underlying codes still work correctly even the formula seems incorrect, I'm not quite sure.