tests

Author: mfherbst

src/response/chi0.jl

@@ -344,11 +344,8 @@ function compute_δψ!(δψ, basis::PlaneWaveBasis{T}, H, ψ, εF, ε, δHψ, ε
         δψk  = δψ[ik]
         tolk = bandtol[ik]
         εk_minus_q = ε_minus_q[ik]
-
-        rnorm_k  = residual_norms[ik]
-        n_iter_k = n_iter[ik]
-        sizehint!(rnorm_k,  length(εk_minus_q))
-        sizehint!(n_iter_k, length(εk_minus_q))
+        sizehint!(residual_norms[ik], length(εk_minus_q))
+        sizehint!(n_iter[ik], length(εk_minus_q))
 
         ψk_extra = ψ_extra[ik]
         Hψk_extra = Hk * ψk_extra
@@ -377,8 +374,8 @@ function compute_δψ!(δψ, basis::PlaneWaveBasis{T}, H, ψ, εF, ε, δHψ, ε
                                     res.tol, res.residual_norm)
 
             δψk[:, n] .+= res.δψkn
-            push!(rnorm_k,  res.residual_norm)
-            push!(n_iter_k, res.n_iter)
+            push!(residual_norms[ik], res.residual_norm)
+            push!(n_iter[ik], res.n_iter)
             converged = converged && res.converged
         end
     end

src/response/hessian.jl

@@ -196,15 +196,15 @@ function (cb::OmegaPlusKDefaultCallback)(info)
                 label_time[1], "  Comment")
         @printf(io, "%s  %s  %s%s  %s  %s%s  %s\n",
                 "-"^4, "-"^7, "-"^6, label_s[2], "-"^10, "-"^7, label_time[2], "-"^15)
-        @printf(io, "%21s%s  %10s  %7.2f%s  %15s\n",
+        @printf(io, "%21s%s  %10s  %7.1f%s  %15s\n",
                 "", label_s[3], "", avgCG, tstr, "Non-interacting")
     elseif info.stage == :iterate
         n_iter  = info.n_iter
         resstr  = format_log8(info.resid_history[n_iter])
         comment = ((n_iter-1) in info.restart_history) ? "Restart" : ""
         sstr    = cb.show_s ? (@sprintf "  %5.2f" info.s) : ""
         restart = length(info.restart_history)
-        @printf(io, "%4i  %7i  %6i%s  %10s  %7.2f%s  %s\n",
+        @printf(io, "%4i  %7i  %6i%s  %10s  %7.1f%s  %s\n",
                 n_iter, restart, info.k, sstr, resstr, avgCG, tstr, comment)
     elseif info.stage == :final
         @printf(io, "%21s%s  %10s  %7.2f%s  %s\n",
@@ -247,7 +247,7 @@ Input parameters:
                                                                eigenvalues; occupation_threshold),
                                     q=zero(Vec3{real(T)}),
                                     maxiter_sternheimer=100,
-                                    maxiter=100, krylovdim=10, s=10.0,
+                                    maxiter=100, krylovdim=20, s=krylovdim,
                                     callback=verbose ? OmegaPlusKDefaultCallback() : identity,
                                     kwargs...) where {T}
     # Using χ04P = -Ω^-1, E extension operator (2P->4P) and R restriction operator:

test/hessian.jl

@@ -120,3 +120,74 @@ end
                         atol=1e-7)
     end
 end
+
+@testitem "Test solve_ΩplusK_split achieves promised accuracies" begin
+    using DFTK
+    using LinearAlgebra
+    using PseudoPotentialData
+
+    function test_solve_ΩplusK(scfres, δVext)
+        # Compute a reference solution
+        δHψ = DFTK.multiply_ψ_by_blochwave(scfres.basis, scfres.ψ, δVext)
+        bandtolalg = DFTK.BandtolGuaranteed(scfres)
+        ref = DFTK.solve_ΩplusK_split(scfres, -δHψ; tol=1e-12, verbose=false, bandtolalg)
+
+        # Test agreement of non-interacting response
+        δρ0 = apply_χ0(scfres, δVext, tol=1e-14).δρ
+        @test maximum(abs, δρ0 - ref.δρ0) < 5e-11
+        @test norm(δρ0 - ref.δρ0)         < 5e-11
+
+        # Check residual is small
+        ε = DFTK.DielectricAdjoint(scfres; bandtolalg=DFTK.BandtolGuaranteed(scfres))
+        εδρ = reshape(DFTK.inexact_mul(ε, ref.δρ; tol=1e-14).Ax, size(δρ0))
+        @test norm(δρ0 - εδρ)         < 1e-12
+        @test maximum(abs, δρ0 - εδρ) < 1e-12
+
+        # TODO δψ agreement not tested
+        @assert length(scfres.basis.kpoints) == 1
+
+        for tol in (1e-2, 1e-4, 1e-6, 1e-8, 1e-10)
+            res = DFTK.solve_ΩplusK_split(scfres, -δHψ; tol, verbose=false)
+            @test norm(res.δρ - ref.δρ) < tol
+
+            # This is probably a bit flaky ... orbitals could rotate
+            @test maximum(abs, res.δψ[1] - ref.δψ[1]) < tol
+        end
+
+        for s in (1, 10, 100, 10^5)
+            tol = 1e-8
+            res = DFTK.solve_ΩplusK_split(scfres, -δHψ; tol, s, verbose=false)
+            @test norm(res.δρ - ref.δρ) < tol
+            @test maximum(abs, res.δψ[1] - ref.δψ[1]) < tol
+        end
+    end
+
+    @testset "Polarisability in Helium" begin
+        # Helium at the centre of the box
+        a = 10.0
+        lattice = a * I(3)  # cube of ``a`` bohrs
+        pseudopotentials = PseudoFamily("cp2k.nc.sr.lda.v0_1.largecore.gth")
+        atoms     = [ElementPsp(:He, pseudopotentials)]
+        positions = [[1/2, 1/2, 1/2]]
+
+        model  = model_DFT(lattice, atoms, positions;
+                           functionals=PBE(), symmetries=false)
+        basis  = PlaneWaveBasis(model; Ecut=30, kgrid=(1, 1, 1))
+        scfres = self_consistent_field(basis; tol=1e-12)
+
+        # δVext is the potential from a uniform field interacting with the dielectric dipole
+        # of the density. Note, that this expression only makes sense for a molecule
+        δVext = [-(r[1] - a/2) for r in r_vectors_cart(basis)]
+        δVext = cat(δVext; dims=4)
+
+        test_solve_ΩplusK(scfres, δVext)
+    end
+
+    @testset "Aluminium atomic displacements" begin
+
+
+        # TODO displace some atoms and solve response
+
+
+    end
+end