There might be a bug in the methods for constructing Hamiltonian.

[ZongYongyue]

For MPSKitModels.jl@v0.4.0,

H = hubbard_model(Float64, U1Irrep, U1Irrep, FiniteChain(2); t=1, U=8);
H.W

gives

2-element Vector{BlockTensorKit.SparseBlockTensorMap{AbstractTensorMap{Float64, GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, Int64}}, 2, 2}, Float64, GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, Int64}}, 2, 2, 4}}:
 1×1×1×3 SparseBlockTensorMap((⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1)) ⊗ ⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1, (1, 1, 1/2)=>1, (1, 1, -1/2)=>1, (0, 2, 0)=>1))) ← (⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1, (1, 1, 1/2)=>1, (1, 1, -1/2)=>1, (0, 2, 0)=>1)) ⊗ (Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1/2)=>1, (1, 1, -1/2)=>1, (1, -1, 1/2)=>1, (1, -1, -1/2)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1)))):
* Block for sector (FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](0)):
[1.0 1.4142135623730945 1.4142135623730945 -0.0]* Block for sector (FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](1/2)):
[1.0 1.4142135623730945 1.414213562373095 -0.0]* Block for sector (FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](-1/2)):
[1.0 1.4142135623730945 -1.414213562373095 -0.0]* Block for sector (FermionParity(0) ⊠ Irrep[U₁](2) ⊠ Irrep[U₁](0)):
[1.0 -1.414213562373095 1.414213562373095 -0.0]
 3×1×1×1 SparseBlockTensorMap(((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1/2)=>1, (1, 1, -1/2)=>1, (1, -1, 1/2)=>1, (1, -1, -1/2)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1)) ⊗ ⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1, (1, 1, 1/2)=>1, (1, 1, -1/2)=>1, (0, 2, 0)=>1))) ← (⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1, (1, 1, 1/2)=>1, (1, 1, -1/2)=>1, (0, 2, 0)=>1)) ⊗ ⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1)))):
* Block for sector (FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](0)):
[-0.0; -0.7071067811865476; -0.7071067811865476; 1.0;;]* Block for sector (FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](1/2)):
[-0.0; -0.7071067811865476; 0.7071067811865475; 1.0;;]* Block for sector (FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](-1/2)):
[-0.0; -0.7071067811865476; -0.7071067811865475; 1.0;;]* Block for sector (FermionParity(0) ⊠ Irrep[U₁](2) ⊠ Irrep[U₁](0)):
[-0.0; 0.7071067811865475; -0.7071067811865475; 1.0;;]

in which the U value disappeared. But when I construct the same model by FiniteMPOHamiltonian

H = FiniteMPOHamiltonian(fill(pspace, len), mpos...)

it gives:

H.W
2-element Vector{BlockTensorKit.SparseBlockTensorMap{AbstractTensorMap{Float64, GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, Int64}}, 2, 2}, Float64, GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, Int64}}, 2, 2, 4}}:
 1×1×1×6 SparseBlockTensorMap((⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1)) ⊗ ⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 1)=>1, (0, 0, -1)=>1, (1, 1, 0)=>1, (1, -1, 0)=>1))) ← (⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 1)=>1, (0, 0, -1)=>1, (1, 1, 0)=>1, (1, -1, 0)=>1)) ⊗ (Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, -1)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 1)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, -1)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1)))):
* Block for sector (FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1)):
[1.0 -1.4142135623730945 -1.4142135623730945 8.0]* Block for sector (FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](-1)):
[1.0 -1.4142135623730945 -1.4142135623730945 0.0]* Block for sector (FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](0)):
[1.0 -1.414213562373095 1.414213562373095 0.0]* Block for sector (FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](0)):
[1.0 -1.414213562373095 1.414213562373095 0.0]
 6×1×1×1 SparseBlockTensorMap(((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, -1)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 1)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, -1)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1)=>1) ⊕ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1)) ⊗ ⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 1)=>1, (0, 0, -1)=>1, (1, 1, 0)=>1, (1, -1, 0)=>1))) ← (⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 1)=>1, (0, 0, -1)=>1, (1, 1, 0)=>1, (1, -1, 0)=>1)) ⊗ ⊕(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1)))):
* Block for sector (FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1)):
[8.0; -0.7071067811865475; -0.7071067811865475; 1.0;;]* Block for sector (FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](-1)):
[0.0; 0.7071067811865475; 0.7071067811865475; 1.0;;]* Block for sector (FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](0)):
[0.0; -0.7071067811865476; -0.7071067811865476; 1.0;;]* Block for sector (FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](0)):
[0.0; 0.7071067811865476; 0.7071067811865476; 1.0;;]

which seems to be correct (it has the U value 8), at least for the ground state energy:

expectation_value(gs, H)
-0.4721359549995791

is same as the value my ED code gives:-0.4721359549995797.
But the former one gives

expectation_value(gs, H)
-1.9999999999999996

which is the exact energy for U=0 but not for U=8.

[lkdvos]

Thanks for reporting this, I should be able to fix this soon.

As a side note, do you happen to have more of these benchmark results for the Hubbard model? I would love to add them as simple tests to avoid future issues.

[ZongYongyue]

I’d be happy to do so. I’m planning to explore the new versions of MPSKit and MPSKitModels soon, and I’ll add some of the test cases under this issue.

[lkdvos]

This issue should now be solved as of MPSKit v0.12.1, so I will close this.
I would still be incredibly happy with some benchmark results to validate our operators.