diff --git a/src/Hamiltonians/FroelichPolaron_nD.jl b/src/Hamiltonians/FroelichPolaron_nD.jl new file mode 100644 index 000000000..abbd0a7d0 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl index 4baf1c052..1fdea9f96 100644 --- a/src/Hamiltonians/Hamiltonians.jl +++ b/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") diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl index ebb1976dc..1e29099b5 100644 --- a/test/Hamiltonians.jl +++ b/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)