feature/HubbardMomSpace

docs/src/hamiltonians.md

@@ -26,6 +26,7 @@ ExtendedHubbardReal1D
 ```@docs
 HubbardMom1D
 HubbardMom1DEP
+HubbardMomSpace
 ExtendedHubbardMom1D
 ```
 

src/ExactDiagonalization/ExactDiagonalization.jl

@@ -26,11 +26,11 @@ using CommonSolve: CommonSolve, solve, init
 using VectorInterface: VectorInterface, add
 using OrderedCollections: freeze
 using NamedTupleTools: delete
-using StaticArrays: setindex
+using StaticArrays: setindex, SMatrix
 import Folds
 
 using Rimu: Rimu, DictVectors, Hamiltonians, Interfaces, BitStringAddresses, replace_keys,
-    clean_and_warn_if_others_present, split_keys
+    clean_and_warn_if_others_present, split_keys, GutzwillerSampling, HubbardMomSpace, ExtendedHubbardMom1D
 using ..Interfaces: AbstractDVec, AbstractHamiltonian, AbstractOperator, AdjointUnknown,
     diagonal_element, offdiagonals, starting_address, LOStructure, IsHermitian, operator_column
 using ..BitStringAddresses: AbstractFockAddress, BoseFS, FermiFS, CompositeFS,

src/ExactDiagonalization/basis_set_representation.jl

@@ -113,6 +113,16 @@ function BasisSetRepresentation(
     return _bsr_ensure_symmetry(LOStructure(hamiltonian), hamiltonian, addr; kwargs...)
 end
 
+function BasisSetRepresentation(
+     hamiltonian::Union{HubbardMomSpace{<:Any,<:Any,<:Any,<:Any,<:SMatrix}, ExtendedHubbardMom1D,
+     GutzwillerSampling{<:Any,<:Any,<:Union{HubbardMomSpace{<:Any,<:Any,<:Any,<:Any,<:SMatrix},ExtendedHubbardMom1D}}}, 
+     addr=starting_address(hamiltonian); kwargs...
+)
+    # For symmetry wrappers it is necessary to explicity symmetrise the matrix to
+    # avoid the loss of matrix symmetry due to floating point rounding errors
+    return _bsr_ensure_symmetry(IsHermitian(), hamiltonian, addr; kwargs...)
+end
+
 # default, does not enforce symmetries
 function _bsr_ensure_symmetry(
     ::LOStructure, hamiltonian::AbstractOperator, addr_or_vec;
@@ -136,6 +146,7 @@ function _bsr_ensure_symmetry(
     return BasisSetRepresentation(sparse_matrix, basis, hamiltonian)
 end
 
+
 """
     fix_approx_hermitian!(A; test_approx_symmetry=true, kwargs...)
 
@@ -181,7 +192,9 @@ Returns boolean `true` is the test is passed and `false` if not.
 Furthermore, the matrix `A` is modified to become exactly equal to `½(A + A')` if the test
 is passed.
 """
-function isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...)
+function isapprox_enforce_hermitian!(
+    A::AbstractSparseMatrixCSC{T}; atol=√eps(T), kwargs...
+) where {T}
     # based on `ishermsym()` from `SparseArrays`; relies on `SparseArrays` internals
     # https://github.com/JuliaSparse/SparseArrays.jl/blob/1bae96dc8f9a8ca8b7879eef4cf71e186598e982/src/sparsematrix.jl#L3793
     m, n = size(A)
@@ -202,7 +215,8 @@ function isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...)
             row = rowval[p]
 
             # Ignore stored zeros
-            if iszero(val)
+            if isapprox(val, zero(T); atol, kwargs...)
+                nzval[p] = zero(T)
                 continue
             end
 
@@ -215,7 +229,7 @@ function isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...)
 
             # Diagonal element
             if row == col
-                if isapprox(val, conj(val); kwargs...)
+                if isapprox(val, conj(val); atol, kwargs...)
                     nzval[p] = real(val)
                 else
                     return false
@@ -239,9 +253,11 @@ function isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...)
                 # We therefore "catch up" here while making sure that
                 # the elements are actually zero.
                 while row2 < col
-                    if !iszero(nzval[offset])
+                    if !isapprox(nzval[offset], zero(T); atol, kwargs...)
                         return false
                     end
+                    nzval[offset] = zero(T)
+
                     offset += 1
                     row2 = rowval[offset]
                     tracker[row] += 1
@@ -254,16 +270,13 @@ function isapprox_enforce_hermitian!(A::AbstractSparseMatrixCSC; kwargs...)
 
                 # A[i,j] and A[j,i] exists
                 if row2 == col
-                    if isapprox(val, conj(nzval[offset]); kwargs...)
+                    if isapprox(val, conj(nzval[offset]); atol, kwargs...)
                         val = 1 / 2 * (val + conj(nzval[offset]))
                         nzval[p] = val
                         nzval[offset] = conj(val)
                     else
                         return false
                     end
-                    # if val != conj(nzval[offset])
-                    #     return false
-                    # end
                     tracker[row] += 1
                 end
             end

src/Hamiltonians/ExtendedHubbardMom1D.jl

@@ -29,28 +29,27 @@ in momentum space.
 * [`HubbardReal1D`](@ref)
 * [`ExtendedHubbardReal1D`](@ref)
 """
-struct ExtendedHubbardMom1D{TT,M,AD<:AbstractFockAddress,U,V,T,BOUNDARY_CONDITION} <: AbstractHamiltonian{TT}
+struct ExtendedHubbardMom1D{M,AD<:AbstractFockAddress,U,V,T,BOUNDARY_CONDITION} <: AbstractHamiltonian{Float64}
     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
+    ks::SVector{M,Float64} # values for k
+    kes::SVector{M,Float64} # 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)
+    ks = SVector{M,Float64}(kr)
+    kes = SVector{M,Float64}(dispersion.(t, kr .+ (boundary_condition/M)))
+    return ExtendedHubbardMom1D{M,typeof(address),Float64(u),Float64(v),t,boundary_condition}(address, ks, kes)
 end
 
 function Base.show(io::IO, h::ExtendedHubbardMom1D)
@@ -64,16 +63,16 @@ end
 
 dimension(::ExtendedHubbardMom1D, address) = number_conserving_dimension(address)
 
-LOStructure(::Type{<:ExtendedHubbardMom1D{<:Real}}) = IsHermitian()
+LOStructure(::Type{<:ExtendedHubbardMom1D}) = 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}, 
+Base.getproperty(h::ExtendedHubbardMom1D{<:Any,<:Any,U}, ::Val{:u}) where {U} = U
+Base.getproperty(h::ExtendedHubbardMom1D{<:Any,<:Any,<:Any,V}, ::Val{:v}) where {V} = V
+Base.getproperty(h::ExtendedHubbardMom1D{<:Any,<:Any,<:Any,<:Any,T}, ::Val{:t}) where {T} = T
+Base.getproperty(h::ExtendedHubbardMom1D{<:Any,<:Any,<:Any,<:Any,<:Any,BOUNDARY_CONDITION}, 
     ::Val{:boundary_condition}) where {BOUNDARY_CONDITION} = BOUNDARY_CONDITION
 
 ks(h::ExtendedHubbardMom1D) = getfield(h, :ks)
@@ -90,26 +89,26 @@ end
     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}
+@inline function diagonal_element(h::ExtendedHubbardMom1D{M,A}, address::A) where {M,A<:SingleComponentFockAddress}
     map = occupied_mode_map(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}
+@inline function diagonal_element(h::ExtendedHubbardMom1D{M,A}, address::A) where {M,A<:FermiFS}
     map = occupied_mode_map(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=occupied_mode_map(address)
+    ham::ExtendedHubbardMom1D{M,A}, address::A, chosen, map=occupied_mode_map(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=occupied_mode_map(address)
+    ham::ExtendedHubbardMom1D{M,A}, address::A, chosen, map=occupied_mode_map(address)
 ) where {M,A<:FermiFS}
     address, onproduct,_,_,q = momentum_transfer_excitation(address, chosen, map)
     return address, -ham.v * onproduct * cos(q * 2π / M) / M

src/Hamiltonians/Hamiltonians.jl

@@ -78,7 +78,8 @@ export dimension, rayleigh_quotient, momentum
 
 export IdentityOperator
 export MatrixHamiltonian
-export HubbardReal1D, HubbardMom1D, ExtendedHubbardReal1D, ExtendedHubbardMom1D, HubbardRealSpace
+export HubbardReal1D, HubbardMom1D, ExtendedHubbardReal1D, ExtendedHubbardMom1D
+export HubbardMomSpace, HubbardRealSpace
 export HubbardReal1DEP, shift_lattice, shift_lattice_inv
 export HubbardMom1DEP
 export GutzwillerSampling, GuidingVectorSampling
@@ -126,6 +127,7 @@ include("ExtendedHubbardMom1D.jl")
 include("HubbardMom1D.jl")
 include("HubbardMom1DEP.jl")
 include("HubbardRealSpace.jl")
+include("HubbardMomSpace.jl")
 include("ExtendedHubbardReal1D.jl")
 
 include("FroehlichPolaron.jl")