Merge branch 'master' into use_pseudo

examples/custom_solvers.jl

@@ -7,7 +7,7 @@ a = 10.26
 lattice = a / 2 * [[0 1 1.];
                    [1 0 1.];
                    [1 1 0.]]
-Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
+Si = ElementPsp(:Si; psp=load_psp("hgh/lda/Si-q4"))
 atoms = [Si, Si]
 positions =  [ones(3)/8, -ones(3)/8]
 
@@ -28,7 +28,7 @@ function my_fp_solver(f, x0, max_iter; tol)
         x = x + mixing_factor * inc
         fx = f(x)
     end
-    (fixpoint=x, converged=norm(fx-x) < tol)
+    (; fixpoint=x, converged=norm(fx-x) < tol)
 end;
 
 # Our eigenvalue solver just forms the dense matrix and diagonalizes
@@ -39,7 +39,7 @@ function my_eig_solver(A, X0; maxiter, tol, kwargs...)
     E = eigen(A)
     λ = E.values[1:n]
     X = E.vectors[:, 1:n]
-    (; λ, X, residual_norms=[], iterations=0, converged=true, n_matvec=0)
+    (; λ, X, residual_norms=[], n_iter=0, converged=true, n_matvec=0)
 end;
 
 # Finally we also define our custom mixing scheme. It will be a mixture

examples/geometry_optimization.jl

@@ -13,7 +13,7 @@ Ecut = 5                # kinetic energy cutoff in Hartree
 tol = 1e-8              # tolerance for the optimization routine
 a = 10                  # lattice constant in Bohr
 lattice = a * I(3)
-H = ElementPsp(:H, psp=load_psp("hgh/lda/h-q1"));
+H = ElementPsp(:H; psp=load_psp("hgh/lda/h-q1"));
 atoms = [H, H];
 
 # We define a Bloch wave and a density to be used as global variables so that we

examples/supercells.jl

@@ -16,7 +16,7 @@ using ASEconvert
 function aluminium_setup(repeat=1; Ecut=7.0, kgrid=[2, 2, 2])
     a = 7.65339
     lattice = a * Matrix(I, 3, 3)
-    Al = ElementPsp(:Al, psp=load_psp("hgh/lda/al-q3"))
+    Al = ElementPsp(:Al; psp=load_psp("hgh/lda/al-q3"))
     atoms     = [Al, Al, Al, Al]
     positions = [[0.0, 0.0, 0.0], [0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]
     unit_cell = periodic_system(lattice, atoms, positions)
@@ -29,7 +29,7 @@ function aluminium_setup(repeat=1; Ecut=7.0, kgrid=[2, 2, 2])
 
     ## Unfortunately right now the conversion to ASE drops the pseudopotential information,
     ## so we need to reattach it:
-    supercell = attach_psp(supercell, Al="hgh/lda/al-q3")
+    supercell = attach_psp(supercell; Al="hgh/lda/al-q3")
 
     ## Construct an LDA model and discretise
     ## Note: We disable symmetries explicitly here. Otherwise the problem sizes

src/eigen/diag.jl

@@ -29,8 +29,8 @@ function diagonalize_all_kblocks(eigensolver, ham::Hamiltonian, nev_per_kpoint::
         # Get ψguessk
         if !isnothing(ψguess)
             if n_Gk != size(ψguess[ik], 1)
-                error("Mismatch in dimension between guess ($(size(ψguess, 1)) and " *
-                      "Hamiltonian $n_Gk")
+                error("Mismatch in dimension between guess ($(size(ψguess[ik], 1)) and " *
+                      "Hamiltonian ($n_Gk)")
             end
             nev_guess = size(ψguess[ik], 2)
             if nev_guess > nev_per_kpoint
@@ -64,10 +64,10 @@ function diagonalize_all_kblocks(eigensolver, ham::Hamiltonian, nev_per_kpoint::
     # Transform results into a nicer datastructure
     # TODO It feels inconsistent to put λ onto the CPU here but none of the other objects.
     #      Better have this handled by the caller of diagonalize_all_kblocks.
-    (λ=[to_cpu(real.(res.λ)) for res in results],  # Always get onto the CPU
+    (; λ=[to_cpu(real.(res.λ)) for res in results],  # Always get onto the CPU
      X=[res.X for res in results],
      residual_norms=[res.residual_norms for res in results],
-     iterations=[res.iterations for res in results],
+     n_iter=[res.n_iter for res in results],
      converged=all(res.converged for res in results),
      n_matvec=sum(res.n_matvec for res in results))
 end
@@ -77,7 +77,7 @@ Function to select a subset of eigenpairs on each ``k``-Point. Works on the
 Tuple returned by `diagonalize_all_kblocks`.
 """
 function select_eigenpairs_all_kblocks(eigres, range)
-    merge(eigres, (λ=[λk[range] for λk in eigres.λ],
+    merge(eigres, (; λ=[λk[range] for λk in eigres.λ],
                    X=[Xk[:, range] for Xk in eigres.X],
                    residual_norms=[resk[range] for resk in eigres.residual_norms]))
 end

src/eigen/diag_full.jl

@@ -6,7 +6,7 @@ function diag_full(A, X0; kwargs...)
     λ = E.values[1:Neig]
     (; λ, X,
      residual_norms=zeros(Neig),
-     iterations=0,
+     n_iter=0,
      converged=true,
      n_matvec=0)
 end

src/eigen/diag_lobpcg_hyper.jl

@@ -10,7 +10,7 @@ function lobpcg_hyper(A, X0; maxiter=100, prec=nothing,
 
     n_conv_check === nothing && (n_conv_check = size(X0, 2))
     converged = maximum(result.residual_norms[1:n_conv_check]) < tol
-    iterations = size(result.residual_history, 2) - 1
+    n_iter = size(result.residual_history, 2) - 1
 
-    merge(result, (; iterations, converged))
+    merge(result, (; n_iter, converged))
 end

src/postprocess/band_structure.jl

@@ -95,7 +95,7 @@ end
     eigres = diagonalize_all_kblocks(eigensolver, ham, n_bands + 3;
                                      n_conv_check=n_bands, tol, show_progress, kwargs...)
     if !eigres.converged
-        @warn "Eigensolver not converged" iterations=band_data.iterations
+        @warn "Eigensolver not converged" n_iter=band_data.n_iter
     end
     merge((; basis=bs_basis), select_eigenpairs_all_kblocks(eigres, 1:n_bands))
 end

src/response/cg.jl

@@ -58,7 +58,7 @@ function cg!(x::AbstractVector{T}, A::LinearMap{T}, b::AbstractVector{T};
         β = γ / γ_prev
         p .= proj(c .+ β .* p)
     end
-    info = (; x, converged, tol, residual_norm, iterations=n_iter, maxiter, stage=:finalize)
+    info = (; x, converged, tol, residual_norm, n_iter, maxiter, stage=:finalize)
     callback(info)
     info
 end

src/response/chi0.jl

@@ -178,7 +178,7 @@ function sternheimer_solver(Hk, ψk, ε, rhs;
     J = LinearMap{eltype(ψk)}(RAR, size(Hk, 1))
     res = cg(J, bb; precon=FunctionPreconditioner(R_ldiv!), tol, proj=R,
              callback=cg_callback)
-    !res.converged && @warn("Sternheimer CG not converged", res.iterations,
+    !res.converged && @warn("Sternheimer CG not converged", res.n_iter,
                             res.tol, res.residual_norm)
     δψknᴿ = res.x
     info = (; basis, kpoint, res)

src/response/hessian.jl

@@ -116,7 +116,7 @@ function solve_ΩplusK(basis::PlaneWaveBasis{T}, ψ, rhs, occupation;
     res = cg(J, rhs_pack; precon=FunctionPreconditioner(f_ldiv!), proj, tol,
              callback)
     (; δψ=unpack(res.x), res.converged, res.tol, res.residual_norm,
-     res.iterations)
+     res.n_iter)
 end
 
 

src/scf/mixing.jl

@@ -248,7 +248,7 @@ is increased by a `factor`, which is then smoothly lowered towards the temperatu
 within the model as the SCF converges. Once the density change is below `above_ρdiff` the
 mixing temperature is equal to the model temperature.
 """
-function IncreaseMixingTemperature(;factor=25, above_ρdiff=1e-2, temperature_max=0.5)
+function IncreaseMixingTemperature(; factor=25, above_ρdiff=1e-2, temperature_max=0.5)
     function callback(temperature; n_iter, ρin=nothing, ρout=nothing, info...)
         if iszero(temperature) || temperature > temperature_max
             return temperature

src/scf/scf_callbacks.jl

@@ -29,7 +29,7 @@ function ScfDefaultCallback(; show_damping=true, show_time=true)
             # Average number of diagonalizations per k-point needed for this SCF step
             # Note: If two Hamiltonian diagonalizations have been used (e.g. adaptive damping),
             # the per k-point values are summed.
-            diagiter = mpi_mean(sum(mean(diag.iterations) for diag in info.diagonalization),
+            diagiter = mpi_mean(sum(mean(diag.n_iter) for diag in info.diagonalization),
                                 info.basis.comm_kpts)
         end
 
@@ -135,7 +135,7 @@ Determine the tolerance used for the next diagonalization. This function takes
 ensuring additionally that the returned value is between `diagtol_min` and `diagtol_max`
 and never increases.
 """
-function ScfDiagtol(;ratio_ρdiff=0.2, diagtol_min=nothing, diagtol_max=0.03)
+function ScfDiagtol(; ratio_ρdiff=0.2, diagtol_min=nothing, diagtol_max=0.03)
     function determine_diagtol(info)
         isnothing(diagtol_min) && (diagtol_min = 100eps(real(eltype(info.ρin))))
         info.n_iter ≤ 1 && return diagtol_max

src/scf/self_consistent_field.jl

@@ -32,7 +32,7 @@ function next_density(ham::Hamiltonian,
 
     eigres = diagonalize_all_kblocks(eigensolver, ham, n_bands_compute;
                                      ψguess=ψ, n_conv_check=n_bands_converge, kwargs...)
-    eigres.converged || (@warn "Eigensolver not converged" iterations=eigres.iterations)
+    eigres.converged || (@warn "Eigensolver not converged" n_iter=eigres.n_iter)
 
     # Check maximal occupation of the unconverged bands is sensible.
     occupation, εF = compute_occupation(ham.basis, eigres.λ, fermialg;

src/terms/ewald.jl

@@ -47,7 +47,7 @@ function energy_forces_ewald(lattice::AbstractArray{T}, charges, positions;
     # TODO should something more clever be done here? For now
     # we assume that we are not interested in the Ewald
     # energy of non-3D systems
-    any(iszero.(eachcol(lattice))) && return zero(T)
+    any(iszero.(eachcol(lattice))) && return (; energy=zero(T), forces=zero(positions))
 
     recip_lattice = compute_recip_lattice(lattice)
     @assert length(charges) == length(positions)

test/compute_jacobian_eigen.jl

@@ -37,7 +37,7 @@ if mpi_nprocs() == 1  # Distributed implementation not yet available
 
         # random starting point on the tangent space to avoid eigenvalue 0
         n_bands = size(ψ[1], 2)
-        x0 = map(1:numval) do n
+        x0 = map(1:numval) do _
             initial = map(basis.kpoints) do kpt
                 n_Gk = length(G_vectors(basis, kpt))
                 randn(Complex{eltype(basis)}, n_Gk, n_bands)
@@ -59,7 +59,7 @@ if mpi_nprocs() == 1  # Distributed implementation not yet available
         J = LinearMap{T}(ΩpK, size(x0, 1))
 
         # compute smallest eigenvalue of Ω with LOBPCG
-        lobpcg_hyper(J, x0, prec=FunctionPreconditioner(f_ldiv!), tol=1e-7)
+        lobpcg_hyper(J, x0; prec=FunctionPreconditioner(f_ldiv!), tol=1e-7)
     end
 
     @testset "Compute eigenvalues" begin

test/lobpcg.jl

@@ -27,14 +27,14 @@ include("./testcases.jl")
 
     @test length(ref_λ) == length(silicon.kcoords)
     @testset "without Preconditioner" begin
-        res = diagonalize_all_kblocks(lobpcg_hyper, ham, nev_per_k, tol=1e-4,
+        res = diagonalize_all_kblocks(lobpcg_hyper, ham, nev_per_k; tol=1e-4,
                                       prec_type=nothing, interpolate_kpoints=false)
 
         @test res.converged
         for ik in 1:length(basis.kpoints)
             @test ref_λ[basis.krange_thisproc[ik]] ≈ res.λ[ik] atol = 1e-4
             @test maximum(res.residual_norms[ik]) < 1e-4
-            @test res.iterations[ik] < 200
+            @test res.n_iter[ik] < 200
         end
     end
 
@@ -46,7 +46,7 @@ include("./testcases.jl")
         for ik in 1:length(basis.kpoints)
             @test ref_λ[basis.krange_thisproc[ik]] ≈ res.λ[ik]
             @test maximum(res.residual_norms[ik]) < 100tol  # TODO Why the 100?
-            @test res.iterations[ik] < 50
+            @test res.n_iter[ik] < 50
         end
     end
 end
@@ -56,13 +56,13 @@ if !isdefined(Main, :FAST_TESTS) || !FAST_TESTS
         Ecut = 25
         fft_size = [33, 33, 33]
 
-        Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/lda/si-q4"))
+        Si = ElementPsp(silicon.atnum; psp=load_psp("hgh/lda/si-q4"))
         model = Model(silicon.lattice, silicon.atoms, silicon.positions;
                       terms=[Kinetic(),AtomicLocal()])
         basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.kweights; fft_size)
         ham = Hamiltonian(basis)
 
-        res = diagonalize_all_kblocks(lobpcg_hyper, ham, 6, tol=1e-8)
+        res = diagonalize_all_kblocks(lobpcg_hyper, ham, 6; tol=1e-8)
 
         ref = [
             [-4.087198659513310, -4.085326314828677, -0.506869382308294,
@@ -84,14 +84,14 @@ end
     Ecut = 10
     fft_size = [21, 21, 21]
 
-    Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/lda/si-q4"))
+    Si = ElementPsp(silicon.atnum; psp=load_psp("hgh/lda/si-q4"))
     model = Model(silicon.lattice, silicon.atoms, silicon.positions;
                   terms=[Kinetic(), AtomicLocal(), AtomicNonlocal()])
 
     basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.kweights; fft_size)
     ham = Hamiltonian(basis)
 
-    res = diagonalize_all_kblocks(lobpcg_hyper, ham, 5, tol=1e-8, interpolate_kpoints=false)
+    res = diagonalize_all_kblocks(lobpcg_hyper, ham, 5; tol=1e-8, interpolate_kpoints=false)
     ref = [
         [0.067955741977536, 0.470244204908046, 0.470244204920801,
          0.470244204998022, 0.578392222232969],
@@ -115,8 +115,8 @@ end
     basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.kweights)
     ham = Hamiltonian(basis; ρ=guess_density(basis))
 
-    res1 = diagonalize_all_kblocks(lobpcg_hyper, ham, 5, tol=1e-8, interpolate_kpoints=false)
-    res2 = diagonalize_all_kblocks(diag_full, ham, 5, tol=1e-8, interpolate_kpoints=false)
+    res1 = diagonalize_all_kblocks(lobpcg_hyper, ham, 5; tol=1e-8, interpolate_kpoints=false)
+    res2 = diagonalize_all_kblocks(diag_full, ham, 5; tol=1e-8, interpolate_kpoints=false)
     for ik in 1:length(basis.kpoints)
         @test res1.λ[ik] ≈ res2.λ[ik] atol=1e-6
     end