From e77e6a0221f42d8136042386a471990575451628 Mon Sep 17 00:00:00 2001 From: Antoine Levitt Date: Tue, 12 Apr 2022 17:56:45 +0200 Subject: [PATCH] Test consistency of local potential in both real and Fourier space (#647) --- Project.toml | 3 ++- src/pseudo/NormConservingPsp.jl | 2 +- src/pseudo/PspHgh.jl | 1 + src/terms/anyonic.jl | 4 ++-- test/PspHgh.jl | 23 +++++++++++++++++++---- test/anyons.jl | 2 +- 6 files changed, 26 insertions(+), 9 deletions(-) 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