From a155418713c211f6a68db878986b01af8b71f0cf Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar Date: Tue, 23 Sep 2025 00:00:42 +1200 Subject: [PATCH 01/25] develop --- src/Hamiltonians/Hamiltonians.jl | 4 +- src/Hamiltonians/HubbardMomSpace.jl | 639 ++++++++++++++++++++++++++++ test/Hamiltonians.jl | 227 ++++++++++ 3 files changed, 869 insertions(+), 1 deletion(-) create mode 100644 src/Hamiltonians/HubbardMomSpace.jl diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl index 4baf1c052..8c299f96e 100644 --- a/src/Hamiltonians/Hamiltonians.jl +++ b/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") diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl new file mode 100644 index 000000000..3aca4638f --- /dev/null +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -0,0 +1,639 @@ +""" + dispersion_MomSpace(ks, geometry, t) + Dispersion relation for a given set of k values, lattice geometry and hopping strengths t. Returns + ``-2(\\sum_{\\bar{k}} \\Re(t_{\\bar{k}}) \\cos(k_{\\bar{k}}) + \\Im(t_{\\bar{k}}) \\sin(k_{\\bar{k}}))`` + where ``\\bar{k}`` runs over all dimensions of the lattice. +""" +function dispersion_MomSpace(ks::SVector{D}, geometry::CubicGrid{D, S}, t::SMatrix) where {D,S} + # Calculate the dispersion relation for a given set of k values and hopping strength t. + C,_ = size(t) + M = prod(S) + kes_mat = zeros(Float64,C,M) + ks_mat = zeros(Float64,D,M) + phase = atan.(imag.(t)./real.(t)) + for j in 1:C + for i in 1:M + mom_val = value_of_mom_mode(M-i+1, ks, geometry) + kes_mat[j,i] = convert(Float64,hub_dis_Mom_Space(t[j,:], mom_val)[1]) + ks_mat[:,i] = mom_val + phase[j,:] + end + end + return SMatrix{C,M,Float64}(kes_mat), SMatrix{D,M,Float64}(ks_mat) +end +function value_of_mom_mode(add_index::Int, ks::SVector{D}, geometry::CubicGrid{D}) where {D} + mom_mode = reverse(geometry[add_index]) + return [ks[i][mode] for (i, mode) in enumerate(mom_mode)] +end + +function hub_dis_Mom_Space(t::SVector, k::Vector) + # Calculate the dispersion relation for a given k value and hopping strength t. + return -2 * (real.(t') * cos.(k) - imag.(t') * sin.(k)) +end + +""" + _mom_hopping(kes, address) + + Calculate the hopping term for a single-component or multi-component Fock state address in momentum space. +""" +@inline function _mom_hopping(kes::SMatrix{C}, address::CompositeFS{C}) where {C} + # Calculate the hopping term for a single component. + comp = address.components + onproduct = 0.0 + for i in 1:C + onproduct += dot(kes[i,:], OccupiedModeMap(comp[i])) + end + return onproduct +end + +@inline _mom_hopping(kes::SMatrix{1}, address::SingleComponentFockAddress) = dot(kes[1,:], OccupiedModeMap(address)) + +""" + mom_transfer_MomSpace(add, chosen, map, g; fold=true) + mom_transfer_MomSpace(add1, add2, chosen, map1, map2, g; fold=true) + +Get the momentum transfer for a given excitation in a same or two different components of a +multi-component Fock state address in momentum space, i.e., `add` or between `add1` and `add2`. +`map`, `map1` and `map2` are the occupied mode maps for the relevant components of the +multi-component Fock state. `chosen` is an integer that determines which excitation and `g` is +a geometry of the lattice. + +""" +@inline function mom_transfer_MomSpace( + add::SingleComponentFockAddress{<:Any, M}, chosen::Int, map::ModeMap, + g::CubicGrid{D,S}; fold=true) where {M,D,S} + # Get the momentum transfer for a given excitation. + singlies = length(map) # number of at least singly occupied modes + + double = chosen - singlies * (singlies - 1) * (M - 2) + + if double > 0 + # Both moves from the same mode. + double, mom_change = fldmod1(double, M - 1) + idx = first(map) # placeholder + for i in map + double -= i.occnum ≥ 2 + if double == 0 + idx = i + break + end + end + src_indices = (idx, idx) + else + # Moves from different modes. + pair, mom_change = fldmod1(chosen, M - 2) + fst, snd = fldmod1(pair, singlies - 1) # where the holes are to be made + if snd < fst # put them in ascending order + f_hole = snd + s_hole = fst + else + f_hole = fst + s_hole = snd + 1 # as we are counting through all singlies + end + src_indices = (map[f_hole], map[s_hole]) + end + src_modes = (src_indices[1].mode, src_indices[2].mode) + src_loc = (g[src_modes[1]], g[src_modes[2]]) + Q = g[mom_change+1] - g[1] + dst_loc = (src_loc[1]+Q, src_loc[2]-Q) + if fold + dst_loc = (mod1.(dst_loc[1], S) , mod1.(dst_loc[2], S)) + if dst_loc == src_loc || reverse(dst_loc) == src_loc + # If the momentum transfer is out of bounds, we return the original address. + Q = g[M] - g[1] + dst_loc = (src_loc[1]+Q, src_loc[2]-Q) + dst_loc = (mod1.(dst_loc[1], S) , mod1.(dst_loc[2], S)) + end + elseif !(all(ones(Int, D) .≤ dst_loc[2]. ≤S) && all(ones(Int, D) .≤ dst_loc[2] .≤ S)) + Q .-= S + dst_loc .= [SRC[1]+Q, SRC[2]-Q] + if !(all(ones(Int, D) .≤ dst_loc[2]. ≤S) && all(ones(Int, D) .≤ dst_loc[2] .≤ S)) + return add, 0.0, src_modes..., -Q + end + end + dst_indices = find_mode(add, (g[dst_loc[1]], g[dst_loc[2]])) + return excitation(add, dst_indices, reverse(src_indices))..., src_modes..., -Q +end + +@inline function mom_transfer_MomSpace( + add1::SingleComponentFockAddress{<:Any, M}, add2::SingleComponentFockAddress{<:Any, M}, + chosen::Int, map1::ModeMap, map2::ModeMap, g::CubicGrid{D,S}; fold=true) where {M,D,S} + # Get the momentum transfer for a given excitation. + singlies = length(map2) + + pair, mom_change = fldmod1(chosen, M - 1) + + f_hole, s_hole = fldmod1(pair, singlies) # where the holes are to be made + src_indices = (map1[f_hole], map2[s_hole]) + src_modes = (src_indices[1].mode, src_indices[2].mode) + src_loc = (g[src_modes[1]], g[src_modes[2]]) + Q = g[mom_change+1] - g[1] + dst_loc = (src_loc[1]+Q, src_loc[2]-Q) + + if fold + dst_loc = (mod1.(dst_loc[1], S) , mod1.(dst_loc[2], S)) + elseif !(all(ones(Int, D) .≤ dst_loc[2]. ≤S) && all(ones(Int, D) .≤ dst_loc[2] .≤ S)) + Q .-= S + dst_loc .= [SRC[1]+Q, SRC[2]-Q] + if !(all(ones(Int, D) .≤ dst_loc[2]. ≤S) && all(ones(Int, D) .≤ dst_loc[2] .≤ S)) + return add, 0.0, src_modes..., -Q + end + end + return excitation(add1, find_mode(add1, (g[dst_loc[1]],)), (src_indices[1],))..., + excitation(add2, find_mode(add2, (g[dst_loc[2]],)), (src_indices[2],))..., src_modes..., -Q +end + +""" + extended_mom_transfer_diagonal(map, g, u, w) + extended_mom_transfer_diagonal(map1, map2) + +Calculate the extended momentum transfer diagonal between a given same component occupied mode map `map` or +between two different component occupied mode maps `map1` and `map2`. `g` is the geometry of the lattice. +`u` and `w` are the on-site and nearest neighbour interaction parameters respectively. + +""" + +@inline function extended_mom_transfer_diag(map::ModeMap, g::CubicGrid{D,S}, u, w) where {D,S} + onproduct = 0 + for i in 1:length(map) + occ_i = map[i].occnum + onproduct += occ_i * (occ_i - 1) * (u/2 + w*D) + for j in 1:i-1 + occ_j = map[j].occnum + q = g[map[i].mode] - g[map[j].mode] + onproduct += 2 * occ_i * occ_j * (u + w*(D + _cosin_sum(q, S))) + end + end + return onproduct +end + +@inline function extended_mom_transfer_diag(map::ModeMap, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} + if iszero(w) + return 0.0 + end + onproduct = 0 + for i in 1:length(map) + occ_i = map[i].occnum + onproduct += occ_i * (occ_i - 1) * (w*D) + for j in 1:i-1 + occ_j = map[j].occnum + q = g[map[i].mode] - g[map[j].mode] + onproduct += 2 * occ_i * occ_j * (D + _cosin_sum(q, S)) + end + end + return onproduct * w +end + +@inline function extended_mom_transfer_diag(map::ModeMap, ::CubicGrid, u, ::Nothing) + if iszero(u) + return 0 + end + 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 += 4 * occ_i * occ_j + end + end + return onproduct * u +end + +@inline function extended_mom_transfer_diag(map::FermiOccupiedModeMap, g::CubicGrid{D,S}, _, w) where {D,S} + if iszero(w) + return 0 + end + 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 + q = g[map[i].mode] - g[map[j].mode] + onproduct += 2 * occ_i * occ_j * (D - _cosin_sum(q, S)) + end + end + return onproduct*w +end + +@inline extended_mom_transfer_diag(::FermiOccupiedModeMap, ::CubicGrid, _, ::Nothing) = 0 + +@inline function extended_mom_transfer_diag(map1::ModeMap, map2::ModeMap) + onproduct = 0 + for i in map1 + occ_i = i.occnum + for j in map2 + occ_j = j.occnum + onproduct += occ_i * occ_j + end + end + return onproduct +end + +@inline extended_mom_transfer_diag(map1::FermiOccupiedModeMap, map2::FermiOccupiedModeMap) = length(map1) * length(map2) + +function _cosin_sum(q::SVector{D}, S::NTuple{D}) where {D} + onproduct = 0.0 + for i in 1:D + onproduct += cos(q[i] * 2π / S[i]) + end + return onproduct +end + +""" + _mom_interactions_dig(component, g) +Calculate the interaction terms for a given component of a multi-component Fock state address in momentum space. +`component` is a tuple of all combination between the pair of relevant components of the multi-component Fock state +for which the interaction term needs to be calculated, + +'''math + \\hat{H}_\\text{int} = \\frac{1}{2}\\sum_{p,q,σ,σ'} V_{σσ'} \\hat{b}^†_{pσ} \\hat{b}^†_{qσ'} \\hat{b}^†_{qσ'} \\hat{b}_{pσ} +''' + +where `V_{σσ}' is the interaction coefficent that depends on `u_{σσ'}' and `w_{σσ'}'. `g` is the geometry of the lattice. + +""" + +function _mom_interactions_dig(component::Tuple, g::CubicGrid{D,S}) where {D,S} + onproduct = 0.0 + for data in component + u = data.u + w = data.w + if !(isnothing(u) && isnothing(w)) + Index = component_index(data) + if Index[1] == Index[2] + # If the occupied modes are the same, we can use the extended mom transfer. + onproduct += extended_mom_transfer_diag(data.occmap1, g, u, w) + else + # Otherwise we need to calculate the interaction between two different occupied modes. + onproduct += extended_mom_transfer_diag(data.occmap1, data.occmap2) * _interaction_parameter_dig(u, w, D) + end + end + end + return onproduct +end + +@inline _interaction_parameter_dig(u::Float64, w::Float64, D::Int) = u + 2 * w * D +@inline _interaction_parameter_dig(::Nothing, w::Float64, D::Int) = 2 * w * D +@inline _interaction_parameter_dig(u::Float64, ::Nothing, _) = u + +""" + HubbardMomSpace(address; geometry=PeriodicBoundaries(M,), t=ones(C, D), u=ones(C, C), w=ones(C, C)) + +Hubbard model in Mom space. Supports single or multi-component Fock state +addresses (with `C` components) and various (rectangular) lattice geometries +in `D` dimensions and of `M` volume. + +```math + \\hat{H} = -\\sum_{k,σ} ϵ_{kσ} n_{kσ} + + \\sum_{p,q,k,σ,σ'} V_{σσ'} a^†_{p+k,σ} a^†_{q-k,σ'} a_{q,σ'} a_{p,σ} +``` +where ``ϵ_{kσ} = -2 (\\sum_{d=1}^{D} \\Re(t_{σ,d}) \\cos(k_d) - \\Im(t_{σ,d}) \\sin(k_d))`` and +``V_{σσ'} = (u_{σσ'}(1- \\frac{δ_{σσ'}}{2}) + w_{σσ'} \\sum_{d=1}^{D} \\cos(q_d))/M`` + +## Address types + +* [`BoseFS`](@ref): Single-component Bose-Hubbard model. +* [`FermiFS`](@ref): Single-component Fermi-Hubbard model. +* [`CompositeFS`](@ref): For multi-component models. + +Note that a single component of fermions cannot interact with itself. A warning +is produced if `address`is incompatible with the interaction parameters `u`. + +## Geometries + +Implemented [`CubicGrid`](@ref)s for keyword `geometry` + +* [`PeriodicBoundaries`](@ref) +* [`HardwallBoundaries`](@ref) +* [`LadderBoundaries`](@ref) + +Default is `geometry=PeriodicBoundaries(M,)`, i.e. a one-dimensional lattice with the +number of sites `M` inferred from the number of modes in `address`. + +## Other parameters + +* `t`: the hopping strengths. Must be a matrix of length `C × D `. The `i`-th and `j`-th element of the + matrix corresponds to the hopping strength of the `i`-th component and `j`-th direction. +* `u`: the on-site interaction parameters. Must be a symmetric matrix. `u[i, j]` + corresponds to the interaction between the `i`-th and `j`-th component. `u[i, i]` + corresponds to the interaction of a component with itself. Note that `u[i,i]` mustadd + be zero for fermionic components. +* `w`: the nearest neighbour interaction parameters. Must be a symmetric matrix. + `w[i, j]` corresponds to the interaction between the `i`-th and `j`-th component. +""" +struct HubbardMomSpace{ + C, # components + D, # dimension + A<:AbstractFockAddress, + G<:CubicGrid, + # The following need to be type params. + KS<:SMatrix{D,<:Any, Float64}, # k values + KES<:SMatrix{C,<:Any,Float64}, + T<:SMatrix{C,D,<:Any}, # hopping strengths + U<:Union{SMatrix{C,C,Float64},Nothing}, + W<:Union{SMatrix{C,C,Float64},Nothing}, +} <: AbstractHamiltonian{Float64} + address::A + ks::KS # k values + kes::KES # kinetic energy values + t::T + u::U # interactions + w::W # nearest neighbour interactions + geometry::G +end + +function HubbardMomSpace( + address::AbstractFockAddress; + geometry::CubicGrid=PeriodicBoundaries((num_modes(address),)), + t=ones(num_components(address), num_dimensions(geometry)), + u=ones(num_components(address), num_components(address)), + w=ones(num_components(address), num_components(address)), +) + C = num_components(address) + D = num_dimensions(geometry) + S = size(geometry) + + # Sanity checks + if prod(size(geometry)) ≠ num_modes(address) + throw(ArgumentError("`geometry` does not have the correct number of sites")) + elseif !(address isa SingleComponentFockAddress || address isa CompositeFS) + throw(ArgumentError( + "unsupported address type detected use `CompositeFS` or `<: SingleComponentFockAddress`" + )) + end + + t_mat = _t_or_v_to_matrix(:t, t, C, D; zero_is_nothing=false) + u_mat = _u_or_w_to_matrix(:u, u, C) + w_mat = _u_or_w_to_matrix(:w, w, C) + + warn_fermi_interaction(address, u_mat) + ks_mat = Array{Float64}[] + for i in eachindex(S) + step = 2π/S[i] + if isodd(S[i]) + start = -π*(1+1/S[i]) + step + else + start = -π + step + end + kr = range(start; step = step, length = S[i]) + if isodd(S[i]) + push!(ks_mat,[j for j in kr]) + else + push!(ks_mat,reverse([j for j in kr])) + end + end + kes, ks = dispersion_MomSpace(SVector{D}(ks_mat), geometry, t_mat) + + return HubbardMomSpace{C,D,typeof(address),typeof(geometry),typeof(ks),typeof(kes),typeof(t_mat), + typeof(u_mat),typeof(w_mat)}( + address, ks, kes, t_mat, u_mat, w_mat, geometry, + ) +end + +LOStructure(::Type{<:HubbardMomSpace}) = IsHermitian() + +function Base.show(io::IO, h::HubbardMomSpace{C}) where {C} + io = IOContext(io, :compact => true) + println(io, "HubbardMomSpace(") + println(io, " ", starting_address(h), ",") + println(io, " geometry = ", h.geometry, ",") + println(io, " t = ", (h.t), ",") + if isnothing(h.u) + println(io, " u = ", zeros(C,C), ",") + else + println(io, " u = ", Float64.(h.u), ",") + end + if isnothing(h.w) + println(io, " w = ", zeros(C,C), ",") + else + println(io, " w = ", Float64.(h.w), ",") + end + print(io, ")") +end + +# Overload equality due to stored potential energy arrays. +Base.:(==)(H::HubbardMomSpace, G::HubbardMomSpace) = all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H))) + +starting_address(h::HubbardMomSpace) = h.address + +dimension(::HubbardMomSpace, address) = number_conserving_dimension(address) + +# offdiaonals =========================================================================== # +# Holds the offdiagonals for a single-component nearest neighbour one-body term. It's +# structured like a matric where the first index determines the occupied site in the adress +# and the second index determines the site the particle will hop to. +struct HubbardMomSpaceComponentData{C,I1,I2,D,G,A,A1,A2,O1,O2} <: AbstractMatrix{Pair{A,Float64}} + geometry::G + parent_address::A + address1::A1 + address2::A2 + u::Union{Float64, Nothing} # interaction strength + w::Union{Float64, Nothing} # nearest neighbour interaction strength + occmap1::O1 + occmap2::O2 + + function HubbardMomSpaceComponentData{C,I1,I2,D}( + geometry::G, + parent::A, + address1::A1, + address2::A2, + u::Union{Float64, Nothing}, + w::Union{Float64, Nothing}, + occmap1::O1=occupied_mode_map(address1), + occmap2::O2=occupied_mode_map(address2), + ) where {C,I1,I2,D,G,A,A1,A2,O1,O2} + return new{C,I1,I2,D,G,A,A1,A2,O1,O2}(geometry, parent, address1, address2, u, w, occmap1, occmap2) + end +end + +function Base.size(data::HubbardMomSpaceComponentData{<:Any,I,I}) where {I} + M= num_modes(data.address1) + s1, d1 = num_singly_doubly_occupied_sites(data.address1) + return s1 * (s1-1) * (M-2) + d1*(M-1) +end + +function Base.size(data::HubbardMomSpaceComponentData) + M= num_modes(data.address1) + s1 = length(data.occmap1) + s2 = length(data.occmap2) + return s1*s2*(M-1) +end + +component_index(::HubbardMomSpaceComponentData{<:Any,I1,I2}) where {I1,I2} = (I1, I2) + +function Base.getindex(data::HubbardMomSpaceComponentData{C,I,I,D}, chosen::Int) where {C,I,D} + geometry = data.geometry + S = size(geometry) + M = prod(S) + map1 = data.occmap1 + new_add, onproduct,_,_,q = mom_transfer_MomSpace(data.address1, chosen, map1, geometry) + if data.parent_address isa CompositeFS + new_parent = BitStringAddresses.update_component( + data.parent_address, new_add, Val(I) + ) + else + new_parent = new_add + end + return new_parent => _interaction_parameter(data.u, data.w, q, S) * onproduct/M +end + +function Base.getindex(data::HubbardMomSpaceComponentData{C,I1,I2,D}, chosen::Int) where {C,I1,I2,D} + geometry = data.geometry + S = size(geometry) + M = prod(S) + new_add1, onproduct1,new_add2, onproduct2,_,_,q = + mom_transfer_MomSpace(data.address1,data.address2,chosen,data.occmap1,data.occmap2,data.geometry) + new_parent = BitStringAddresses.update_component( + data.parent_address, new_add1, Val(I1) + ) + new_parent = BitStringAddresses.update_component( + new_parent, new_add2, Val(I2) + ) + return new_parent => 2 * _interaction_parameter(data.u, data.w, q, S) * onproduct1*onproduct2/M +end + +@inline _interaction_parameter(u::Float64, w::Float64, q::SVector, S::NTuple) = u/2 + w * _cosin_sum(q, S) +@inline _interaction_parameter(::Nothing, w::Float64, q::SVector, S::NTuple) = w * _cosin_sum(q, S) +@inline _interaction_parameter(u::Float64, ::Nothing, ::SVector, ::NTuple) = u/2 + +# column ================================================================================= # +struct HubbardMomSpaceColumn{H,G,A,C<:Tuple} <: AbstractOperatorColumn{A,Float64,H} + hamiltonian::H + geometry::G + address::A + components::C + num_offdiagonals::Int +end + +parent_operator(column::HubbardMomSpaceColumn) = column.hamiltonian +starting_address(column::HubbardMomSpaceColumn) = column.address + +function diagonal_element(col::HubbardMomSpaceColumn) + M = num_modes(col.address) + return _mom_hopping(col.hamiltonian.kes, col.address) + _mom_interactions_dig(col.components, col.geometry)/M +end + +function operator_column(h::HubbardMomSpace{<:Any,<:Any,A,G}, address) where {A,G} + components = _column_components(h, address) + return HubbardMomSpaceColumn{typeof(h),G,A,typeof(components)}( + h, h.geometry, address, components, sum(length, components) + ) +end + +# Collect HubbardMomSpaceComponentData for each component of the address. +@inline function _column_components(h::HubbardMomSpace{1,D}, address::SingleComponentFockAddress) where {D} + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, h.u[1,1], h.w[1,1]),) +end +@inline function _column_components(h::HubbardMomSpace{1,D}, address::FermiFS) where {D} + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, h.u[1,1], h.w[1,1]),) +end +@inline function _column_components(h::HubbardMomSpace, address::CompositeFS) + return _column_components(h, address, address.components, h.u, h.w, Val(1)) +end +@inline function _column_components(::HubbardMomSpace, _, ::Tuple{}, ::Union{SMatrix{0,0},Nothing}, + ::Union{SMatrix{0,0},Nothing}, ::Val) + return () +end +@inline function _column_components( + h::HubbardMomSpace{C,D}, address, (a, as...), + m::Union{SMatrix{N,N},Nothing}, σ::Union{SMatrix{N,N},Nothing}, + ::Val{I1}) where {C,D,N,I1} + # Split the matrix into the column we need now, and the rest. + (u, u_column...) = isnothing(m) ? (nothing, nothing) : Tuple(m[:, 1]) + (w, w_column...) = isnothing(σ) ? (nothing, nothing) : Tuple(σ[:, 1]) + # Type-stable way to subset SMatrix: + m_rest = isnothing(m) ? nothing : SMatrix{N-1,N-1}(view(m, 2:N, 2:N)) + σ_rest = isnothing(σ) ? nothing : SMatrix{N-1,N-1}(view(σ, 2:N, 2:N)) + + return (HubbardMomSpaceComponentData{C,I1,I1,D}(h.geometry, address, a, a, u, w), _mom_interactions_col(h,address,a,as,u_column,w_column,Val(I1),Val(I1+1))..., + _column_components(h, address, as, m_rest, σ_rest, Val(I1+1) )...,) +end + +@inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, + ::Tuple{},::Union{Tuple{},Tuple{Nothing}}, ::Union{Tuple{},Tuple{Nothing}}, ::Val, ::Val) = () + +@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, a::SingleComponentFockAddress, (b,as...)::NTuple{N}, + (u, us...)::NTuple{N}, (w, ws...)::NTuple{N}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} + return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, u, w), _mom_interactions_col(h,address,a, as, us, ws, Val(I1), Val(I2+1))...) +end +@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, a::SingleComponentFockAddress, (b,as...)::NTuple{N}, + m::NTuple{N}, σ::Tuple{Nothing}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} + return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), _mom_interactions_col(h,address,a, as, m[2:N],σ,Val(I1),Val(I2+1))...) +end +@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, a::SingleComponentFockAddress, (b,as...)::Tuple{N}, + m::Tuple{Nothing}, σ::NTuple{N}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} + return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), _mom_interactions_col(h,address,a, as, m, σ[2:N],Val(I1),Val(I2+1))...) +end + +# Split one-dimensional array `index` that indexes over many components simultaneously into +# a two-dimensional one. The dimension of the new index picks the component, while the +# second picks the offdiagonal within the component. + + +function random_offdiagonal(column::HubbardMomSpaceColumn) + random_number = rand(1:column.num_offdiagonals) + component, remainder = _split_component_from_index(column, random_number) + + addr, val = index_apply(getindex, column.components, component, remainder) + return addr, 1/column.num_offdiagonals, val +end + +struct HubbardMomSpaceColumnOffdiagonals{A,G,C<:Tuple} <: AbstractVector{Pair{A,Float64}} + address::A + geometry::G + components::C + num_offdiagonals::Int +end + +function offdiagonals(column::HubbardMomSpaceColumn{<:Any,G,A,C}) where {G,A,C} + return HubbardMomSpaceColumnOffdiagonals{A,G,C}( + column.address, + column.hamiltonian.geometry, + column.components, + column.num_offdiagonals, + ) +end + +@inline function Base.iterate(ods::HubbardMomSpaceColumnOffdiagonals, state=(1,1)) + component_index, chosen = state + if chosen > index_apply(size, ods.components, component_index) + chosen = 1 + component_index += 1 + end + if component_index > length(ods.components) + return nothing + else + result = index_apply( + getindex, ods.components, component_index, chosen + ) + return result, (component_index, chosen + 1) + end +end +Base.size(ods::HubbardMomSpaceColumnOffdiagonals) = (ods.num_offdiagonals,) +Base.eltype(::HubbardMomSpaceColumnOffdiagonals{TT,A}) where {TT,A} = Pair{A,TT} + +function Base.getindex(column::HubbardMomSpaceColumnOffdiagonals, index) + component_index, inner_index = _split_component_from_index(column, index) + return index_apply(getindex, column.components, component_index, inner_index) +end + +########################################################################################################## +# Momentum operator in momentum space +########################################################################################################## + +struct MomentumMomSpace{T,C,D,H<:AbstractHamiltonian{T}} <: AbstractHamiltonian{Vector{T}} + ham::H +end +LOStructure(::Type{MomentumMomSpace{H,T}}) where {H,T <: Real} = IsDiagonal() +num_offdiagonals(ham::MomentumMomSpace, _) = 0 +function diagonal_element(mom::MomentumMomSpace{<:Any,1,D}, address::SingleComponentFockAddress) where D + return [dot(mom.ham.ks[i,:], OccupiedModeMap(address)) for i in 1:D] +end +function diagonal_element(mom::MomentumMomSpace{<:Any,C,D}, address::CompositeFS) where {C,D} + return [sum(dot(mom.ham.ks[i,:], OccupiedModeMap(c)) for c in address.components) for i in 1:D] +end +# fold into (-π, π] +starting_address(mom::MomentumMomSpace) = starting_address(mom.ham) + +momentum(ham::HubbardMomSpace{C,D}) where {C,D} = MomentumMomSpace{Float64,C,D,typeof(ham)}(ham) \ No newline at end of file diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl index 68c278728..c3a1b4e5b 100644 --- a/test/Hamiltonians.jl +++ b/test/Hamiltonians.jl @@ -47,6 +47,15 @@ end FermiFS((1, 1, 1, 1, 0, 0, 0, 0)), ); t=[1, 2], u=[0 3; 3 0], w=[1 0.5; 0.5 1] ), + HubbardMomSpace(BoseFS((1, 2, 3)); u=[1], t=[3], w=[1]), + HubbardMomSpace(FermiFS((1, 1, 1, 1, 1, 0, 0, 0)); u=[0], t=[3]), + HubbardMomSpace(FermiFS((1, 1, 1, 1, 1, 0, 0, 0)); u=[0], t=[3*im]), + HubbardMomSpace( + CompositeFS( + FermiFS((1, 1, 1, 1, 1, 0, 0, 0)), + FermiFS((1, 1, 1, 1, 0, 0, 0, 0)), + ); t=[1, 2], u=[0 3; 3 0], w=[1 0.5; 0.5 1] + ), GutzwillerSampling(HubbardReal1D(BoseFS((1, 2, 3)); u=6 + 2im); g=0.3), GutzwillerSampling(Transcorrelated1D(FermiFS2C((0, 0, 1, 1), (0, 1, 1, 0))); g=0.1), @@ -553,6 +562,224 @@ end end end +@testset "HubbardMomSpace" begin + @testset "Constructor" begin + bose = BoseFS((1, 2, 3, 4, 5, 6)) + @test_throws MethodError HubbardMomSpace(BoseFS{10,10}) + @test_throws ArgumentError HubbardMomSpace(bose; geometry=PeriodicBoundaries(3,3)) + @test_throws ArgumentError HubbardMomSpace( + bose; geometry=PeriodicBoundaries(3,2), t=[1, 2], + ) + @test_throws ArgumentError HubbardMomSpace( + bose; geometry=PeriodicBoundaries(3,2), u=[1 1; 1 1], + ) + @test_throws InexactError HubbardMomSpace( + bose; geometry=PeriodicBoundaries(3,2), u=[1.0im], t=[1.0im] + ) + + comp = CompositeFS(bose, bose) + @test_throws ArgumentError HubbardMomSpace( + comp; geometry=PeriodicBoundaries(3,2), t=[1, 2], u=[1 2; 3 4], + ) + @test_throws ArgumentError HubbardMomSpace( + comp; geometry=PeriodicBoundaries(3,2), t=[1, 2], w=[1 2; 3 4], + ) + @test_throws ArgumentError HubbardMomSpace( + comp; geometry=PeriodicBoundaries(3,2), t=[1, 2], u=[2 2; 2 2; 2 2], + ) + @test_throws ArgumentError HubbardMomSpace( + comp; geometry=PeriodicBoundaries(3,2), v=[1 1; 1 1; 1 1], + ) + @test_throws ArgumentError HubbardMomSpace( + comp; t=[1 2] + ) + + @test_logs (:warn,) HubbardMomSpace(FermiFS((1,0)), u=[2]) + @test_logs (:warn,) HubbardMomSpace( + CompositeFS(BoseFS((1,1)), FermiFS((1,0))); u=[2 2; 2 2] + ) + + H = HubbardMomSpace(comp, t=[1,2], u=[1 2; 2 3]) + @test eval(Meta.parse(repr(H))) == H + end + @testset "1D Bosons (single)" begin + H1 = HubbardMom1D(BoseFS((0, 0, 5, 0, 0, 0)); u=2, t=3) + H2 = HubbardMomSpace(BoseFS((0, 0, 5, 0, 0, 0)); u=[2], t=[3]) + + @test exact_energy(H1) == exact_energy(H2) + end + @testset "1D Bosons (2-component)" begin + add2 = CompositeFS( + BoseFS((0, 0, 3, 0, 0, 0)), + BoseFS((0, 0, 1, 0, 0, 0)), + ) + H2 = HubbardMomSpace(add2, t=[1,4], u=[2 3; 3 0]) + + add3 = CompositeFS( + BoseFS((0, 0, 3, 0, 0, 0)), + FermiFS((0, 0, 1, 0, 0, 0)), + ) + H3 = HubbardMomSpace(add3, t=[1,4], u=[2 3; 3 0]) + + add4 = CompositeFS( + BoseFS((0, 0, 1, 0, 0, 0)), + BoseFS((0, 0, 3, 0, 0, 0)), + ) + H4 = HubbardMomSpace(add4, t=[4,1], u=[0 3; 3 2]) + + add5 = CompositeFS( + FermiFS((1, 0, 0, 0, 0, 0)), + BoseFS((1, 1, 1, 0, 0, 0)), + ) + H5 = HubbardMomSpace(add5, t=[4,1], u=[0 3; 3 2]) + + E2 = exact_energy(H2) + E3 = exact_energy(H3) + E4 = exact_energy(H4) + E5 = exact_energy(H5) + + @test E2 ≈ E3 rtol=0.0001 + @test E3 ≈ E4 rtol=0.0001 + @test E4 ≈ E5 rtol=0.0001 + end + @testset "1D Fermions" begin + H1 = HubbardMomSpace(FermiFS((0, 1, 1, 1, 0, 0)), t=[3.5]) + + # Kinetic energies [+1, -1, -2, -1, +1, +2] can be multiplied by t to get the exact + # energy. + @test exact_energy(H1) ≈ -14 rtol=0.0001 + + # Not interacting, we can sum the parts together. + H2 = HubbardMomSpace( + CompositeFS(FermiFS((0, 1, 1, 1, 1, 0)), FermiFS((0, 0, 1, 1, 0, 0))), + t=[1, 2], u=[0 0; 0 0], + ) + + @test exact_energy(H2) ≈ -3 + -6 rtol=0.0001 + + # Repulsive interactions increase energy. + H3 = HubbardMomSpace( + CompositeFS(FermiFS((0, 1, 1, 1, 1, 0)), FermiFS((0, 0, 1, 1, 0, 0))), + t=[1, 2], u=[0 1; 1 0], + ) + @test exact_energy(H3) > -9 + + # Attractive interactions reduce energy. + H4 = HubbardMomSpace( + CompositeFS(FermiFS((0, 1, 1, 1, 1, 0)), FermiFS((0, 0, 1, 1, 0, 0))), + t=[1, 2], u=[0 -1; -1 0], + ) + @test exact_energy(H4) < -9 + end + @testset "2D Fermions" begin + @testset "2 × 2" begin + p22 = PeriodicBoundaries(2, 2) + @test exact_energy( + HubbardMomSpace(FermiFS(1,0,0,0), geometry=p22, t=[2]) + ) ≈ -8 rtol=0.001 + @test exact_energy( + HubbardMomSpace(FermiFS(1,1,0,0), geometry=p22, t=[2]) + ) ≈ -8 rtol=0.001 + @test exact_energy( + HubbardMomSpace(FermiFS(1,1,1,0), geometry=p22, t=[2]) + ) ≈ -8 rtol=0.001 + @test exact_energy( + HubbardMomSpace(FermiFS(1,1,1,1), geometry=p22, t=[2]) + ) ≈ 0 rtol=0.001 + end + @testset "4 × 4" begin + p44 = PeriodicBoundaries(4, 4) + @test exact_energy( + HubbardMomSpace(FermiFS(0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0); geometry=p44) + ) ≈ -4 rtol=0.001 + @test exact_energy( + HubbardMomSpace(FermiFS(0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0); geometry=p44) + ) ≈ -6 rtol=0.001 + @test exact_energy( + HubbardMomSpace(FermiFS(0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0); geometry=p44) + ) ≈ -8 rtol=0.001 + end + @testset "Two-component" begin + addr = CompositeFS( + FermiFS{3,9}(0,0,0,1,1,1,0,0,0), + FermiFS{2,9}(0,0,0,1,1,0,0,0,0), + ) + H1 = HubbardMomSpace( + addr; + t=[1,2], + u=[0 0; 0 0], + geometry=PeriodicBoundaries(3, 3), + ) + @test exact_energy(H1) ≈ -16 rtol=0.001 + + H2 = HubbardMomSpace( + addr; + t=[1,2], + u=[0 1; 1 0], + geometry=PeriodicBoundaries(3, 3), + ) + @test exact_energy(H2) > -16 + + H3 = HubbardMomSpace( + addr; + t=[1,2], + u=[0 -1; -1 0], + geometry=PeriodicBoundaries(3, 3), + ) + @test exact_energy(H3) < -16 + + H4 = HubbardRealSpace( + addr; + t=[1,2], + u=[0 1; 1 0], + w= [-1 0.5; 0.5 -2], + geometry=PeriodicBoundaries(3, 3), + ) + + H5 = HubbardMomSpace( + addr; + t=[1,2], + u=[0 1; 1 0], + w= [-1 0.5; 0.5 -2], + geometry=PeriodicBoundaries(3, 3), + ) + + eig1 = eigsolve(BasisSetRep(H4; sizelim=1e12).sparse_matrix, 1, :SR)[1][1] + eig2 = eigsolve(BasisSetRep(H5; sizelim=1e12).sparse_matrix, 1, :SR)[1][1] + @test round(real(eig1); digits=10) == round(eig2; digits=10) + end + end + @testset "Complex hopping" begin + address = FermiFS(0, 1, 1, 0) + for H in ( + HubbardMomSpace(address; t=[2.0 + 3im], geometry=CubicGrid(2, 2)), + HubbardMomSpace(address; w=[6], t=[2.0 + 3im], geometry=CubicGrid(4)), + HubbardMomSpace(address; t=[im 2im], geometry=CubicGrid(2, 2)), + ) + @test eltype(H) ≡ Float64 + @test LOStructure(H) ≡ IsHermitian() + test_hamiltonian_structure(H) + end + end + @testset "Nearest neighbour interaction" begin + addr = BoseFS(0,4,0,0) + H1 = HubbardMomSpace(addr; geometry=PeriodicBoundaries(4), w=[2.0]) + H2 = ExtendedHubbardMom1D(addr; v=2.0) + @test Matrix(H1) == Matrix(H2) + addr = FermiFS{2,4}(0,1,1,0) + H1 = HubbardMomSpace(addr; geometry=PeriodicBoundaries(4), w=[-1.0]) + H2 = ExtendedHubbardMom1D(addr; v=-1.0) + @test Matrix(H1) == Matrix(H2) + + addr = BoseFS(0,0,0, 0,2,0, 0,0,0) + H1 = HubbardMomSpace(addr; geometry=PeriodicBoundaries(3, 3), w=[2]) + H2 = HubbardRealSpace(addr; geometry=PeriodicBoundaries(3, 3), w=[2]) + eig1 = eigsolve(BasisSetRep(H1; sizelim=1e12).sparse_matrix, 1, :SR)[1][1] + eig2 = eigsolve(BasisSetRep(H2; sizelim=1e12).sparse_matrix, 1, :SR)[1][1] + @test round(real(eig1); digits=10) == round(eig2; digits=10) + end +end + @testset "Importance sampling" begin @testset "Gutzwiller" begin @testset "Gutzwiller transformation" begin From 0f321a1a0e70bd7b017b54a2215318c8c03e8baf Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 23 Sep 2025 12:50:48 +1200 Subject: [PATCH 02/25] Update HubbardMomSpace.jl --- src/Hamiltonians/HubbardMomSpace.jl | 37 ++++++++++++++--------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl index 3aca4638f..376e952d1 100644 --- a/src/Hamiltonians/HubbardMomSpace.jl +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -143,8 +143,8 @@ end end """ - extended_mom_transfer_diagonal(map, g, u, w) - extended_mom_transfer_diagonal(map1, map2) + extended_mom_transfer_diag(map, g, u, w) + extended_mom_transfer_diag(map1, map2) Calculate the extended momentum transfer diagonal between a given same component occupied mode map `map` or between two different component occupied mode maps `map1` and `map2`. `g` is the geometry of the lattice. @@ -152,7 +152,10 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the """ -@inline function extended_mom_transfer_diag(map::ModeMap, g::CubicGrid{D,S}, u, w) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u, w::Float64) where {D,S} + if iszero(u) && iszero(w) + return 0 + end onproduct = 0 for i in 1:length(map) occ_i = map[i].occnum @@ -166,9 +169,9 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the return onproduct end -@inline function extended_mom_transfer_diag(map::ModeMap, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} if iszero(w) - return 0.0 + return 0 end onproduct = 0 for i in 1:length(map) @@ -183,7 +186,7 @@ end return onproduct * w end -@inline function extended_mom_transfer_diag(map::ModeMap, ::CubicGrid, u, ::Nothing) +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, ::CubicGrid, u, ::Nothing) if iszero(u) return 0 end @@ -255,18 +258,15 @@ where `V_{σσ}' is the interaction coefficent that depends on `u_{σσ'}' and """ function _mom_interactions_dig(component::Tuple, g::CubicGrid{D,S}) where {D,S} - onproduct = 0.0 + onproduct = 0 for data in component - u = data.u - w = data.w - if !(isnothing(u) && isnothing(w)) - Index = component_index(data) - if Index[1] == Index[2] + if !(isnothing(data.u) && isnothing(data.w)) + if component_index(data)[3] # If the occupied modes are the same, we can use the extended mom transfer. - onproduct += extended_mom_transfer_diag(data.occmap1, g, u, w) + onproduct += extended_mom_transfer_diag(data.occmap1, g, data.u, data.w) else # Otherwise we need to calculate the interaction between two different occupied modes. - onproduct += extended_mom_transfer_diag(data.occmap1, data.occmap2) * _interaction_parameter_dig(u, w, D) + onproduct += extended_mom_transfer_diag(data.occmap1, data.occmap2) * _interaction_parameter_dig(data.u, data.w, D) end end end @@ -348,7 +348,7 @@ function HubbardMomSpace( geometry::CubicGrid=PeriodicBoundaries((num_modes(address),)), t=ones(num_components(address), num_dimensions(geometry)), u=ones(num_components(address), num_components(address)), - w=ones(num_components(address), num_components(address)), + w=zeros(num_components(address), num_components(address)), ) C = num_components(address) D = num_dimensions(geometry) @@ -460,7 +460,7 @@ function Base.size(data::HubbardMomSpaceComponentData) return s1*s2*(M-1) end -component_index(::HubbardMomSpaceComponentData{<:Any,I1,I2}) where {I1,I2} = (I1, I2) +component_index(::HubbardMomSpaceComponentData{<:Any,I1,I2}) where {I1,I2} = (I1, I2, I1 == I2) function Base.getindex(data::HubbardMomSpaceComponentData{C,I,I,D}, chosen::Int) where {C,I,D} geometry = data.geometry @@ -510,8 +510,7 @@ parent_operator(column::HubbardMomSpaceColumn) = column.hamiltonian starting_address(column::HubbardMomSpaceColumn) = column.address function diagonal_element(col::HubbardMomSpaceColumn) - M = num_modes(col.address) - return _mom_hopping(col.hamiltonian.kes, col.address) + _mom_interactions_dig(col.components, col.geometry)/M + return _mom_hopping(col.hamiltonian.kes, col.address) + _mom_interactions_dig(col.components, col.geometry)/num_modes(col.address) end function operator_column(h::HubbardMomSpace{<:Any,<:Any,A,G}, address) where {A,G} @@ -636,4 +635,4 @@ end # fold into (-π, π] starting_address(mom::MomentumMomSpace) = starting_address(mom.ham) -momentum(ham::HubbardMomSpace{C,D}) where {C,D} = MomentumMomSpace{Float64,C,D,typeof(ham)}(ham) \ No newline at end of file +momentum(ham::HubbardMomSpace{C,D}) where {C,D} = MomentumMomSpace{Float64,C,D,typeof(ham)}(ham) From 0e9ee7a2fe6cd8317b12fb836096c018c78a8701 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 23 Sep 2025 13:08:18 +1200 Subject: [PATCH 03/25] Update Hamiltonians.jl --- test/Hamiltonians.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl index c3a1b4e5b..362097237 100644 --- a/test/Hamiltonians.jl +++ b/test/Hamiltonians.jl @@ -628,8 +628,8 @@ end H4 = HubbardMomSpace(add4, t=[4,1], u=[0 3; 3 2]) add5 = CompositeFS( - FermiFS((1, 0, 0, 0, 0, 0)), - BoseFS((1, 1, 1, 0, 0, 0)), + FermiFS((0, 0, 1, 0, 0, 0)), + BoseFS((0, 0, 3, 0, 0, 0)), ) H5 = HubbardMomSpace(add5, t=[4,1], u=[0 3; 3 2]) From 0807f25c92e586ae994f90e6a0dc7116412f559e Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 23 Sep 2025 13:09:17 +1200 Subject: [PATCH 04/25] Update HubbardMomSpace.jl --- src/Hamiltonians/HubbardMomSpace.jl | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl index 376e952d1..470d4b03a 100644 --- a/src/Hamiltonians/HubbardMomSpace.jl +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -152,7 +152,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the """ -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u, w::Float64) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u, w) where {D,S} if iszero(u) && iszero(w) return 0 end @@ -522,10 +522,22 @@ end # Collect HubbardMomSpaceComponentData for each component of the address. @inline function _column_components(h::HubbardMomSpace{1,D}, address::SingleComponentFockAddress) where {D} - return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, h.u[1,1], h.w[1,1]),) + if isnothing(h.w) && isnothing(h.u) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, nothing, nothing),) + elseif isnothing(h.w) && !isnothing(h.u) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, h.u[1], nothing),) + elseif !isnothing(h.w) && isnothing(h.u) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, nothing, h.w[1]),) + else + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, h.u[1], h.w[1]),) + end end @inline function _column_components(h::HubbardMomSpace{1,D}, address::FermiFS) where {D} - return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, h.u[1,1], h.w[1,1]),) + if !isnothing(h.w) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, nothing, h.w[1,1]),) + else + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, nothing, nothing),) + end end @inline function _column_components(h::HubbardMomSpace, address::CompositeFS) return _column_components(h, address, address.components, h.u, h.w, Val(1)) From 711f75e66f203096eb83d13081c7c8b1f73394bb Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 23 Sep 2025 16:01:51 +1200 Subject: [PATCH 05/25] Update HubbardMomSpace.jl --- src/Hamiltonians/HubbardMomSpace.jl | 163 +++++++++++++++++----------- 1 file changed, 98 insertions(+), 65 deletions(-) diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl index 470d4b03a..c9141d648 100644 --- a/src/Hamiltonians/HubbardMomSpace.jl +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -14,7 +14,7 @@ function dispersion_MomSpace(ks::SVector{D}, geometry::CubicGrid{D, S}, t::SMatr for j in 1:C for i in 1:M mom_val = value_of_mom_mode(M-i+1, ks, geometry) - kes_mat[j,i] = convert(Float64,hub_dis_Mom_Space(t[j,:], mom_val)[1]) + kes_mat[j,i] = convert(Float64,hub_dis_mom_space(t[j,:], mom_val)[1]) ks_mat[:,i] = mom_val + phase[j,:] end end @@ -25,7 +25,7 @@ function value_of_mom_mode(add_index::Int, ks::SVector{D}, geometry::CubicGrid{D return [ks[i][mode] for (i, mode) in enumerate(mom_mode)] end -function hub_dis_Mom_Space(t::SVector, k::Vector) +function hub_dis_mom_space(t::SVector, k::Vector) # Calculate the dispersion relation for a given k value and hopping strength t. return -2 * (real.(t') * cos.(k) - imag.(t') * sin.(k)) end @@ -48,8 +48,8 @@ end @inline _mom_hopping(kes::SMatrix{1}, address::SingleComponentFockAddress) = dot(kes[1,:], OccupiedModeMap(address)) """ - mom_transfer_MomSpace(add, chosen, map, g; fold=true) - mom_transfer_MomSpace(add1, add2, chosen, map1, map2, g; fold=true) + mom_transfer_mom_space(add, chosen, map, g; fold=true) + mom_transfer_mom_space(add1, add2, chosen, map1, map2, g; fold=true) Get the momentum transfer for a given excitation in a same or two different components of a multi-component Fock state address in momentum space, i.e., `add` or between `add1` and `add2`. @@ -58,7 +58,7 @@ multi-component Fock state. `chosen` is an integer that determines which excitat a geometry of the lattice. """ -@inline function mom_transfer_MomSpace( +@inline function mom_transfer_mom_space( add::SingleComponentFockAddress{<:Any, M}, chosen::Int, map::ModeMap, g::CubicGrid{D,S}; fold=true) where {M,D,S} # Get the momentum transfer for a given excitation. @@ -114,7 +114,7 @@ a geometry of the lattice. return excitation(add, dst_indices, reverse(src_indices))..., src_modes..., -Q end -@inline function mom_transfer_MomSpace( +@inline function mom_transfer_mom_space( add1::SingleComponentFockAddress{<:Any, M}, add2::SingleComponentFockAddress{<:Any, M}, chosen::Int, map1::ModeMap, map2::ModeMap, g::CubicGrid{D,S}; fold=true) where {M,D,S} # Get the momentum transfer for a given excitation. @@ -131,10 +131,10 @@ end if fold dst_loc = (mod1.(dst_loc[1], S) , mod1.(dst_loc[2], S)) - elseif !(all(ones(Int, D) .≤ dst_loc[2]. ≤S) && all(ones(Int, D) .≤ dst_loc[2] .≤ S)) + elseif !(all(x -> 1 ≤ x ≤ S, dst_loc[1]) && all(x -> 1 ≤ x ≤ S, dst_loc[2])) Q .-= S dst_loc .= [SRC[1]+Q, SRC[2]-Q] - if !(all(ones(Int, D) .≤ dst_loc[2]. ≤S) && all(ones(Int, D) .≤ dst_loc[2] .≤ S)) + if !(all(x -> 1 ≤ x ≤ S, dst_loc[1]) && all(x -> 1 ≤ x ≤ S, dst_loc[2])) return add, 0.0, src_modes..., -Q end end @@ -144,7 +144,7 @@ end """ extended_mom_transfer_diag(map, g, u, w) - extended_mom_transfer_diag(map1, map2) + extended_mom_transfer_diag(map1, map2, g, u, w) Calculate the extended momentum transfer diagonal between a given same component occupied mode map `map` or between two different component occupied mode maps `map1` and `map2`. `g` is the geometry of the lattice. @@ -152,7 +152,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the """ -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u, w) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u::Float64, w::Float64) where {D,S} if iszero(u) && iszero(w) return 0 end @@ -169,7 +169,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the return onproduct end -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, ::Nothing, w::Float64) where {D,S} if iszero(w) return 0 end @@ -186,7 +186,7 @@ end return onproduct * w end -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, ::CubicGrid, u, ::Nothing) +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, ::CubicGrid, u::Float64, ::Nothing) if iszero(u) return 0 end @@ -202,7 +202,7 @@ end return onproduct * u end -@inline function extended_mom_transfer_diag(map::FermiOccupiedModeMap, g::CubicGrid{D,S}, _, w) where {D,S} +@inline function extended_mom_transfer_diag(map::FermiOccupiedModeMap, g::CubicGrid{D,S}, _, w::Float64) where {D,S} if iszero(w) return 0 end @@ -221,7 +221,7 @@ end @inline extended_mom_transfer_diag(::FermiOccupiedModeMap, ::CubicGrid, _, ::Nothing) = 0 -@inline function extended_mom_transfer_diag(map1::ModeMap, map2::ModeMap) +@inline function extended_mom_transfer_diag(map1::ModeMap, map2::ModeMap, ::CubicGrid{D}, u, w) where D onproduct = 0 for i in map1 occ_i = i.occnum @@ -230,7 +230,7 @@ end onproduct += occ_i * occ_j end end - return onproduct + return onproduct * _interaction_parameter_diag(u, w, D) end @inline extended_mom_transfer_diag(map1::FermiOccupiedModeMap, map2::FermiOccupiedModeMap) = length(map1) * length(map2) @@ -244,7 +244,7 @@ function _cosin_sum(q::SVector{D}, S::NTuple{D}) where {D} end """ - _mom_interactions_dig(component, g) + _mom_interactions_diag(component, g) Calculate the interaction terms for a given component of a multi-component Fock state address in momentum space. `component` is a tuple of all combination between the pair of relevant components of the multi-component Fock state for which the interaction term needs to be calculated, @@ -257,7 +257,7 @@ where `V_{σσ}' is the interaction coefficent that depends on `u_{σσ'}' and """ -function _mom_interactions_dig(component::Tuple, g::CubicGrid{D,S}) where {D,S} +function _mom_interactions_diag(component::Tuple, g::CubicGrid) onproduct = 0 for data in component if !(isnothing(data.u) && isnothing(data.w)) @@ -266,16 +266,16 @@ function _mom_interactions_dig(component::Tuple, g::CubicGrid{D,S}) where {D,S} onproduct += extended_mom_transfer_diag(data.occmap1, g, data.u, data.w) else # Otherwise we need to calculate the interaction between two different occupied modes. - onproduct += extended_mom_transfer_diag(data.occmap1, data.occmap2) * _interaction_parameter_dig(data.u, data.w, D) + onproduct += extended_mom_transfer_diag(data.occmap1, data.occmap2, g, data.u, data.w) end end end return onproduct end -@inline _interaction_parameter_dig(u::Float64, w::Float64, D::Int) = u + 2 * w * D -@inline _interaction_parameter_dig(::Nothing, w::Float64, D::Int) = 2 * w * D -@inline _interaction_parameter_dig(u::Float64, ::Nothing, _) = u +@inline _interaction_parameter_diag(u::Float64, w::Float64, D::Int) = u + 2 * w * D +@inline _interaction_parameter_diag(::Nothing, w::Float64, D::Int) = 2 * w * D +@inline _interaction_parameter_diag(u::Float64, ::Nothing, _) = u """ HubbardMomSpace(address; geometry=PeriodicBoundaries(M,), t=ones(C, D), u=ones(C, C), w=ones(C, C)) @@ -359,7 +359,8 @@ function HubbardMomSpace( throw(ArgumentError("`geometry` does not have the correct number of sites")) elseif !(address isa SingleComponentFockAddress || address isa CompositeFS) throw(ArgumentError( - "unsupported address type detected use `CompositeFS` or `<: SingleComponentFockAddress`" + "unsupported address type detected use `CompositeFS` or + `<: SingleComponentFockAddress`" )) end @@ -385,8 +386,8 @@ function HubbardMomSpace( end kes, ks = dispersion_MomSpace(SVector{D}(ks_mat), geometry, t_mat) - return HubbardMomSpace{C,D,typeof(address),typeof(geometry),typeof(ks),typeof(kes),typeof(t_mat), - typeof(u_mat),typeof(w_mat)}( + return HubbardMomSpace{C,D,typeof(address),typeof(geometry),typeof(ks),typeof(kes), + typeof(t_mat),typeof(u_mat),typeof(w_mat)}( address, ks, kes, t_mat, u_mat, w_mat, geometry, ) end @@ -413,7 +414,8 @@ function Base.show(io::IO, h::HubbardMomSpace{C}) where {C} end # Overload equality due to stored potential energy arrays. -Base.:(==)(H::HubbardMomSpace, G::HubbardMomSpace) = all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H))) +Base.:(==)(H::HubbardMomSpace, G::HubbardMomSpace) = + all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H))) starting_address(h::HubbardMomSpace) = h.address @@ -423,41 +425,52 @@ dimension(::HubbardMomSpace, address) = number_conserving_dimension(address) # Holds the offdiagonals for a single-component nearest neighbour one-body term. It's # structured like a matric where the first index determines the occupied site in the adress # and the second index determines the site the particle will hop to. -struct HubbardMomSpaceComponentData{C,I1,I2,D,G,A,A1,A2,O1,O2} <: AbstractMatrix{Pair{A,Float64}} +struct HubbardMomSpaceComponentData{ + C,I1,I2,D,G,A,A1,A2,U<:Union{Float64,Nothing},W<:Union{Float64,Nothing},O1,O2 +} <: AbstractMatrix{Pair{A,Float64}} geometry::G parent_address::A address1::A1 address2::A2 - u::Union{Float64, Nothing} # interaction strength - w::Union{Float64, Nothing} # nearest neighbour interaction strength + u::U # interaction strength + w::W # nearest neighbour interaction strength occmap1::O1 occmap2::O2 - function HubbardMomSpaceComponentData{C,I1,I2,D}( geometry::G, parent::A, address1::A1, address2::A2, - u::Union{Float64, Nothing}, - w::Union{Float64, Nothing}, + u::U, + w::W, occmap1::O1=occupied_mode_map(address1), occmap2::O2=occupied_mode_map(address2), - ) where {C,I1,I2,D,G,A,A1,A2,O1,O2} - return new{C,I1,I2,D,G,A,A1,A2,O1,O2}(geometry, parent, address1, address2, u, w, occmap1, occmap2) + ) where {C,I1,I2,D,G,A,A1,A2,U,W,O1,O2} + return new{C,I1,I2,D,G,A,A1,A2,U,W,O1,O2}( + geometry, parent, address1, address2, u, w, occmap1, occmap2 + ) end end function Base.size(data::HubbardMomSpaceComponentData{<:Any,I,I}) where {I} - M= num_modes(data.address1) - s1, d1 = num_singly_doubly_occupied_sites(data.address1) - return s1 * (s1-1) * (M-2) + d1*(M-1) + if isnothing(data.u) && isnothing(data.w) + return 0 + else + M= num_modes(data.address1) + s1, d1 = num_singly_doubly_occupied_sites(data.address1) + return s1 * (s1 - 1) * (M - 2) + d1*(M - 1) + end end function Base.size(data::HubbardMomSpaceComponentData) - M= num_modes(data.address1) - s1 = length(data.occmap1) - s2 = length(data.occmap2) - return s1*s2*(M-1) + if isnothing(data.u) && isnothing(data.w) + return 0 + else + M= num_modes(data.address1) + s1 = length(data.occmap1) + s2 = length(data.occmap2) + return s1 * s2 * (M - 1) + end end component_index(::HubbardMomSpaceComponentData{<:Any,I1,I2}) where {I1,I2} = (I1, I2, I1 == I2) @@ -467,7 +480,7 @@ function Base.getindex(data::HubbardMomSpaceComponentData{C,I,I,D}, chosen::Int) S = size(geometry) M = prod(S) map1 = data.occmap1 - new_add, onproduct,_,_,q = mom_transfer_MomSpace(data.address1, chosen, map1, geometry) + new_add, onproduct,_,_,q = mom_transfer_mom_space(data.address1, chosen, map1, geometry) if data.parent_address isa CompositeFS new_parent = BitStringAddresses.update_component( data.parent_address, new_add, Val(I) @@ -483,7 +496,7 @@ function Base.getindex(data::HubbardMomSpaceComponentData{C,I1,I2,D}, chosen::In S = size(geometry) M = prod(S) new_add1, onproduct1,new_add2, onproduct2,_,_,q = - mom_transfer_MomSpace(data.address1,data.address2,chosen,data.occmap1,data.occmap2,data.geometry) + mom_transfer_mom_space(data.address1,data.address2,chosen,data.occmap1,data.occmap2,data.geometry) new_parent = BitStringAddresses.update_component( data.parent_address, new_add1, Val(I1) ) @@ -510,7 +523,8 @@ parent_operator(column::HubbardMomSpaceColumn) = column.hamiltonian starting_address(column::HubbardMomSpaceColumn) = column.address function diagonal_element(col::HubbardMomSpaceColumn) - return _mom_hopping(col.hamiltonian.kes, col.address) + _mom_interactions_dig(col.components, col.geometry)/num_modes(col.address) + return _mom_hopping(col.hamiltonian.kes, col.address) + + _mom_interactions_diag(col.components, col.geometry)/num_modes(col.address) end function operator_column(h::HubbardMomSpace{<:Any,<:Any,A,G}, address) where {A,G} @@ -521,22 +535,29 @@ function operator_column(h::HubbardMomSpace{<:Any,<:Any,A,G}, address) where {A, end # Collect HubbardMomSpaceComponentData for each component of the address. -@inline function _column_components(h::HubbardMomSpace{1,D}, address::SingleComponentFockAddress) where {D} +@inline function _column_components(h::HubbardMomSpace{1,D}, + address::SingleComponentFockAddress) where {D} if isnothing(h.w) && isnothing(h.u) - return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, nothing, nothing),) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, + address, nothing, nothing),) elseif isnothing(h.w) && !isnothing(h.u) - return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, h.u[1], nothing),) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, + address, h.u[1], nothing),) elseif !isnothing(h.w) && isnothing(h.u) - return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, nothing, h.w[1]),) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, + address, nothing, h.w[1]),) else - return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, h.u[1], h.w[1]),) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, + address, h.u[1], h.w[1]),) end end @inline function _column_components(h::HubbardMomSpace{1,D}, address::FermiFS) where {D} if !isnothing(h.w) - return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, nothing, h.w[1,1]),) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, + address, nothing, h.w[1,1]),) else - return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, address, nothing, nothing),) + return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, + address, nothing, nothing),) end end @inline function _column_components(h::HubbardMomSpace, address::CompositeFS) @@ -557,24 +578,36 @@ end m_rest = isnothing(m) ? nothing : SMatrix{N-1,N-1}(view(m, 2:N, 2:N)) σ_rest = isnothing(σ) ? nothing : SMatrix{N-1,N-1}(view(σ, 2:N, 2:N)) - return (HubbardMomSpaceComponentData{C,I1,I1,D}(h.geometry, address, a, a, u, w), _mom_interactions_col(h,address,a,as,u_column,w_column,Val(I1),Val(I1+1))..., + return (HubbardMomSpaceComponentData{C,I1,I1,D}(h.geometry, address, a, a, u, w), + _mom_interactions_col(h,address,a,as,u_column,w_column,Val(I1),Val(I1+1))..., _column_components(h, address, as, m_rest, σ_rest, Val(I1+1) )...,) end -@inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, - ::Tuple{},::Union{Tuple{},Tuple{Nothing}}, ::Union{Tuple{},Tuple{Nothing}}, ::Val, ::Val) = () - -@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, a::SingleComponentFockAddress, (b,as...)::NTuple{N}, - (u, us...)::NTuple{N}, (w, ws...)::NTuple{N}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} - return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, u, w), _mom_interactions_col(h,address,a, as, us, ws, Val(I1), Val(I2+1))...) +@inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, + ::Tuple{},::Tuple{}, ::Tuple{}, ::Val, ::Val) = () +@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, + a::SingleComponentFockAddress, (b,as...)::NTuple{N}, (u, us...)::NTuple{N}, (w, ws...)::NTuple{N}, + ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} + return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, u, w), + _mom_interactions_col(h,address,a, as, us, ws, Val(I1), Val(I2+1))...) end -@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, a::SingleComponentFockAddress, (b,as...)::NTuple{N}, - m::NTuple{N}, σ::Tuple{Nothing}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} - return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), _mom_interactions_col(h,address,a, as, m[2:N],σ,Val(I1),Val(I2+1))...) + +@inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, + ::Tuple{},::Tuple{}, ::Tuple{Nothing}, ::Val, ::Val) = () +@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, + a::SingleComponentFockAddress, (b,as...)::NTuple{N}, m::NTuple{N}, σ::Tuple{Nothing}, ::Val{I1}, + ::Val{I2}) where {C,D,N,I1,I2} + return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), + _mom_interactions_col(h,address,a, as, m[2:N],σ,Val(I1),Val(I2+1))...) end -@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, a::SingleComponentFockAddress, (b,as...)::Tuple{N}, - m::Tuple{Nothing}, σ::NTuple{N}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} - return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), _mom_interactions_col(h,address,a, as, m, σ[2:N],Val(I1),Val(I2+1))...) + +@inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, + ::Tuple{},::Tuple{Nothing}, ::Tuple{}, ::Val, ::Val) = () +@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, + a::SingleComponentFockAddress, (b,as...)::Tuple{N}, m::Tuple{Nothing}, σ::NTuple{N}, ::Val{I1}, + ::Val{I2}) where {C,D,N,I1,I2} + return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), + _mom_interactions_col(h,address,a, as, m, σ[2:N],Val(I1),Val(I2+1))...) end # Split one-dimensional array `index` that indexes over many components simultaneously into @@ -629,9 +662,9 @@ function Base.getindex(column::HubbardMomSpaceColumnOffdiagonals, index) return index_apply(getindex, column.components, component_index, inner_index) end -########################################################################################################## +######################################################################################################## # Momentum operator in momentum space -########################################################################################################## +######################################################################################################## struct MomentumMomSpace{T,C,D,H<:AbstractHamiltonian{T}} <: AbstractHamiltonian{Vector{T}} ham::H From 8897302e57e6b377b0863ffc1675395ecb75d755 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 23 Sep 2025 16:14:07 +1200 Subject: [PATCH 06/25] Update hamiltonians.md --- docs/src/hamiltonians.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/hamiltonians.md b/docs/src/hamiltonians.md index 673304563..f9c91ea70 100644 --- a/docs/src/hamiltonians.md +++ b/docs/src/hamiltonians.md @@ -26,6 +26,7 @@ ExtendedHubbardReal1D ```@docs HubbardMom1D HubbardMom1DEP +HubbardMomSpace ExtendedHubbardMom1D ``` From db80bfb2487d970f43cf9f5ac2208e47c55d61fd Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 23 Sep 2025 18:34:39 +1200 Subject: [PATCH 07/25] Update HubbardMomSpace.jl --- src/Hamiltonians/HubbardMomSpace.jl | 68 ++++++++++++++++++++--------- 1 file changed, 47 insertions(+), 21 deletions(-) diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl index c9141d648..261a40a75 100644 --- a/src/Hamiltonians/HubbardMomSpace.jl +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -1,10 +1,11 @@ """ - dispersion_MomSpace(ks, geometry, t) + dispersion_mom_space(ks, geometry, t) + Dispersion relation for a given set of k values, lattice geometry and hopping strengths t. Returns ``-2(\\sum_{\\bar{k}} \\Re(t_{\\bar{k}}) \\cos(k_{\\bar{k}}) + \\Im(t_{\\bar{k}}) \\sin(k_{\\bar{k}}))`` where ``\\bar{k}`` runs over all dimensions of the lattice. """ -function dispersion_MomSpace(ks::SVector{D}, geometry::CubicGrid{D, S}, t::SMatrix) where {D,S} +function dispersion_mom_space(ks::SVector{D}, geometry::CubicGrid{D, S}, t::SMatrix) where {D,S} # Calculate the dispersion relation for a given set of k values and hopping strength t. C,_ = size(t) M = prod(S) @@ -152,7 +153,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the """ -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u::Float64, w::Float64) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u, w) where {D,S} if iszero(u) && iszero(w) return 0 end @@ -163,13 +164,14 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the for j in 1:i-1 occ_j = map[j].occnum q = g[map[i].mode] - g[map[j].mode] - onproduct += 2 * occ_i * occ_j * (u + w*(D + _cosin_sum(q, S))) + onproduct += 2 * occ_i * occ_j * (u + w * (D + _cosin_sum(q, S))) end end return onproduct end -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, ::Nothing, w::Float64) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} + if iszero(w) return 0 end @@ -186,7 +188,7 @@ end return onproduct * w end -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, ::CubicGrid, u::Float64, ::Nothing) +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, ::CubicGrid, u, ::Nothing) if iszero(u) return 0 end @@ -202,7 +204,7 @@ end return onproduct * u end -@inline function extended_mom_transfer_diag(map::FermiOccupiedModeMap, g::CubicGrid{D,S}, _, w::Float64) where {D,S} +@inline function extended_mom_transfer_diag(map::FermiOccupiedModeMap, g::CubicGrid{D,S}, _, w) where {D,S} if iszero(w) return 0 end @@ -233,7 +235,10 @@ end return onproduct * _interaction_parameter_diag(u, w, D) end -@inline extended_mom_transfer_diag(map1::FermiOccupiedModeMap, map2::FermiOccupiedModeMap) = length(map1) * length(map2) +@inline function extended_mom_transfer_diag(map1::FermiOccupiedModeMap, map2::FermiOccupiedModeMap, + ::CubicGrid{D}, u, w) where D + return length(map1) * length(map2) * _interaction_parameter_diag(u, w, D) +end function _cosin_sum(q::SVector{D}, S::NTuple{D}) where {D} onproduct = 0.0 @@ -278,7 +283,7 @@ end @inline _interaction_parameter_diag(u::Float64, ::Nothing, _) = u """ - HubbardMomSpace(address; geometry=PeriodicBoundaries(M,), t=ones(C, D), u=ones(C, C), w=ones(C, C)) + HubbardMomSpace(address; geometry=PeriodicBoundaries(M,), t=ones(C, D), u=ones(C, C), w=zeros(C, C)) Hubbard model in Mom space. Supports single or multi-component Fock state addresses (with `C` components) and various (rectangular) lattice geometries @@ -384,7 +389,7 @@ function HubbardMomSpace( push!(ks_mat,reverse([j for j in kr])) end end - kes, ks = dispersion_MomSpace(SVector{D}(ks_mat), geometry, t_mat) + kes, ks = dispersion_mom_space(SVector{D}(ks_mat), geometry, t_mat) return HubbardMomSpace{C,D,typeof(address),typeof(geometry),typeof(ks),typeof(kes), typeof(t_mat),typeof(u_mat),typeof(w_mat)}( @@ -456,9 +461,9 @@ function Base.size(data::HubbardMomSpaceComponentData{<:Any,I,I}) where {I} if isnothing(data.u) && isnothing(data.w) return 0 else - M= num_modes(data.address1) + M = num_modes(data.address1) s1, d1 = num_singly_doubly_occupied_sites(data.address1) - return s1 * (s1 - 1) * (M - 2) + d1*(M - 1) + return s1 * (s1 - 1) * (M - 2) + d1 * (M - 1) end end @@ -466,7 +471,7 @@ function Base.size(data::HubbardMomSpaceComponentData) if isnothing(data.u) && isnothing(data.w) return 0 else - M= num_modes(data.address1) + M = num_modes(data.address1) s1 = length(data.occmap1) s2 = length(data.occmap2) return s1 * s2 * (M - 1) @@ -488,7 +493,7 @@ function Base.getindex(data::HubbardMomSpaceComponentData{C,I,I,D}, chosen::Int) else new_parent = new_add end - return new_parent => _interaction_parameter(data.u, data.w, q, S) * onproduct/M + return new_parent => _interaction_parameter(data.u, data.w, q, S) * onproduct / M end function Base.getindex(data::HubbardMomSpaceComponentData{C,I1,I2,D}, chosen::Int) where {C,I1,I2,D} @@ -503,7 +508,7 @@ function Base.getindex(data::HubbardMomSpaceComponentData{C,I1,I2,D}, chosen::In new_parent = BitStringAddresses.update_component( new_parent, new_add2, Val(I2) ) - return new_parent => 2 * _interaction_parameter(data.u, data.w, q, S) * onproduct1*onproduct2/M + return new_parent => 2 * _interaction_parameter(data.u, data.w, q, S) * onproduct1 * onproduct2 / M end @inline _interaction_parameter(u::Float64, w::Float64, q::SVector, S::NTuple) = u/2 + w * _cosin_sum(q, S) @@ -551,6 +556,7 @@ end address, h.u[1], h.w[1]),) end end + @inline function _column_components(h::HubbardMomSpace{1,D}, address::FermiFS) where {D} if !isnothing(h.w) return (HubbardMomSpaceComponentData{1,1,1,D}(h.geometry, address, address, @@ -560,9 +566,11 @@ end address, nothing, nothing),) end end + @inline function _column_components(h::HubbardMomSpace, address::CompositeFS) return _column_components(h, address, address.components, h.u, h.w, Val(1)) end + @inline function _column_components(::HubbardMomSpace, _, ::Tuple{}, ::Union{SMatrix{0,0},Nothing}, ::Union{SMatrix{0,0},Nothing}, ::Val) return () @@ -584,21 +592,21 @@ end end @inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, - ::Tuple{},::Tuple{}, ::Tuple{}, ::Val, ::Val) = () + ::Tuple{}, ::Tuple{}, ::Tuple{}, ::Val, ::Val) = () @inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, a::SingleComponentFockAddress, (b,as...)::NTuple{N}, (u, us...)::NTuple{N}, (w, ws...)::NTuple{N}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, u, w), - _mom_interactions_col(h,address,a, as, us, ws, Val(I1), Val(I2+1))...) + _mom_interactions_col(h, address, a, as, us, ws, Val(I1), Val(I2+1))...) end @inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, - ::Tuple{},::Tuple{}, ::Tuple{Nothing}, ::Val, ::Val) = () + ::Tuple{}, ::Tuple{}, ::Tuple{Nothing}, ::Val, ::Val) = () @inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, a::SingleComponentFockAddress, (b,as...)::NTuple{N}, m::NTuple{N}, σ::Tuple{Nothing}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), - _mom_interactions_col(h,address,a, as, m[2:N],σ,Val(I1),Val(I2+1))...) + _mom_interactions_col(h, address, a, as, m[2:N], σ, Val(I1), Val(I2+1))...) end @inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, @@ -607,14 +615,22 @@ end a::SingleComponentFockAddress, (b,as...)::Tuple{N}, m::Tuple{Nothing}, σ::NTuple{N}, ::Val{I1}, ::Val{I2}) where {C,D,N,I1,I2} return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), - _mom_interactions_col(h,address,a, as, m, σ[2:N],Val(I1),Val(I2+1))...) + _mom_interactions_col(h, address, a, as, m, σ[2:N], Val(I1), Val(I2+1))...) +end + +@inline _mom_interactions_col(::HubbardMomSpace, ::AbstractFockAddress, ::SingleComponentFockAddress, + ::Tuple{},::Tuple{Nothing}, ::Tuple{Nothing}, ::Val, ::Val) = () +@inline function _mom_interactions_col(h::HubbardMomSpace{C,D}, address::AbstractFockAddress, + a::SingleComponentFockAddress, (b,as...)::Tuple{N}, m::Tuple{Nothing}, σ::Tuple{Nothing}, ::Val{I1}, + ::Val{I2}) where {C,D,N,I1,I2} + return (HubbardMomSpaceComponentData{C,I1,I2,D}(h.geometry, address, a, b, m[1], σ[1]), + _mom_interactions_col(h, address, a, as, m, σ, Val(I1), Val(I2+1))...) end # Split one-dimensional array `index` that indexes over many components simultaneously into # a two-dimensional one. The dimension of the new index picks the component, while the # second picks the offdiagonal within the component. - function random_offdiagonal(column::HubbardMomSpaceColumn) random_number = rand(1:column.num_offdiagonals) component, remainder = _split_component_from_index(column, random_number) @@ -641,6 +657,7 @@ end @inline function Base.iterate(ods::HubbardMomSpaceColumnOffdiagonals, state=(1,1)) component_index, chosen = state + component_index = _test_interaction_coeff(ods, component_index) if chosen > index_apply(size, ods.components, component_index) chosen = 1 component_index += 1 @@ -654,6 +671,15 @@ end return result, (component_index, chosen + 1) end end + +@inline function _test_interaction_coeff(ods::HubbardMomSpaceColumnOffdiagonals, component_index::Int) + if index_apply(size, ods.components, component_index) == 0 && component_index < length(ods.components) + return _test_interaction_coeff(ods, component_index + 1) + else + return component_index + end +end + Base.size(ods::HubbardMomSpaceColumnOffdiagonals) = (ods.num_offdiagonals,) Base.eltype(::HubbardMomSpaceColumnOffdiagonals{TT,A}) where {TT,A} = Pair{A,TT} From aa7aefd903fb89065e66fe6877d5c7c08e4db56e Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Thu, 25 Sep 2025 11:08:27 +1200 Subject: [PATCH 08/25] Update HubbardMomSpace.jl --- src/Hamiltonians/HubbardMomSpace.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl index 261a40a75..3cbea1aff 100644 --- a/src/Hamiltonians/HubbardMomSpace.jl +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -41,12 +41,12 @@ end comp = address.components onproduct = 0.0 for i in 1:C - onproduct += dot(kes[i,:], OccupiedModeMap(comp[i])) + onproduct += dot(kes[i,:], occupied_mode_map(comp[i])) end return onproduct end -@inline _mom_hopping(kes::SMatrix{1}, address::SingleComponentFockAddress) = dot(kes[1,:], OccupiedModeMap(address)) +@inline _mom_hopping(kes::SMatrix{1}, address::SingleComponentFockAddress) = dot(kes[1,:], occupied_mode_map(address)) """ mom_transfer_mom_space(add, chosen, map, g; fold=true) @@ -153,7 +153,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the """ -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u, w) where {D,S} +@inline function extended_mom_transfer_diag(map::Boseoccupied_mode_map, g::CubicGrid{D,S}, u, w) where {D,S} if iszero(u) && iszero(w) return 0 end @@ -170,7 +170,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the return onproduct end -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} +@inline function extended_mom_transfer_diag(map::Boseoccupied_mode_map, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} if iszero(w) return 0 @@ -188,7 +188,7 @@ end return onproduct * w end -@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, ::CubicGrid, u, ::Nothing) +@inline function extended_mom_transfer_diag(map::Boseoccupied_mode_map, ::CubicGrid, u, ::Nothing) if iszero(u) return 0 end @@ -204,7 +204,7 @@ end return onproduct * u end -@inline function extended_mom_transfer_diag(map::FermiOccupiedModeMap, g::CubicGrid{D,S}, _, w) where {D,S} +@inline function extended_mom_transfer_diag(map::Fermioccupied_mode_map, g::CubicGrid{D,S}, _, w) where {D,S} if iszero(w) return 0 end @@ -221,7 +221,7 @@ end return onproduct*w end -@inline extended_mom_transfer_diag(::FermiOccupiedModeMap, ::CubicGrid, _, ::Nothing) = 0 +@inline extended_mom_transfer_diag(::Fermioccupied_mode_map, ::CubicGrid, _, ::Nothing) = 0 @inline function extended_mom_transfer_diag(map1::ModeMap, map2::ModeMap, ::CubicGrid{D}, u, w) where D onproduct = 0 @@ -235,7 +235,7 @@ end return onproduct * _interaction_parameter_diag(u, w, D) end -@inline function extended_mom_transfer_diag(map1::FermiOccupiedModeMap, map2::FermiOccupiedModeMap, +@inline function extended_mom_transfer_diag(map1::Fermioccupied_mode_map, map2::Fermioccupied_mode_map, ::CubicGrid{D}, u, w) where D return length(map1) * length(map2) * _interaction_parameter_diag(u, w, D) end @@ -698,10 +698,10 @@ end LOStructure(::Type{MomentumMomSpace{H,T}}) where {H,T <: Real} = IsDiagonal() num_offdiagonals(ham::MomentumMomSpace, _) = 0 function diagonal_element(mom::MomentumMomSpace{<:Any,1,D}, address::SingleComponentFockAddress) where D - return [dot(mom.ham.ks[i,:], OccupiedModeMap(address)) for i in 1:D] + return [dot(mom.ham.ks[i,:], occupied_mode_map(address)) for i in 1:D] end function diagonal_element(mom::MomentumMomSpace{<:Any,C,D}, address::CompositeFS) where {C,D} - return [sum(dot(mom.ham.ks[i,:], OccupiedModeMap(c)) for c in address.components) for i in 1:D] + return [sum(dot(mom.ham.ks[i,:], occupied_mode_map(c)) for c in address.components) for i in 1:D] end # fold into (-π, π] starting_address(mom::MomentumMomSpace) = starting_address(mom.ham) From 5233865ff75d6f06d6f1bb412446a939da4f17b6 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Thu, 25 Sep 2025 12:14:09 +1200 Subject: [PATCH 09/25] Update HubbardMomSpace.jl --- src/Hamiltonians/HubbardMomSpace.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl index 3cbea1aff..aa1b84c5e 100644 --- a/src/Hamiltonians/HubbardMomSpace.jl +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -153,7 +153,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the """ -@inline function extended_mom_transfer_diag(map::Boseoccupied_mode_map, g::CubicGrid{D,S}, u, w) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u, w) where {D,S} if iszero(u) && iszero(w) return 0 end @@ -170,7 +170,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the return onproduct end -@inline function extended_mom_transfer_diag(map::Boseoccupied_mode_map, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} if iszero(w) return 0 @@ -188,7 +188,7 @@ end return onproduct * w end -@inline function extended_mom_transfer_diag(map::Boseoccupied_mode_map, ::CubicGrid, u, ::Nothing) +@inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, ::CubicGrid, u, ::Nothing) if iszero(u) return 0 end @@ -204,7 +204,7 @@ end return onproduct * u end -@inline function extended_mom_transfer_diag(map::Fermioccupied_mode_map, g::CubicGrid{D,S}, _, w) where {D,S} +@inline function extended_mom_transfer_diag(map::FermiOccupiedModeMap, g::CubicGrid{D,S}, _, w) where {D,S} if iszero(w) return 0 end @@ -221,7 +221,7 @@ end return onproduct*w end -@inline extended_mom_transfer_diag(::Fermioccupied_mode_map, ::CubicGrid, _, ::Nothing) = 0 +@inline extended_mom_transfer_diag(::FermiOccupiedModeMap, ::CubicGrid, _, ::Nothing) = 0 @inline function extended_mom_transfer_diag(map1::ModeMap, map2::ModeMap, ::CubicGrid{D}, u, w) where D onproduct = 0 @@ -235,7 +235,7 @@ end return onproduct * _interaction_parameter_diag(u, w, D) end -@inline function extended_mom_transfer_diag(map1::Fermioccupied_mode_map, map2::Fermioccupied_mode_map, +@inline function extended_mom_transfer_diag(map1::FermiOccupiedModeMap, map2::FermiOccupiedModeMap, ::CubicGrid{D}, u, w) where D return length(map1) * length(map2) * _interaction_parameter_diag(u, w, D) end From 541831918c2beb5e403dd92d2d5d2a60994671cd Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Thu, 25 Sep 2025 13:07:39 +1200 Subject: [PATCH 10/25] Update ExactDiagonalization.jl --- src/ExactDiagonalization/ExactDiagonalization.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ExactDiagonalization/ExactDiagonalization.jl b/src/ExactDiagonalization/ExactDiagonalization.jl index 79ecdee8b..e3f667dea 100644 --- a/src/ExactDiagonalization/ExactDiagonalization.jl +++ b/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, HubbardMomSpace using ..Interfaces: AbstractDVec, AbstractHamiltonian, AbstractOperator, AdjointUnknown, diagonal_element, offdiagonals, starting_address, LOStructure, IsHermitian, operator_column using ..BitStringAddresses: AbstractFockAddress, BoseFS, FermiFS, CompositeFS, From 65b86c16b1fa8a906b8469c84662b7c3d1e5ed0b Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Thu, 25 Sep 2025 13:11:40 +1200 Subject: [PATCH 11/25] Update basis_set_representation.jl - enforcing Hermitian property in ``BasisSetRepresentation(HubbardMomSpace)`` which seems to be lost due to the use of trigonometric functions during the build of the matrix. --- src/ExactDiagonalization/basis_set_representation.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/ExactDiagonalization/basis_set_representation.jl b/src/ExactDiagonalization/basis_set_representation.jl index beb1e882e..6a297d5b5 100644 --- a/src/ExactDiagonalization/basis_set_representation.jl +++ b/src/ExactDiagonalization/basis_set_representation.jl @@ -113,6 +113,15 @@ function BasisSetRepresentation( return _bsr_ensure_symmetry(LOStructure(hamiltonian), hamiltonian, addr; kwargs...) end +function BasisSetRepresentation( + hamiltonian::HubbardMomSpace{<:Any,<:Any,<:Any,<:Any,<:SMatrix}, 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(LOStructure(hamiltonian), hamiltonian, addr; kwargs...) +end + # default, does not enforce symmetries function _bsr_ensure_symmetry( ::LOStructure, hamiltonian::AbstractOperator, addr_or_vec; @@ -136,6 +145,7 @@ function _bsr_ensure_symmetry( return BasisSetRepresentation(sparse_matrix, basis, hamiltonian) end + """ fix_approx_hermitian!(A; test_approx_symmetry=true, kwargs...) From d4f8cee44e03820e9e3589a184807ff487361d1e Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Fri, 26 Sep 2025 12:36:49 +1200 Subject: [PATCH 12/25] Update basis_set_representation.jl --- .../basis_set_representation.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/ExactDiagonalization/basis_set_representation.jl b/src/ExactDiagonalization/basis_set_representation.jl index 6a297d5b5..0c07cee0f 100644 --- a/src/ExactDiagonalization/basis_set_representation.jl +++ b/src/ExactDiagonalization/basis_set_representation.jl @@ -191,7 +191,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) @@ -212,7 +214,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 @@ -225,7 +228,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 @@ -249,9 +252,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 @@ -264,16 +269,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 From 89a52945c8c749804a71dde65daf93f350c26fc5 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Fri, 26 Sep 2025 12:38:58 +1200 Subject: [PATCH 13/25] Update ExactDiagonalization.jl --- test/ExactDiagonalization.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/ExactDiagonalization.jl b/test/ExactDiagonalization.jl index d2c872832..8f03083d1 100644 --- a/test/ExactDiagonalization.jl +++ b/test/ExactDiagonalization.jl @@ -134,6 +134,22 @@ using SparseArrays @test even_sm ≈ even_m # still approximately the same! end + @testset "isapprox_enforce_hermitian!" begin + matrix = sprand(100, 100, 0.2) + matrix .+= matrix' + for _ in 1:1000 + matrix[rand(1:100), rand(1:100)] += 1e-9 + end + + matrix1 = copy(matrix) + @test !ishermitian(matrix1) + @test isapprox_enforce_hermitian!(matrix1) + @test ishermitian(matrix1) + + matrix2 = copy(matrix) + @test !isapprox_enforce_hermitian!(matrix2; atol=1e-12) + end + @testset "basis-only" begin m = 5 n = 5 From 064a7dde110b28bed06fc8598dfc7c994d4fad7c Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Fri, 26 Sep 2025 13:28:17 +1200 Subject: [PATCH 14/25] Update ExactDiagonalization.jl --- src/ExactDiagonalization/ExactDiagonalization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ExactDiagonalization/ExactDiagonalization.jl b/src/ExactDiagonalization/ExactDiagonalization.jl index e3f667dea..fad85ef10 100644 --- a/src/ExactDiagonalization/ExactDiagonalization.jl +++ b/src/ExactDiagonalization/ExactDiagonalization.jl @@ -30,7 +30,7 @@ using StaticArrays: setindex, SMatrix import Folds using Rimu: Rimu, DictVectors, Hamiltonians, Interfaces, BitStringAddresses, replace_keys, - clean_and_warn_if_others_present, split_keys, HubbardMomSpace + clean_and_warn_if_others_present, split_keys, HubbardMomSpace, ExtendedHubbardReal1D using ..Interfaces: AbstractDVec, AbstractHamiltonian, AbstractOperator, AdjointUnknown, diagonal_element, offdiagonals, starting_address, LOStructure, IsHermitian, operator_column using ..BitStringAddresses: AbstractFockAddress, BoseFS, FermiFS, CompositeFS, From cf1113a6c7bd9f33e6610a86843e9792b73e5b92 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Fri, 26 Sep 2025 13:30:45 +1200 Subject: [PATCH 15/25] Update basis_set_representation.jl --- src/ExactDiagonalization/basis_set_representation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ExactDiagonalization/basis_set_representation.jl b/src/ExactDiagonalization/basis_set_representation.jl index 0c07cee0f..dd82180b8 100644 --- a/src/ExactDiagonalization/basis_set_representation.jl +++ b/src/ExactDiagonalization/basis_set_representation.jl @@ -114,8 +114,8 @@ function BasisSetRepresentation( end function BasisSetRepresentation( - hamiltonian::HubbardMomSpace{<:Any,<:Any,<:Any,<:Any,<:SMatrix}, addr=starting_address(hamiltonian); - kwargs... + hamiltonian::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 From e1b0fa23dcfc2806eb0092afb6044c409007cb6e Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Fri, 26 Sep 2025 13:31:13 +1200 Subject: [PATCH 16/25] Update ExactDiagonalization.jl --- src/ExactDiagonalization/ExactDiagonalization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ExactDiagonalization/ExactDiagonalization.jl b/src/ExactDiagonalization/ExactDiagonalization.jl index fad85ef10..9cdce3854 100644 --- a/src/ExactDiagonalization/ExactDiagonalization.jl +++ b/src/ExactDiagonalization/ExactDiagonalization.jl @@ -30,7 +30,7 @@ using StaticArrays: setindex, SMatrix import Folds using Rimu: Rimu, DictVectors, Hamiltonians, Interfaces, BitStringAddresses, replace_keys, - clean_and_warn_if_others_present, split_keys, HubbardMomSpace, ExtendedHubbardReal1D + clean_and_warn_if_others_present, split_keys, HubbardMomSpace, ExtendedHubbardMom1D using ..Interfaces: AbstractDVec, AbstractHamiltonian, AbstractOperator, AdjointUnknown, diagonal_element, offdiagonals, starting_address, LOStructure, IsHermitian, operator_column using ..BitStringAddresses: AbstractFockAddress, BoseFS, FermiFS, CompositeFS, From 2f81c7252fdf453870cdf9fad6ad2a5dd8dbb944 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Fri, 26 Sep 2025 13:36:53 +1200 Subject: [PATCH 17/25] Update Hamiltonians.jl --- test/Hamiltonians.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl index 362097237..bbfc289ba 100644 --- a/test/Hamiltonians.jl +++ b/test/Hamiltonians.jl @@ -587,9 +587,6 @@ end @test_throws ArgumentError HubbardMomSpace( comp; geometry=PeriodicBoundaries(3,2), t=[1, 2], u=[2 2; 2 2; 2 2], ) - @test_throws ArgumentError HubbardMomSpace( - comp; geometry=PeriodicBoundaries(3,2), v=[1 1; 1 1; 1 1], - ) @test_throws ArgumentError HubbardMomSpace( comp; t=[1 2] ) From 6f17d1d5504a701b12cf7bc6bb1e95eb63fbc34d Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Sat, 4 Oct 2025 20:17:28 +1300 Subject: [PATCH 18/25] Update basis_set_representation.jl --- src/ExactDiagonalization/basis_set_representation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ExactDiagonalization/basis_set_representation.jl b/src/ExactDiagonalization/basis_set_representation.jl index dd82180b8..c24b6da6f 100644 --- a/src/ExactDiagonalization/basis_set_representation.jl +++ b/src/ExactDiagonalization/basis_set_representation.jl @@ -119,7 +119,7 @@ function BasisSetRepresentation( ) # 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(LOStructure(hamiltonian), hamiltonian, addr; kwargs...) + return _bsr_ensure_symmetry(IsHermitian(), hamiltonian, addr; kwargs...) end # default, does not enforce symmetries From 1f5777dba83c8df8a073e98a5444dd602bc5a957 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Sun, 5 Oct 2025 12:15:47 +1300 Subject: [PATCH 19/25] Update ExtendedHubbardMom1D.jl # Fixing the bug - ``ExtendedHubbardMom1D`` is always real. Previously, it was also allowed to be typed as ``AbstractHamiltonian{TT}`` where ``TT`` was allowed to be `ComplexF64`, which was not required. --- src/Hamiltonians/ExtendedHubbardMom1D.jl | 29 ++++++++++++------------ 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/Hamiltonians/ExtendedHubbardMom1D.jl b/src/Hamiltonians/ExtendedHubbardMom1D.jl index 2af72ad62..eab5421a9 100644 --- a/src/Hamiltonians/ExtendedHubbardMom1D.jl +++ b/src/Hamiltonians/ExtendedHubbardMom1D.jl @@ -29,10 +29,10 @@ 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( @@ -40,7 +40,6 @@ function ExtendedHubbardMom1D( 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 @@ -48,9 +47,9 @@ function ExtendedHubbardMom1D( 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) @@ -70,10 +69,10 @@ 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 From 6c2370d200a654874da2cdc79b0fa8d9f7781afd Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 7 Oct 2025 15:48:57 +1300 Subject: [PATCH 20/25] Update ExactDiagonalization.jl --- src/ExactDiagonalization/ExactDiagonalization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ExactDiagonalization/ExactDiagonalization.jl b/src/ExactDiagonalization/ExactDiagonalization.jl index 9cdce3854..52533f3aa 100644 --- a/src/ExactDiagonalization/ExactDiagonalization.jl +++ b/src/ExactDiagonalization/ExactDiagonalization.jl @@ -30,7 +30,7 @@ using StaticArrays: setindex, SMatrix import Folds using Rimu: Rimu, DictVectors, Hamiltonians, Interfaces, BitStringAddresses, replace_keys, - clean_and_warn_if_others_present, split_keys, HubbardMomSpace, ExtendedHubbardMom1D + 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, From 358d345cced0e1d97ff5e579cc536ac6b7cea0be Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 7 Oct 2025 15:49:55 +1300 Subject: [PATCH 21/25] Update basis_set_representation.jl --- src/ExactDiagonalization/basis_set_representation.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ExactDiagonalization/basis_set_representation.jl b/src/ExactDiagonalization/basis_set_representation.jl index c24b6da6f..5628d2089 100644 --- a/src/ExactDiagonalization/basis_set_representation.jl +++ b/src/ExactDiagonalization/basis_set_representation.jl @@ -114,7 +114,8 @@ function BasisSetRepresentation( end function BasisSetRepresentation( - hamiltonian::Union{HubbardMomSpace{<:Any,<:Any,<:Any,<:Any,<:SMatrix}, ExtendedHubbardMom1D}, + 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 From 19656be6141c15de42f8bf2346af070e3e17e45b Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 7 Oct 2025 16:38:12 +1300 Subject: [PATCH 22/25] Update ExtendedHubbardMom1D.jl --- src/Hamiltonians/ExtendedHubbardMom1D.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Hamiltonians/ExtendedHubbardMom1D.jl b/src/Hamiltonians/ExtendedHubbardMom1D.jl index eab5421a9..7fe8f688c 100644 --- a/src/Hamiltonians/ExtendedHubbardMom1D.jl +++ b/src/Hamiltonians/ExtendedHubbardMom1D.jl @@ -63,7 +63,7 @@ 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) From 55850507bd397771bc2d56fb9f93f8521ae0d682 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 7 Oct 2025 16:57:47 +1300 Subject: [PATCH 23/25] Update HubbardMomSpace.jl --- src/Hamiltonians/HubbardMomSpace.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl index aa1b84c5e..5235890b7 100644 --- a/src/Hamiltonians/HubbardMomSpace.jl +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -195,10 +195,10 @@ end onproduct = 0 for i in 1:length(map) occ_i = map[i].occnum - onproduct += occ_i * (occ_i - 1) + onproduct += occ_i * (occ_i - 1) / 2 for j in 1:i-1 occ_j = map[j].occnum - onproduct += 4 * occ_i * occ_j + onproduct += 2 * occ_i * occ_j end end return onproduct * u From ab5db4161d55548d799439546eae36eccb0ea32d Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 7 Oct 2025 17:20:33 +1300 Subject: [PATCH 24/25] Update Hamiltonians.jl --- test/Hamiltonians.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl index bbfc289ba..82b074a43 100644 --- a/test/Hamiltonians.jl +++ b/test/Hamiltonians.jl @@ -28,7 +28,7 @@ end HubbardReal1DEP(BoseFS(1, 2, 3, 4); u=1.0im), HubbardMom1D(BoseFS((6, 0, 0, 4)); t=1.0, u=0.5), HubbardMom1D(OccupationNumberFS(6, 0, 0, 4); t=1.0, u=0.5), - HubbardMom1D(BoseFS((6, 0, 0, 4)); t=1.0, u=0.5 + im), + HubbardMom1D(BoseFS((6, 0, 0, 4)); t=1.0, u=0.5), 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), From 820798109a61c1aa44208a03b2bc45d70a4adf53 Mon Sep 17 00:00:00 2001 From: Satyanand Kuwar <168614969+Skuwar1@users.noreply.github.com> Date: Tue, 7 Oct 2025 17:35:01 +1300 Subject: [PATCH 25/25] Update HubbardMomSpace.jl --- src/Hamiltonians/HubbardMomSpace.jl | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/src/Hamiltonians/HubbardMomSpace.jl b/src/Hamiltonians/HubbardMomSpace.jl index 5235890b7..a1345e56b 100644 --- a/src/Hamiltonians/HubbardMomSpace.jl +++ b/src/Hamiltonians/HubbardMomSpace.jl @@ -154,9 +154,7 @@ between two different component occupied mode maps `map1` and `map2`. `g` is the """ @inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, u, w) where {D,S} - if iszero(u) && iszero(w) - return 0 - end + onproduct = 0 for i in 1:length(map) occ_i = map[i].occnum @@ -172,9 +170,6 @@ end @inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, g::CubicGrid{D,S}, ::Nothing, w) where {D,S} - if iszero(w) - return 0 - end onproduct = 0 for i in 1:length(map) occ_i = map[i].occnum @@ -189,9 +184,7 @@ end end @inline function extended_mom_transfer_diag(map::BoseOccupiedModeMap, ::CubicGrid, u, ::Nothing) - if iszero(u) - return 0 - end + onproduct = 0 for i in 1:length(map) occ_i = map[i].occnum @@ -205,9 +198,7 @@ end end @inline function extended_mom_transfer_diag(map::FermiOccupiedModeMap, g::CubicGrid{D,S}, _, w) where {D,S} - if iszero(w) - return 0 - end + onproduct = 0 for i in 1:length(map) occ_i = map[i].occnum