Update version to 0.11.1 and add FroehlichPolaron Hamiltonian

Project.toml

@@ -1,7 +1,7 @@
 name = "Rimu"
 uuid = "c53c40cc-bd84-11e9-2cf4-a9fde2b9386e"
 authors = ["Joachim Brand <j.brand@massey.ac.nz>"]
-version = "0.11.0"
+version = "0.11.1"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

docs/src/hamiltonians.md

@@ -73,6 +73,7 @@ HOCartesianCentralImpurity
 ```@docs
 MatrixHamiltonian
 Transcorrelated1D
+FroehlichPolaron
 ```
 
 ### Convenience functions

src/Hamiltonians/FroehlichPolaron.jl

@@ -0,0 +1,107 @@
+"""
+    FroehlichPolaron(
+        address;
+        alpha = 1.0,
+        total_mom = 0.0,
+        geometry = nothing,
+        num_dimensions = 1
+    )
+
+The Froehlich polaron Hamiltonian is given by
+
+```math
+H = (P̂_f - P)^2 + N̂ + Σₖ νₖ (âₖ^† + â_{-k})
+```
+
+where ``P`` is the total momentum `total_mom`, ``P̂_f = Σ_k k âₖ^† âₖ`` is the momentum
+operator for the bosons, ``N̂ = Σ_k âₖ^† âₖ`` is the number operator for the bosons, and
+``νₖ = \\sqrt{α}/|k|`` is the coupling strength determined by the coupling parameter
+`α == alpha`.
+
+The optional `geometry` argument specifies the geometry of the lattice for ``k``-space and
+should be of the type [`PeriodicBoundaries`](@ref). A simplified way of specifying the
+geometry is to provide the number of dimensions `num_dimensions`. In this case the
+[`num_modes(address)`](@ref) of `address` must be a square number for `num_dimensions = 2`,
+or a cube number for `num_dimensions = 3`.
+
+The `address` must be of type [`OccupationNumberFS`](@ref).
+"""
+struct FroehlichPolaron{
+    P, # total momentum
+    T, # eltype
+    A<:OccupationNumberFS, # address type
+    G # lattice type
+} <: AbstractHamiltonian{T}
+    addr::A
+    alpha::T
+    geometry::G
+    # nu_k::Vector{T}
+    # p_squared_k::Vector{T}
+end
+
+function FroehlichPolaron(
+    addr::OccupationNumberFS;
+    geometry=nothing,
+    alpha=1.0,
+    total_mom=0.0,
+    num_dimensions=1
+)
+    M = num_modes(addr) # this is compile-time information
+    if isnothing(geometry)
+        if num_dimensions == 1
+            geometry = PeriodicBoundaries(M)
+        elseif num_dimensions == 2
+            sm = round(Int, sqrt(M))
+            if sm^2 == M
+                geometry = PeriodicBoundaries(sm, sm)
+            else
+                throw(ArgumentError("Number of modes $M is not a square. Specify geometry."))
+            end
+        elseif num_dimensions == 3
+            cm = round(Int, cbrt(M))
+            if cm^3 == M
+                geometry = PeriodicBoundaries(cm, cm, cm)
+            else
+                throw(ArgumentError("Number of modes $M is not a cube. Specify geometry."))
+            end
+        else
+            throw(ArgumentError("Invalid number of dimensions: $num_dimensions"))
+        end
+    end
+    geometry isa PeriodicBoundaries || throw(ArgumentError("Invalid geometry: $geometry"))
+    alpha, P = promote(float(alpha), float(total_mom))
+    return FroehlichPolaron{P,typeof(P),typeof(addr),typeof(geometry)}(addr, alpha, geometry)
+end
+
+function Base.show(io::IO, h::FroehlichPolaron)
+    print(io, "FroehlichPolaron($(h.addr); alpha=$(h.alpha), total_mom=$(h.total_mom), ")
+    print(io, "geometry=$(h.geometry))")
+end
+
+function starting_address(h::FroehlichPolaron)
+    return getfield(h, :addr)
+end
+
+LOStructure(::Type{<:FroehlichPolaron{<:Real}}) = IsHermitian()
+
+Base.getproperty(h::FroehlichPolaron, s::Symbol) = getproperty(h, Val(s))
+Base.getproperty(h::FroehlichPolaron, ::Val{:addr}) = getfield(h, :addr)
+Base.getproperty(h::FroehlichPolaron, ::Val{:alpha}) = getfield(h, :alpha)
+Base.getproperty(::FroehlichPolaron{P}, ::Val{:total_mom}) where P = P
+Base.getproperty(h::FroehlichPolaron, ::Val{:geometry}) = getfield(h, :geometry)
+
+num_dimensions(h::FroehlichPolaron) = num_dimensions(h.geometry)
+
+# function diagonal_element(h::FroehlichPolaron, addr::OccupationNumberFS)
+#     return (
+#         (-h.total_mom)^2 +
+#         num_particles(addr) +
+#         sum(νk(h, k) for k in momentum(h, addr))
+#     )
+# end
+
+# TODO: Sort out scales for momentum and energy; additional parameters?
+# TODO: pre-compute and store grid for p_squared_k (k.e.) and nu_k (coupling strength)
+# TODO: compute diagonal and off-diagonal elements
+# TODO: rest of AbstractHamiltonian interface
+# TODO: write unit tests for FroehlichPolaron

src/Hamiltonians/Hamiltonians.jl

@@ -24,6 +24,7 @@ Harmonic oscillator models
 - [`HOCartesianCentralImpurity`](@ref)
 
 Other
+- [`FroehlichPolaron`](@ref)
 - [`MatrixHamiltonian`](@ref)
 - [`Transcorrelated1D`](@ref)
 
@@ -73,6 +74,7 @@ export TimeReversalSymmetry
 export Stoquastic
 export Transcorrelated1D
 export hubbard_dispersion, continuum_dispersion
+export FroehlichPolaron
 
 export G2MomCorrelator, G2RealCorrelator, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum
 export StringCorrelator
@@ -102,6 +104,7 @@ include("ExtendedHubbardReal1D.jl")
 
 include("BoseHubbardReal1D2C.jl")
 include("BoseHubbardMom1D2C.jl")
+include("FroehlichPolaron.jl")
 
 include("GutzwillerSampling.jl")
 include("GuidingVectorSampling.jl")

test/Hamiltonians.jl

@@ -880,7 +880,7 @@ end
     m = 6
     n1 = 4
     n2 = m
-    
+
     # 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)
@@ -893,7 +893,7 @@ end
     S2 = StringCorrelator(2)
 
     @test num_offdiagonals(S0, localised_state) == 0
-    
+
     # 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)
@@ -918,7 +918,7 @@ end
     d = 5
     output = @capture_out print(StringCorrelator(d))
     @test output == "StringCorrelator($d)"
-    
+
 end
 
 @testset "Momentum" begin
@@ -1587,4 +1587,15 @@ end
         null_addr = BoseFS(prod(S), 1=>0)
         @test isempty(fock_to_cart(null_addr, S))
     end
+
+    @testset "FroehlichPolaron" begin
+        addr = OccupationNumberFS(0, 0, 0, 0, 1, 0, 0, 0)
+        ham = FroehlichPolaron(addr; total_mom=3, alpha=6, num_dimensions=3)
+        @test num_dimensions(ham) == num_dimensions(ham.geometry) == 3
+        @test ham.geometry == PeriodicBoundaries(2, 2, 2)
+        @test eval(Meta.parse(repr(ham))) == ham
+        @test starting_address(ham) == ham.addr == addr
+
+        # TODO: test the rest of the interface (add `FroehlichPolaron` to interface tests)
+    end
 end