diff --git a/src/scf/chi0models.jl b/src/scf/chi0models.jl index d7c097f347..a4e56f802d 100644 --- a/src/scf/chi0models.jl +++ b/src/scf/chi0models.jl @@ -16,6 +16,9 @@ For details see [Herbst, Levitt 2020](https://arxiv.org/abs/2009.01665). """ @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 @@ function (χ0::LdosModel)(basis::PlaneWaveBasis{T}; eigenvalues, ψ, εF, kwargs 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 diff --git a/src/scf/mixing.jl b/src/scf/mixing.jl index b1ed7d67fe..03d9e40e9a 100644 --- a/src/scf/mixing.jl +++ b/src/scf/mixing.jl @@ -188,8 +188,9 @@ Important `kwargs` passed on to [`χ0Mixing`](@ref) - `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 δρ 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 @@ within the model as the SCF converges. Once the density change is below `above_ 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) diff --git a/src/scf/potential_mixing.jl b/src/scf/potential_mixing.jl index 89a5335ee1..88eedf456b 100644 --- a/src/scf/potential_mixing.jl +++ b/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,6 +198,8 @@ 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 @@ -209,11 +207,10 @@ Simple SCF algorithm using potential mixing. Parameters are largely the same as α_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 )) diff --git a/src/terms/Hamiltonian.jl b/src/terms/Hamiltonian.jl index 712ac0fb48..c2c1ac94ad 100644 --- a/src/terms/Hamiltonian.jl +++ b/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 diff --git a/test/scf_compare.jl b/test/scf_compare.jl index 022f8ecc7e..211a94a003 100644 --- a/test/scf_compare.jl +++ b/test/scf_compare.jl @@ -140,7 +140,7 @@ 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 @@ -148,8 +148,14 @@ 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!)