[WIP] Additional exx kernels: Voxel averaged and Wigner-Seitz truncation by toschaefer · Pull Request #1282 · JuliaMolSim/DFTK.jl

toschaefer · 2026-03-11T08:03:44Z

Two additional useful strategies for the Coulomb kernel for electron repulsion integrals (e.g. for exx).

Wigner-Seitz truncation
R. Sundararaman, T. A. Arias. Phys. Rev. B 87, 165122 (2013)

Voxel averaged (mean potential)
T. Schaefer et al, J. Chem. Phys. 160, 051101 (2024)
This requires Gauss-Legendre quadrature points, hence I added the package FastGaussQuadrature to the dependencies.

Sync with upstream

antoine-levitt · 2026-04-27T16:17:37Z

Nice! Particularly interested in the Wigner-Seitz truncation - is there any reason to do anything else? It seems to me this is the superior approach? Do we really want to support the voxel averaged? Also maybe some opportunities to share code with ewald.jl, but I haven't looked too closely.

It'd be nice to have a generalized API for all these functions that return centered functions (we have a couple of those: densities guesses, local pseudos, projectors, coulomb kernels, etc.). I think long-term we want the ability to get their real-space and partial Fourier transforms (wrt 1 or 2 dimensions), so we can do eg 2D materials. But maybe it's not worth unifying since they have quite different requirements (numerical psp vs analytic but singular coulomb really are quite different computationally)

CC @denizunel

mfherbst · 2026-04-28T06:45:01Z

Also maybe some opportunities to share code with ewald.jl,

I think so, but I'd keep that for later, same as the generalized API. I'm currently looking at this on the side to streamline the code a bit, will hopefully push a review in the next week.

toschaefer · 2026-05-16T15:47:05Z

Nice! Particularly interested in the Wigner-Seitz truncation - is there any reason to do anything else? It seems to me this is the superior approach?

For EXX of gapped systems Wigner-Seitz truncation is probably the best approach to reach the thermodynamic llimit.
But for beyond-EXX stuff like many-electron correlation methods (MP2, CC, ...) it may not be quite as superior. In particular when we consider relatively small cells where a truncation is simply too rough (no matter if spehrically or Wigner-Seitz).

Do we really want to support the voxel averaged?

We need it for surface science (strongly anisotropic cells) with coupled cluster theory. Since our plan is to make DFTK a standard workflow, voxel averaged would be essential.
It would also be nice to have it for tests with the new (planned) mixers for SCF calculations with EXX. In particular for metals, voxel averaged could be advantageous.

toschaefer · 2026-05-21T06:46:08Z

rebased to make merge easier

antoine-levitt · 2026-05-22T13:24:01Z

Is this still WIP or ready for review?

toschaefer · 2026-05-22T13:56:34Z

from my end, ready for review

antoine-levitt

Very nice!

antoine-levitt · 2026-05-27T09:31:50Z

 - [`SphericallyTruncatedCoulomb`](@ref): θ(R-r)/r
+- [`WignerSeitzTruncatedCoulomb`](@ref): χ(r)/r where χ(r)=1 inside Wigner-Seitz cell, otherwise 0.

 If an interaction model features a singularity, that requires some special treatment,


no comma after singularity (germanism)

features a singularity -> is long-range

antoine-levitt · 2026-05-27T09:35:48Z

-#      where Γ is BZ volume) is used.
+#      (where Γ is BZ volume) is used.
 #   _compute_kernel_fourier(::InteractionKernel, basis, qpt, q)
 #      The single q-point version of compute_kernel_fourier


there's compute_ and eval_, should all these be merged?

(we really need to have a unified API for localized functions that can be evaluated in real or reciprocal space)

also should we have a eval_kernel_real? This is often useful for debug purposes (can be left as TODO if you don't feel like doing it)

compute_ sets up and returns the full vector with all the complexity of how the sampling points (including the singularity) are defined for each G+q while eval_ is more the mathematical evaluation of the kernel at a given point G+q.

I added a TODO for eval_kernel_real.

I don't get this. Both compute the same thing and should be consistent. If the point is that _compute_kernel_fourier(::InteractionKernel, basis, qpt, q) takes multiple kpoints then that should be a separate method of the same function. To be clear, I'd suggest

# computes khat(p). Not guaranteed to be fast eval_kernel_fourier(::Coulomb, p) = 1/p^2 # eval_kernel_fourier computes khat(G+pshift) for all G in basis using a possibly faster algorithm than just looping over all G+pshift # default implementation function eval_kernel_fourier(k::InteractionKernel, basis, pshift) map(G_vectors(basis)) do G eval_kernel_fourier(G+pshift) end end function eval_kernel_fourier(k::WignerSeitz, p) r = eval_kernel_fourier(ShortRange(k.ω)) # add to r the long-range by "slow Fourier transform", just looping over all r_vectors end function eval_kernel_fourier(k::WignerSeitz, pshift) # FFT based algorithm end # add a test that this overload is correct: for all kernels, the result eval_kernel_fourier(k,basis,pshift)

I see. I think there was a misleading comment. I updated the doc string for InteractionKernel to explain why we have eval_kernel_fourier and compute_kernel_fourier for the splitting: (i) formula and (ii) discretization/singularity handling.

@doc raw""" Abstract type for different interaction models. ### Architecture Computing interaction kernels is split into two parts: the mathematical formula (e.g. 4\pi/G^2) and the grid discretization. This split is primarily driven by the need to handle singularities in long-range kernels. 1. **InteractionKernel:** Defines the pure mathematical formula (via `eval_kernel_fourier`). 2. **regularization:** Necessary for long-range kernels (like `Coulomb` and `LongRangeCoulomb`) diverge as ``G+q \to 0``. Evaluating them on a periodic grid requires a specific strategy to handle this divergence. Because of this divergence, long-range `InteractionKernel`s contain a `regularization` field to dictate how the ``G +q=0`` component is built via `_compute_kernel_fourier`. Short-range kernels have a finite limit at `G+q \to 0``` a nd don't need a regularizatin. ### Available models: - [`Coulomb`](@ref): 1/r - [`ShortRangeCoulomb`](@ref): erfc(μr)/r - [`LongRangeCoulomb`](@ref): erf(μr)/r - [`SphericallyTruncatedCoulomb`](@ref): θ(R-r)/r - [`WignerSeitzTruncatedCoulomb`](@ref): χ(r)/r (1 inside Wigner-Seitz cell, 0 otherwise) ### Available singularity corrections (regularizations): - [`ProbeCharge`](@ref): Gygi-Baldereschi probe charge method - [`ReplaceSingularity`](@ref): Set the G+q=0 component to a specific value - [`VoxelAveraged`](@ref): Average the continuous kernel over the Brillouin zone voxel See also: [`compute_kernel_fourier`](@ref) """ abstract type InteractionKernel end Base.Broadcast.broadcastable(k::InteractionKernel) = Ref(k) # TODO: should we have a eval_kernel_real? # TODO: rename "k" in _compute_kernel_fourier(k... # TODO: change notation: p instead of G, G+q, ...

antoine-levitt · 2026-05-27T09:38:33Z


 See also: [`compute_kernel_fourier`](@ref)
 """
 abstract type InteractionKernel end


AbstractCoulombKernel? Interaction seems a bit generic (these are all Coulomb-like)

After some back-and-forth we called it InteractionKernel in #1223 (comment)

Hm OK, fine.

antoine-levitt · 2026-05-27T09:39:09Z

I think we use relatively consistently (or at least we wanted to) p for the Fourier-space variable in \R^3 (so usually things like p=G+k, and q=k-k'). Why mention G at all in this file? Can't you just use p?

Principally, I am happy to follow your stylistic preferences. I don't think we've mentioned this yet during the EXX development, and if I shall change it, I would prefer to postpone it for now, since I'm already working on the next PR (which also includes some changes in src/coulomb.jl).

OK, maybe leave a TODO? I think generally this file should be about just providing interaction(p), regardless of where p is coming from (G+q, G...)

antoine-levitt · 2026-05-27T09:41:50Z

give the equation defining this. Gsq -> p? Or ps?)

Gsq refers to G squared so psq? But I am not sure if I understand correctly...

antoine-levitt · 2026-05-27T10:15:13Z

+            d_min = min(d_min, d)
+        end
+        if d_min > sqrt(eps(T))
+            V_lr_real[idx] = erf(ω * d_min) / d_min


as above, define a erf(x)/x function?

Hm I think we should have a general way of doing f(x)/x stably and AD-friendly, this is coming up pretty frequently (also with @Technici4n in the spherical harmonics stuff)

Ah but of course it's my old friend the divided difference (high order in the case of f(x)/x^l). We have a general solution in https://github.com/xuequan818/MatrixFuns.jl but it might be a bit overkill...?

(the reason I want to get this right is otherwise it will come up to bite us in the ass with AD)

I want to do things like

function expm1_over_x(x) if abs(x) < 1000floatmin(Float64) x = 1000floatmin(Float64) end expm1(x) / x end

but of course this screws up AD at 0. We can have a cutoff point and switch to a finite order taylor expansion. This is ugly and getting good precision of higher derivatives is annoying, but I'm not sure what else to do.

Ok, I will extend the TODO comments by a general statement, that we need a clever and AD-friendly solution for a possible f(x)/x machinery.

Has this been already done ?

yes, below abstract type InteractionKernel end.

antoine-levitt · 2026-05-27T10:24:23Z

+    V_lr_real = zeros(Complex{T}, basis.fft_size...)
+    for idx in CartesianIndices(V_lr_real)
+        r_frac = r_vectors[idx]
+        r_centered = r_frac .- round.(r_frac) # MIC


write it out explicitly (I had to think to get what it was referring to)

antoine-levitt · 2026-05-27T10:27:11Z

+        r_centered = r_frac .- round.(r_frac) # MIC
+        r_cart = model.lattice * r_centered
+        d_min = norm(r_cart)
+        for dx in -nx:nx, dy in -ny:ny, dz in -nz:nz # Check neighbors for non-orthorhombic cells


I have to say I don't understand what's going on here, why isn't this just erf(r)/r for r in r_vectors?, maybe factor out the geometry stuff to another function?

Took me a while to understand what I did here. Added better comments to make it understandable.

antoine-levitt · 2026-05-27T10:28:31Z

+        Gnorm2 = sum(abs2, G_cart)
+        found_singularity = (iG==1 && iszero(q))
+        Rcut = cbrt(basis.model.unit_cell_volume*3/4/π)
+        if !found_singularity


why don't you just write the condition here? Also why not iszero(G_cart)?

(also same comment as above, hopefully we can avoid this special casing)

antoine-levitt · 2026-05-27T10:30:18Z

    end
    kernel_fourier
 end
+


Running out of time today, note to self to resume the review from here.

toschaefer · 2026-06-01T13:05:40Z

Thanks for the review!
I have another branch for the EXX k-points feature waiting on this one, so I would love to get this merged soon to keep the rebase simple. Does everything look ready to go on your end?

antoine-levitt · 2026-06-01T14:58:11Z

I'll try to get to it soon but it's pretty hectic at the moment, sorry... Since you're very active on this I don't want to be a blocker, but at the same time I think one of the key differentiator of DFTK is that we have relatively clean code (not really in the sense that we follow best practices, but rather in the sense that the code is relatively minimal and understandable), and code review really helps with that I find (in the sense that the code we let in has been understood at least once by at least two people).

I think the best way to do that is that you fork DFTK and make your changes there. Then, when we merge a PR on the main DFTK repo (or even along the way when you make commits to the PR), you merge it on to your fork. Since the git history is preserved I think the whole process should be relatively painless (in the sense that you shouldn't have to deal with much conflicts) and you don't have to wait on us to merge PRs to progress. (also AI made everything git much smoother, I find) In exchange I promise I'll get to it eventually in a reasonable timeframe. @mfherbst might have a better idea?

toschaefer · 2026-06-01T18:06:20Z

I totally agree with all of that. That’s why DFTK is my choice ;)

Regarding the fork strategy, normally I do it exactly like that. In this case, I shot myself in the foot because I refactored the coulomb.jl files in both branches in a somewhat non-orthogonal way. This is why I put the k-points branch on hold until this one merged.

That’s the only reason I brought it up, but please don't stress. I appreciate you taking the time to review it properly when you find time.

antoine-levitt

OK, I'm fine with merging the current code with lots of todo, to be adressed (by you, me or Claude :-p) at some point in the future. To me the API is the most pressing concern - it'd be nice to merge eval_kernel_fourier and _compute_kernel_fourier and clarify the relationship between the two.

antoine-levitt · 2026-06-03T11:20:42Z


 See also: [`compute_kernel_fourier`](@ref)
 """
 abstract type InteractionKernel end


Hm OK, fine.

antoine-levitt · 2026-06-03T11:22:39Z

OK, maybe leave a TODO? I think generally this file should be about just providing interaction(p), regardless of where p is coming from (G+q, G...)

antoine-levitt · 2026-06-03T11:39:02Z

-#      where Γ is BZ volume) is used.
+#      (where Γ is BZ volume) is used.
 #   _compute_kernel_fourier(::InteractionKernel, basis, qpt, q)
 #      The single q-point version of compute_kernel_fourier


I don't get this. Both compute the same thing and should be consistent. If the point is that _compute_kernel_fourier(::InteractionKernel, basis, qpt, q) takes multiple kpoints then that should be a separate method of the same function. To be clear, I'd suggest

# computes khat(p). Not guaranteed to be fast eval_kernel_fourier(::Coulomb, p) = 1/p^2 # eval_kernel_fourier computes khat(G+pshift) for all G in basis using a possibly faster algorithm than just looping over all G+pshift # default implementation function eval_kernel_fourier(k::InteractionKernel, basis, pshift) map(G_vectors(basis)) do G eval_kernel_fourier(G+pshift) end end function eval_kernel_fourier(k::WignerSeitz, p) r = eval_kernel_fourier(ShortRange(k.ω)) # add to r the long-range by "slow Fourier transform", just looping over all r_vectors end function eval_kernel_fourier(k::WignerSeitz, pshift) # FFT based algorithm end # add a test that this overload is correct: for all kernels, the result eval_kernel_fourier(k,basis,pshift)

antoine-levitt · 2026-06-03T11:40:22Z

 end
 ShortRangeCoulomb(; μ=0.2/u"Å") = ShortRangeCoulomb(austrip(μ))
 ShortRangeCoulomb(μ::Quantity) = ShortRangeCoulomb(austrip(μ))
 function eval_kernel_fourier(k::ShortRangeCoulomb, Gsq::T) where {T}


sorry I really should have reviewed the previous PR...

antoine-levitt · 2026-06-03T11:40:59Z

 end
 function LongRangeCoulomb(; μ=0.2/u"Å", regularization=ProbeCharge())
    LongRangeCoulomb(austrip(μ), regularization)
 end


@mfherbst I think we should keep the unit stuff away from the code as much as possible

mfherbst

Some minor comments. Main point is still the missing test. I've added some ideas what one could test. Feel free to object and defer if you don't have the time now, but I think it will help in the long run to have such consistency tests.

mfherbst · 2026-06-15T07:56:26Z

+
+### Architecture
+
+Computing interaction kernels is split into two parts: the mathematical formula (e.g. 4\pi/G^2) and the grid discretization. This split is primarily driven by the need to handle singularities in long-range kernels.


Refer to the new issue here to make it clear this is not yet finalised.

mfherbst · 2026-06-15T07:56:38Z

-# Each InteractionKernel should support the following functions:
-#   eval_kernel_fourier(::InteractionKernel, Gsq)
-#   eval_probe_charge_integral(::InteractionKernel, α)
-#      Should return ∫_{BZ}  kernel(q) * e^(-α * q^2) dq
-#      This is needed for the ProbeCharge regularisation. Note, that no factor 1/Γ
-#      where Γ is BZ volume) is used.
-#   _compute_kernel_fourier(::InteractionKernel, basis, qpt, q)
-#      The single q-point version of compute_kernel_fourier


I would not remove this, it's a pretty concise and compact overview of what needs to be done.

mfherbst · 2026-06-15T07:57:10Z

+# TODO: introduce a clever and AD-friendly way to deal with f(x)/x for x->0. E.g. intoduce phi(x) = iszero(x) ? one(x) : expm1(x) / x 
+# Also very useful for other InteractionKernels.


This feels a bit out of place here.

moved it to the other TODOs below InteractionKernel

mfherbst · 2026-06-15T07:58:31Z

+    # TODO: This is a bit hackish as the parameter needs to be re-computed every kernel
+    #       evaluation. Cleaner would be to move this further up in the call hierarchy,
+    #       such that compute_kernel_fourier is never called without Rcut being set to
+    #       not nothing


I agree with the first sentence, but I don't agree with the second (in the sense that this would be the solution. Maybe just have the first sentence and let us consider how to solve it when it actually becomes an issue. It does not cost that much compute after all ...

mfherbst · 2026-06-15T07:59:44Z

+            d_min = min(d_min, d)
+        end
+        if d_min > sqrt(eps(T))
+            V_lr_real[idx] = erf(ω * d_min) / d_min


Has this been already done ?

mfherbst · 2026-06-15T07:59:55Z

+
+    # Use ReplaceSingularity regularisation to explicitly set as the G==0
+    # component the exact limit of the kernel for G->0
+    _compute_kernel_fourier(kRcut, ReplaceSingularity(2π*Rcut^2), basis, qpt, q)


Has this been already done ?

mfherbst · 2026-06-15T08:06:57Z

+        k_wstrunc = compute_kernel_fourier(WignerSeitzTruncatedCoulomb(), basis)
+        E_wstrunc = exx_energy_only(basis, kpt, k_wstrunc, ψk_real, occk)
+        E_ref = -2.3456813523805415
+        @test abs(E_ref - E_wstrunc) < 1e-6


Is it possible to add a test that this is somewhat related to SphericallyTruncatedCoulomb, e.g. test that we get similar entries in k_wstrunc versus k_strunc or that some behaviour between the two agrees ? One may need to set Rcut for the SphericallyTruncatedCoulomb to a special value for this (e.g. the actual size of the cell etc.), but that sounds fair to me.

A good way to test this could be to consider evaluating the kernel on cells of the form a * Diagonal([1, 1, 1]) with larger and larger values of a, I would expect the two implementations to agree more and more as we make a larger.

I'll also add a few tests on non-cubic cells, just that we have something where Wigner-Seitz should make a real difference.

Thanks for the non-cubic test.

I am not sure if the sampling points of the kernel of SphericallyTruncatedCoulomb and WignerSeitzTruncatedCoulomb will get more and more similar when you increase the volume of a cubic cell. When I look at the spherically truncated case, the cos function will oscillate wildly with increasing R: v(G) = 4π/G² * (1 - cos(GR)).
Wigner-Seitz will also oscillate like crazy but probably not in the same way.

In the limit, the physics will be the same, but the individual sampling points probably not.

In the limit, the physics will be the same, but the individual sampling points probably not. But does this not imply that at least the energy should get similar ?

Well if it does not really work, than let's not get too hung up about it.

Yes, the energy should get similar, but in the limit this should hold for all kernels. Of course, converging at different rates...
So what you mean is to define a large a and then calculate exx_energy_only. The energies will agree to some extend and we simply define a test like @test E_wstrunc ≈ E_strunc atol=1e-3 with some apropriate atol?

Yes, essentially this captures we get the right limit and that things converge at the right order (if we take two different as ... at least in theory

toschaefer added 3 commits March 6, 2026 08:25

Merge pull request #8 from JuliaMolSim/master

2db40bc

Sync with upstream

Merge branch 'JuliaMolSim:master' into master

2d8936e

Merge branch 'JuliaMolSim:master' into master

e2135a3

toschaefer added 2 commits May 19, 2026 17:10

Merge branch 'JuliaMolSim:master' into master

aba4d00

reintroduce WignerSeitzTruncatedCoulomb and VoxelAveraged

e92b7d8

toschaefer force-pushed the exx_kernels branch from ec14222 to e92b7d8 Compare May 21, 2026 06:44

antoine-levitt reviewed May 27, 2026

View reviewed changes

little changes based on PR review

ffc8984

antoine-levitt reviewed Jun 3, 2026

View reviewed changes

further changes based on PR review

74f5ab0

mfherbst reviewed Jun 15, 2026

View reviewed changes

mfherbst and others added 4 commits June 15, 2026 10:20

Add non-spherical tests

4fe0f07

Split into extmodule

1a59c59

Fix non-spherical cell coulomb test

fa42d03

changes based on PR review

c82a863


		### Architecture

		Computing interaction kernels is split into two parts: the mathematical formula (e.g. 4\pi/G^2) and the grid discretization. This split is primarily driven by the need to handle singularities in long-range kernels.

		# TODO: introduce a clever and AD-friendly way to deal with f(x)/x for x->0. E.g. intoduce phi(x) = iszero(x) ? one(x) : expm1(x) / x
		# Also very useful for other InteractionKernels.

Conversation

toschaefer commented Mar 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

antoine-levitt commented Apr 27, 2026

Uh oh!

mfherbst commented Apr 28, 2026

Uh oh!

toschaefer commented May 16, 2026

Uh oh!

toschaefer commented May 21, 2026

Uh oh!

antoine-levitt commented May 22, 2026

Uh oh!

toschaefer commented May 22, 2026

Uh oh!

antoine-levitt 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!

toschaefer May 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

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!

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

toschaefer commented Mar 11, 2026 •

edited

Loading

toschaefer May 29, 2026 •

edited

Loading

toschaefer commented Jun 1, 2026 •

edited

Loading