From 98d48934cd64f53d821ccaa97282478ef5f6e3a9 Mon Sep 17 00:00:00 2001 From: mtsch Date: Mon, 6 May 2024 12:15:15 +1200 Subject: [PATCH] cleanup --- src/Hamiltonians/Hamiltonians.jl | 1 - src/Hamiltonians/correlation_functions.jl | 259 ++++++--------- src/Hamiltonians/geometry.jl | 45 ++- test/Hamiltonians.jl | 373 ++++++++++++---------- 4 files changed, 334 insertions(+), 344 deletions(-) diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl index 5e705b253..b9fddc487 100644 --- a/src/Hamiltonians/Hamiltonians.jl +++ b/src/Hamiltonians/Hamiltonians.jl @@ -88,7 +88,6 @@ export G2MomCorrelator, G2RealCorrelator, G2RealSpace, SuperfluidCorrelator, Den export StringCorrelator export Geometry, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries -export num_neighbours, neighbour_site, num_dimensions export sparse # from SparseArrays diff --git a/src/Hamiltonians/correlation_functions.jl b/src/Hamiltonians/correlation_functions.jl index 4c9fbae50..cd22cf4f7 100644 --- a/src/Hamiltonians/correlation_functions.jl +++ b/src/Hamiltonians/correlation_functions.jl @@ -1,5 +1,36 @@ -using StaticArrays -# TO-DO: add geometry for higher dimensions. +""" + circshift_dot(arr1, arr2, inds) + +Fast, non-allocation version of + +```julia +dot(arr1, circshift(arr2, inds)) +``` +""" +function circshift_dot(arr1, arr2, inds) + _circshift_dot!(arr1, (), arr2, (), axes(arr2), inds) +end + +# The following is taken from Julia's implementation of circshift. +@inline function _circshift_dot!( + dst, rdest, src, rsrc, + inds::Tuple{AbstractUnitRange,Vararg{Any}}, + shiftamt::Tuple{Integer,Vararg{Any}} +) + ind1, d = inds[1], shiftamt[1] + s = mod(d, length(ind1)) + sf, sl = first(ind1)+s, last(ind1)-s + r1, r2 = first(ind1):sf-1, sf:last(ind1) + r3, r4 = first(ind1):sl, sl+1:last(ind1) + tinds, tshiftamt = Base.tail(inds), Base.tail(shiftamt) + + return _circshift_dot!(dst, (rdest..., r1), src, (rsrc..., r4), tinds, tshiftamt) + + _circshift_dot!(dst, (rdest..., r2), src, (rsrc..., r3), tinds, tshiftamt) +end +@inline function _circshift_dot!(dst, rdest, src, rsrc, inds, shiftamt) + return dot(view(dst, rdest...), view(src, rsrc...)) +end + """ G2RealCorrelator(d::Int) <: AbstractHamiltonian{Float64} @@ -53,11 +84,7 @@ function diagonal_element(::G2RealCorrelator{D}, add::SingleComponentFockAddress M = num_modes(add) d = mod(D, M) v = onr(add) - result = 0 - for i in eachindex(v) - result += v[i] * v[mod1(i + d, M)] - end - return result / M + return circshift_dot(v, v, (d,)) end function diagonal_element(::G2RealCorrelator{0}, add::CompositeFS) @@ -69,22 +96,79 @@ function diagonal_element(::G2RealCorrelator{D}, add::CompositeFS) where {D} M = num_modes(add) d = mod(D, M) v = sum(map(onr, add.components)) - result = 0 - for i in eachindex(v) - result += v[i] * v[mod1(i + d, M)] - end - return result / M + return circshift_dot(v, v, (d,)) end num_offdiagonals(::G2RealCorrelator, ::SingleComponentFockAddress) = 0 num_offdiagonals(::G2RealCorrelator, ::CompositeFS) = 0 -# not needed: -# get_offdiagonal(::G2RealCorrelator, add) -# starting_address(::G2RealCorrelator) +""" + G2RealSpace(g::Geometry) <: AbstractHamiltonian{SArray} + +Two-body operator for density-density correlation for all displacements. + +```math + \\hat{G}^{(2)}(d) = \\frac{1}{M} ∑_i^M∑_v \\hat{n}_i (\\hat{n}_{i+v} - \\delta_{0d}). +``` + +# See also + +* [`HubbardReal1D`](@ref) +* [`G2MomCorrelator`](@ref) +* [`G2RealCorrelator`](@ref) +* [`AbstractHamiltonian`](@ref) +* [`AllOverlaps`](@ref) +""" +struct G2RealSpace{A,B,G<:Geometry,S} <: AbstractHamiltonian{S} + geometry::G + init::S +end +function G2RealSpace(geometry::Geometry, source=1, target=source) + init = zeros(SArray{Tuple{size(geometry)...}}) + return G2RealSpace{source,target,typeof(geometry),typeof(init)}(geometry, init) +end + +function Base.show(io::IO, g2::G2RealSpace{A,B}) where {A,B} + print(io, "G2RealSpace($(g2.geometry), $A,$B)") +end + +LOStructure(::Type{<:G2RealSpace}) = IsDiagonal() + +num_offdiagonals(g2::G2RealSpace, _) = 0 + +function diagonal_element( + g2::G2RealSpace{A,B}, addr1::SingleComponentFockAddress, addr2::SingleComponentFockAddress +) where {A, B} + geo = g2.geometry + result = g2.init + if addr1 ≡ addr2 + v1 = v2 = onr(addr1, geo) + else + v1 = onr(addr1, geo) + v2 = onr(addr2, geo) + end + + @inbounds for i in eachindex(result) + res_i = 0.0 + δ_vec = Offsets(geo)[i] + + # Case of n_i(n_i - 1) on the same component + if A == B && all(==(0), δ_vec) + v2_offset = max.(v2 .- 1, 0) + result = setindex(result, dot(v1, v2_offset), i) + else + result = setindex(result, circshift_dot(v1, v2, δ_vec), i) + end + end + return result ./ length(geo) +end +function diagonal_element(g2::G2RealSpace{A,A}, addr::SingleComponentFockAddress) where {A} + return diagonal_element(g2, addr, addr) +end +function diagonal_element(g2::G2RealSpace{A,B}, addr::CompositeFS) where {A,B} + return diagonal_element(g2, addr.components[A], addr.components[B]) +end -# Methods that need to be implemented: -# • starting_address(::AbstractHamiltonian) - not needed """ G2MomCorrelator(d::Int,c=:cross) <: AbstractHamiltonian{ComplexF64} @@ -166,7 +250,6 @@ function num_offdiagonals(g::G2MomCorrelator{3}, add::BoseFS2C) sa = num_occupied_modes(add.bsa) sb = num_occupied_modes(add.bsb) return sa*(m-1)*sb - # number of excitations that can be made end function num_offdiagonals(g::G2MomCorrelator, add::SingleComponentFockAddress) @@ -377,143 +460,3 @@ function _string_diagonal_real(d, add) return result / M end - -""" - G2RealSpace(g::Geometry) <: AbstractHamiltonian{SArray} - -Two-body operator for density-density correlation for all displacements. - -```math - \\hat{G}^{(2)}(d) = \\frac{1}{M} ∑_i^M∑_v \\hat{n}_i (\\hat{n}_{i+v} - \\delta_{0d}). -``` - -# See also - -* [`HubbardReal1D`](@ref) -* [`G2MomCorrelator`](@ref) -* [`G2RealCorrelator`](@ref) -* [`AbstractHamiltonian`](@ref) -* [`AllOverlaps`](@ref) -""" -struct G2RealSpace{A,B,G<:Geometry,S} <: AbstractHamiltonian{S} - geometry::G - init::S -end -function G2RealSpace(geometry::Geometry, source=1, target=source) - init = zeros(SArray{Tuple{size(geometry)...}}) - return G2RealSpace{source,target,typeof(geometry),typeof(init)}(geometry, init) -end - -LOStructure(::Type{<:G2RealSpace}) = IsDiagonal() - -num_offdiagonals(g2::G2RealSpace, _) = 0 - -#= -function diagonal_element(g2::G2RealSpace{1,1,D}, addr::SingleComponentFockAddress) where {D} - geo = g2.geometry - result = g2.init - v = onr(addr) - - for i in eachindex(result) - res_i = 0.0 - δ_vec = Offsets(geo)[i] - for p in eachindex(v) - p_vec = geo[p] - q = geo[add(p_vec, δ_vec)] - if q ≠ 0 - res_i += v[p] * v[q] - (p == q) - end - end - result = setindex(result, res_i, i) - end - return result ./ length(geo) -end -=# - -function diagonal_element2( - g2::G2RealSpace{A,B}, addr1::SingleComponentFockAddress, addr2::SingleComponentFockAddress -) where {A,B} - geo = g2.geometry - result = g2.init - if addr1 ≡ addr2 - v1 = v2 = onr(addr1) - else - v1 = onr(addr1) - v2 = onr(addr2) - end - - for i in eachindex(result) - res_i = 0.0 - δ_vec = Offsets(geo)[i] - for p in eachindex(v1) - p_vec = geo[p] - q = geo[add(p_vec, δ_vec)] - if q ≠ 0 - # A = B implies addr1 and addr2 are the same address, in which case we need - # to subtract δ_{p,q} - res_i += v1[p] * v2[q] - (A == B) * (p == q) - end - end - result = setindex(result, res_i, i) - end - return result ./ length(geo) -end - -function diagonal_element( - g2::G2RealSpace{A,B}, addr1::SingleComponentFockAddress, addr2::SingleComponentFockAddress -) where {A, B} - geo = g2.geometry - result = g2.init - if addr1 ≡ addr2 - v1 = v2 = onr(addr1, geo) - else - v1 = onr(geo, addr1) - v2 = onr(geo, addr2) - end - - @inbounds for i in eachindex(result) - res_i = 0.0 - δ_vec = Offsets(geo)[i] - - # Case of n_i(n_i - 1) - if A == B && all(==(0), δ_vec) - v2_offset = max.(v2 .- 1, 0) - result = setindex(result, dot(v1, v2_offset), i) - else - #circshift!(v2_offset, v2, δ_vec) - result = setindex(result, csh(v2, δ_vec), i) - end - end - return result ./ length(geo) -end - -function diagonal_element(g2::G2RealSpace{A,A}, addr::SingleComponentFockAddress) where {A} - return diagonal_element(g2, addr, addr) -end -function diagonal_element(g2::G2RealSpace{A,B}, addr::CompositeFS) where {A,B} - return diagonal_element(g2, addr.components[A], addr.components[B]) -end - - -# TODO: clean up! -function csh(arr, inds) - _circshift_dot!((), arr, (), axes(arr), inds) -end - - -@inline function _circshift_dot!(rdest, src, rsrc, - inds::Tuple{AbstractUnitRange,Vararg{Any}}, - shiftamt::Tuple{Integer,Vararg{Any}}) - ind1, d = inds[1], shiftamt[1] - s = mod(d, length(ind1)) - sf, sl = first(ind1)+s, last(ind1)-s - r1, r2 = first(ind1):sf-1, sf:last(ind1) - r3, r4 = first(ind1):sl, sl+1:last(ind1) - tinds, tshiftamt = Base.tail(inds), Base.tail(shiftamt) - _circshift_dot!((rdest..., r1), src, (rsrc..., r4), tinds, tshiftamt) + - _circshift_dot!((rdest..., r2), src, (rsrc..., r3), tinds, tshiftamt) -end -@inline function _circshift_dot!(rdest, src, rsrc, inds, shiftamt) - dot(view(src, rdest...), view(src, rsrc...)) - #copyto!(CartesianIndices(rdest), src, CartesianIndices(rsrc)) -end diff --git a/src/Hamiltonians/geometry.jl b/src/Hamiltonians/geometry.jl index b5f214431..3bcb3a992 100644 --- a/src/Hamiltonians/geometry.jl +++ b/src/Hamiltonians/geometry.jl @@ -44,15 +44,24 @@ struct Geometry{D,Dims,Fold} end end -function PeriodicBoundaries(dims::Vararg{Int,D}) where {D} +function PeriodicBoundaries(dims::NTuple{D,Int}) where {D} return Geometry(dims, ntuple(Returns(true), Val(D))) end -function HardwallBoundaries(dims::Vararg{Int,D}) where {D} +PeriodicBoundaries(dims::Vararg{Int}) = PeriodicBoundaries(dims) + +function HardwallBoundaries(dims::NTuple{D,Int}) where {D} return Geometry(dims, ntuple(Returns(false), Val(D))) end -function LadderBoundaries(dims::Vararg{Int,D}) where {D} +HardwallBoundaries(dims::Vararg{Int}) = HardwallBoundaries(dims) + +function LadderBoundaries(dims::NTuple{D,Int}) where {D} return Geometry(dims, ntuple(>(1), Val(D))) end +LadderBoundaries(dims::Vararg{Int}) = LadderBoundaries(dims) + +function Base.show(io::IO, g::Geometry{<:Any,Dims,Fold}) where {Dims,Fold} + print(io, "Geometry($Dims, $Fold)") +end Base.size(g::Geometry{<:Any,Dims}) where {Dims} = Dims Base.size(g::Geometry{<:Any,Dims}, i) where {Dims} = Dims[i] @@ -96,7 +105,8 @@ end Base.getindex(g::Geometry, i::Int) = Tuple(CartesianIndices(size(g))[i]) """ - UnitVectors{D} <: AbstractVector{NTuple{D,Int}} + UnitVectors(D) <: AbstractVector{NTuple{D,Int}} + UnitVectors(geometry::Geometry) <: AbstractVector{NTuple{D,Int}} Iterate over unit vectors in `D` dimensions. @@ -114,10 +124,12 @@ julia> UnitVectors(3) struct UnitVectors{D} <: AbstractVector{NTuple{D,Int}} end UnitVectors(D) = UnitVectors{D}() +UnitVectors(::Geometry{D}) where {D} = UnitVectors{D}() Base.size(::UnitVectors{D}) where {D} = (2D,) + function Base.getindex(uv::UnitVectors{D}, i) where {D} - @boundscheck 0 < i ≤ 2D || throw(BoundsError(uv, i)) + @boundscheck 0 < i ≤ length(uv) || throw(BoundsError(uv, i)) if i ≤ D return _unit_vec(Val(D), i, 1) else @@ -131,26 +143,43 @@ end return (_unit_vec(Val(I-1), i, x)..., val) end +""" + Offsets(geometry::Geometry) <: AbstractVector{NTuple{D,Int}} + +```jldoctest +julia> geometry = Geometry((3,4)); + +julia> reshape(Offsets(geometry), (3,4)) +3×4 reshape(::Offsets{2}, 3, 4) with eltype Tuple{Int64, Int64}: + (-1, -1) (-1, 0) (-1, 1) (-1, 2) + (0, -1) (0, 0) (0, 1) (0, 2) + (1, -1) (1, 0) (1, 1) (1, 2) + +``` +""" struct Offsets{D} <: AbstractVector{NTuple{D,Int}} geometry::Geometry{D} end + Base.size(off::Offsets) = (length(off.geometry),) + @inline function Base.getindex(off::Offsets{D}, i) where {D} + @boundscheck 0 < i ≤ length(off) || throw(BoundsError(off, i)) geo = off.geometry vec = geo[i] return add(vec, ntuple(i -> -cld(size(geo, i), 2), Val(D))) end """ - neighbour_site(geom::Geometry, site, i) + neighbor_site(geom::Geometry, site, i) -Find the `i`-th neighbour of `site` in the geometry. If the move is illegal, return 0. +Find the `i`-th neighbor of `site` in the geometry. If the move is illegal, return 0. """ function neighbor_site(g::Geometry{D}, mode, chosen) where {D} return g[add(g[mode], UnitVectors(D)[chosen])] end -function BitStringAddresses.onr(add, geom::Geometry{D,S}) where {D,S} +function BitStringAddresses.onr(add, geom::Geometry{<:Any,S}) where {S} return SArray{Tuple{S...}}(onr(add)) end function BitStringAddresses.onr(add::CompositeFS, geom::Geometry) diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl index f64e62280..fe1e33897 100644 --- a/test/Hamiltonians.jl +++ b/test/Hamiltonians.jl @@ -773,205 +773,224 @@ end @test f.shift ≈ a.shift end -@testset "G2MomCorrelator" begin - # v0 is the exact ground state from BoseHubbardMom1D2C(aIni;ua=0,ub=0,v=0.1) - bfs1=BoseFS([0,2,0]) - bfs2=BoseFS([0,1,0]) - aIni = BoseFS2C(bfs1,bfs2) - v0 = DVec( - BoseFS2C((0, 2, 0), (0, 1, 0)) => 0.9999389545691221, - BoseFS2C((1, 1, 0), (0, 0, 1)) => -0.007812695959057453, - BoseFS2C((0, 1, 1), (1, 0, 0)) => -0.007812695959057453, - BoseFS2C((2, 0, 0), (1, 0, 0)) => 4.046694762039993e-5, - BoseFS2C((0, 0, 2), (0, 0, 1)) => 4.046694762039993e-5, - BoseFS2C((1, 0, 1), (0, 1, 0)) => 8.616127793651117e-5, - ) - g0 = G2MomCorrelator(0) - g1 = G2MomCorrelator(1) - g2 = G2MomCorrelator(2) - g3 = G2MomCorrelator(3) - @test imag(dot(v0,g0,v0)) == 0 # should be strictly real - @test abs(imag(dot(v0,g3,v0))) < 1e-10 - @test dot(v0,g0,v0) ≈ 0.65 rtol=0.01 - @test dot(v0,g1,v0) ≈ 0.67 rtol=0.01 - @test dot(v0,g2,v0) ≈ 0.67 rtol=0.01 - @test dot(v0,g3,v0) ≈ 0.65 rtol=0.01 - @test num_offdiagonals(g0,aIni) == 2 - - # on first component - g0f = G2MomCorrelator(0,:first) - g1f = G2MomCorrelator(1,:first) - @test imag(dot(v0,g0f,v0)) == 0 # should be strictly real - @test dot(v0,g0f,v0) ≈ 1.33 rtol=0.01 - @test dot(v0,g1f,v0) ≈ 1.33 + 7.08e-5im rtol=0.01 - # on second component - g0s = G2MomCorrelator(0,:second) - g1s = G2MomCorrelator(1,:second) - #@test_throws ErrorException("invalid ONR") get_offdiagonal(g0s,aIni,1) # should fail due to invalid ONR - @test dot(v0,g0s,v0) ≈ 1/3 - @test dot(v0,g1s,v0) ≈ 1/3 - # test against BoseFS - ham1 = HubbardMom1D(bfs1) - ham2 = HubbardMom1D(bfs2) - @test num_offdiagonals(g0f,aIni) == num_offdiagonals(ham1,bfs1) - @test num_offdiagonals(g0s,aIni) == num_offdiagonals(ham2,bfs2) - aIni = BoseFS2C(bfs2,bfs1) # flip bfs1 and bfs2 - @test get_offdiagonal(g0s,aIni,1) == (BoseFS2C(BoseFS{1,3}((0, 1, 0)),BoseFS{2,3}((1, 0, 1))), 0.47140452079103173) - # test on BoseFS - @test diagonal_element(g0s,bfs1) == 4/3 - @test diagonal_element(g0s,bfs2) == 1/3 -end +using Rimu.Hamiltonians: circshift_dot + +@testset "Correlation functions" begin + @testset "circhshift_dot" begin + for i in 1:10 + A = rand(3, 4, 5) + B = rand(3, 4, 5) + inds = (rand(0:3), rand(0:4), rand(0:5)) + @test circshift_dot(A, B, inds) ≈ dot(A, circshift(B, inds)) + end + end -@testset "G2RealCorrelator" begin - m = 6 - n1 = 4 - n2 = m - add1 = BoseFS((n1,0,0,0,0,0)) - add2 = near_uniform(BoseFS{n2,m}) - - # localised state - @test diagonal_element(G2RealCorrelator(0), add1) == n1 * (n1 - 1) / m - @test diagonal_element(G2RealCorrelator(1), add1) == 0. - - # constant density state - @test diagonal_element(G2RealCorrelator(0), add2) == (n2/m) * ((n2/m) - 1) - @test diagonal_element(G2RealCorrelator(1), add2) == (n2/m)^2 - - # local-local - comp = CompositeFS(add1,add1) - @test diagonal_element(G2RealCorrelator(0), comp) == 2n1 * (2n1 - 1) / m - @test diagonal_element(G2RealCorrelator(1), comp) == 0. - - # local-uniform (assuming unit filling) - comp = CompositeFS(add1,add2) - @test diagonal_element(G2RealCorrelator(0), comp) == (n1 + 1) * n1 / m - @test diagonal_element(G2RealCorrelator(1), comp) == (2 * (n1 + 1) + (m - 2)) / m - - # uniform-uniform - comp = CompositeFS(add2,add2) - @test diagonal_element(G2RealCorrelator(0), comp) == (2n2 / m) * (2 * (n2 / m) - 1) - @test diagonal_element(G2RealCorrelator(1), comp) == (2n2 / m)^2 - - # offdiagonals - @test num_offdiagonals(G2RealCorrelator(0), add1) == 0 - @test num_offdiagonals(G2RealCorrelator(0), comp) == 0 - - # Test show method - d = 5 - output = @capture_out print(G2RealCorrelator(d)) - @test output == "G2RealCorrelator($d)" -end + @testset "G2RealCorrelator" begin + m = 6 + n1 = 4 + n2 = m + add1 = BoseFS((n1,0,0,0,0,0)) + add2 = near_uniform(BoseFS{n2,m}) + + # localised state + @test diagonal_element(G2RealCorrelator(0), add1) == n1 * (n1 - 1) / m + @test diagonal_element(G2RealCorrelator(1), add1) == 0.0 + + # constant density state + @test diagonal_element(G2RealCorrelator(0), add2) == (n2/m) * ((n2/m) - 1) + @test diagonal_element(G2RealCorrelator(1), add2) == (n2/m)^2 + + # local-local + comp = CompositeFS(add1,add1) + @test diagonal_element(G2RealCorrelator(0), comp) == 2n1 * (2n1 - 1) / m + @test diagonal_element(G2RealCorrelator(1), comp) == 0.0 + + # local-uniform (assuming unit filling) + comp = CompositeFS(add1,add2) + @test diagonal_element(G2RealCorrelator(0), comp) == (n1 + 1) * n1 / m + @test diagonal_element(G2RealCorrelator(1), comp) == (2 * (n1 + 1) + (m - 2)) / m + + # uniform-uniform + comp = CompositeFS(add2,add2) + @test diagonal_element(G2RealCorrelator(0), comp) == (2n2 / m) * (2 * (n2 / m) - 1) + @test diagonal_element(G2RealCorrelator(1), comp) == (2n2 / m)^2 + + # offdiagonals + @test num_offdiagonals(G2RealCorrelator(0), add1) == 0 + @test num_offdiagonals(G2RealCorrelator(0), comp) == 0 + + # Test show method + d = 5 + output = @capture_out print(G2RealCorrelator(d)) + @test output == "G2RealCorrelator($d)" + end -@testset "SuperfluidCorrelator" begin - m = 6 - n1 = 4 - n2 = m - add1 = BoseFS((n1,0,0,0,0,0)) - add2 = near_uniform(BoseFS{n2,m}) - - # localised state - @test @inferred diagonal_element(SuperfluidCorrelator(0), add1) == n1/m - @test @inferred diagonal_element(SuperfluidCorrelator(1), add1) == 0. - - # constant density state - @test diagonal_element(SuperfluidCorrelator(0), add2) == n2/m - @test diagonal_element(SuperfluidCorrelator(1), add2) == 0. - - # offdiagonals - @test num_offdiagonals(SuperfluidCorrelator(0), add1) == 1 - @test num_offdiagonals(SuperfluidCorrelator(0), add2) == 6 - - # get_offdiagonal - @test get_offdiagonal(SuperfluidCorrelator(0), add1, 1) == (add1, n1/m) - @test get_offdiagonal(SuperfluidCorrelator(1), add1, 1) == (BoseFS((3,1,0,0,0,0)), sqrt(n1)/m) - @test get_offdiagonal(SuperfluidCorrelator(0), add2, 1) == (add2, 1/m) - @test get_offdiagonal(SuperfluidCorrelator(1), add2, 1) == (BoseFS((0,2,1,1,1,1)), sqrt(2)/m) - - # Test show method - d = 5 - output = @capture_out print(SuperfluidCorrelator(d)) - @test output == "SuperfluidCorrelator($d)" -end + @testset "G2RealSpace" begin + @testset "1D" begin + + end + end -@testset "StringCorrelator" begin - m = 6 - n1 = 4 - n2 = m + @testset "G2MomCorrelator" begin + # v0 is the exact ground state from BoseHubbardMom1D2C(aIni;ua=0,ub=0,v=0.1) + bfs1=BoseFS([0,2,0]) + bfs2=BoseFS([0,1,0]) + aIni = BoseFS2C(bfs1,bfs2) + v0 = DVec( + BoseFS2C((0, 2, 0), (0, 1, 0)) => 0.9999389545691221, + BoseFS2C((1, 1, 0), (0, 0, 1)) => -0.007812695959057453, + BoseFS2C((0, 1, 1), (1, 0, 0)) => -0.007812695959057453, + BoseFS2C((2, 0, 0), (1, 0, 0)) => 4.046694762039993e-5, + BoseFS2C((0, 0, 2), (0, 0, 1)) => 4.046694762039993e-5, + BoseFS2C((1, 0, 1), (0, 1, 0)) => 8.616127793651117e-5, + ) + g0 = G2MomCorrelator(0) + g1 = G2MomCorrelator(1) + g2 = G2MomCorrelator(2) + g3 = G2MomCorrelator(3) + @test imag(dot(v0,g0,v0)) == 0 # should be strictly real + @test abs(imag(dot(v0,g3,v0))) < 1e-10 + @test dot(v0,g0,v0) ≈ 0.65 rtol=0.01 + @test dot(v0,g1,v0) ≈ 0.67 rtol=0.01 + @test dot(v0,g2,v0) ≈ 0.67 rtol=0.01 + @test dot(v0,g3,v0) ≈ 0.65 rtol=0.01 + @test num_offdiagonals(g0,aIni) == 2 + + # on first component + g0f = G2MomCorrelator(0,:first) + g1f = G2MomCorrelator(1,:first) + @test imag(dot(v0,g0f,v0)) == 0 # should be strictly real + @test dot(v0,g0f,v0) ≈ 1.33 rtol=0.01 + @test dot(v0,g1f,v0) ≈ 1.33 + 7.08e-5im rtol=0.01 + # on second component + g0s = G2MomCorrelator(0,:second) + g1s = G2MomCorrelator(1,:second) + #@test_throws ErrorException("invalid ONR") get_offdiagonal(g0s,aIni,1) # should fail due to invalid ONR + @test dot(v0,g0s,v0) ≈ 1/3 + @test dot(v0,g1s,v0) ≈ 1/3 + # test against BoseFS + ham1 = HubbardMom1D(bfs1) + ham2 = HubbardMom1D(bfs2) + @test num_offdiagonals(g0f,aIni) == num_offdiagonals(ham1,bfs1) + @test num_offdiagonals(g0s,aIni) == num_offdiagonals(ham2,bfs2) + aIni = BoseFS2C(bfs2,bfs1) # flip bfs1 and bfs2 + @test get_offdiagonal(g0s,aIni,1) == (BoseFS2C(BoseFS{1,3}((0, 1, 0)),BoseFS{2,3}((1, 0, 1))), 0.47140452079103173) + # test on BoseFS + @test diagonal_element(g0s,bfs1) == 4/3 + @test diagonal_element(g0s,bfs2) == 1/3 + end - # unital refers to n̄=1 - non_unital_localised_state = BoseFS((n1,0,0,0,0,0)) - non_unital_uniform_state = near_uniform(non_unital_localised_state) + @testset "SuperfluidCorrelator" begin + m = 6 + n1 = 4 + n2 = m + add1 = BoseFS((n1,0,0,0,0,0)) + add2 = near_uniform(BoseFS{n2,m}) + + # localised state + @test @inferred diagonal_element(SuperfluidCorrelator(0), add1) == n1/m + @test @inferred diagonal_element(SuperfluidCorrelator(1), add1) == 0. + + # constant density state + @test diagonal_element(SuperfluidCorrelator(0), add2) == n2/m + @test diagonal_element(SuperfluidCorrelator(1), add2) == 0. + + # offdiagonals + @test num_offdiagonals(SuperfluidCorrelator(0), add1) == 1 + @test num_offdiagonals(SuperfluidCorrelator(0), add2) == 6 + + # get_offdiagonal + @test get_offdiagonal(SuperfluidCorrelator(0), add1, 1) == (add1, n1/m) + @test get_offdiagonal(SuperfluidCorrelator(1), add1, 1) == (BoseFS((3,1,0,0,0,0)), sqrt(n1)/m) + @test get_offdiagonal(SuperfluidCorrelator(0), add2, 1) == (add2, 1/m) + @test get_offdiagonal(SuperfluidCorrelator(1), add2, 1) == (BoseFS((0,2,1,1,1,1)), sqrt(2)/m) + + # Test show method + d = 5 + output = @capture_out print(SuperfluidCorrelator(d)) + @test output == "SuperfluidCorrelator($d)" + end - localised_state = BoseFS((n2,0,0,0,0,0)) - uniform_state = near_uniform(BoseFS{n2,m}) + @testset "StringCorrelator" begin + m = 6 + n1 = 4 + n2 = m - S0 = StringCorrelator(0) - S1 = StringCorrelator(1) - S2 = StringCorrelator(2) + # unital refers to n̄=1 + non_unital_localised_state = BoseFS((n1,0,0,0,0,0)) + non_unital_uniform_state = near_uniform(non_unital_localised_state) - @test num_offdiagonals(S0, localised_state) == 0 + localised_state = BoseFS((n2,0,0,0,0,0)) + uniform_state = near_uniform(BoseFS{n2,m}) - # non unital localised state - @test @inferred diagonal_element(S0, non_unital_localised_state) ≈ 20/9 - @test @inferred diagonal_element(S1, non_unital_localised_state) ≈ (-4/9)*exp(im * -2pi/3) + S0 = StringCorrelator(0) + S1 = StringCorrelator(1) + S2 = StringCorrelator(2) - # non unital near uniform state - @test @inferred diagonal_element(S0, non_unital_uniform_state) ≈ 2/9 + @test num_offdiagonals(S0, localised_state) == 0 - # constant density localised state - @test @inferred diagonal_element(S0, localised_state) == 5. - @test @inferred diagonal_element(S1, localised_state) ≈ 1 - @test @inferred diagonal_element(S2, localised_state) ≈ -1 + # non unital localised state + @test @inferred diagonal_element(S0, non_unital_localised_state) ≈ 20/9 + @test @inferred diagonal_element(S1, non_unital_localised_state) ≈ (-4/9)*exp(im * -2pi/3) - # constant density uniform state - @test @inferred diagonal_element(S0, uniform_state) == 0 - @test @inferred diagonal_element(S2, uniform_state) == 0 + # non unital near uniform state + @test @inferred diagonal_element(S0, non_unital_uniform_state) ≈ 2/9 - # Test return type for integer, and non-integer filling - @test @inferred diagonal_element(S0, localised_state) isa Float64 - @test @inferred diagonal_element(S1, non_unital_localised_state) isa ComplexF64 + # constant density localised state + @test @inferred diagonal_element(S0, localised_state) == 5. + @test @inferred diagonal_element(S1, localised_state) ≈ 1 + @test @inferred diagonal_element(S2, localised_state) ≈ -1 - # Test show method - d = 5 - output = @capture_out print(StringCorrelator(d)) - @test output == "StringCorrelator($d)" + # constant density uniform state + @test @inferred diagonal_element(S0, uniform_state) == 0 + @test @inferred diagonal_element(S2, uniform_state) == 0 -end + # Test return type for integer, and non-integer filling + @test @inferred diagonal_element(S0, localised_state) isa Float64 + @test @inferred diagonal_element(S1, non_unital_localised_state) isa ComplexF64 -@testset "Momentum" begin - @test diagonal_element(Momentum(), BoseFS((0,0,2,1,3))) ≡ 2.0 - @test diagonal_element(Momentum(fold=false), BoseFS((0,0,2,1,3))) ≡ 7.0 - @test diagonal_element(Momentum(1), BoseFS((1,0,0,0))) ≡ -1.0 - @test_throws MethodError diagonal_element(Momentum(2), BoseFS((0,1,0))) + # Test show method + d = 5 + output = @capture_out print(StringCorrelator(d)) + @test output == "StringCorrelator($d)" - for add in (BoseFS2C((0,1,2,3,0), (1,2,3,4,5)), FermiFS2C((1,0,0,1), (0,0,1,0))) - @test diagonal_element(Momentum(1), add) + diagonal_element(Momentum(2), add) ≡ - diagonal_element(Momentum(0), add) end - @test num_offdiagonals(Momentum(), BoseFS((0,1,0))) == 0 - @test LOStructure(Momentum(2; fold=true)) == IsDiagonal() - @test Momentum(1)' === Momentum(1) -end - -@testset "DensityMatrixDiagonal" begin - @test diagonal_element(DensityMatrixDiagonal(5), FermiFS((0,1,0,1,0,1,0))) == 0 - @test diagonal_element(DensityMatrixDiagonal(2; component=1), BoseFS((1,5,1,0))) == 5 + @testset "Momentum" begin + @test diagonal_element(Momentum(), BoseFS((0,0,2,1,3))) ≡ 2.0 + @test diagonal_element(Momentum(fold=false), BoseFS((0,0,2,1,3))) ≡ 7.0 + @test diagonal_element(Momentum(1), BoseFS((1,0,0,0))) ≡ -1.0 + @test_throws MethodError diagonal_element(Momentum(2), BoseFS((0,1,0))) - for add in ( - CompositeFS(BoseFS((1,2,3,4,5)), BoseFS((5,4,3,2,1))), - BoseFS2C((1,2,3,4,5), (5,4,3,2,1)) - ) - for i in 1:5 - @test diagonal_element(DensityMatrixDiagonal(i, component=1), add) == i - @test diagonal_element(DensityMatrixDiagonal(i, component=2), add) == 6 - i - @test diagonal_element(DensityMatrixDiagonal(i), add) == 6 + for add in (BoseFS2C((0,1,2,3,0), (1,2,3,4,5)), FermiFS2C((1,0,0,1), (0,0,1,0))) + @test diagonal_element(Momentum(1), add) + diagonal_element(Momentum(2), add) ≡ + diagonal_element(Momentum(0), add) end + + @test num_offdiagonals(Momentum(), BoseFS((0,1,0))) == 0 + @test LOStructure(Momentum(2; fold=true)) == IsDiagonal() + @test Momentum(1)' === Momentum(1) end - @test num_offdiagonals(DensityMatrixDiagonal(1), BoseFS((0,1,0))) == 0 - @test LOStructure(DensityMatrixDiagonal(2)) == IsDiagonal() - @test DensityMatrixDiagonal(15)' === DensityMatrixDiagonal(15) + @testset "DensityMatrixDiagonal" begin + @test diagonal_element(DensityMatrixDiagonal(5), FermiFS((0,1,0,1,0,1,0))) == 0 + @test diagonal_element(DensityMatrixDiagonal(2; component=1), BoseFS((1,5,1,0))) == 5 + + for add in ( + CompositeFS(BoseFS((1,2,3,4,5)), BoseFS((5,4,3,2,1))), + BoseFS2C((1,2,3,4,5), (5,4,3,2,1)) + ) + for i in 1:5 + @test diagonal_element(DensityMatrixDiagonal(i, component=1), add) == i + @test diagonal_element(DensityMatrixDiagonal(i, component=2), add) == 6 - i + @test diagonal_element(DensityMatrixDiagonal(i), add) == 6 + end + end + + @test num_offdiagonals(DensityMatrixDiagonal(1), BoseFS((0,1,0))) == 0 + @test LOStructure(DensityMatrixDiagonal(2)) == IsDiagonal() + @test DensityMatrixDiagonal(15)' === DensityMatrixDiagonal(15) + end end @testset "HubbardReal1DEP" begin