LdosMixing for potential-based SCF

src/scf/chi0models.jl

@@ -16,6 +16,9 @@
 """
 @kwdef struct LdosModel <: χ0Model
     adjust_temperature = IncreaseMixingTemperature()
+    characteristic_length = 0  # in Bohr (0 means low-pass filtering is disabled)
+                               # Good values for filtering are between 2 and 0.1 Bohr
+                               # Fluctuations of the LDOS below this length are smoothened.
 end
 function (χ0::LdosModel)(basis::PlaneWaveBasis{T}; eigenvalues, ψ, εF, kwargs...) where {T}
     n_spin = basis.model.n_spin_components
@@ -26,6 +29,16 @@
     ldos = compute_ldos(εF, basis, eigenvalues, ψ; temperature)
     maximum(abs, ldos) < sqrt(eps(T)) && return nothing
 
+    if χ0.characteristic_length > 0
+        wavelength = 2π / χ0.characteristic_length
+        G₀ = 2wavelength / 3
+        Gslope = wavelength / 6
+        f_lowpass(G) = erfc((G - G₀) / Gslope) / 2
+
+        lowpass = [f_lowpass(norm(G)) for G in G_vectors_cart(basis)]
+        ldos = irfft(basis, fft(basis, ldos) .* lowpass)
+    end
+
     tdos = sum(sum, ldos) * basis.dvol  # Integrate LDOS to form total DOS
     function apply!(δρ, δV, α=1)
         δεF = dot(ldos, δV) .* basis.dvol

src/scf/mixing.jl

@@ -188,8 +188,9 @@
 - `verbose`: Run the GMRES in verbose mode.
 - `reltol`: Relative tolerance for GMRES
 """
-function LdosMixing(; adjust_temperature=IncreaseMixingTemperature(), kwargs...)
-    χ0Mixing(; χ0terms=[LdosModel(;adjust_temperature)], kwargs...)
+function LdosMixing(; adjust_temperature=IncreaseMixingTemperature(),
+                    characteristic_length=0, kwargs...)
+    χ0Mixing(; χ0terms=[LdosModel(;adjust_temperature, characteristic_length)], kwargs...)
 end
 
 
@@ -238,8 +239,49 @@
     δρ
 end
 
-@timing "χ0Mixing" function mix_potential(mixing::Mixing, basis::χ0Mixing, δF::AbstractArray; kwargs...)
-    error("Not yet implemented.")
+@timing "χ0Mixing" function mix_potential(mixing::χ0Mixing, basis, δF::AbstractArray; ρin, kwargs...)
+    # Initialise χ0terms and remove nothings (terms that don't yield a contribution)
+    χ0applies = filter(!isnothing, [χ0(basis; ρin, kwargs...) for χ0 in mixing.χ0terms])
+
+    # If no applies left, do not bother running GMRES and directly do simple mixing
+    isempty(χ0applies) && return mix_potential(SimpleMixing(), basis, δF)
+
+    # Note: Since in the potential-mixing version the χ₀ model is directly applied to δF
+    #       (instead of first being "low-pass" filtered by the 1/G² in the Hartree kernel
+    #       like in the density-mixing version), mix_potential is much more susceptible
+    #       to having a good model. For example when using LdosMixing this means one needs
+    #       to choose a good enough k-Point sampling / high enough smearing temperature.
+    #       I also tried experimenting with some low-pass filtering in the LdosModel, but
+    #       so far without a fully satisfactory result. As of now LdosMixing should be avoided
+    #       with potential mixing.
+    @warn("LdosMixing / χ0Mixing not yet fine-tuned for potential mixing. You're on your own. " *
+          "Make sure to use sufficient k-Point sampling and maybe low-pass filtering.", maxlog=1)
+
+    # Solve ε δV = δF with ε = (1 - vc χ₀) and χ₀ given as the sum of the χ0terms
+    devec(x) = reshape(x, size(δF))
+    function dielectric(δF)
+        δF = devec(δF)
+
+        δρ = zero(δF)
+        for apply_term! in χ0applies
+            apply_term!(δρ, δF)  # δρ .+= χ₀ * δF
+        end
+        δρ .-= mean(δρ)
+
+        # Maybe needed ??
+        δρ = symmetrize_ρ(basis, δρ; do_lowpass=true)
+
+        εδF = δF .- apply_kernel(basis, δρ; ρ=ρin, RPA=mixing.RPA)
+        εδF .-= mean(εδF)
+        vec(εδF)
+    end
+
+    DC_δF = mean(δF)
+    δF .-= DC_δF
+    ε  = LinearMap(dielectric, length(δF))
+    δV = devec(gmres(ε, vec(δF), verbose=mixing.verbose, reltol=mixing.reltol))
+    δV .+= DC_δF  # Set DC from δF
+    δV
 end
 
 
@@ -250,19 +292,20 @@
 mixing temperature is equal to the model temperature.
 """
 function IncreaseMixingTemperature(; factor=25, above_ρdiff=1e-2, temperature_max=0.5)
-    function callback(temperature; n_iter, ρin=nothing, ρout=nothing, info...)
+    function callback(temperature; n_iter, history_Δρ=nothing, info...)
         if iszero(temperature) || temperature > temperature_max
             return temperature
-        elseif isnothing(ρin) || isnothing(ρout)
+        elseif isnothing(history_Δρ)
             return temperature
         elseif n_iter ≤ 1
             return factor * temperature
         end
 
         # Continuous piecewise linear function on a logarithmic scale
-        # In [log(above_ρdiff), log(above_ρdiff) + 1] it switches from 1 to factor
-        ρdiff = norm(ρout .- ρin)
-        enhancement = clamp(1 + (factor - 1) * log10(ρdiff / above_ρdiff), 1, factor)
+        # In [log(above_ρdiff), log(above_ρdiff) + switch_slope] it switches from 1 to factor
+        switch_slope = 1
+        ρdiff = last(history_Δρ)
+        enhancement = clamp(1 + (factor - 1) / switch_slope * log10(ρdiff / above_ρdiff), 1, factor)
 
         # Between SCF iterations temperature may never grow
         temperature = clamp(enhancement * temperature, temperature, temperature_max)

src/scf/potential_mixing.jl

@@ -176,10 +176,6 @@ Simple SCF algorithm using potential mixing. Parameters are largely the same as
     accept_step=ScfAcceptStepAll(),
     max_backtracks=3,  # Maximal number of backtracking line searches
 )
-    # TODO Test other mixings and lift this
-    @assert (   mixing isa SimpleMixing
-             || mixing isa KerkerMixing
-             || mixing isa KerkerDosMixing)
     damping isa Number && (damping = FixedDamping(damping))
 
     if !isnothing(ψ)
@@ -202,18 +198,19 @@ Simple SCF algorithm using potential mixing. Parameters are largely the same as
          Vin, Vout=total_local_potential(new_ham), res_V...)
     end
 
+    history_Etot = eltype(ρ)[]
+    history_Δρ   = eltype(ρ)[]
     n_iter    = 1
     converged = false
     ΔEdown    = 0.0
     start_ns  = time_ns()
     α_trial   = trial_damping(damping)
     diagtol   = determine_diagtol(diagtolalg, (; ρin=ρ, Vin=V, n_iter))
     info      = EVρ(V; diagtol, ψ)
-    Pinv_δV   = mix_potential(mixing, basis, info.Vout - info.Vin; n_iter, info...)
+    Pinv_δV   = mix_potential(mixing, basis, info.Vout - info.Vin;
+                              ρin=ρ, n_iter, history_Δρ, info...)
     info      = merge(info, (; α=NaN, diagonalization=[info.diagonalization], ρin=ρ,
                              n_iter, Pinv_δV))
-    history_Etot = eltype(ρ)[]
-    history_Δρ   = eltype(ρ)[]
 
     while n_iter < maxiter
         push!(history_Etot, info.energies.total)
@@ -247,7 +244,7 @@ Simple SCF algorithm using potential mixing. Parameters are largely the same as
 
             info_next    = EVρ(Vnext; ψ=guess, diagtol, info.eigenvalues, info.occupation)
             Pinv_δV_next = mix_potential(mixing, basis, info_next.Vout - info_next.Vin;
-                                         n_iter, info_next...)
+                                         ρin=info_next.ρout, n_iter, history_Δρ, info_next...)
             push!(diagonalization, info_next.diagonalization)
             info_next = merge(info_next, (; α, diagonalization, ρin=info.ρout, n_iter,
                                           Pinv_δV=Pinv_δV_next, history_Δρ, history_Etot ))

src/terms/Hamiltonian.jl

@@ -164,7 +164,6 @@ end
         if Threads.threadid() == 1
             merge!(DFTK.timer, to; tree_point=[t.name for t in DFTK.timer.timer_stack])
         end
-
         synchronize_device(H.basis.architecture)
     end
 

test/scf_compare.jl

@@ -140,16 +140,22 @@ end
 
     for mixing_str in ("KerkerMixing()", "KerkerDosMixing()", "DielectricMixing(; εr=10)",
                        "HybridMixing(; εr=10)", "χ0Mixing(; χ0terms=[Applyχ0Model()], RPA=false)")
-        @testset "Testing $mixing_str" begin
+        @testset "Testing $mixing_str (density)" begin
             mixing = eval(Meta.parse(mixing_str))
             scfres = self_consistent_field(basis; ρ=ρ0, mixing, tol, damping=0.8)
             @test maximum(abs, scfres.ρ - ρ_ref) < 10tol
         end
     end
 
     # Potential mixing
-    scfres = DFTK.scf_potential_mixing(basis; mixing=KerkerMixing(), tol, ρ=ρ0)
-    @test maximum(abs, scfres.ρ - ρ_ref) < 10tol
+    for mixing_str in ("KerkerMixing()", "KerkerDosMixing()", "DielectricMixing(; εr=10)",
+                       "HybridMixing(; εr=10)", "χ0Mixing(; χ0terms=[Applyχ0Model()], RPA=false)")
+        @testset "Testing $mixing_str (potential)" begin
+            mixing = eval(Meta.parse(mixing_str))
+            scfres = DFTK.scf_potential_mixing(basis; ρ=ρ0, mixing, tol, damping=0.8)
+            @test maximum(abs.(scfres.ρ - ρ_ref)) < 10tol
+        end
+    end
 
     # Adaptive potential mixing (started deliberately with the very bad damping
     #          of 1.5 to provoke backtrack steps ... don't do this in production runs!)