diff --git a/Project.toml b/Project.toml
index 5a6d583409..21307ebd57 100644
--- a/Project.toml
+++ b/Project.toml
@@ -88,10 +88,11 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 wannier90_jll = "c5400fa0-8d08-52c2-913f-1e3f656c1ce9"
 
 [targets]
-test = ["Test", "Aqua", "DoubleFloats", "FiniteDiff", "GenericLinearAlgebra", "IntervalArithmetic", "Plots", "Random", "KrylovKit", "Logging", "JLD2", "WriteVTK", "wannier90_jll"]
+test = ["Test", "Aqua", "DoubleFloats", "FiniteDiff", "GenericLinearAlgebra", "IntervalArithmetic", "Plots", "Random", "KrylovKit", "Logging", "JLD2", "WriteVTK", "wannier90_jll", "QuadGK"]
diff --git a/src/pseudo/NormConservingPsp.jl b/src/pseudo/NormConservingPsp.jl
index 5064ae59db..727f86cbe1 100644
--- a/src/pseudo/NormConservingPsp.jl
+++ b/src/pseudo/NormConservingPsp.jl
@@ -51,7 +51,7 @@ eval_psp_local_real(psp::NormConservingPsp, r::AbstractVector) =
 
 Evaluate the local part of the pseudopotential in reciprocal space:
 V(q) = ∫_R^3 Vloc(r) e^{-iqr} dr
-     = 4π ∫_{R+} sin(qr)/q r e^{-iqr} dr
+     = 4π ∫_{R+} Vloc(r) sin(qr)/q r dr
 """
 eval_psp_local_fourier(psp::NormConservingPsp, q::Real) =
     error("Not implemented")
diff --git a/src/pseudo/PspHgh.jl b/src/pseudo/PspHgh.jl
index 9e71394e41..b75cbf56f5 100644
--- a/src/pseudo/PspHgh.jl
+++ b/src/pseudo/PspHgh.jl
@@ -157,6 +157,7 @@ end
 
 # [GTH98] (1)
 function eval_psp_local_real(psp::PspHgh, r::T) where {T <: Real}
+    r == 0 && return eval_psp_local_real(psp, eps(T)) # quick hack for the division by zero below
     cloc = psp.cloc
     rr = r / psp.rloc
     convert(T,
diff --git a/src/terms/anyonic.jl b/src/terms/anyonic.jl
index 69152f87b5..693c793f62 100644
--- a/src/terms/anyonic.jl
+++ b/src/terms/anyonic.jl
@@ -16,7 +16,7 @@
 # This way, ρ-ρref has zero mass, and A_SR is shorter range.
 # The long-range part is computed explicitly
 
-function ρref_real(x::T, y::T, M, σ) where {T <: Real}
+function ρref_real_2D(x::T, y::T, M, σ) where {T <: Real}
     r = hypot(x, y)
     # gaussian normalized to have integral M
     M * exp(-T(1)/2 * (r/σ)^2) / (σ^2 * (2T(π))^(2/2))
@@ -89,7 +89,7 @@ function TermAnyonic(basis::PlaneWaveBasis{T}, hbar, β) where {T}
     σ = 2
     for (ir, r) in enumerate(r_vectors(basis))
         rcart = basis.model.lattice * (r - @SVector[.5, .5, .0])
-        ρref[ir] = ρref_real(rcart[1], rcart[2], M, σ)
+        ρref[ir] = ρref_real_2D(rcart[1], rcart[2], M, σ)
         Aref[1][ir], Aref[2][ir] = magnetic_potential_produced_by_ρref(rcart[1], rcart[2], M, σ)
     end
     # Aref is not divergence-free in the finite basis, so we explicitly project it
diff --git a/test/PspHgh.jl b/test/PspHgh.jl
index 780e13d963..70e61497cb 100644
--- a/test/PspHgh.jl
+++ b/test/PspHgh.jl
@@ -3,6 +3,7 @@ using DFTK: load_psp, eval_psp_projector_fourier, eval_psp_local_fourier
 using DFTK: eval_psp_projector_real, psp_local_polynomial
 using DFTK: psp_projector_polynomial, qcut_psp_projector, qcut_psp_local
 using SpecialFunctions: besselj
+using QuadGK
 
 
 @testset "Check reading 'C-lda-q4'" begin
@@ -117,7 +118,7 @@ end
     end
 end
 
-@testset "Numerical integration to obtain fourier-space projectors" begin
+@testset "Projectors are consistent in real and Fourier space" begin
     # The spherical bessel function of the first kind in terms of ordinary bessels:
     function j(n, x::T) where {T}
         x == 0 ? zero(T) : sqrt(π/2x) * besselj(n+1/2, x)
@@ -129,16 +130,30 @@ end
         4π * x^2 * eval_psp_projector_real(psp, i, l, x) * j(l, q*x)
     end
 
-    dx, xmax = 0.01, 10
     for pspfile in ["Au-q11", "Ba-q10"]
         psp = load_psp("hgh/lda/" * pspfile)
         for (l, i) in [(0, 1), (0, 2), (0, 3), (1, 1), (1, 2), (1, 3),
                        (2, 1), (2, 2), (3, 1)]
             l > length(psp.rp) - 1 && continue  # Overshooting available AM
-            for q in [0.01, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 100]
-                reference = sum(x -> integrand(psp, i, l, q, x) * dx, 0:dx:xmax)
+            for q in [0.01, 0.1, 0.2, 0.5, 1, 2, 5, 10]
+                reference = quadgk(r -> integrand(psp, i, l, q, r), 0, Inf)[1]
                 @test reference ≈ eval_psp_projector_fourier(psp, i, l, q) atol=5e-15 rtol=1e-8
             end
         end
     end
 end
+
+@testset "Potentials are consistent in real and Fourier space" begin
+    reg_param = 0.001  # divergent integral, needs regularization
+    function integrand(psp, q, r)
+        4π*DFTK.eval_psp_local_real(psp, r) * exp(-reg_param * r) * sin(q*r) / q * r
+    end
+
+    for pspfile in ["Au-q11", "Ba-q10"]
+        psp = load_psp("hgh/lda/" * pspfile)
+        for q in [0.01, 0.2, 1, 1.3]
+            reference = quadgk(r -> integrand(psp, q, r), 0, Inf)[1]
+            @test reference ≈ eval_psp_local_fourier(psp, q) rtol=.1 atol = .1
+        end
+    end
+end
diff --git a/test/anyons.jl b/test/anyons.jl
index 3dbe71a5d0..366d5dd712 100644
--- a/test/anyons.jl
+++ b/test/anyons.jl
@@ -15,6 +15,6 @@ using Test
             - DFTK.magnetic_potential_produced_by_ρref(x, y,     M, σ)) / ε
     curlA = dAdx[2] - dAdy[1]
     divA = dAdx[1] + dAdy[2]
-    @test norm(curlA - 2π*DFTK.ρref_real(x, y, M, σ)) < 1e-4
+    @test norm(curlA - 2π*DFTK.ρref_real_2D(x, y, M, σ)) < 1e-4
     @test abs(divA) < 1e-6
 end