Add FroehlichPolaron_nD structure and implementation #355

yosrabi · 2025-09-26T01:35:31Z

Hi, After some testing I think this is working as intended. All feedback and improvement ideas are welcome.

mtsch

Looks great! I haven't tested the code yet, but I didn't notice any potential issues, so I only have some style comments for now. Generally, we try to adhere to something close to this style guide https://github.com/JuliaDiff/BlueStyle

The PR also needs some tests. Add a few examples of interface tests (see first test set in test/Hamiltonians.jl). Other things you could test for is that this code is equivalent to FroelichPolaron for 1D and some general property tests (you can use the FroelichPolaron tests for inspiration there).

mtsch · 2025-09-26T10:02:36Z

src/Hamiltonians/FroelichPolaron_nD.jl

+function FroehlichPolaron_nD(
+    address::OccupationNumberFS{M,AT};
+    D::Int = 2,
+    geometry::CubicGrid = PeriodicBoundaries(ntuple(_ -> Int(round(num_modes(address)^(1/D))), D)),


Suggested change

geometry::CubicGrid = PeriodicBoundaries(ntuple(_ -> Int(round(num_modes(address)^(1/D))), D)),

geometry::CubicGrid = PeriodicBoundaries(ntuple(Returns(round(Int, M^(1/D))), D)),

A bit more compact, but equivalent

mtsch · 2025-09-26T10:03:15Z

src/Hamiltonians/FroelichPolaron_nD.jl

+) where {M, AT}
+
+    # Check positive l
+    if length(l) != D || any(x -> x <= 0, l)


Feel free to ignore this comment, but I would write it like this

Suggested change

if length(l) != D || any(x -> x <= 0, l)

if length(l) != D || any(<=(0), l)

mtsch · 2025-09-26T10:03:59Z

src/Hamiltonians/FroelichPolaron_nD.jl

+    end
+
+    #check m_lin is possible
+    if abs(M^(1/D) - round(M^(1/D))) >0.01  


Suggested change

if abs(M^(1/D) - round(M^(1/D))) >0.01

if abs(M^(1/D) - round(M^(1/D))) > 0.01

Try to put spaces on both sides of operators, except in very long expressions where having some terms close together makes it more readable.

mtsch · 2025-09-26T10:07:31Z

src/Hamiltonians/FroelichPolaron_nD.jl

+    for idx_mode in 1:M
+    idx = Tuple(geometry[idx_mode])
+    kv = ntuple(d -> begin
+        m = idx[d] - 1 - div(round(M^(1/D)),2)   
+        (2π / l[d]) * m
+    end, D)
+    ks_tmp[idx_mode] = SVector{D,Float64}(kv)
+    end


Suggested change

for idx_mode in 1:M

idx = Tuple(geometry[idx_mode])

kv = ntuple(d -> begin

m = idx[d] - 1 - div(round(M^(1/D)),2)

(2π / l[d]) * m

end, D)

ks_tmp[idx_mode] = SVector{D,Float64}(kv)

end

for idx_mode in 1:M

idx = Tuple(geometry[idx_mode])

ntuple(D) do d

m = idx[d] - 1 - div(round(M^(1/D)), 2)

(2π / l[d]) * m

end

ks_tmp[idx_mode] = SVector{D,Float64}(kv)

end

Please make sure the code is indented properly! The do notation is also helpful here.

mtsch · 2025-09-26T10:09:53Z

src/Hamiltonians/FroelichPolaron_nD.jl

+    Pph = zero(h.ks[1])
+    Nph = 0


It's not clear what these are. Try to use more descriptive variable names.

mtsch · 2025-09-26T10:12:21Z

src/Hamiltonians/FroelichPolaron_nD.jl

+using SpecialFunctions
+
+#structure 
+struct FroehlichPolaron_nD{


We usually avoid underscores in struct names. Maybe FroelichPolaronND?

Since this supports 1D as well, do we still need the regular FroelichPolaron? If this is just as efficient and can be used as a drop-in replacement, I think the old Hamiltonian should be removed.

Thanks, @yosrabi for the PR! The code looks nice and I'm looking forward to having the multi-dimensional Froehlich model represented in Rimu.

For now I have just a few comments supporting what @mtsch has already said:

It is common in Julia to use 'CamelCase' for structs and 'snake_case' for functions, and we follow this convention (you can google the terms).

I would also strongly prefer to not have two different Froehlich polaron Hamiltonians in Rimu, in particular if they have intersecting functionality. What I had in mind for the extension to 2 and 3 dimensions was an implementation that just adds new features to the existing FroehlichPolaron by adding new keyword arguments (e.g. for dimension with default set to 1) and reusing the old code, user interface, and tests where possible and extending or replacing them only where necessary. Please consider whether this is still feasible! Otherwise we could go with Matija's suggestion to let this new implementation supersede the old one. This would make sense if the new code was faster/more efficient/more elegant or in some other way superior to the old one. As a minimum it should not be any worse.

Ideally, the call signature should be compatible such that replacing the old with the new Hamiltonian would be a non-breaking change (i.e. all previous scripts making use of the documented features of the old FroehlichPolaron would still run without errors). If there is a good reason to change the behaviour or call signature such that a breaking change is required, we can discuss that as well, of course.

Adding tests for the new code and features that run automatically is important (and a hard requirement for any new code to be added to the codebase). These test are all in the test/ folder. The entry point is test/runtests.jl but your test should be added to test/Hamiltonians.jl. I suggest that you look into that file and search for FroelichPolaron to see what tests are already being done for the existing implementation. Ideally, the tests should run on the smallest units (functions) of the new code, check for the correctness of the results, and cover all (feasible) use cases. The minimum requirement is that all (meaningful) parts of the new code are run at least once during the tests.

Finally, I see that GitHub has detected conflicts of your branch with the main develop branch. The reason for this to happen is likely that we recently merged a few other PRs that may be conflicting with yours. The best way to deal with this is to (first fetch and then) merge develop into your branch locally on your computer, and resolve any merge conflicts with a code editor (e.g. vscode) before committing the changes and pushing the local branch back to the repo. There should not be any serious merge conficts as you've mostly edited a new file. You'll simplify your life a lot if you do this before acting on any other comments from this code review.

mtsch · 2025-09-26T10:15:47Z

src/Hamiltonians/FroelichPolaron_nD.jl

+function operator_column(h::FroehlichPolaron_nD, address)
+    M = num_modes(address)
+    return FroehlichPolaron_nDColumn(h, address, 2M)
+end


I would move this function next to the struct definition since it acts like a constructor.

mtsch · 2025-09-26T10:17:47Z

src/Hamiltonians/FroelichPolaron_nD.jl

+offdiagonals(col::FroehlichPolaron_nDColumn) = FroehlichPolaron_nDOffdiagonals(col.address, col.hamiltonian, col.num_offdiagonals)
+
+Base.size(ods::FroehlichPolaron_nDOffdiagonals) = (ods.num_offdiagonals,)
+Base.eltype(::FroehlichPolaron_nDOffdiagonals{A}) where {A} = Pair{A,Float64}


This is not necessary, a default definition is created when you make something a subtype of AbstractVector.

mtsch · 2025-09-26T10:19:43Z

src/Hamiltonians/FroelichPolaron_nD.jl

+end
+
+
+#diagonal elements calculator


Try avoiding comments that simply state what the code does. They should be used to explain how/why something is done in a particular way or to break up long chunks of code.

mtsch · 2025-09-26T10:20:15Z

src/Hamiltonians/FroelichPolaron_nD.jl

+#calculate V_k for a specific k. Edge case of D = 1.
+
+


Here, a docstring would be more useful than a comment.

…tions and tests

yosrabi · 2025-10-14T09:20:10Z

Sorry i think I accidentally tried to merge the branch. The requested changes have been implemented and the tests updated.

jamie-tay

This is looking great! Just a few small comments.

jamie-tay · 2025-10-15T23:05:08Z

src/Hamiltonians/FroelichPolaron_nD.jl

+function operator_column(h::FroehlichPolaronND, address)
+    M = num_modes(address)
+    return FroehlichPolaronNDColumn(h, address, 2M)
+end


This should be under the definition of the FroehlichPolaronNDColumn struct.

jamie-tay · 2025-10-15T23:37:28Z

src/Hamiltonians/FroelichPolaron_nD.jl

+function Base.iterate(ods::FroehlichPolaronNDOffdiagonals, state::Int=1)
+    state > ods.num_offdiagonals && return nothing
+    return ods[state], state + 1
+end


This is not needed, since the offdiagonals struct is an AbstractVector.

jamie-tay · 2025-10-15T23:45:26Z

src/Hamiltonians/FroelichPolaron_nD.jl

+    end
+
+
+    if h.momentum_cutoff != nothing


Suggested change

if h.momentum_cutoff != nothing

if !isnothing(h.momentum_cutoff)

This is slightly more efficient.

jamie-tay · 2025-10-16T00:14:36Z

src/Hamiltonians/FroelichPolaron_nD.jl

+@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 (-1) * sqrt((gamma(((h.d-1)/2)) * h.alpha * 2^(h.d-(3/2)) * pi^((h.d-1)/2) * h.omega^2) / (sqrt(h.mass * h.omega) * (knorm^(h.d-1) * vol)))
+        end
+    end
+end


Most of this doesn't need to be calculated over and over again; I'd suggest moving the constant part into the Hamiltonian as h.v, and then this function simplifies to something like

@inline function calc_vk(h::FroehlichPolaronND, kidx::Int) if h.d == 1 return h.v else knorm = sqrt(dot(h.ks[kidx], h.ks[kidx])) if knorm == 0.0 return 0.0 else return h.v / sqrt(knorm^(h.d - 1)) end end end

jamie-tay · 2025-10-16T00:26:19Z

src/Hamiltonians/FroelichPolaron_nD.jl

+    if abs(M^(1/D) - round(M^(1/D))) > 0.01 
+        throw(ArgumentError("M = $M is not a perfect $D th power"))
+    end


Suggested change

if abs(M^(1/D) - round(M^(1/D))) > 0.01

throw(ArgumentError("M = $M is not a perfect $D th power"))

end

if D > 1 && !isinteger(log(D, M))

throw(ArgumentError("M = $M is not a perfect power of $D"))

end

jamie-tay · 2025-10-16T02:13:30Z

test/Hamiltonians.jl

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


Suggested change

FroehlichPolaronND(OccupationNumberFS(0,0,0,0))

FroehlichPolaronND(OccupationNumberFS(0,0,0,0)),

jamie-tay · 2025-10-16T02:19:18Z

test/Hamiltonians.jl

    @test get_offdiagonal(f6, addr5, 1) == get_offdiagonal(f7, addr5, 1)
 end

+@testset "FroehlichPolaronND" begin


These tests look good, we should probably have some checking that it works as expected for higher dimensions as well.

jamie-tay · 2025-10-16T02:24:38Z

src/Hamiltonians/FroelichPolaron_nD.jl

+end
+
+
+function Base.show(io::IO, h::FroehlichPolaronND)


The result of this should be able to be pasted back into the REPL.

jamie-tay · 2025-10-16T02:30:59Z

src/Hamiltonians/FroelichPolaron_nD.jl

@@ -0,0 +1,235 @@
+using SpecialFunctions


This is already done in Hamiltonians.jl, so we don't need it here.

jamie-tay · 2025-10-16T02:34:58Z

src/Hamiltonians/FroelichPolaron_nD.jl

+using SpecialFunctions
+
+
+struct FroehlichPolaronND{


This should have a docstring; see the other Hamiltonians for examples.

yosrabi requested review from ElkePahl, jamie-tay, joachimbrand and mtsch September 26, 2025 01:35

Add FroehlichPolaron_nD structure and implementation

0eb6437

mtsch requested changes Sep 26, 2025

View reviewed changes

yosrabi and others added 2 commits October 14, 2025 12:48

Updates and corrections to FroelichPolaron_nD and update related func…

2ea9f9a

…tions and tests

Merge branch 'develop' into feature-Froelich-hamiltonian-nD

92acdb9

yosrabi requested a review from mtsch October 14, 2025 09:16

jamie-tay requested changes Oct 16, 2025

View reviewed changes

	geometry::CubicGrid = PeriodicBoundaries(ntuple(_ -> Int(round(num_modes(address)^(1/D))), D)),
	geometry::CubicGrid = PeriodicBoundaries(ntuple(Returns(round(Int, M^(1/D))), D)),

	if length(l) != D \|\| any(x -> x <= 0, l)
	if length(l) != D \|\| any(<=(0), l)

	if abs(M^(1/D) - round(M^(1/D))) >0.01
	if abs(M^(1/D) - round(M^(1/D))) > 0.01

	if h.momentum_cutoff != nothing
	if !isnothing(h.momentum_cutoff)

	FroehlichPolaronND(OccupationNumberFS(0,0,0,0))
	FroehlichPolaronND(OccupationNumberFS(0,0,0,0)),

Add FroehlichPolaron_nD structure and implementation #355

Are you sure you want to change the base?

Add FroehlichPolaron_nD structure and implementation #355

Uh oh!

Conversation

yosrabi commented Sep 26, 2025

Uh oh!

mtsch left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yosrabi commented Oct 14, 2025

Uh oh!

jamie-tay left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants