Skip to content
Permalink
Browse files

Merge pull request #126 from JuliaMolSim/flexible-convergence

First go at flexible SCF convergence criteria
  • Loading branch information
mfherbst committed Feb 14, 2020
2 parents d2a8438 + e0fb004 commit 0826dafc93f75750ae1fbbadf54128af1d38549e
Showing with 46 additions and 17 deletions.
  1. +2 −1 src/densities.jl
  2. +6 −6 src/scf/scf_solvers.jl
  3. +33 −5 src/scf/self_consistent_field.jl
  4. +5 −5 test/scf_compare.jl
@@ -117,7 +117,8 @@ This interpolation uses a very basic real-space algorithm, and makes
a DWIM-y attempt to take into account the fact that b_out can be a supercell of b_in
"""
function interpolate_density(ρ_in::RealFourierArray, b_out::PlaneWaveBasis)
ρ_out = interpolate_density(ρ_in.real, ρ_in.basis.fft_size, b_out.fft_size, ρ_in.basis.lattice, b_out.lattice)
ρ_out = interpolate_density(ρ_in.real, ρ_in.basis.fft_size, b_out.fft_size,
ρ_in.basis.model.lattice, b_out.model.lattice)
from_real(b_out, ρ_out)
end

@@ -14,7 +14,7 @@ fixed-point scheme, keeping `m` steps for Anderson acceleration. See the
NLSolve documentation for details about the other parameters and methods.
"""
function scf_nlsolve_solver(m=5, method=:anderson; kwargs...)
function fp_solver(f, x0, tol, max_iter)
function fp_solver(f, x0, max_iter; tol=1e-6)
res = nlsolve(x -> f(x) - x, x0; method=method, m=m, xtol=tol,
ftol=0.0, show_trace=false, iterations=max_iter, kwargs...)
(fixpoint=res.zero, converged=converged(res))
@@ -27,13 +27,13 @@ Create a damped SCF solver updating the density as
`x = β * x_new + (1 - β) * x`
"""
function scf_damping_solver=0.2)
function fp_solver(f, x0, tol, max_iter)
function fp_solver(f, x0, max_iter; tol=1e-6)
converged = false
x = copy(x0)
for i in 1:max_iter
x_new = f(x)

if 20 * norm(x_new - x) < tol
if norm(x_new - x) < tol
x = x_new
converged = true
break
@@ -56,7 +56,7 @@ function anderson(f, x0, m::Int, max_iter::Int, tol::Real, warming=0)

# Cheat support for multidimensional arrays
if length(size(x0)) != 1
x, conv= anderson(x -> vec(f(reshape(x, size(x0)...))), vec(x0), m, max_iter, tol, warming)
x, conv = anderson(x -> vec(f(reshape(x, size(x0)...))), vec(x0), m, max_iter, tol, warming)
return (fixpoint=reshape(x, size(x0)...), converged=conv)
end
N = size(x0, 1)
@@ -94,7 +94,7 @@ function anderson(f, x0, m::Int, max_iter::Int, tol::Real, warming=0)
(fixpoint=xs[:,1], converged=err < tol)
end
function scf_anderson_solver(m=5)
(f, x0, tol, max_iter) -> anderson(x -> f(x) - x, x0, m, max_iter, tol)
(f, x0, max_iter; tol=1e-6) -> anderson(x -> f(x) - x, x0, m, max_iter, tol)
end

"""
@@ -157,4 +157,4 @@ function CROP(f, x0, m::Int, max_iter::Int, tol::Real, warming=0)
end
(fixpoint=xs[:, 1], converged=err < tol)
end
scf_CROP_solver(m=5) = (f, x0, tol, max_iter) -> CROP(x -> f(x) - x, x0, m, max_iter, tol)
scf_CROP_solver(m=5) = (f, x0, max_iter; tol=1e-6) -> CROP(x -> f(x) - x, x0, m, max_iter, tol)
@@ -35,6 +35,28 @@ function scf_default_callback(info)
@printf "%3d %-15.12f %E\n" neval E res
end

"""
Flag convergence as soon as total energy change drops below tolerance
"""
function scf_convergence_energy_difference(tolerance)
energy_total = NaN

function is_converged(info)
etot_old = energy_total
energy_total = sum(values(info.energies))
abs(energy_total - etot_old) < tolerance
end
return is_converged
end

"""
Flag convergence by using the L2Norm of the change between
input density and unpreconditioned output density (ρout)
"""
function scf_convergence_density_difference(tolerance)
info -> norm(info.ρout.fourier - info.ρin.fourier) < tolerance
end


"""
Solve the KS equations with a SCF algorithm, starting at `ham`
@@ -43,7 +65,8 @@ function self_consistent_field(ham::Hamiltonian, n_bands;
Psi=nothing, tol=1e-6, max_iter=100,
solver=scf_nlsolve_solver(),
eigensolver=lobpcg_hyper, n_ep_extra=3, diagtol=tol / 10,
mixing=SimpleMixing(), callback=scf_default_callback)
mixing=SimpleMixing(), callback=scf_default_callback,
is_converged=scf_convergence_energy_difference(tol))
T = real(eltype(ham.density))
basis = ham.basis
model = basis.model
@@ -77,14 +100,19 @@ function self_consistent_field(ham::Hamiltonian, n_bands;
# mix it with ρin to get a proposal step
ρnext = mix(mixing, basis, ρin, ρout)
neval += 1
callback((ham=ham, energies=energies, ρin=ρin, ρout=ρout, ρnext=ρnext,
orben=orben, occupation=occupation, εF=εF, neval=neval))

info = (ham=ham, energies=energies, ρin=ρin, ρout=ρout, ρnext=ρnext,
orben=orben, occupation=occupation, εF=εF, neval=neval)
callback(info)
is_converged(info) && return x

ρnext.real
end

fpres = solver(fixpoint_map, ρout.real, tol, max_iter)
# We do not use the return value of fpres but rather the one that got updated by fixpoint_map
fpres = solver(fixpoint_map, ρout.real, max_iter; tol=10eps(T))
# Tolerance is only dummy here: Convergence is flagged by is_converged
# inside the fixpoint_map. Also we do not use the return value of fpres but rather the
# one that got updated by fixpoint_map

energies = update_energies(ham, Psi, occupation, ρout)

@@ -7,7 +7,7 @@ include("testcases.jl")
Ecut = 3
n_bands = 6
fft_size = [9, 9, 9]
tol = 1e-6
tol = 1e-7

Si = Element(silicon.atnum, psp=load_psp(silicon.psp))
model = model_dft(silicon.lattice, :lda_xc_teter93, [Si => silicon.positions])
@@ -20,9 +20,9 @@ include("testcases.jl")

# Run DM
println("\nTesting direct minimization")
dmres = direct_minimization(basis; g_tol=1e-8)
dmres = direct_minimization(basis; g_tol=tol)
ρ_dm = dmres.ρ.fourier
@test maximum(abs.(ρ_dm - ρ_nl)) < 30tol
@test maximum(abs.(ρ_dm - ρ_nl)) < sqrt(tol) / 10

# Run other SCFs with SAD guess
ρ0 = guess_density(basis, [Si => silicon.positions])
@@ -31,14 +31,14 @@ include("testcases.jl")
println("\nTesting $solver")
scfres = self_consistent_field(Hamiltonian(basis, ρ0), n_bands, tol=tol, solver=solver())
ρ_alg = scfres.ρ.fourier
@test maximum(abs.(ρ_alg - ρ_nl)) < 30tol
@test maximum(abs.(ρ_alg - ρ_nl)) < sqrt(tol) / 10
end

# Run other mixing with nlsolve (the others are too slow...)
for mixing in (KerkerMixing(), SimpleMixing(), SimpleMixing(.5))
println("\n Testing $mixing")
scfres = self_consistent_field(Hamiltonian(basis, ρ0), n_bands, tol=tol, solver=scf_nlsolve_solver(), mixing=mixing)
ρ_alg = scfres.ρ.fourier
@test maximum(abs.(ρ_alg - ρ_nl)) < 30tol
@test maximum(abs.(ρ_alg - ρ_nl)) < sqrt(tol) / 10
end
end

0 comments on commit 0826daf

Please sign in to comment.
You can’t perform that action at this time.