diff --git a/.gitignore b/.gitignore index 0dba073..8cccf88 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ /Manifest.toml .DS_Store /test/Manifest.toml +.vscode diff --git a/Project.toml b/Project.toml index 7c5d44f..99185e1 100644 --- a/Project.toml +++ b/Project.toml @@ -7,7 +7,7 @@ version = "0.1.0" Rimu = "c53c40cc-bd84-11e9-2cf4-a9fde2b9386e" [compat] -Rimu = "0.14" +Rimu = "0.14, 0.15" julia = "1.9" [extras] diff --git a/src/BoseHubbardReal1D2C.jl b/src/BoseHubbardReal1D2C.jl index 905c134..b9c03bc 100644 --- a/src/BoseHubbardReal1D2C.jl +++ b/src/BoseHubbardReal1D2C.jl @@ -48,6 +48,27 @@ function Rimu.starting_address(h::BoseHubbardReal1D2C) return BoseFS2C(starting_address(h.ha), starting_address(h.hb)) end +struct BoseHubbardReal1D2CColumn{ + A,T,O<:BoseHubbardReal1D2C{T} +} <: Rimu.AbstractOperatorColumn{A,T,O} + address::A + hamiltonian::O + diagonal_element::T +end + +function operator_column(h::BoseHubbardReal1D2C, address::BoseFS2C) + diagonal_element = Rimu.diagonal_element(h, address) + return BoseHubbardReal1D2CColumn(address, h, diagonal_element) +end + +function parent_operator(c::BoseHubbardReal1D2CColumn) + return c.hamiltonian +end +function starting_address(c::BoseHubbardReal1D2CColumn) + return c.address +end + + Rimu.LOStructure(::Type{<:BoseHubbardReal1D2C{<:Real}}) = IsHermitian() Base.getproperty(h::BoseHubbardReal1D2C, s::Symbol) = getproperty(h, Val(s)) @@ -56,9 +77,17 @@ Base.getproperty(h::BoseHubbardReal1D2C, ::Val{:hb}) = getfield(h, :hb) Base.getproperty(h::BoseHubbardReal1D2C{<:Any,<:Any,<:Any,V}, ::Val{:v}) where {V} = V # number of excitations that can be made -function Rimu.num_offdiagonals(::BoseHubbardReal1D2C, address) +function Rimu.num_offdiagonals(c::BoseHubbardReal1D2CColumn) + address = starting_address(c) return 2*(num_occupied_modes(address.bsa) + num_occupied_modes(address.bsb)) end +function Rimu.num_offdiagonals(h::BoseHubbardReal1D2C) + address = starting_address(h) + num_particles_a = Rimu.num_particles(address.bsa) + num_particles_b = Rimu.num_particles(address.bsb) + return 2 * (min(num_particles_a, num_modes(h)) + + min(num_particles_b, num_modes(h))) +end """ bose_hubbard_2c_interaction(::BoseFS2C) @@ -79,27 +108,29 @@ end function Rimu.diagonal_element(ham::BoseHubbardReal1D2C, address::BoseFS2C) return ( - diagonal_element(ham.ha, address.bsa) + - diagonal_element(ham.hb, address.bsb) + + diagonal_element(ham.ha, address.bsa) + # still makes use of the old interface + diagonal_element(ham.hb, address.bsb) + # old interface ham.v * bose_hubbard_2c_interaction(address) ) end - -function Rimu.get_offdiagonal(ham::BoseHubbardReal1D2C, address, chosen) - nhops = num_offdiagonals(ham,address) - nhops_a = 2 * num_occupied_modes(address.bsa) - if chosen ≤ nhops_a - naddress_from_bsa, onproduct = hopnextneighbour(address.bsa, chosen) - elem = - ham.ha.t * onproduct - return BoseFS2C(naddress_from_bsa,address.bsb), elem - else - chosen -= nhops_a - naddress_from_bsb, onproduct = hopnextneighbour(address.bsb, chosen) - elem = -ham.hb.t * onproduct - return BoseFS2C(address.bsa,naddress_from_bsb), elem - end - # return new address and matrix element -end +Rimu.diagonal_element(c::BoseHubbardReal1D2CColumn) = c.diagonal_element + + +# function Rimu.get_offdiagonal(ham::BoseHubbardReal1D2C, address, chosen) +# nhops = num_offdiagonals(ham,address) +# nhops_a = 2 * num_occupied_modes(address.bsa) +# if chosen ≤ nhops_a +# naddress_from_bsa, onproduct = hopnextneighbour(address.bsa, chosen) +# elem = - ham.ha.t * onproduct +# return BoseFS2C(naddress_from_bsa,address.bsb), elem +# else +# chosen -= nhops_a +# naddress_from_bsb, onproduct = hopnextneighbour(address.bsb, chosen) +# elem = -ham.hb.t * onproduct +# return BoseFS2C(address.bsa,naddress_from_bsb), elem +# end +# # return new address and matrix element +# end struct OffdiagonalsBoseReal1D2C{ A<:BoseFS2C,T,H<:TwoComponentHamiltonian{T} @@ -110,9 +141,11 @@ struct OffdiagonalsBoseReal1D2C{ num_hops_a::Int end -function Rimu.offdiagonals(h::BoseHubbardReal1D2C, a::BoseFS2C) - hops_a = num_offdiagonals(h.ha, a.bsa) - hops_b = num_offdiagonals(h.hb, a.bsb) +function Rimu.offdiagonals(c::BoseHubbardReal1D2CColumn) + h = c.hamiltonian + a = c.address + hops_a = num_offdiagonals(h.ha, a.bsa) # old interface + hops_b = num_offdiagonals(h.hb, a.bsb) # old interface length = hops_a + hops_b return OffdiagonalsBoseReal1D2C(h, a, length, hops_a) @@ -122,13 +155,26 @@ function Base.getindex(s::OffdiagonalsBoseReal1D2C{A,T}, i)::Tuple{A,T} where {A @boundscheck 1 ≤ i ≤ s.length || throw(BoundsError(s, i)) if i ≤ s.num_hops_a new_a, matrix_element = get_offdiagonal(s.hamiltonian.ha, s.address.bsa, i) + # old interface new_add = A(new_a, s.address.bsb) else i -= s.num_hops_a new_b, matrix_element = get_offdiagonal(s.hamiltonian.hb, s.address.bsb, i) + # old interface new_add = A(s.address.bsa, new_b) end return new_add, matrix_element end Base.size(s::OffdiagonalsBoseReal1D2C) = (s.length,) + +function random_offdiagonal(c::BoseHubbardReal1D2CColumn) + return random_from_offdiagonal(offdiagonals(c)) +end + +function random_from_offdiagonal(offdiags::AbstractOffdiagonals) + nl = length(offdiags) # check how many sites we could reach + chosen = rand(1:nl) # choose one of them + naddress, melem = offdiags[chosen] + return naddress, 1.0 / nl, melem +end diff --git a/src/RimuLegacyHamiltonians.jl b/src/RimuLegacyHamiltonians.jl index e8c462e..2e34169 100644 --- a/src/RimuLegacyHamiltonians.jl +++ b/src/RimuLegacyHamiltonians.jl @@ -14,6 +14,9 @@ using Rimu.Hamiltonians: Hamiltonians, dimension, HubbardMom1D, starting_address export BoseFS2C export BoseHubbardMom1D2C, BoseHubbardReal1D2C, G2MomCorrelator +import Rimu.Interfaces: operator_column, parent_operator, starting_address, + random_offdiagonal + include("BoseFS2C.jl") include("BoseHubbardMom1D2C.jl") include("BoseHubbardReal1D2C.jl") diff --git a/test/runtests.jl b/test/runtests.jl index ca832e3..4c3d99a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -90,12 +90,12 @@ using RimuLegacyHamiltonians: bose_hubbard_2c_interaction @test hamA == Ĥ2cReal.ha @test hamB == Ĥ2cReal.hb - @test num_offdiagonals(Ĥ2cReal, aIni2cReal) == 16 - @test num_offdiagonals(Ĥ2cReal, aIni2cReal) == num_offdiagonals(Ĥ2cReal.ha, aIni2cReal.bsa) + num_offdiagonals(Ĥ2cReal.hb, aIni2cReal.bsb) + @test num_offdiagonals(operator_column(Ĥ2cReal, aIni2cReal)) == 16 + @test 16 == num_offdiagonals(Ĥ2cReal.ha, aIni2cReal.bsa) + num_offdiagonals(Ĥ2cReal.hb, aIni2cReal.bsb) @test dimension(Ĥ2cReal) == 1225 @test dimension(Float64, Ĥ2cReal) == 1225.0 - hp2c = offdiagonals(Ĥ2cReal, aIni2cReal) + hp2c = offdiagonals(operator_column(Ĥ2cReal, aIni2cReal)) @test length(hp2c) == 16 @test hp2c[1][1] == BoseFS2C(BoseFS((0, 2, 1, 1)), BoseFS((1, 1, 1, 1))) @test hp2c[1][2] ≈ -1.4142135623730951 @@ -138,28 +138,26 @@ using RimuLegacyHamiltonians: bose_hubbard_2c_interaction addr1 = near_uniform(BoseFS2C{1,100,20}) addr2 = near_uniform(BoseFS2C{100,1,20}) - for Hamiltonian in (BoseHubbardReal1D2C, BoseHubbardMom1D2C) - @testset "$Hamiltonian" begin - H1 = BoseHubbardReal1D2C(addr1; ta=1.0, tb=2.0, ua=0.5, ub=0.7, v=0.2) - H2 = BoseHubbardReal1D2C(addr2; ta=2.0, tb=1.0, ua=0.7, ub=0.5, v=0.2) - @test starting_address(H1) == addr1 - @test LOStructure(H1) == IsHermitian() - - hops1 = collect(offdiagonals(H1, addr1)) - hops2 = collect(offdiagonals(H2, addr2)) - sort!(hops1, by=a -> first(a).bsa) - sort!(hops2, by=a -> first(a).bsb) - - addrs1 = first.(hops1) - addrs2 = flip.(first.(hops2)) - values1 = last.(hops1) - values2 = last.(hops1) - @test addrs1 == addrs2 - @test values1 == values2 - - @test eval(Meta.parse(repr(H1))) == H1 - @test eval(Meta.parse(repr(H2))) == H2 - end + @testset "Hamiltonian" begin + H1 = BoseHubbardReal1D2C(addr1; ta=1.0, tb=2.0, ua=0.5, ub=0.7, v=0.2) + H2 = BoseHubbardReal1D2C(addr2; ta=2.0, tb=1.0, ua=0.7, ub=0.5, v=0.2) + @test starting_address(H1) == addr1 + @test LOStructure(H1) == IsHermitian() + + hops1 = collect(offdiagonals(operator_column(H1, addr1))) + hops2 = collect(offdiagonals(operator_column(H2, addr2))) + sort!(hops1, by=a -> first(a).bsa) + sort!(hops2, by=a -> first(a).bsb) + + addrs1 = first.(hops1) + addrs2 = flip.(first.(hops2)) + values1 = last.(hops1) + values2 = last.(hops1) + @test addrs1 == addrs2 + @test values1 == values2 + + @test eval(Meta.parse(repr(H1))) == H1 + @test eval(Meta.parse(repr(H2))) == H2 end end @testset "1D Bosons (2-component)" begin