Add FroehlichPolaron_nD structure and implementation

src/Hamiltonians/FroelichPolaron_nD.jl

@@ -0,0 +1,238 @@
+
+
+struct FroehlichPolaronND{
+    T,
+    M,
+    D,
+    A<:OccupationNumberFS{M}, 
+    MC,
+    G<:CubicGrid
+    } <: AbstractHamiltonian{T}
+    address::A
+    d:: Int64
+    geometry::G
+    alpha::T
+    mass::T
+    omega::T
+    l::SVector{D, Float64}
+    p::SVector{D, Float64}
+    ks::SVector{M, SVector{D, T}}
+    momentum_cutoff::Union{MC, Nothing}
+    mode_cutoff::Union{Int, Nothing}
+    vkc::T
+end
+
+
+function FroehlichPolaronND(
+    address::OccupationNumberFS{M,AT};
+    D::Int = 1,
+    geometry::CubicGrid = PeriodicBoundaries(ntuple(Returns(round(Int, M^(1/D))), D)),
+    alpha = 1,
+    mass = 1,
+    omega = 1,
+    l = ones(D),
+    p = zeros(D),
+    momentum_cutoff = nothing,
+    mode_cutoff = nothing,
+    vk_constant = (-1) * sqrt((gamma(((D-1)/2)) * alpha * 2^(D-(3/2)) * pi^((D-1)/2) * omega^2) / (sqrt(mass * omega))),
+) where {M, AT}
+
+
+    if length(l) != D || any(x -> x <= 0, l)
+        throw(ArgumentError("`l` must be a positive vector of length $D"))
+    end
+
+
+    if abs(M^(1/D) - round(M^(1/D))) > 0.01 
+        throw(ArgumentError("M = $M is not a perfect $D th power"))
+    end
+
+    alpha, mass, omega = promote(float(alpha), float(mass), float(omega))
+    l = SVector{D,Float64}(float.(l))
+    p = SVector{D,Float64}(float.(p))
+
+    ks_tmp = Vector{SVector{D,Float64}}(undef, M)
+
+
+    for idx_mode in 1:M
+        idx = Tuple(geometry[idx_mode])
+
+        kv = ntuple(d -> begin
+            if isodd(round(M^(1/D)))
+
+                m = idx[d] - 1 - div(round(M^(1/D)), 2)
+            else
+
+                m = idx[d] - div(round(M^(1/D)) ,2)
+            end
+            (2π / l[d]) * m
+        end, D)
+        ks_tmp[idx_mode] = SVector{D,Float64}(kv)
+    end
+
+    ks = SVector{M,SVector{D,Float64}}(Tuple(ks_tmp))
+
+
+    if !isnothing(momentum_cutoff)
+        momentum_cutoff = typeof(alpha)(momentum_cutoff)
+        momentum = dot(ks,onr(address))
+        if abs(momentum) > momentum_cutoff
+            throw(ArgumentError("Starting address has momentum $momentum which cannot exceed momentum_cutoff $momentum_cutoff"))
+        end
+    end
+
+
+    if isnothing(mode_cutoff)
+        mode_cutoff = Int(typemax(AT))
+    end
+    mode_cutoff = floor(Int, mode_cutoff)::Int
+    if _exceed_mode_cutoff(mode_cutoff, address)
+        throw(ArgumentError("Starting address cannot have occupations that exceed mode_cutoff"))
+    end
+
+
+
+
+
+    return FroehlichPolaronND{typeof(alpha), M, D, typeof(address), typeof(momentum_cutoff),typeof(geometry)}(
+        address,D, geometry, alpha, mass, omega, l, p, ks, momentum_cutoff, mode_cutoff,vk_constant)
+end
+
+
+function Base.show(io::IO, h::FroehlichPolaronND)
+    io = IOContext(io, :compact => true)
+    println(io, "FroehlichPolaronND(")
+    println(io, "  ", starting_address(h), ",")
+    println(io, "  geometry = ", h.geometry, ",")
+    println(io, "  alpha = ", h.alpha, ", mass = ", h.mass, ", omega = ", h.omega, ", D = ", h.d, ",vk_constant = ",h.vkc,",")
+    println(io, "  l = ", Float64.(h.l), ", p = ", Float64.(h.p),  ",")
+    !isnothing(h.momentum_cutoff) && println(io, "  momentum_cutoff = ", h.momentum_cutoff, ",")
+    println(io, "  mode_cutoff = ", isnothing(h.mode_cutoff) ? "nothing" : string(h.mode_cutoff))
+    print(io, ")")
+end
+
+LOStructure(::Type{<:FroehlichPolaronND}) = IsHermitian()
+
+starting_address(h::FroehlichPolaronND) = h.address
+
+function dimension(h::FroehlichPolaronND, address)
+    # takes into account `mode_cutoff` but not `momentum_cutoff`
+    M = num_modes(address)
+    n = h.mode_cutoff
+    return BigInt(n + 1)^BigInt(M)
+end
+
+struct FroehlichPolaronNDColumn{H,A} <: AbstractOperatorColumn{A,Float64,H}
+    hamiltonian::H
+    address::A
+    num_offdiagonals::Int
+end
+
+function operator_column(h::FroehlichPolaronND, address)
+    M = num_modes(address)
+    return FroehlichPolaronNDColumn(h, address, 2M)
+end
+
+
+function diagonal_element(col::FroehlichPolaronNDColumn)
+    h = col.hamiltonian
+    occ = onr(col.address)
+
+    Pphonon = zero(h.ks[1])
+    Nphonon = 0
+    for m in 1:num_modes(col.address)
+        nm = occ[m]
+        Nphonon += nm
+        Pphonon += h.ks[m] * nm
+    end
+    ek = dot(h.p - Pphonon, h.p - Pphonon) / (h.mass)
+    return h.omega * Nphonon + ek
+end
+
+
+struct FroehlichPolaronNDOffdiagonals{A,H} <: AbstractVector{Pair{A,Float64}}
+    address::A
+    h::H
+    num_offdiagonals::Int
+end
+
+offdiagonals(col::FroehlichPolaronNDColumn) = FroehlichPolaronNDOffdiagonals(col.address, col.hamiltonian, col.num_offdiagonals)
+
+
+Base.size(ods::FroehlichPolaronNDOffdiagonals) = (ods.num_offdiagonals,)
+
+
+@inline function calc_vk(h::FroehlichPolaronND, kidx::Int)
+    knorm = sqrt(dot(h.ks[kidx], h.ks[kidx]))
+    vol = prod(h.l)
+    if h.d == 1
+        return (2 * h.alpha /(h.l[1]))^0.5
+    else
+        if knorm == 0.0
+            return 0.0
+        else
+            return h.vkc * sqrt(1/ (knorm^(h.d-1) * vol))
+        end
+    end
+end
+
+@inline function phonon_op(h::FroehlichPolaronND, addr, chosen)
+    M = num_modes(addr)
+    if chosen ≤ M
+        if !isnothing(h.mode_cutoff) && onr(addr)[chosen] ≥ h.mode_cutoff
+            return addr => 0.0
+        end
+        new_addr, val = excitation(addr, (chosen,), ())
+        vk = calc_vk(h, chosen)
+        amp = -vk * val
+
+    else
+        k = chosen - M
+        new_addr, val = excitation(addr, (), (k,))
+        vk = calc_vk(h, k)
+        amp = -vk * val
+
+    end
+
+
+    if !isnothing(h.momentum_cutoff)
+        occ = onr(new_addr)
+        phononmom = zeros(h.d)
+        for m in 1:M
+
+            phononmom += h.ks[m] * occ[m]
+
+        end
+        if norm(phononmom) ≤ h.momentum_cutoff
+            return new_addr => amp
+        else
+            return addr => 0.0
+        end
+    else
+        return new_addr => amp
+    end
+
+end
+
+
+function Base.getindex(ods::FroehlichPolaronNDOffdiagonals, i::Int)
+    return phonon_op(ods.h, ods.address, i)
+end
+
+
+function Base.iterate(ods::FroehlichPolaronNDOffdiagonals, state::Int=1)
+    state > ods.num_offdiagonals && return nothing
+    return ods[state], state + 1
+end
+
+function random_offdiagonal(col::FroehlichPolaronNDColumn)
+    M2 = col.num_offdiagonals
+    i = rand(1:M2)
+    addr_val = phonon_op(col.hamiltonian, col.address, i)
+    return first(addr_val), 1/M2, last(addr_val)
+end
+
+
+parent_operator(col::FroehlichPolaronNDColumn) = col.hamiltonian
+starting_address(col::FroehlichPolaronNDColumn) = col.address
+num_offdiagonals(col::FroehlichPolaronNDColumn) = col.num_offdiagonals

src/Hamiltonians/Hamiltonians.jl

@@ -23,6 +23,9 @@ Other
 - [`FroehlichPolaron`](@ref)
 - [`MatrixHamiltonian`](@ref)
 - [`Transcorrelated1D`](@ref)
+- [`HamiltonianProduct`](@ref)
+- ['FroelichPolaronND'](@ref)
+
 - [`MolecularHamiltonian`](@ref)
 
 ## [Wrappers](#Hamiltonian-wrappers)
@@ -88,6 +91,7 @@ export Stoquastic
 export Transcorrelated1D
 export hubbard_dispersion, continuum_dispersion
 export FroehlichPolaron
+export FroehlichPolaronND
 export ParticleNumberOperator
 
 export MolecularHamiltonian
@@ -129,6 +133,8 @@ include("HubbardRealSpace.jl")
 include("ExtendedHubbardReal1D.jl")
 
 include("FroehlichPolaron.jl")
+include("FroelichPolaron_nD.jl")
+
 
 include("Transcorrelated1D.jl")
 

test/Hamiltonians.jl

@@ -74,8 +74,10 @@ end
         momentum(HubbardMom1D(BoseFS(0, 1, 5, 1, 0))),
         Rimu.FirstOrderTransitionOperator(HubbardRealSpace(BoseFS(1,1,1,1)), -5.0, 0.01),
         HubbardReal1D(BoseFS(2,0,0); u=1.0im) * ExtendedHubbardReal1D(BoseFS(2,0,0)),
+        FroehlichPolaronND(OccupationNumberFS(0,0,0,0)),        
         HubbardReal1D(BoseFS(2,0,0); u=1.0im) + ExtendedHubbardReal1D(BoseFS(2,0,0)),
         2*HubbardReal1D(BoseFS(2,0,0); u=1.0im)
+
     ]
         test_hamiltonian_interface(H)
         # Check that the result of show can be pasted into the REPL. Does not work with
@@ -1207,6 +1209,7 @@ end
     # test offdiagonal element
     f2_offdiag = (OccupationNumberFS(1,3,3), -f2.v*sqrt(3))
     @test get_offdiagonal(f2, addr2, 2) == f2_offdiag
+
 
     f3_offdiag = (OccupationNumberFS(1,2,3,3), -f3.v*sqrt(4))
     @test get_offdiagonal(f3, addr3, 8) == f3_offdiag
@@ -1228,7 +1231,7 @@ end
     f5 = FroehlichPolaron(addr5; l, mode_cutoff=1, momentum_cutoff)
     basis5 = build_basis(f5)
     mom_vec = map(o -> dot(o, f5.ks), onr.(basis5))
-    @test all(abs.(mom_vec) .≤ momentum_cutoff) == true
+    @test all(abs.(mom_vec) .≤ momentum_cutoff) ==  true
     @test length(basis5) == 20
 
     # with and without momentum cutoff
@@ -1237,6 +1240,81 @@ end
     @test get_offdiagonal(f6, addr5, 1) == get_offdiagonal(f7, addr5, 1)
 end
 
+@testset "FroehlichPolaronND" begin
+    addr1 = OccupationNumberFS(1,1,1)
+
+    # test momentum_cutoff, mode_cutoff and dimention when initialising
+    addr2 = OccupationNumberFS(1,2,3)
+    @test_throws ArgumentError FroehlichPolaronND(addr2; mode_cutoff=1.0)
+    @test_throws ArgumentError FroehlichPolaronND(addr2; momentum_cutoff=10.0)
+    @test_throws ArgumentError FroehlichPolaronND(OccupationNumberFS(3,2,1); momentum_cutoff=10.0)
+    @test_throws ArgumentError FroehlichPolaronND(OccupationNumberFS(3,2,1); D=2)
+
+    addr3 = OccupationNumberFS(1,2,3,4)
+    f2 = FroehlichPolaronND(addr2)
+    f3 = FroehlichPolaronND(addr3; mode_cutoff=20.0)
+
+    @test starting_address(f2) == f2.address == addr2
+
+    # test ks vector
+
+    step = (2π/3)
+    ks2 = [(3/1) * [k] for k in range(-π*(1+1/3) + step; step=step, length=3)]
+    @test Vector(f2.ks) == ks2
+    step = (2π/4)
+    ks3 = [(4/1)* [k] for k in range(-π+step; step=step, length=4)]
+    @test Vector(f3.ks) == ks3
+
+    # test num_offdiagonals
+    @test num_offdiagonals(operator_column(f2, addr1)) == 2*3
+
+    # test diagonal_element
+    f2_diag = f2.omega*6 + (1/f2.mass) * (norm(f2.p) - dot(f2.ks, onr(addr2)))^2
+    @test diagonal_element(f2*addr2) == f2_diag
+
+    # test offdiagonal element
+    offd = offdiagonals(operator_column(f2,addr2))
+    @test (offd[2]) == (OccupationNumberFS(1,3,3) => -(2* f2.alpha /(f2.l[1]))^0.5 *sqrt(3))
+
+
+    f3_offdiag = (OccupationNumberFS(1,2,3,3) => -(2* f2.alpha /(f2.l[1]))^0.5*sqrt(4))
+    offdf3 = offdiagonals(operator_column(f3,addr3))
+    @test (offdf3[8]) == f3_offdiag
+
+    # test mode_cutoff
+    offd =offdiagonals(operator_column(f2,OccupationNumberFS(10,3,4)))
+    @test offd[1][2] ≠ 0.0
+    offdf3 =offdiagonals(operator_column(f3,OccupationNumberFS(1,3,20,10)))
+    @test offdf3[3][2] == 0.0
+
+    # test momentum_cutoff
+    # addr2 has momentum 12.56
+    addr4 = OccupationNumberFS(1,2,1)
+    f4 = FroehlichPolaronND(addr4; momentum_cutoff=10.0)
+    offdf4 = offdiagonals(operator_column(f4,addr2))
+    @test offdf4[3][2] == 0.0
+
+
+    m = 5; l = 6
+    addr5 = OccupationNumberFS{5}()
+    mom_unit = 2π/l
+    momentum_cutoff = 1.5 * mom_unit
+    f5 = FroehlichPolaronND(addr5; l, mode_cutoff=1, momentum_cutoff)
+    basis5 = build_basis(f5)
+    mom_vec = map(o -> dot(o, f5.ks), onr.(basis5))
+    @test all(abs.(mom_vec) .≤ momentum_cutoff) == true
+    @test length(basis5) == 20
+
+    # with and without momentum cutoff
+    f6 = FroehlichPolaronND(addr5;  mode_cutoff=1)
+    f7 = FroehlichPolaronND(addr5;  mode_cutoff=1, momentum_cutoff=100)
+    offdf6 = offdiagonals(operator_column(f6, addr5))
+    offdf7 = offdiagonals(operator_column(f7, addr5))
+    @test offdf6[1] == offdf7[1]
+
+
+end
+
 """
     compare_to_bethe(g, nf, m)