Feature/ExtendedHubbardMom1D

docs/src/hamiltonians.md

@@ -31,6 +31,7 @@ ExtendedHubbardReal1D
 HubbardMom1D
 BoseHubbardMom1D2C
 HubbardMom1DEP
+ExtendedHubbardMom1D
 ```
 
 ### Harmonic oscillator models

src/Hamiltonians/ExtendedHubbardMom1D.jl

@@ -0,0 +1,117 @@
+"""
+    ExtendedHubbardMom1D(
+        address; 
+        u=1.0, t=1.0, v=1.0, dispersion=hubbard_dispersion, boundary_condition = 0.0
+    )
+
+Implements a one-dimensional extended Hubbard chain, also known as the ``t - V`` model, 
+in momentum space.
+
+```math
+\\hat{H} =  \\sum_{k} ϵ_k n_k + \\frac{1}{2M} \\sum_{kpqr} (u + 2v \\cos(q-p)) a^†_{r} a^†_{q} a_p a_k δ_{r+q,p+k}
+```
+
+# Arguments
+
+* `address`: the starting address, defines number of particles and sites.
+* `u`: the interaction parameter.
+* `t`: the hopping strength.
+* `boundary_condition`: `θ <: Number`: hopping over the boundary incurs a
+    factor ``\\exp(iθ)`` for a hop to the right and ``\\exp(−iθ)`` for a hop to the left.
+
+* `dispersion`: defines ``ϵ_k =``` dispersion(t, k + θ)`
+    - [`hubbard_dispersion`](@ref): ``ϵ_k = -2 (\\Re(t) \\cos(k + θ) + \\Im(t) \\sin(k + θ))``
+    - [`continuum_dispersion`](@ref): ``ϵ_k = \\Re(t) (k + θ)^2 - 2 \\Im(t) (k + θ)``
+
+# See also
+
+* [`HubbardMom1D`](@ref)
+* [`HubbardReal1D`](@ref)
+* [`ExtendedHubbardReal1D`](@ref)
+"""
+struct ExtendedHubbardMom1D{TT,M,AD<:AbstractFockAddress,U,V,T,BOUNDARY_CONDITION} <: AbstractHamiltonian{TT}
+    address::AD # default starting address, should have N particles and M modes
+    ks::SVector{M,TT} # values for k
+    kes::SVector{M,TT} # values for kinetic energy
+end
+
+function ExtendedHubbardMom1D(
+    address::SingleComponentFockAddress;
+    u=1.0, v=1.0, t=1.0, dispersion = hubbard_dispersion, boundary_condition = 0.0
+)
+    M = num_modes(address)
+    U, V, T= promote(float(u), float(v), float(t))
+    step = 2π/M
+    if isodd(M)
+        start = -π*(1+1/M) + step
+    else
+        start = -π + step
+    end
+    kr = range(start; step = step, length = M)
+    ks = SVector{M}(kr)
+    kes = SVector{M}(dispersion.(T , kr .+ (boundary_condition/M)))
+    return ExtendedHubbardMom1D{typeof(U),M,typeof(address),U,V,T,boundary_condition}(address, ks, kes)
+end
+
+function Base.show(io::IO, h::ExtendedHubbardMom1D)
+    compact_addr = repr(h.address, context=:compact => true) # compact print address
+    print(io, "ExtendedHubbardMom1D($(compact_addr); u=$(h.u), v=$(h.v), t=$(h.t), boundary_condition=$(h.boundary_condition))")
+end
+
+function starting_address(h::ExtendedHubbardMom1D)
+    return h.address
+end
+
+dimension(::ExtendedHubbardMom1D, address) = number_conserving_dimension(address)
+
+LOStructure(::Type{<:ExtendedHubbardMom1D{<:Real}}) = IsHermitian()
+
+Base.getproperty(h::ExtendedHubbardMom1D, s::Symbol) = getproperty(h, Val(s))
+Base.getproperty(h::ExtendedHubbardMom1D, ::Val{:ks}) = getfield(h, :ks)
+Base.getproperty(h::ExtendedHubbardMom1D, ::Val{:kes}) = getfield(h, :kes)
+Base.getproperty(h::ExtendedHubbardMom1D, ::Val{:address}) = getfield(h, :address)
+Base.getproperty(h::ExtendedHubbardMom1D{<:Any,<:Any,<:Any,U}, ::Val{:u}) where {U} = U
+Base.getproperty(h::ExtendedHubbardMom1D{<:Any,<:Any,<:Any,<:Any,V}, ::Val{:v}) where {V} = V
+Base.getproperty(h::ExtendedHubbardMom1D{<:Any,<:Any,<:Any,<:Any,<:Any,T}, ::Val{:t}) where {T} = T
+Base.getproperty(h::ExtendedHubbardMom1D{<:Any,<:Any,<:Any,<:Any,<:Any,<:Any,BOUNDARY_CONDITION}, 
+    ::Val{:boundary_condition}) where {BOUNDARY_CONDITION} = BOUNDARY_CONDITION
+
+ks(h::ExtendedHubbardMom1D) = getfield(h, :ks)
+
+# standard interface function
+function num_offdiagonals(ham::ExtendedHubbardMom1D, address::SingleComponentFockAddress)
+    singlies, doublies = num_singly_doubly_occupied_sites(address)
+    return num_offdiagonals(ham, address, singlies, doublies)
+end
+
+# 4-argument version
+@inline function num_offdiagonals(ham::ExtendedHubbardMom1D, ::SingleComponentFockAddress, singlies, doublies)
+    M = num_modes(ham)
+    return singlies * (singlies - 1) * (M - 2) + doublies * (M - 1)
+end
+
+@inline function diagonal_element(h::ExtendedHubbardMom1D{<:Any,M,A}, address::A) where {M,A<:SingleComponentFockAddress}
+    map = OccupiedModeMap(address)
+    return (dot(h.kes, map) + (h.u/ 2M) * momentum_transfer_diagonal(map) 
+        + (h.v/ M) * extended_momentum_transfer_diagonal(map, 2π / M))
+end
+
+@inline function diagonal_element(h::ExtendedHubbardMom1D{<:Any,M,A}, address::A) where {M,A<:FermiFS}
+    map = OccupiedModeMap(address)
+    return dot(h.kes, map) + (h.v/ M) * extended_momentum_transfer_diagonal(map, 2π / M)
+end
+
+@inline function get_offdiagonal(
+    ham::ExtendedHubbardMom1D{<:Any,M,A}, address::A, chosen, map=OccupiedModeMap(address)
+) where {M,A<:SingleComponentFockAddress}
+    address, onproduct,_,_,q = momentum_transfer_excitation(address, chosen, map)
+    return address, ham.u * onproduct / 2M + ham.v * cos(q * 2π / M) * onproduct / M
+end
+
+@inline function get_offdiagonal(
+    ham::ExtendedHubbardMom1D{<:Any,M,A}, address::A, chosen, map=OccupiedModeMap(address)
+) where {M,A<:FermiFS}
+    address, onproduct,_,_,q = momentum_transfer_excitation(address, chosen, map)
+    return address, -ham.v * onproduct * cos(q * 2π / M) / M
+end
+momentum(ham::ExtendedHubbardMom1D) = MomentumMom1D(ham)

src/Hamiltonians/Hamiltonians.jl

@@ -71,7 +71,7 @@ import ..Interfaces: diagonal_element, num_offdiagonals, get_offdiagonal, starti
 export dimension, rayleigh_quotient, momentum
 
 export MatrixHamiltonian
-export HubbardReal1D, HubbardMom1D, ExtendedHubbardReal1D, HubbardRealSpace
+export HubbardReal1D, HubbardMom1D, ExtendedHubbardReal1D, ExtendedHubbardMom1D, HubbardRealSpace
 export HubbardReal1DEP, shift_lattice, shift_lattice_inv
 export HubbardMom1DEP
 export BoseHubbardMom1D2C, BoseHubbardReal1D2C
@@ -103,6 +103,7 @@ include("MatrixHamiltonian.jl")
 
 include("HubbardReal1D.jl")
 include("HubbardReal1DEP.jl")
+include("ExtendedHubbardMom1D.jl")
 include("HubbardMom1D.jl")
 include("HubbardMom1DEP.jl")
 include("HubbardRealSpace.jl")

src/Hamiltonians/HubbardMom1D.jl

@@ -132,7 +132,7 @@ function num_offdiagonals(ham::HubbardMom1D, address::SingleComponentFockAddress
     singlies, doublies = num_singly_doubly_occupied_sites(address)
     return num_offdiagonals(ham, address, singlies, doublies)
 end
-
+num_offdiagonals(ham::HubbardMom1D, address::FermiFS) = 0
 # 4-argument version
 @inline function num_offdiagonals(ham::HubbardMom1D, ::SingleComponentFockAddress, singlies, doublies)
     M = num_modes(ham)
@@ -175,6 +175,10 @@ end
     map = OccupiedModeMap(address)
     return dot(h.kes, map) + momentum_transfer_diagonal(h, map)
 end
+@inline function diagonal_element(h::HubbardMom1D, address::FermiFS)
+    map = OccupiedModeMap(address)
+    return dot(h.kes, map)
+end
 @inline function diagonal_element(h::HubbardMom1D, address::FermiFS2C)
     map_a = OccupiedModeMap(address.components[1])
     map_b = OccupiedModeMap(address.components[2])

src/Hamiltonians/HubbardMom1DEP.jl

@@ -57,8 +57,8 @@ is an external harmonic potential in momentum space,
 * `address`: the starting address, defines number of particles and sites.
 * `u`: the interaction parameter.
 * `t`: the hopping strength.
-* `dispersion`: defines ``ϵ_k =``` t*dispersion(k)`
-    - [`hubbard_dispersion`](@ref): ``ϵ_k = -2( \\Re(t) cos(k) + \\Im(t) sin(k))``
+* `dispersion`: defines ``ϵ_k =``` dispersion(t, k)`
+    - [`hubbard_dispersion`](@ref): ``ϵ_k = -2[\\Re(t) \\cos(k) + \\Im(t) \\sin(k)]``
     - [`continuum_dispersion`](@ref): ``ϵ_k = \\Re(t) k^2 - 2 \\Im(t) k``
 * `v_ho`: strength of the external harmonic oscillator potential ``v_\\mathrm{ho}``.
 

src/Hamiltonians/excitations.jl

@@ -79,10 +79,6 @@ See [`excitation`](@ref), [`OccupiedModeMap`](@ref).
     return excitation(add, dst_indices, src_indices)..., src_modes..., -mom_change
 end
 
-function momentum_transfer_excitation(add::FermiFS, chosen::Integer, map; fold=true)
-    return add, 0.0, 0, 0, 0
-end
-
 @inline function momentum_transfer_excitation(
     add_a, add_b, chosen, map_a, map_b; fold=true
 )
@@ -125,7 +121,7 @@ end
 """
     momentum_transfer_diagonal(map)
 
-The diagonal part of [`momentum_transfer_excitation`](@ref).
+The diagonal part of onsite [`momentum_transfer_excitation`](@ref).
 """
 function momentum_transfer_diagonal(map::BoseOccupiedModeMap)
     onproduct = 0
@@ -139,6 +135,39 @@ function momentum_transfer_diagonal(map::BoseOccupiedModeMap)
     end
     return float(onproduct)
 end
+
+"""
+    extended_momentum_transfer_diagonal(map, step)
+
+The diagonal part of nearest neighbour term [`momentum_transfer_excitation`](@ref) in [`ExtendedHubbardMom1D`](@ref).
+Where `step` is the separation of single-particle momenta in the momentum grid.
+"""
+function extended_momentum_transfer_diagonal(map::OccupiedModeMap, step)
+    onproduct = 0
+    for i in 1:length(map)
+        occ_i = map[i].occnum
+        onproduct += occ_i * (occ_i - 1)
+        for j in 1:i-1
+            occ_j = map[j].occnum
+            onproduct += 2 * occ_i * occ_j * (1 + cos((map[j].mode - map[i].mode)*step))
+        end
+    end
+    return float(onproduct)
+end
+
+function extended_momentum_transfer_diagonal(map::FermiOccupiedModeMap,step)
+    onproduct = 0
+    for i in 1:length(map)
+        occ_i = map[i].occnum
+        onproduct += occ_i * (occ_i - 1)
+        for j in 1:i-1
+            occ_j = map[j].occnum
+            onproduct += 2*occ_i * occ_j * (1 - cos((map[j].mode - map[i].mode)*step))
+        end
+    end
+    return float(onproduct)
+end
+
 function momentum_transfer_diagonal(
     map_a::FermiOccupiedModeMap, map_b::FermiOccupiedModeMap
 )

test/Hamiltonians.jl

@@ -173,6 +173,11 @@ end
         ExtendedHubbardReal1D(BoseFS((1, 0, 0, 0, 1)); u=1.0, v=2.0, t=3.0),
         ExtendedHubbardReal1D(BoseFS(1, 0, 2, 1); u=1 + 0.5im),
         ExtendedHubbardReal1D(BoseFS(1, 0, 2, 1); t=1 + 0.5im),
+        ExtendedHubbardMom1D(BoseFS((1, 0, 0, 0, 1)); u=1.0, v=2.0, t=3.0),
+        ExtendedHubbardMom1D(BoseFS(1, 0, 2, 1); u=1 + 0.5im),
+        ExtendedHubbardMom1D(BoseFS(1, 0, 2, 1); t=1 + 0.5im),
+        ExtendedHubbardMom1D(OccupationNumberFS(1,2,0,0); u=1.0, v=2.0, t=3.0),
+        ExtendedHubbardMom1D(FermiFS(1,1,0,0); u=1.0, v=2.0, t=3.0),
         HubbardRealSpace(BoseFS((1, 2, 3)); u=[1], t=[3]),
         HubbardRealSpace(FermiFS((1, 1, 1, 1, 1, 0, 0, 0)); u=[0], t=[3]),
         HubbardRealSpace(
@@ -672,6 +677,7 @@ end
             for H in (
                 HubbardMom1D(BoseFS((2,2,2)), u=6),
                 ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
+                ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
                 BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0),
             )
                 # GutzwillerSampling with parameter zero is exactly equal to the original H
@@ -702,6 +708,7 @@ end
                 HubbardReal1D(BoseFS((2,2,2)), u=6),
                 HubbardMom1D(BoseFS((2,2,2)), u=6),
                 ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
+                ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
                 # BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0), # multicomponent not implemented for G2RealCorrelator
             )
                 # energy
@@ -785,6 +792,7 @@ end
                 HubbardReal1D(BoseFS((2,2,2)), u=6),
                 HubbardMom1D(BoseFS((2,2,2)), u=6),
                 ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
+                ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
                 # BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0), # multicomponent not implemented for G2RealCorrelator
             )
                 # energy
@@ -854,6 +862,7 @@ end
         for H in (
             HubbardMom1D(BoseFS((2,2,2)), u=6),
             ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
+            ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
             BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0),
         )
             @test_throws ArgumentError Rimu.Hamiltonians.TransformUndoer(H)
@@ -1871,3 +1880,16 @@ end
     @test LOStructure(h) isa IsDiagonal
     @test_throws ArgumentError adjoint(h2)
 end
+
+@testset "Comparison of ExtendedHubbardMom1D with ExtendedHubbardReal1D" begin
+    addr_f = FermiFS{3,6}(0,1,1,1,0,0)
+    addr_b = BoseFS{3,6}(0,0,3,0,0,0)
+    for boundary_condition in ([i*π for i in 0.0:0.2:1.0]...,)
+        HR_f = ExtendedHubbardReal1D(addr_f; boundary_condition)
+        HM_f = ExtendedHubbardMom1D(addr_f; boundary_condition)
+        HR_b = ExtendedHubbardReal1D(addr_b; boundary_condition)
+        HM_b = ExtendedHubbardMom1D(addr_b; boundary_condition)
+        @test round.(eigvals(Matrix(HM_f)), digits=8) ⊆ round.(eigvals(Matrix(HR_f)), digits=8)
+        @test round.(eigvals(Matrix(HM_b)), digits=8) ⊆ round.(eigvals(Matrix(HR_b)), digits=8)
+    end
+end