New Hamiltonian Interface

.gitignore

@@ -1,3 +1,4 @@
 /Manifest.toml
 .DS_Store
 /test/Manifest.toml
+.vscode

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]

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

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")

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