diff --git a/.github/workflows/documentation.yaml b/.github/workflows/documentation.yaml
index 790ddfd1da..adb8afadfc 100644
--- a/.github/workflows/documentation.yaml
+++ b/.github/workflows/documentation.yaml
@@ -44,3 +44,9 @@ jobs:
             using DFTK
             DocMeta.setdocmeta!(DFTK, :DocTestSetup, :(using DFTK); recursive=true)
             doctest(DFTK)'
+      - name: Upload Docs
+        uses: actions/upload-artifact@v3
+        with:
+          name: documentation
+          path: docs/build
+          retention-days: 1
diff --git a/docs/make.jl b/docs/make.jl
index 476f816d2d..3254a79900 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -90,7 +90,7 @@ EXAMPLE_ASSETS = ["examples/Fe_afm.pwi", "examples/Si.extxyz"]
 #
 # Configuration and setup
 #
-DEBUG = false  # Set to true to disable some checks and cleanup
+DEBUG = false  # set to `true` to disable some checks and cleanup
 
 import LibGit2
 import Pkg
@@ -107,7 +107,7 @@ DFTKREPO   = DFTKGH * ".git"
 # Setup julia dependencies for docs generation if not yet done
 Pkg.activate(@__DIR__)
 if !isfile(joinpath(@__DIR__, "Manifest.toml"))
-    Pkg.develop(Pkg.PackageSpec(path=ROOTPATH))
+    Pkg.develop(Pkg.PackageSpec(; path=ROOTPATH))
     Pkg.instantiate()
 end
 
@@ -152,9 +152,9 @@ end
 # The examples go to docs/literate_build/examples, the .jl files stay where they are
 literate_files = map(filter!(endswith(".jl"), extract_paths(PAGES))) do file
     if startswith(file, "examples/")
-        (src=joinpath(ROOTPATH, file), dest=joinpath(SRCPATH, "examples"), example=true)
+        (; src=joinpath(ROOTPATH, file), dest=joinpath(SRCPATH, "examples"), example=true)
     else
-        (src=joinpath(SRCPATH, file), dest=joinpath(SRCPATH, dirname(file)), example=false)
+        (; src=joinpath(SRCPATH, file), dest=joinpath(SRCPATH, dirname(file)), example=false)
     end
 end
 
@@ -186,6 +186,7 @@ for file in literate_files
 end
 
 # Generate the docs in BUILDPATH
+remote_args = CONTINUOUS_INTEGRATION ? (; ) : (; remotes=nothing)
 makedocs(;
     modules=[DFTK],
     format=Documenter.HTML(
@@ -208,6 +209,7 @@ makedocs(;
     pages=transform_to_md(PAGES),
     checkdocs=:exports,
     warnonly=DEBUG,
+    remote_args...,
 )
 
 # Dump files for managing dependencies in binder
diff --git a/docs/src/assets/0_pregenerate.jl b/docs/src/assets/0_pregenerate.jl
index 9810a9f7c8..082d73dda8 100644
--- a/docs/src/assets/0_pregenerate.jl
+++ b/docs/src/assets/0_pregenerate.jl
@@ -1,7 +1,7 @@
 using MKL
 using DFTK
 using LinearAlgebra
-setup_threading(n_blas=2)
+setup_threading(; n_blas=2)
 
 let
     include("../../../../examples/convergence_study.jl")
@@ -61,7 +61,7 @@ let
     conv_hgh = converge_Ecut(Ecuts, psp_hgh, tol)
     println("HGH: $(conv_hgh.Ecut_conv)")
 
-    plt = plot(yaxis=:log10, xlabel="Ecut [Eh]", ylabel="Error [Eh]")
+    plt = plot(; yaxis=:log10, xlabel="Ecut [Eh]", ylabel="Error [Eh]")
     plot!(plt, conv_hgh.Ecuts, conv_hgh.errors, label="HGH",
           markers=true, linewidth=3)
     plot!(plt, conv_upf.Ecuts, conv_upf.errors, label="PseudoDojo NC SR LDA UPF",
diff --git a/examples/anyons.jl b/examples/anyons.jl
index a149d21d5e..8c4fd33aea 100644
--- a/examples/anyons.jl
+++ b/examples/anyons.jl
@@ -29,5 +29,5 @@ scfres = direct_minimization(basis, tol=1e-14)  # Reduce tol for production
 E = scfres.energies.total
 s = 2
 E11 = π/2 * (2(s+1)/s)^((s+2)/s) * (s/(s+2))^(2(s+1)/s) * E^((s+2)/s) / β
-println("e(1,1) / (2π)= ", E11 / (2π))
+println("e(1,1) / (2π) = ", E11 / (2π))
 heatmap(scfres.ρ[:, :, 1, 1], c=:blues)
diff --git a/examples/custom_potential.jl b/examples/custom_potential.jl
index 7f8456a840..e33c0fe8a7 100644
--- a/examples/custom_potential.jl
+++ b/examples/custom_potential.jl
@@ -74,6 +74,8 @@ tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dime
 # Extract other quantities before plotting them
 ρ = scfres.ρ[:, 1, 1, 1]        # converged density, first spin component
 ψ_fourier = scfres.ψ[1][:, 1]   # first k-point, all G components, first eigenvector
+
+# Transform the wave function to real space and fix the phase:
 ψ = ifft(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]
 ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));
 
diff --git a/examples/error_estimates_forces.jl b/examples/error_estimates_forces.jl
index 4f643bb312..7e7000ca2b 100644
--- a/examples/error_estimates_forces.jl
+++ b/examples/error_estimates_forces.jl
@@ -54,8 +54,7 @@ tol = 1e-5;
 # We compute the reference solution ``P_*`` from which we will compute the
 # references forces.
 scfres_ref = self_consistent_field(basis_ref; tol, callback=identity)
-ψ_ref, _ = DFTK.select_occupied_orbitals(basis_ref, scfres_ref.ψ,
-                                         scfres_ref.occupation);
+ψ_ref = DFTK.select_occupied_orbitals(basis_ref, scfres_ref.ψ, scfres_ref.occupation).ψ;
 
 # We compute a variational approximation of the reference solution with
 # smaller `Ecut`. `ψr`, `ρr` and `Er` are the quantities computed with `Ecut`
@@ -79,7 +78,7 @@ Er, hamr = energy_hamiltonian(basis_ref, ψr, scfres.occupation; ρ=ρr);
 #   in [`src/scf/newton.jl`](https://github.com/JuliaMolSim/DFTK.jl/blob/fedc720dab2d194b30d468501acd0f04bd4dd3d6/src/scf/newton.jl#L121).
 res = DFTK.compute_projected_gradient(basis_ref, ψr, scfres.occupation)
 res, occ = DFTK.select_occupied_orbitals(basis_ref, res, scfres.occupation)
-ψr, _ = DFTK.select_occupied_orbitals(basis_ref, ψr, scfres.occupation);
+ψr = DFTK.select_occupied_orbitals(basis_ref, ψr, scfres.occupation).ψ;
 
 # - Compute the error ``P-P_*`` on the associated orbitals ``ϕ-ψ`` after aligning
 #   them: this is done by solving ``\min |ϕ - ψU|`` for ``U`` unitary matrix of
@@ -94,7 +93,7 @@ function compute_error(basis, ϕ, ψ)
 end
 err = compute_error(basis_ref, ψr, ψ_ref);
 
-# - Compute ``{\bm M}^{-1}R(P)`` with ``{\bm M}^{-1}`` defined in [^CDKL2021]:
+# - Compute ``{\boldsymbol M}^{-1}R(P)`` with ``{\boldsymbol M}^{-1}`` defined in [^CDKL2021]:
 P = [PreconditionerTPA(basis_ref, kpt) for kpt in basis_ref.kpoints]
 map(zip(P, ψr)) do (Pk, ψk)
     DFTK.precondprep!(Pk, ψk)
@@ -115,7 +114,7 @@ function apply_inv_M(φk, Pk, δφnk, n)
         DFTK.proj_tangent_kpt!(x, φk)
     end
     J = LinearMap{eltype(φk)}(op, size(δφnk, 1))
-    δφnk = cg(J, δφnk, Pl=DFTK.FunctionPreconditioner(f_ldiv!),
+    δφnk = cg(J, δφnk; Pl=DFTK.FunctionPreconditioner(f_ldiv!),
               verbose=false, reltol=0, abstol=1e-15)
     DFTK.proj_tangent_kpt!(δφnk, φk)
 end
@@ -136,8 +135,8 @@ Mres = apply_metric(ψr, P, res, apply_inv_M);
 #
 # ```math
 # \begin{bmatrix}
-# (\bm Ω + \bm K)_{11} & (\bm Ω + \bm K)_{12} \\
-# 0 & {\bm M}_{22}
+# (\boldsymbol Ω + \boldsymbol K)_{11} & (\boldsymbol Ω + \boldsymbol K)_{12} \\
+# 0 & {\boldsymbol M}_{22}
 # \end{bmatrix}
 # \begin{bmatrix}
 # P_{1} - P_{*1} \\ P_{2}-P_{*2}
@@ -151,7 +150,7 @@ Mres = apply_metric(ψr, P, res, apply_inv_M);
 resLF = DFTK.transfer_blochwave(res, basis_ref, basis)
 resHF = res - DFTK.transfer_blochwave(resLF, basis, basis_ref);
 
-# - Compute ``{\bm M}^{-1}_{22}R_2(P)``:
+# - Compute ``{\boldsymbol M}^{-1}_{22}R_2(P)``:
 e2 = apply_metric(ψr, P, resHF, apply_inv_M);
 
 # - Compute the right hand side of the Schur system:
@@ -166,9 +165,9 @@ end
 rhs = resLF - ΩpKe2;
 
 # - Solve the Schur system to compute ``R_{\rm Schur}(P)``: this is the most
-#   costly step, but inverting ``\bm{Ω} + \bm{K}`` on the small space has
+#   costly step, but inverting ``\boldsymbol{Ω} + \boldsymbol{K}`` on the small space has
 #   the same cost than the full SCF cycle on the small grid.
-ψ, _ = DFTK.select_occupied_orbitals(basis, scfres.ψ, scfres.occupation)
+(; ψ) = DFTK.select_occupied_orbitals(basis, scfres.ψ, scfres.occupation)
 e1 = DFTK.solve_ΩplusK(basis, ψ, rhs, occ; tol).δψ
 e1 = DFTK.transfer_blochwave(e1, basis, basis_ref)
 res_schur = e1 + Mres;
diff --git a/examples/geometry_optimization.jl b/examples/geometry_optimization.jl
index 0820d6912b..8f27b6464b 100644
--- a/examples/geometry_optimization.jl
+++ b/examples/geometry_optimization.jl
@@ -68,7 +68,7 @@ end;
 
 x0 = vcat(lattice \ [0., 0., 0.], lattice \ [1.4, 0., 0.])
 xres = optimize(Optim.only_fg!(fg!), x0, LBFGS(),
-                Optim.Options(show_trace=true, f_tol=tol))
+                Optim.Options(; show_trace=true, f_tol=tol))
 xmin = Optim.minimizer(xres)
 dmin = norm(lattice*xmin[1:3] - lattice*xmin[4:6])
 @printf "\nOptimal bond length for Ecut=%.2f: %.3f Bohr\n" Ecut dmin
diff --git a/examples/gross_pitaevskii.jl b/examples/gross_pitaevskii.jl
index 33dd8c1709..1f298bb0ab 100644
--- a/examples/gross_pitaevskii.jl
+++ b/examples/gross_pitaevskii.jl
@@ -94,8 +94,8 @@ H = ham.blocks[1];
 
 # `H` can be used as a linear operator (efficiently using FFTs), or converted to a dense matrix:
 ψ11 = scfres.ψ[1][:, 1] # first k-point, first eigenvector
-Hmat = Array(H) # This is now just a plain Julia matrix,
-##                which we can compute and store in this simple 1D example
+Hmat = Array(H)  # This is now just a plain Julia matrix,
+##                  which we can compute and store in this simple 1D example
 @assert norm(Hmat * ψ11 - H * ψ11) < 1e-10
 
 # Let's check that ψ11 is indeed an eigenstate:
diff --git a/examples/scf_callbacks.jl b/examples/scf_callbacks.jl
index 78129825f7..d7c1e194d5 100644
--- a/examples/scf_callbacks.jl
+++ b/examples/scf_callbacks.jl
@@ -39,7 +39,7 @@ basis  = PlaneWaveBasis(model; Ecut=5, kgrid=[3, 3, 3]);
 # has finished. For this we first define the empty plot canvas
 # and an empty container for all the density differences:
 using Plots
-p = plot(yaxis=:log)
+p = plot(; yaxis=:log)
 density_differences = Float64[];
 
 # The callback function itself gets passed a named tuple
diff --git a/ext/DFTKWriteVTK.jl b/ext/DFTKWriteVTK.jl
index 6cc85364f2..044ccfdf9b 100644
--- a/ext/DFTKWriteVTK.jl
+++ b/ext/DFTKWriteVTK.jl
@@ -19,8 +19,8 @@ function DFTK.save_scfres_master(filename::AbstractString, scfres::NamedTuple, :
 
     # Storing the bloch waves
     if save_ψ
-        for ik in 1:length(basis.kpoints)
-            for iband in 1:size(scfres.ψ[ik])[2]
+        for ik = 1:length(basis.kpoints)
+            for iband = 1:size(scfres.ψ[ik])[2]
                 ψnk_real = ifft(basis, basis.kpoints[ik], scfres.ψ[ik][:, iband])
                 vtkfile["ψ_k$(ik)_band$(iband)_real"] = real.(ψnk_real)
                 vtkfile["ψ_k$(ik)_band$(iband)_imag"] = imag.(ψnk_real)
diff --git a/src/PlaneWaveBasis.jl b/src/PlaneWaveBasis.jl
index 4cc2bd5dec..682366c04a 100644
--- a/src/PlaneWaveBasis.jl
+++ b/src/PlaneWaveBasis.jl
@@ -124,7 +124,7 @@ Base.eltype(::PlaneWaveBasis{T}) where {T} = T
 @timing function build_kpoints(model::Model{T}, fft_size, kcoords, Ecut;
                                variational=true,
                                architecture::AbstractArchitecture) where {T}
-    kpoints_per_spin = [Kpoint[] for _ in 1:model.n_spin_components]
+    kpoints_per_spin = [Kpoint[] for _ = 1:model.n_spin_components]
     for k in kcoords
         k = Vec3{T}(k)  # rationals are sloooow
         mapping = Int[]
@@ -367,7 +367,7 @@ function G_vectors(fft_size::Union{Tuple,AbstractVector})
     # than this implementation, hence the code duplication.
     start = .- cld.(fft_size .- 1, 2)
     stop  = fld.(fft_size .- 1, 2)
-    axes  = [[collect(0:stop[i]); collect(start[i]:-1)] for i in 1:3]
+    axes  = [[collect(0:stop[i]); collect(start[i]:-1)] for i = 1:3]
     [Vec3{Int}(i, j, k) for i in axes[1], j in axes[2], k in axes[3]]
 end
 
@@ -376,7 +376,7 @@ function G_vectors_generator(fft_size::Union{Tuple,AbstractVector})
     # accumulate_over_symmetries!, which are 100-fold slower with G_vector(fft_size).
     start = .- cld.(fft_size .- 1, 2)
     stop  = fld.(fft_size .- 1, 2)
-    axes = [[collect(0:stop[i]); collect(start[i]:-1)] for i in 1:3]
+    axes = [[collect(0:stop[i]); collect(start[i]:-1)] for i = 1:3]
     (Vec3{Int}(i, j, k) for i in axes[1], j in axes[2], k in axes[3])
 end
 
@@ -542,7 +542,7 @@ function gather_kpts(data::AbstractArray, basis::PlaneWaveBasis)
     if MPI.Comm_rank(basis.comm_kpts) == master
         allk_data = similar(data, n_kpts)
         allk_data[basis.krange_allprocs[1]] = data
-        for rank in 1:mpi_nprocs(basis.comm_kpts) - 1  # Note: MPI ranks are 0-based
+        for rank = 1:mpi_nprocs(basis.comm_kpts) - 1  # Note: MPI ranks are 0-based
             rk_data, status = MPI.recv(basis.comm_kpts, MPI.Status; source=rank, tag=tag)
             @assert MPI.Get_error(status) == 0  # all went well
             allk_data[basis.krange_allprocs[rank + 1]] = rk_data
diff --git a/src/common/quadrature.jl b/src/common/quadrature.jl
index c588a8c33c..52e418acb4 100644
--- a/src/common/quadrature.jl
+++ b/src/common/quadrature.jl
@@ -14,7 +14,7 @@ trapezoidal
     n == length(y) || error("vectors `x` and `y` must have the same number of elements")
     n == 1 && return zero(promote_type(eltype(x), eltype(y)))
     I = (x[2] - x[1]) * y[1]
-    @turbo for i in 2:(n-1)
+    @turbo for i = 2:(n-1)
         # dx[i] + dx[i - 1] = (x[i + 1] - x[i]) + (x[i] - x[i - 1])
         #                   = x[i + 1] - x[i - 1]
         I += (x[i + 1] - x[i - 1]) * y[i]
@@ -46,10 +46,10 @@ end
     istop = isodd(n_intervals) ? n - 1 : n - 2
 
     I = 1 / 3 * dx * y[1]
-    @turbo for i in 2:2:istop
+    @turbo for i = 2:2:istop
         I += 4 / 3 * dx * y[i]
     end
-    @turbo for i in 3:2:istop
+    @turbo for i = 3:2:istop
         I += 2 / 3 * dx * y[i]
     end
 
@@ -70,7 +70,7 @@ end
 
     I = zero(promote_type(eltype(x), eltype(y)))
     # This breaks when @turbo'd
-    @simd for i in 1:2:istop
+    @simd for i = 1:2:istop
         dx0 = x[i + 1] - x[i]
         dx1 = x[i + 2] - x[i + 1]
         c = (dx0 + dx1) / 6
diff --git a/src/common/threading.jl b/src/common/threading.jl
index f52a6b0a8e..b052affc3a 100644
--- a/src/common/threading.jl
+++ b/src/common/threading.jl
@@ -1,7 +1,7 @@
 import FFTW
 using LinearAlgebra
 
-function setup_threading(;n_fft=1, n_blas=Threads.nthreads())
+function setup_threading(; n_fft=1, n_blas=Threads.nthreads())
     n_julia = Threads.nthreads()
     FFTW.set_num_threads(n_fft)
     BLAS.set_num_threads(n_blas)
diff --git a/src/densities.jl b/src/densities.jl
index 7fd187b50d..7c1b48a59b 100644
--- a/src/densities.jl
+++ b/src/densities.jl
@@ -77,7 +77,7 @@ end
     τ .= 0
     dαψnk_real = zeros(complex(eltype(basis)), basis.fft_size)
     for (ik, kpt) in enumerate(basis.kpoints)
-        G_plus_k = [[Gk[α] for Gk in Gplusk_vectors_cart(basis, kpt)] for α in 1:3]
+        G_plus_k = [[Gk[α] for Gk in Gplusk_vectors_cart(basis, kpt)] for α = 1:3]
         for n = 1:size(ψ[ik], 2), α = 1:3
             ifft!(dαψnk_real, basis, kpt, im .* G_plus_k[α] .* ψ[ik][:, n])
             @. τ[:, :, :, kpt.spin] += occupation[ik][n] * basis.kweights[ik] / 2 * abs2(dαψnk_real)
diff --git a/src/eigen/lobpcg_hyper_impl.jl b/src/eigen/lobpcg_hyper_impl.jl
index 8e20cc3f19..b1543358f3 100644
--- a/src/eigen/lobpcg_hyper_impl.jl
+++ b/src/eigen/lobpcg_hyper_impl.jl
@@ -302,7 +302,7 @@ function final_retval(X, AX, resid_history, niter, n_matvec)
     λ = real(diag(X' * AX))
     residuals = AX .- X*Diagonal(λ)
     (; λ, X,
-     residual_norms=[norm(residuals[:, i]) for i in 1:size(residuals, 2)],
+     residual_norms=[norm(residuals[:, i]) for i = 1:size(residuals, 2)],
      residual_history=resid_history[:, 1:niter+1], n_matvec)
 end
 
diff --git a/src/elements.jl b/src/elements.jl
index 839ad131e2..1e2728088f 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -142,17 +142,17 @@ or an element name (e.g. `"silicon"`)
 function ElementCohenBergstresser(key; lattice_constant=nothing)
     # Form factors from Cohen-Bergstresser paper Table 2
     # Lattice constants from Table 1
-    data = Dict(:Si => (form_factors=Dict( 3 => -0.21u"Ry",
-                                           8 =>  0.04u"Ry",
-                                          11 =>  0.08u"Ry"),
+    data = Dict(:Si => (; form_factors=Dict( 3 => -0.21u"Ry",
+                                             8 =>  0.04u"Ry",
+                                            11 =>  0.08u"Ry"),
                         lattice_constant=5.43u"Å"),
-                :Ge => (form_factors=Dict( 3 => -0.23u"Ry",
-                                           8 =>  0.01u"Ry",
-                                          11 =>  0.06u"Ry"),
+                :Ge => (; form_factors=Dict( 3 => -0.23u"Ry",
+                                             8 =>  0.01u"Ry",
+                                            11 =>  0.06u"Ry"),
                         lattice_constant=5.66u"Å"),
-                :Sn => (form_factors=Dict( 3 => -0.20u"Ry",
-                                           8 =>  0.00u"Ry",
-                                          11 =>  0.04u"Ry"),
+                :Sn => (; form_factors=Dict( 3 => -0.20u"Ry",
+                                             8 =>  0.00u"Ry",
+                                            11 =>  0.04u"Ry"),
                         lattice_constant=6.49u"Å"),
             )
 
diff --git a/src/external/wannier90.jl b/src/external/wannier90.jl
index 40defa3020..9ce0472211 100644
--- a/src/external/wannier90.jl
+++ b/src/external/wannier90.jl
@@ -48,7 +48,7 @@ function write_w90_win(fileprefix::String, basis::PlaneWaveBasis;
             path = kpath.paths[1]
 
             println(fp, "begin kpoint_path")
-            for i in 1:length(path)-1
+            for i = 1:length(path)-1
                 A, B = path[i:i+1]  # write segment A -> B
                 @printf(fp, "%s %10.6f %10.6f %10.6f  ", A, round.(kpath.points[A], digits=5)...)
                 @printf(fp, "%s %10.6f %10.6f %10.6f\n", B, round.(kpath.points[B], digits=5)...)
@@ -97,7 +97,7 @@ function read_w90_nnkp(fileprefix::String)
         # 1st: Index of k-point
         # 2nd: Index of periodic image of k+b k-point
         # 3rd: Shift vector to get k-point of ikpb to the actual k+b point required
-        (ik=splitted[1], ikpb=splitted[2], G_shift=splitted[3:5])
+        (; ik=splitted[1], ikpb=splitted[2], G_shift=splitted[3:5])
     end
     (; nntot, nnkpts)
 end
@@ -124,9 +124,9 @@ end
     for ik in krange_spin(basis, spin)
         open(dirname(fileprefix) * (@sprintf "/UNK%05i.%i" ik spin), "w") do fp
             println(fp, "$(fft_size[1]) $(fft_size[2]) $(fft_size[3]) $ik $n_bands")
-            for iband in 1:n_bands
+            for iband = 1:n_bands
                 ψnk_real = ifft(basis, basis.kpoints[ik], @view ψ[ik][:, iband])
-                for iz in 1:fft_size[3], iy in 1:fft_size[2], ix in 1:fft_size[1]
+                for iz = 1:fft_size[3], iy = 1:fft_size[2], ix = 1:fft_size[1]
                     println(fp, real(ψnk_real[ix, iy, iz]), " ", imag(ψnk_real[ix, iy, iz]))
                 end
             end
@@ -159,8 +159,8 @@ We use here that:
     # Compute overlaps
     # TODO This should be improved ...
     Mkb = zeros(ComplexF64, (n_bands, n_bands))
-    for n in 1:n_bands
-        for m in 1:n_bands
+    for n = 1:n_bands
+        for m = 1:n_bands
             # Select the coefficient in right order
             Mkb[m, n] = dot(ψ[ik][iGk, m], ψ[ikpb][iGkpb, n])
         end
@@ -208,13 +208,13 @@ function compute_Ak_gaussian_guess(basis::PlaneWaveBasis, ψk, kpt, centers, n_b
     Ak = zeros(eltype(ψk), (n_bands, n_wannier))
 
     # Compute Ak
-    for n in 1:n_wannier
+    for n = 1:n_wannier
         # Functions are l^2 normalized in Fourier, in DFTK conventions.
         norm_gn_per = norm(fourier_gn(centers[n], qs), 2)
         # Fourier coeffs of gn_per for k+G in common with ψk
         coeffs_gn_per = fourier_gn(centers[n], qs[kpt.mapping]) ./ norm_gn_per
         # Compute overlap
-        for m in 1:n_bands
+        for m = 1:n_bands
             # TODO Check the ordering of m and n here!
             Ak[m, n] = dot(ψk[:, m], coeffs_gn_per)
         end
@@ -230,8 +230,8 @@ end
 
         for (ik, (ψk, kpt)) in enumerate(zip(ψ, basis.kpoints))
             Ak = compute_Ak_gaussian_guess(basis, ψk, kpt, centers, n_bands)
-            for n in 1:size(Ak, 2)
-                for m in 1:size(Ak, 1)
+            for n = 1:size(Ak, 2)
+                for m = 1:size(Ak, 1)
                     @printf(fp, "%3i %3i %3i  %22.18f %22.18f \n",
                             m, n, ik, real(Ak[m, n]), imag(Ak[m, n]))
                 end
@@ -245,7 +245,7 @@ end
 Default random Gaussian guess for maximally-localised wannier functions
 generated in reduced coordinates.
 """
-default_wannier_centres(n_wannier) = [rand(1, 3) for _ in 1:n_wannier]
+default_wannier_centres(n_wannier) = [rand(1, 3) for _ = 1:n_wannier]
 
 @timing function run_wannier90(scfres;
                                n_bands=scfres.n_bands_converge,
diff --git a/src/interpolation.jl b/src/interpolation.jl
index b983d26fea..1a219dc69a 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -96,7 +96,7 @@ function interpolate_kpoint(data_in::AbstractVecOrMat,
     n_Gk_out = length(G_vectors(basis_out, kpoint_out))
     data_out = similar(data_in, n_Gk_out, n_bands) .= 0
     # TODO: use a map, or this will not be GPU compatible (scalar indexing)
-    for iin in 1:size(data_in, 1)
+    for iin = 1:size(data_in, 1)
         idx_fft = kpoint_in.mapping[iin]
         idx_fft in keys(kpoint_out.mapping_inv) || continue
         iout = kpoint_out.mapping_inv[idx_fft]
diff --git a/src/plotting.jl b/src/plotting.jl
index a58b95be92..2d82efed09 100644
--- a/src/plotting.jl
+++ b/src/plotting.jl
@@ -1,14 +1,14 @@
 # This is needed to flag that the plots-dependent code has been loaded
 const PLOTS_LOADED = true
 
-function ScfPlotTrace(plt=Plots.plot(yaxis=:log); kwargs...)
+function ScfPlotTrace(plt=Plots.plot(; yaxis=:log); kwargs...)
     energies = nothing
     function callback(info)
         if info.stage == :finalize
             minenergy = minimum(energies[max(1, end-5):end])
             error = abs.(energies .- minenergy)
             error[error .== 0] .= NaN
-            extra = ifelse(:mark in keys(kwargs), (), (mark=:x, ))
+            extra = ifelse(:mark in keys(kwargs), (), (; mark=:x))
             Plots.plot!(plt, error; extra..., kwargs...)
             display(plt)
         elseif info.n_iter == 1
@@ -42,9 +42,9 @@ function plot_band_data(kpath::KPathInterpolant, band_data;
     to_unit = ustrip(auconvert(unit, 1.0))
 
     # Plot all bands, spins and errors
-    p = Plots.plot(xlabel="wave vector")
+    p = Plots.plot(; xlabel="wave vector")
     margs = length(kpath) < 70 ? (; markersize=2, markershape=:circle) : (; )
-    for σ in 1:data.n_spin, iband = 1:data.n_bands, branch in data.kbranches
+    for σ = 1:data.n_spin, iband = 1:data.n_bands, branch in data.kbranches
         yerror = nothing
         if hasproperty(data, :λerror)
             yerror = data.λerror[:, iband, σ][branch] .* to_unit
@@ -90,7 +90,7 @@ function plot_dos(basis, eigenvalues; εF=nothing, unit=u"hartree",
     spinlabels = spin_components(basis.model)
     colors = [:blue, :red]
     Dεs = compute_dos.(εs, Ref(basis), Ref(eigenvalues))
-    for σ in 1:n_spin
+    for σ = 1:n_spin
         D = [Dσ[σ] for Dσ in Dεs]
         label = n_spin > 1 ? "DOS $(spinlabels[σ]) spin" : "DOS"
         Plots.plot!(p, (εs .- eshift) .* to_unit, D; label, color=colors[σ])
diff --git a/src/postprocess/band_structure.jl b/src/postprocess/band_structure.jl
index 550a7af272..464f85870f 100644
--- a/src/postprocess/band_structure.jl
+++ b/src/postprocess/band_structure.jl
@@ -133,7 +133,7 @@ function kdistances_and_ticks(kcoords, klabels::Dict, kbranches)
             end
         end
     end
-    ticks = (distances=tick_distances, labels=tick_labels)
+    ticks = (; distances=tick_distances, labels=tick_labels)
     (; kdistances, ticks)
 end
 
@@ -167,7 +167,7 @@ function data_for_plotting(kpath::KPathInterpolant, band_data; datakeys=[:λ, :
         hasproperty(band_data, key) || continue
         n_bands = length(band_data[key][1])
         data_per_kσ = similar(band_data[key][1], n_kcoord, n_bands, n_spin)
-        for σ in 1:n_spin, (ito, ik) in enumerate(krange_spin(basis, σ))
+        for σ = 1:n_spin, (ito, ik) in enumerate(krange_spin(basis, σ))
             data_per_kσ[ito, :, σ] = band_data[key][ik]
         end
         data[key] = data_per_kσ
@@ -187,7 +187,7 @@ i.e. where bands are cut by the Fermi level.
 """
 function is_metal(eigenvalues, εF; tol=1e-4)
     n_bands = length(eigenvalues[1])
-    for ib in 1:n_bands
+    for ib = 1:n_bands
         some_larger  = any(εk[ib] > εF + tol for εk in eigenvalues)
         some_smaller = any(εk[ib] < εF - tol for εk in eigenvalues)
         some_larger && some_smaller && return true
diff --git a/src/postprocess/dos.jl b/src/postprocess/dos.jl
index 05b6780116..141dd1379d 100644
--- a/src/postprocess/dos.jl
+++ b/src/postprocess/dos.jl
@@ -19,7 +19,7 @@ function compute_dos(ε, basis, eigenvalues; smearing=basis.model.smearing,
     filled_occ = filled_occupation(basis.model)
 
     D = zeros(typeof(ε), basis.model.n_spin_components)
-    for σ in 1:basis.model.n_spin_components, ik = krange_spin(basis, σ)
+    for σ = 1:basis.model.n_spin_components, ik = krange_spin(basis, σ)
         for (iband, εnk) in enumerate(eigenvalues[ik])
             enred = (εnk - ε) / temperature
             D[σ] -= (filled_occ * basis.kweights[ik] / temperature
diff --git a/src/postprocess/phonon.jl b/src/postprocess/phonon.jl
index b4671dca3e..07c592ebbf 100644
--- a/src/postprocess/phonon.jl
+++ b/src/postprocess/phonon.jl
@@ -13,7 +13,7 @@ function dynmat_red_to_cart(model::Model, dynamical_matrix)
     #   ⇒ D_cart = lattice⁻ᵀ · D_red · lattice⁻¹.
 
     cart_mat = zero.(dynamical_matrix)
-    for τ in 1:size(cart_mat, 2), η in 1:size(cart_mat, 4)
+    for τ = 1:size(cart_mat, 2), η = 1:size(cart_mat, 4)
         cart_mat[:, η, :, τ] = inv_lattice' * dynamical_matrix[:, η, :, τ] * inv_lattice
     end
     cart_mat
diff --git a/src/pseudo/NormConservingPsp.jl b/src/pseudo/NormConservingPsp.jl
index c41d38fb18..c7b1dd628e 100644
--- a/src/pseudo/NormConservingPsp.jl
+++ b/src/pseudo/NormConservingPsp.jl
@@ -149,7 +149,7 @@ import Base.Broadcast.broadcastable
 Base.Broadcast.broadcastable(psp::NormConservingPsp) = Ref(psp)
 
 function projector_indices(psp::NormConservingPsp)
-    ((i, l, m) for l in 0:psp.lmax for i in 1:size(psp.h[l+1], 1) for m = -l:l)
+    ((i, l, m) for l = 0:psp.lmax for i = 1:size(psp.h[l+1], 1) for m = -l:l)
 end
 
 """
diff --git a/src/pseudo/PspHgh.jl b/src/pseudo/PspHgh.jl
index d76a78cee8..8bf8746234 100644
--- a/src/pseudo/PspHgh.jl
+++ b/src/pseudo/PspHgh.jl
@@ -48,7 +48,7 @@ function PspHgh(path; identifier=path, kwargs...)
     h = Vector{Matrix{Float64}}(undef, lmax + 1)
     cur = 5  # Current line to parse
 
-    for l in 0:lmax
+    for l = 0:lmax
         # loop over all AM channels and extract projectors,
         # these are given in blocks like
         #
@@ -77,7 +77,7 @@ function PspHgh(path; identifier=path, kwargs...)
         # Else we have to parse the extra parts of the hcoeff matrix.
         # This is done here.
         hcoeff = [parse(Float64, part) for part in split(m[3])]
-        for i in 1:nproj
+        for i = 1:nproj
             for j in i:nproj
                 h[l + 1][j, i] = h[l + 1][i, j] = hcoeff[j - i + 1]
             end
diff --git a/src/pseudo/PspUpf.jl b/src/pseudo/PspUpf.jl
index 21b08b5dfb..4d005ab165 100644
--- a/src/pseudo/PspUpf.jl
+++ b/src/pseudo/PspUpf.jl
@@ -119,11 +119,11 @@ function PspUpf(path; identifier=path, rcut=nothing)
         count += nproj_l
     end
 
-    r2_pswfcs = [Vector{Float64}[] for _ in 0:lmax-1]
-    pswfc_occs = [Float64[] for _ in 0:lmax-1]
-    pswfc_energies = [Float64[] for _ in 0:lmax-1]
-    pswfc_labels = [String[] for _ in 0:lmax-1]
-    for l in 0:lmax-1
+    r2_pswfcs = [Vector{Float64}[] for _ = 0:lmax-1]
+    pswfc_occs = [Float64[] for _ = 0:lmax-1]
+    pswfc_energies = [Float64[] for _ = 0:lmax-1]
+    pswfc_labels = [String[] for _ = 0:lmax-1]
+    for l = 0:lmax-1
         pswfcs_l = filter(χ -> χ["angular_momentum"] == l, pseudo["atomic_wave_functions"])
         for pswfc_li in pswfcs_l
             # Ry -> Ha, rχ -> r²χ
diff --git a/src/pseudo/list_psp.jl b/src/pseudo/list_psp.jl
index 8b1553d2d6..ede35cf724 100644
--- a/src/pseudo/list_psp.jl
+++ b/src/pseudo/list_psp.jl
@@ -8,11 +8,11 @@ to restrict the displayed files.
 
 # Examples
 ```julia-repl
-julia> list_psp(family="hgh")
+julia> list_psp(; family="hgh")
 ```
 will list all HGH-type pseudopotentials and
 ```julia-repl
-julia> list_psp(family="hgh", functional="lda")
+julia> list_psp(; family="hgh", functional="lda")
 ```
 will only list those for LDA (also known as Pade in this context)
 and
@@ -47,7 +47,7 @@ function list_psp(element=nothing; family=nothing, functional=nothing, core=noth
 
             f_identifier = joinpath(root, file)
             Sys.iswindows() && (f_identifier = replace(f_identifier, "\\" => "/"))
-            push!(psplist, (identifier=f_identifier, family=f_family,
+            push!(psplist, (; identifier=f_identifier, family=f_family,
                   functional=f_functional, element=f_element,
                   n_elec_valence=parse(Int, f_nvalence[2:end]),
                   path=joinpath(datadir_psp(), root, file)))
@@ -59,10 +59,10 @@ function list_psp(element=nothing; family=nothing, functional=nothing, core=noth
     psp_per_element = map(per_elem) do elgroup
         @assert length(elgroup) > 0
         if length(elgroup) == 1
-            cores = [(core=:fullcore, )]
+            cores = [(; core=:fullcore)]
         else
-            cores = append!(fill((core=:other, ), length(elgroup) - 2),
-                            [(core=:fullcore, ), (core=:semicore, )])
+            cores = append!(fill((; core=:other), length(elgroup) - 2),
+                            [(; core=:fullcore), (; core=:semicore)])
         end
         merge.(sort(elgroup, by=psp -> psp.n_elec_valence), cores)
     end
diff --git a/src/response/chi0.jl b/src/response/chi0.jl
index 32f11440db..6f3610ac43 100644
--- a/src/response/chi0.jl
+++ b/src/response/chi0.jl
@@ -391,8 +391,8 @@ function apply_χ0(ham, ψ, occupation, εF, eigenvalues, δV::AbstractArray{T};
 
     δHψ = [RealSpaceMultiplication(basis, kpt, @views δV[:, :, :, kpt.spin]) * ψ[ik]
            for (ik, kpt) in enumerate(basis.kpoints)]
-    δψ, δoccupation, δεF = apply_χ0_4P(ham, ψ, occupation, εF, eigenvalues, δHψ;
-                                       occupation_threshold, kwargs_sternheimer...)
+    (; δψ, δoccupation) = apply_χ0_4P(ham, ψ, occupation, εF, eigenvalues, δHψ;
+                                      occupation_threshold, kwargs_sternheimer...)
     δρ = compute_δρ(basis, ψ, δψ, occupation, δoccupation; occupation_threshold)
     δρ * normδV
 end
diff --git a/src/scf/direct_minimization.jl b/src/scf/direct_minimization.jl
index 76b55e2fb2..99c9e34454 100644
--- a/src/scf/direct_minimization.jl
+++ b/src/scf/direct_minimization.jl
@@ -17,9 +17,7 @@ function Optim.project_tangent!(m::DMManifold, g, x)
     g_unpack = m.unpack(g)
     x_unpack = m.unpack(x)
     for ik = 1:m.Nk
-        Optim.project_tangent!(Optim.Stiefel(),
-                               g_unpack[ik],
-                               x_unpack[ik])
+        Optim.project_tangent!(Optim.Stiefel(), g_unpack[ik], x_unpack[ik])
     end
     g
 end
@@ -84,7 +82,7 @@ function direct_minimization(basis::PlaneWaveBasis{T}, ψ0;
     if ψ0 === nothing
         ψ0 = [random_orbitals(basis, kpt, n_bands) for kpt in basis.kpoints]
     end
-    occupation = [filled_occ * ones(T, n_bands) for ik = 1:Nk]
+    occupation = [filled_occ * ones(T, n_bands) for _ = 1:Nk]
 
     # we need to copy the reinterpret array here to not raise errors in Optim.jl
     # TODO raise this issue in Optim.jl
diff --git a/src/scf/mixing.jl b/src/scf/mixing.jl
index 845065dd55..b2b9f08f8b 100644
--- a/src/scf/mixing.jl
+++ b/src/scf/mixing.jl
@@ -163,7 +163,7 @@ Important `kwargs` passed on to [`χ0Mixing`](@ref)
 - `verbose`: Run the GMRES in verbose mode.
 - `reltol`: Relative tolerance for GMRES
 """
-function HybridMixing(;εr=1.0, kTF=0.8, localization=identity,
+function HybridMixing(; εr=1.0, kTF=0.8, localization=identity,
                       adjust_temperature=IncreaseMixingTemperature(), kwargs...)
     χ0terms = [DielectricModel(; εr, kTF, localization),
                LdosModel(;adjust_temperature)]
diff --git a/src/scf/newton.jl b/src/scf/newton.jl
index e3835bc95a..f8b1f1e713 100644
--- a/src/scf/newton.jl
+++ b/src/scf/newton.jl
@@ -97,7 +97,7 @@ function newton(basis::PlaneWaveBasis{T}, ψ0;
 
     # number of kpoints and occupation
     Nk = length(basis.kpoints)
-    occupation = [filled_occ * ones(T, n_bands) for ik = 1:Nk]
+    occupation = [filled_occ * ones(T, n_bands) for _ = 1:Nk]
 
     # iterators
     n_iter = 0
@@ -116,7 +116,7 @@ function newton(basis::PlaneWaveBasis{T}, ψ0;
         res = compute_projected_gradient(basis, ψ, occupation)
         # solve (Ω+K) δψ = -res so that the Newton step is ψ <- ψ + δψ
         δψ = solve_ΩplusK(basis, ψ, -res, occupation; tol=tol_cg).δψ
-        ψ  = [ortho_qr(ψ[ik] + δψ[ik]) for ik in 1:Nk]
+        ψ  = [ortho_qr(ψ[ik] + δψ[ik]) for ik = 1:Nk]
 
         ρ_next = compute_density(basis, ψ, occupation)
         energies, H = energy_hamiltonian(basis, ψ, occupation; ρ=ρ_next)
diff --git a/src/scf/scf_solvers.jl b/src/scf/scf_solvers.jl
index 534df8a517..04b6b1aa5d 100644
--- a/src/scf/scf_solvers.jl
+++ b/src/scf/scf_solvers.jl
@@ -15,7 +15,7 @@ function scf_damping_solver(β=0.2)
         β = convert(eltype(x0), β)
         converged = false
         x = copy(x0)
-        for i in 1:max_iter
+        for i = 1:max_iter
             x_new = f(x)
 
             if norm(x_new - x) < tol
@@ -71,7 +71,7 @@ function CROP(f, x0, m::Int, max_iter::Int, tol::Real, warming=0)
     # Cheat support for multidimensional arrays
     if length(size(x0)) != 1
         x, conv= CROP(x -> vec(f(reshape(x, size(x0)...))), vec(x0), m, max_iter, tol, warming)
-        return (fixpoint=reshape(x, size(x0)...), converged=conv)
+        return (; fixpoint=reshape(x, size(x0)...), converged=conv)
     end
     N = size(x0,1)
     T = eltype(x0)
@@ -110,6 +110,6 @@ function CROP(f, x0, m::Int, max_iter::Int, tol::Real, warming=0)
         fs[:,1] = ftnp1
         # fs[:,1] = f(xs[:,1])
     end
-    (fixpoint=xs[:, 1], converged=err < tol)
+    (; fixpoint=xs[:, 1], converged=err < tol)
 end
 scf_CROP_solver(m=10) = (f, x0, max_iter; tol=1e-6) -> CROP(x -> f(x) - x, x0, m, max_iter, tol)
diff --git a/src/scf/self_consistent_field.jl b/src/scf/self_consistent_field.jl
index 17d52447da..fcd1990ee3 100644
--- a/src/scf/self_consistent_field.jl
+++ b/src/scf/self_consistent_field.jl
@@ -128,7 +128,7 @@ Overview of parameters:
         # Diagonalize `ham` to get the new state
         nextstate = next_density(ham, nbandsalg, fermialg; eigensolver, ψ, eigenvalues,
                                  occupation, miniter=1, tol=determine_diagtol(info))
-        ψ, eigenvalues, occupation, εF, ρout = nextstate
+        (; ψ, eigenvalues, occupation, εF, ρout) = nextstate
 
         # Update info with results gathered so far
         info = (; ham, basis, converged, stage=:iterate, algorithm="SCF",
diff --git a/src/structure.jl b/src/structure.jl
index 45afe9b2f6..949625506e 100644
--- a/src/structure.jl
+++ b/src/structure.jl
@@ -51,7 +51,7 @@ function estimate_integer_lattice_bounds(M::AbstractMatrix{T}, δ, shift=zeros(3
     # As a general statement, with M a lattice matrix, then if ||Mx|| <= δ,
     # then xi = <ei, M^-1 Mx> = <M^-T ei, Mx> <= ||M^-T ei|| δ.
     inv_lattice_t = compute_inverse_lattice(M')
-    xlims = [norm(inv_lattice_t[:, i]) * δ + shift[i] for i in 1:3]
+    xlims = [norm(inv_lattice_t[:, i]) * δ + shift[i] for i = 1:3]
 
     # Round up, unless exactly zero (in which case keep it zero in
     # order to just have one x vector for 1D or 2D systems)
diff --git a/src/supercell.jl b/src/supercell.jl
index c10dabd215..8bc33b6efb 100644
--- a/src/supercell.jl
+++ b/src/supercell.jl
@@ -12,8 +12,8 @@ function create_supercell(lattice, atoms, positions, supercell_size)
 
     for (atom, position) in zip(atoms, positions)
         append!(positions_supercell, [(position .+ [i;j;k]) ./ [nx, ny, nz]
-                                      for i in 0:nx-1, j in 0:ny-1, k in 0:nz-1])
-        append!(atoms_supercell, vcat([atom for _ in 1:nx*ny*nz]...))
+                                      for i = 0:nx-1, j = 0:ny-1, k = 0:nz-1])
+        append!(atoms_supercell, vcat([atom for _ = 1:nx*ny*nz]...))
     end
 
     (; lattice=lattice_supercell, atoms=atoms_supercell, positions=positions_supercell)
diff --git a/src/symmetry.jl b/src/symmetry.jl
index 1c79a24757..24af8250ef 100644
--- a/src/symmetry.jl
+++ b/src/symmetry.jl
@@ -257,7 +257,7 @@ function apply_symop(symop::SymOp, basis, kpoint, ψk::AbstractVecOrMat)
     invS = Mat3{Int}(inv(S))
     Gs_full = [G + kshift for G in G_vectors(basis, Skpoint)]
     ψSk = zero(ψk)
-    for iband in 1:size(ψk, 2)
+    for iband = 1:size(ψk, 2)
         for (ig, G_full) in enumerate(Gs_full)
             igired = index_G_vectors(basis, kpoint, invS * G_full)
             @assert igired !== nothing
@@ -420,7 +420,7 @@ function unfold_array_(basis_irred, basis_unfolded, data, is_ψ)
         error("Brillouin zone symmetry unfolding not supported with MPI yet")
     end
     data_unfolded = similar(data, length(basis_unfolded.kpoints))
-    for ik_unfolded in 1:length(basis_unfolded.kpoints)
+    for ik_unfolded = 1:length(basis_unfolded.kpoints)
         kpt_unfolded = basis_unfolded.kpoints[ik_unfolded]
         ik_irred, symop = unfold_mapping(basis_irred, kpt_unfolded)
         if is_ψ
diff --git a/src/terms/Hamiltonian.jl b/src/terms/Hamiltonian.jl
index 23e2ecbf54..37f8247c6b 100644
--- a/src/terms/Hamiltonian.jl
+++ b/src/terms/Hamiltonian.jl
@@ -107,9 +107,9 @@ Base.:*(H::Hamiltonian, ψ) = mul!(deepcopy(ψ), H, ψ)
         ifft!(ψ_real, H.basis, H.kpoint, ψ[:, iband])
         for op in H.optimized_operators
             @timing "$(nameof(typeof(op)))" begin
-                apply!((fourier=Hψ_fourier, real=Hψ_real),
+                apply!((; fourier=Hψ_fourier, real=Hψ_real),
                        op,
-                       (fourier=ψ[:, iband], real=ψ_real))
+                       (; fourier=ψ[:, iband], real=ψ_real))
             end
         end
         Hψ[:, iband] .= Hψ_fourier
@@ -147,9 +147,9 @@ end
 
             if have_divAgrad
                 @timeit to "divAgrad" begin
-                    apply!((fourier=Hψ[:, iband], real=nothing),
+                    apply!((; fourier=Hψ[:, iband], real=nothing),
                            H.divAgrad_op,
-                           (fourier=ψ[:, iband], real=nothing),
+                           (; fourier=ψ[:, iband], real=nothing),
                            ψ_real)  # ψ_real used as scratch
                 end
             end
@@ -165,7 +165,9 @@ end
     # Apply the nonlocal operator
     if !isnothing(H.nonlocal_op)
         @timing "nonlocal" begin
-            apply!((fourier=Hψ, real=nothing), H.nonlocal_op, (fourier=ψ, real=nothing))
+            apply!((; fourier=Hψ, real=nothing),
+                   H.nonlocal_op,
+                   (; fourier=ψ, real=nothing))
         end
     end
 
diff --git a/src/terms/ewald.jl b/src/terms/ewald.jl
index 73aebbd824..3d23d166e8 100644
--- a/src/terms/ewald.jl
+++ b/src/terms/ewald.jl
@@ -112,7 +112,7 @@ function energy_forces_ewald(S, lattice::AbstractArray{T}, charges, positions, q
     # In the real-space term we have erfc(η ||A(rj - rk - R)||),
     # where A is the real-space lattice, rj and rk are atomic positions and
     # thus use the bound  ||A(rj - rk - R)|| * η ≤ max_erfc_arg
-    poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i in 1:3]
+    poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i = 1:3]
     Rlims = estimate_integer_lattice_bounds(lattice, max_erfc_arg / η, poslims)
 
     #
@@ -238,8 +238,8 @@ function compute_dynmat(ewald::TermEwald, basis::PlaneWaveBasis{T}, ψ, occupati
 
     dynmat = zeros(complex(T), n_dim, n_atoms, n_dim, n_atoms)
     # Real part
-    for τ in 1:n_atoms
-        for γ in 1:n_dim
+    for τ = 1:n_atoms
+        for γ = 1:n_dim
             displacement = zero.(model.positions)
             displacement[τ] = setindex(displacement[τ], one(T), γ)
             real_part = -ForwardDiff.derivative(zero(T)) do ε
@@ -253,8 +253,8 @@ function compute_dynmat(ewald::TermEwald, basis::PlaneWaveBasis{T}, ψ, occupati
         end
     end
     # Reciprocal part
-    for τ in 1:n_atoms
-        for σ in 1:n_atoms
+    for τ = 1:n_atoms
+        for σ = 1:n_atoms
             dynmat[:, σ, :, τ] += dynmat_ewald_recip(model, σ, τ; ewald.η, q)
         end
     end
diff --git a/src/terms/local.jl b/src/terms/local.jl
index d37df2c589..3ee08cb40b 100644
--- a/src/terms/local.jl
+++ b/src/terms/local.jl
@@ -131,7 +131,7 @@ end
 
     # energy = sum of form_factor(G) * struct_factor(G) * rho(G)
     # where struct_factor(G) = e^{-i G·r}
-    forces = [zero(Vec3{T}) for _ in 1:length(model.positions)]
+    forces = [zero(Vec3{T}) for _ = 1:length(model.positions)]
     for group in model.atom_groups
         element = model.atoms[first(group)]
         form_factors = [Complex{T}(local_potential_fourier(element, norm(G)))
diff --git a/src/terms/nonlocal.jl b/src/terms/nonlocal.jl
index 3bfb4d0cdc..e7dfedf985 100644
--- a/src/terms/nonlocal.jl
+++ b/src/terms/nonlocal.jl
@@ -56,7 +56,7 @@ end
     isempty(psp_groups) && return nothing
 
     # energy terms are of the form <psi, P C P' psi>, where P(G) = form_factor(G) * structure_factor(G)
-    forces = [zero(Vec3{T}) for _ in 1:length(model.positions)]
+    forces = [zero(Vec3{T}) for _ = 1:length(model.positions)]
     for group in psp_groups
         element = model.atoms[first(group)]
 
@@ -119,7 +119,7 @@ function build_projection_coefficients_(T, psp::NormConservingPsp)
     n_proj = count_n_proj(psp)
     proj_coeffs = zeros(T, n_proj, n_proj)
     count = 0
-    for l in 0:psp.lmax, m in -l:l
+    for l = 0:psp.lmax, m in -l:l
         n_proj_l = count_n_proj_radial(psp, l)  # Number of i's
         range = count .+ (1:n_proj_l)
         proj_coeffs[range, range] = psp.h[l + 1]
@@ -205,7 +205,7 @@ function build_form_factors(psp, qs::Array)
         q_norm = norm(q)
         if !haskey(radials, q_norm)
             radials_q = Matrix{T}(undef, n_proj_max, psp.lmax + 1)
-            for l in 0:psp.lmax, iproj_l in 1:count_n_proj_radial(psp, l)
+            for l = 0:psp.lmax, iproj_l = 1:count_n_proj_radial(psp, l)
                 radials_q[iproj_l, l+1] = eval_psp_projector_fourier(psp, iproj_l, l, q_norm)
             end
             radials[q_norm] = radials_q
@@ -216,10 +216,10 @@ function build_form_factors(psp, qs::Array)
     for (iq, q) in enumerate(qs)
         radials_q = radials[norm(q)]
         count = 1
-        for l in 0:psp.lmax, m in -l:l
+        for l = 0:psp.lmax, m in -l:l
             # see "Fourier transforms of centered functions" in the docs for the formula
             angular = (-im)^l * ylm_real(l, m, q)
-            for iproj_l in 1:count_n_proj_radial(psp, l)
+            for iproj_l = 1:count_n_proj_radial(psp, l)
                 form_factors[iq, count] = radials_q[iproj_l, l+1] * angular
                 count += 1
             end
diff --git a/src/terms/operators.jl b/src/terms/operators.jl
index f552ed333b..1602e0db26 100644
--- a/src/terms/operators.jl
+++ b/src/terms/operators.jl
@@ -18,7 +18,9 @@ function LinearAlgebra.mul!(Hψ::AbstractVector, op::RealFourierOperator, ψ::Ab
     Hψ_real = similar(ψ_real)
     Hψ_fourier .= 0
     Hψ_real .= 0
-    apply!((real=Hψ_real, fourier=Hψ_fourier), op, (real=ψ_real, fourier=ψ))
+    apply!((; real=Hψ_real, fourier=Hψ_fourier),
+           op,
+           (; real=ψ_real, fourier=ψ))
     Hψ .= Hψ_fourier .+ fft(op.basis, op.kpoint, Hψ_real)
     Hψ
 end
@@ -140,7 +142,7 @@ function apply!(Hψ, op::DivAgradOperator, ψ,
                 ψ_scratch=zeros(complex(eltype(op.basis)), op.basis.fft_size...))
     # TODO: Performance improvements: Unscaled plans, avoid remaining allocations
     #       (which are only on the small k-point-specific Fourier grid
-    G_plus_k = [[Gk[α] for Gk in Gplusk_vectors_cart(op.basis, op.kpoint)] for α in 1:3]
+    G_plus_k = [[Gk[α] for Gk in Gplusk_vectors_cart(op.basis, op.kpoint)] for α = 1:3]
     for α = 1:3
         ∂αψ_real = ifft!(ψ_scratch, op.basis, op.kpoint, im .* G_plus_k[α] .* ψ.fourier)
         A∇ψ      = fft(op.basis, op.kpoint, ∂αψ_real .* op.A)
diff --git a/src/terms/pairwise.jl b/src/terms/pairwise.jl
index fd34e04725..55f597459c 100644
--- a/src/terms/pairwise.jl
+++ b/src/terms/pairwise.jl
@@ -34,7 +34,7 @@ struct TermPairwisePotential{TV, Tparams, T} <:Term
 end
 
 function ene_ops(term::TermPairwisePotential, basis::PlaneWaveBasis, ψ, occupation; kwargs...)
-    (E=term.energy, ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
+    (; E=term.energy, ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
 end
 compute_forces(term::TermPairwisePotential, ::PlaneWaveBasis, ψ, occ; kwargs...) = term.forces
 
@@ -94,7 +94,7 @@ function energy_forces_pairwise(S, lattice::AbstractArray{T}, symbols, positions
     # The potential V(dist) decays very quickly with dist = ||A (rj - rk - R)||,
     # so we cut off at some point. We use the bound  ||A (rj - rk - R)|| ≤ max_radius
     # where A is the real-space lattice, rj and rk are atomic positions.
-    poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i in 1:3]
+    poslims = [maximum(rj[i] - rk[i] for rj in positions for rk in positions) for i = 1:3]
     Rlims = estimate_integer_lattice_bounds(lattice, max_radius, poslims)
 
     # Check if some coordinates are not used.
@@ -148,8 +148,8 @@ function compute_dynmat(term::TermPairwisePotential, basis::PlaneWaveBasis{T}, 
     symbols = Symbol.(atomic_symbol.(model.atoms))
 
     dynmat = zeros(complex(T), 3, n_atoms, 3, n_atoms)
-    for τ in 1:n_atoms
-        for γ in 1:n_dim
+    for τ = 1:n_atoms
+        for γ = 1:n_dim
             displacement = zero.(model.positions)
             displacement[τ] = setindex(displacement[τ], one(T), γ)
             dynmat[:, :, γ, τ] = -ForwardDiff.derivative(zero(T)) do ε
diff --git a/src/terms/psp_correction.jl b/src/terms/psp_correction.jl
index 2c0e30903b..b3d029f21a 100644
--- a/src/terms/psp_correction.jl
+++ b/src/terms/psp_correction.jl
@@ -16,7 +16,7 @@ function TermPspCorrection(basis::PlaneWaveBasis)
 end
 
 function ene_ops(term::TermPspCorrection, basis::PlaneWaveBasis, ψ, occupation; kwargs...)
-    (E=term.energy, ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
+    (; E=term.energy, ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
 end
 
 """
diff --git a/src/terms/terms.jl b/src/terms/terms.jl
index 70eb920002..b4bf7e4c52 100644
--- a/src/terms/terms.jl
+++ b/src/terms/terms.jl
@@ -1,11 +1,11 @@
 include("operators.jl")
 
 ### Terms
-# - A Term is something that, given a state, returns a named tuple (E, hams) with an energy
+# - A Term is something that, given a state, returns a named tuple (; E, hams) with an energy
 #   and a list of RealFourierOperator (for each kpoint).
 # - Each term must overload
 #     `ene_ops(term, basis, ψ, occupation; kwargs...)`
-#         -> (E::Real, ops::Vector{RealFourierOperator}).
+#         -> (; E::Real, ops::Vector{RealFourierOperator}).
 # - Note that terms are allowed to hold on to references to ψ (eg Fock term),
 #   so ψ should not mutated after ene_ops
 
@@ -27,7 +27,7 @@ A term with a constant zero energy.
 """
 struct TermNoop <: Term end
 function ene_ops(term::TermNoop, basis::PlaneWaveBasis{T}, ψ, occupation; kwargs...) where {T}
-    (E=zero(eltype(T)), ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
+    (; E=zero(eltype(T)), ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
 end
 
 include("Hamiltonian.jl")
diff --git a/src/terms/xc.jl b/src/terms/xc.jl
index 5b3a726a38..537149fb85 100644
--- a/src/terms/xc.jl
+++ b/src/terms/xc.jl
@@ -95,7 +95,7 @@ function xc_potential_real(term::TermXc, basis::PlaneWaveBasis{T}, ψ, occupatio
 
     # Potential contributions Vρ -2 ∇⋅(Vσ ∇ρ) + ΔVl
     potential = zero(ρ)
-    @views for s in 1:n_spin
+    @views for s = 1:n_spin
         Vρ = to_device(basis.architecture, reshape(terms.Vρ, n_spin, basis.fft_size...))
 
         potential[:, :, :, s] .+= Vρ[s, :, :, :]
@@ -109,7 +109,7 @@ function xc_potential_real(term::TermXc, basis::PlaneWaveBasis{T}, ψ, occupatio
                 # in the energy expression. See comment block below on spin-polarised XC.
                 sum((s == t ? one(T) : one(T)/2)
                     .* Vσ[tσ(s, t), :, :, :] .* density.∇ρ_real[t, :, :, :, α]
-                    for t in 1:n_spin)
+                    for t = 1:n_spin)
             end
         end
         if haskey(terms, :Vl) && any(x -> abs(x) > potential_threshold, terms.Vl)
@@ -158,7 +158,7 @@ end
     # early return if nlcc is disabled / no elements have model core charges
     isnothing(term.ρcore) && return nothing
 
-    _, Vxc_real, _ = xc_potential_real(term, basis, ψ, occupation; ρ, τ)
+    Vxc_real = xc_potential_real(term, basis, ψ, occupation; ρ, τ).potential
     # TODO: the factor of 2 here should be associated with the density, not the potential
     if basis.model.spin_polarization in (:none, :spinless)
         Vxc_fourier = fft(basis, Vxc_real[:,:,:,1])
@@ -172,7 +172,7 @@ end
                    if has_core_density(model.atoms[first(group)])]
     @assert !isnothing(nlcc_groups)
 
-    forces = [zero(Vec3{T}) for _ in 1:length(model.positions)]
+    forces = [zero(Vec3{T}) for _ = 1:length(model.positions)]
     for (igroup, group) in nlcc_groups
         for iatom in group
             r = model.positions[iatom]
@@ -309,7 +309,7 @@ function LibxcDensities(basis, max_derivative::Integer, ρ, τ)
 
         tσ = DftFunctionals.spinindex_σ  # Spin index transformation (s, t) => st as expected by Libxc
         σ_real .= 0
-        @views for α in 1:3
+        @views for α = 1:3
             σ_real[tσ(1, 1), :, :, :] .+= ∇ρ_real[1, :, :, :, α] .* ∇ρ_real[1, :, :, :, α]
             if n_spin > 1
                 σ_real[tσ(1, 2), :, :, :] .+= ∇ρ_real[1, :, :, :, α] .* ∇ρ_real[2, :, :, :, α]
@@ -376,7 +376,7 @@ function apply_kernel(term::TermXc, basis::PlaneWaveBasis{T}, δρ; ρ, kwargs..
     cross_derivatives = Dict{Symbol, Any}()
     if max_ρ_derivs > 0
         cross_derivatives[:δσ] = [
-            @views 2sum(∇ρ[I[1], :, :, :, α] .* ∇δρ[I[2], :, :, :, α] for α in 1:3)
+            @views 2sum(∇ρ[I[1], :, :, :, α] .* ∇δρ[I[2], :, :, :, α] for α = 1:3)
             for I in CartesianIndices((n_spin, n_spin))
         ]
     end
@@ -386,7 +386,7 @@ function apply_kernel(term::TermXc, basis::PlaneWaveBasis{T}, δρ; ρ, kwargs..
     δV = zero(ρ)  # [ix, iy, iz, iσ]
 
     Vρρ = to_device(basis.architecture, reshape(terms.Vρρ, n_spin, n_spin, basis.fft_size...))
-    @views for s in 1:n_spin, t in 1:n_spin  # LDA term
+    @views for s = 1:n_spin, t = 1:n_spin  # LDA term
         δV[:, :, :, s] .+= Vρρ[s, t, :, :, :] .* δρ[t, :, :, :]
     end
     if haskey(terms, :Vρσ)  # GGA term
@@ -433,8 +433,8 @@ function add_kernel_gradient_correction!(δV, terms, density, perturbation, cros
     tσ  = DftFunctionals.spinindex_σ
 
     # Note: δV[ix, iy, iz, iσ] unlike the other quantities ...
-    @views for s in 1:n_spin
-        for t in 1:n_spin, u in 1:n_spin
+    @views for s = 1:n_spin
+        for t = 1:n_spin, u = 1:n_spin
             spinfac_tu = (t == u ? one(T) : one(T)/2)
             @. δV[:, :, :, s] += spinfac_tu * Vρσ[s, tσ(t, u), :, :, :] * δσ[t, u][:, :, :]
         end
@@ -444,16 +444,16 @@ function add_kernel_gradient_correction!(δV, terms, density, perturbation, cros
         δV[:, :, :, s] .+= divergence_real(density.basis) do α
             ret_α = similar(density.ρ_real, basis.fft_size...)
             ret_α .= 0
-            for t in 1:n_spin
+            for t = 1:n_spin
                 spinfac_st = (t == s ? one(T) : one(T)/2)
                 ret_α .+= -2spinfac_st .* Vσ[tσ(s, t), :, :, :] .* ∇δρ[t, :, :, :, α]
 
-                for u in 1:n_spin
+                for u = 1:n_spin
                     spinfac_su = (s == u ? one(T) : one(T)/2)
                     ret_α .+= (-2spinfac_su .* Vρσ[t, tσ(s, u), :, :, :]
                                .* ∇ρ[u, :, :, :, α] .* δρ[t, :, :, :])
 
-                    for v in 1:n_spin
+                    for v = 1:n_spin
                         spinfac_uv = (u == v ? one(T) : one(T)/2)
                         ret_α .+= (-2spinfac_uv .* spinfac_st
                                    .* Vσσ[tσ(s, t), tσ(u, v), :, :, :]
@@ -496,7 +496,7 @@ for fun in (:potential_terms, :kernel_terms)
         function DftFunctionals.$fun(xcs::Vector{Functional}, density::LibxcDensities)
             isempty(xcs) && return NamedTuple()
             result = $fun(xcs[1], density)
-            for i in 2:length(xcs)
+            for i = 2:length(xcs)
                 result = mergesum(result, $fun(xcs[i], density))
             end
             result
diff --git a/src/transfer.jl b/src/transfer.jl
index 8a53126380..a9c80d9bac 100644
--- a/src/transfer.jl
+++ b/src/transfer.jl
@@ -84,7 +84,7 @@ function compute_transfer_matrix(basis_in::PlaneWaveBasis{T}, basis_out::PlaneWa
     @assert basis_in.model.lattice == basis_out.model.lattice
     @assert length(basis_in.kpoints) == length(basis_out.kpoints)
     @assert all(basis_in.kpoints[ik].coordinate == basis_out.kpoints[ik].coordinate
-                for ik in 1:length(basis_in.kpoints))
+                for ik = 1:length(basis_in.kpoints))
     [compute_transfer_matrix(basis_in, kpt_in, basis_out, kpt_out)
      for (kpt_in, kpt_out) in zip(basis_in.kpoints, basis_out.kpoints)]
 end
@@ -134,7 +134,7 @@ function transfer_blochwave(ψ_in, basis_in::PlaneWaveBasis{T},
     @assert basis_in.model.lattice == basis_out.model.lattice
     @assert length(basis_in.kpoints) == length(basis_out.kpoints)
     @assert all(basis_in.kpoints[ik].coordinate == basis_out.kpoints[ik].coordinate
-                for ik in 1:length(basis_in.kpoints))
+                for ik = 1:length(basis_in.kpoints))
 
     # If, for some kpt ik, basis_in has less vectors than basis_out, then idcs_out[ik] is
     # the array of the indices of the G_vectors from basis_in in basis_out.
diff --git a/src/workarounds/fft_generic.jl b/src/workarounds/fft_generic.jl
index a9b8c4ff0d..492d6a9169 100644
--- a/src/workarounds/fft_generic.jl
+++ b/src/workarounds/fft_generic.jl
@@ -50,13 +50,13 @@ end
 function generic_apply(p::GenericPlan, X::AbstractArray)
     pl1, pl2, pl3 = p.subplans
     ret = similar(X)
-    for i in 1:size(X, 1), j in 1:size(X, 2)
+    for i = 1:size(X, 1), j = 1:size(X, 2)
         @views ret[i, j, :] .= pl3 * X[i, j, :]
     end
-    for i in 1:size(X, 1), k in 1:size(X, 3)
+    for i = 1:size(X, 1), k = 1:size(X, 3)
         @views ret[i, :, k] .= pl2 * ret[i, :, k]
     end
-    for j in 1:size(X, 2), k in 1:size(X, 3)
+    for j = 1:size(X, 2), k = 1:size(X, 3)
         @views ret[:, j, k] .= pl1 * ret[:, j, k]
     end
     p.factor .* ret
diff --git a/src/workarounds/intervals.jl b/src/workarounds/intervals.jl
index 7acbc2f926..9a1c9e3dee 100644
--- a/src/workarounds/intervals.jl
+++ b/src/workarounds/intervals.jl
@@ -54,6 +54,6 @@ function estimate_integer_lattice_bounds(M::AbstractMatrix{<:Interval}, δ, shif
     # As a general statement, with M a lattice matrix, then if ||Mx|| <= δ,
     # then xi = <ei, M^-1 Mx> = <M^-T ei, Mx> <= ||M^-T ei|| δ.
     # Below code does not support non-3D systems.
-    xlims = [norm(inv(M')[:, i]) * δ + shift[i] for i in 1:3]
+    xlims = [norm(inv(M')[:, i]) * δ + shift[i] for i = 1:3]
     map(x -> ceil(Int, x.hi), xlims)
 end
diff --git a/test/PlaneWaveBasis.jl b/test/PlaneWaveBasis.jl
index 1ea3d315a7..0a117ddee9 100644
--- a/test/PlaneWaveBasis.jl
+++ b/test/PlaneWaveBasis.jl
@@ -61,7 +61,7 @@ end
     basis = PlaneWaveBasis(model; Ecut=3, fft_size=(15, 15, 15), kgrid=(1, 1, 1))
     g_all = collect(G_vectors(basis))
 
-    for i in 1:15, j in 1:15, k in 1:15
+    for i = 1:15, j = 1:15, k = 1:15
         @test index_G_vectors(basis, g_all[i, j, k]) == CartesianIndex(i, j, k)
     end
     @test index_G_vectors(basis, [15, 1, 1]) === nothing
@@ -105,7 +105,7 @@ end
     for kpt in basis.kpoints
         Gs_basis = collect(G_vectors(basis))
         Gs_kpt   = collect(G_vectors(basis, kpt))
-        for i in 1:length(kpt.mapping)
+        for i = 1:length(kpt.mapping)
             @test Gs_basis[kpt.mapping[i]] == Gs_kpt[i]
         end
         for i in keys(kpt.mapping_inv)
diff --git a/test/PspHgh.jl b/test/PspHgh.jl
index 8ac17ae62d..bdf828f560 100644
--- a/test/PspHgh.jl
+++ b/test/PspHgh.jl
@@ -95,7 +95,7 @@ end
         @test (res[1] < res[2]) == (res[3] < res[2])
     end
 
-    for i in 1:2, l in 0:2
+    for i = 1:2, l = 0:2
         qcut = qcut_psp_projector(psp, i, l)
         res = eval_psp_projector_fourier.(psp, i, l, [qcut - ε, qcut, qcut + ε])
         @test (res[1] < res[2]) == (res[3] < res[2])
diff --git a/test/PspUpf.jl b/test/PspUpf.jl
index 5a8d4371ea..6542369553 100644
--- a/test/PspUpf.jl
+++ b/test/PspUpf.jl
@@ -17,8 +17,8 @@ upf_pseudos = Dict(
     :Cr => load_psp(artifact"pd_nc_sr_pbe_standard_0.4.1_upf", "Cu.upf"; rcut=12.0)
 )
 hgh_pseudos = [
-    (hgh=load_psp("hgh/pbe/si-q4.hgh"), upf=upf_pseudos[:Si]),
-    (hgh=load_psp("hgh/pbe/tl-q13.hgh"), upf=upf_pseudos[:Tl])
+    (; hgh=load_psp("hgh/pbe/si-q4.hgh"), upf=upf_pseudos[:Si]),
+    (; hgh=load_psp("hgh/pbe/tl-q13.hgh"), upf=upf_pseudos[:Tl])
 ]
 end
 
@@ -76,11 +76,11 @@ end
         hgh = psp_pair.hgh
 
         @test upf.lmax == hgh.lmax
-        for l in 0:upf.lmax
+        for l = 0:upf.lmax
             @test count_n_proj_radial(upf, l) == count_n_proj_radial(hgh, l)
         end
 
-        for l in 0:upf.lmax, i in count_n_proj_radial(upf, l)
+        for l = 0:upf.lmax, i in count_n_proj_radial(upf, l)
             ircut = length(upf.r2_projs[l+1][i])
             for q in (0.01, 0.1, 0.2, 0.5, 1., 2., 5., 10.)
                 reference_hgh = eval_psp_projector_fourier(hgh, i, l, q)
@@ -140,7 +140,7 @@ end
 
     for psp in values(mPspUpf.upf_pseudos)
         ir_start = iszero(psp.rgrid[1]) ? 2 : 1
-        for l in 0:psp.lmax, i in count_n_proj_radial(psp, l)
+        for l = 0:psp.lmax, i in count_n_proj_radial(psp, l)
             ir_cut = min(psp.ircut, length(psp.r2_projs[l+1][i]))
             for q in (0.01, 0.1, 0.2, 0.5, 1., 2., 5., 10.)
                 reference = quadgk(r -> integrand(psp, i, l, q, r),
diff --git a/test/aqua.jl b/test/aqua.jl
index 5b6d6c1d3b..0447f970e3 100644
--- a/test/aqua.jl
+++ b/test/aqua.jl
@@ -7,6 +7,6 @@
     Aqua.test_all(DFTK;
                   ambiguities=false,
                   piracies=false,
-                  deps_compat=(check_extras=false, ),
-                  stale_deps=(ignore=[:Primes, ], ))
+                  deps_compat=(; check_extras=false),
+                  stale_deps=(; ignore=[:Primes, ]))
 end
diff --git a/test/bzmesh.jl b/test/bzmesh.jl
index 89ff6a1387..613f5aa32b 100644
--- a/test/bzmesh.jl
+++ b/test/bzmesh.jl
@@ -13,7 +13,7 @@
         kcoords_spglib = normalize_kpoint_coordinate.(Spglib.eachpoint(spg_mesh))
         sort!(kcoords_spglib)
 
-        kcoords, _ = bzmesh_uniform(kgrid_size; kshift)
+        (; kcoords) = bzmesh_uniform(kgrid_size; kshift)
         sort!(kcoords)
 
         @test kcoords == kcoords_spglib
@@ -43,16 +43,16 @@ end
             system = pyconvert(AbstractSystem, ase_atoms * pytuple(supercell))
         end
 
-        red_kcoords, _ = bzmesh_uniform(kgrid_size; kshift)
+        red_kcoords = bzmesh_uniform(kgrid_size; kshift).kcoords
         symmetries = symmetry_operations(system)
-        irred_kcoords, _ = bzmesh_ir_wedge(kgrid_size, symmetries; kshift)
+        irred_kcoords = bzmesh_ir_wedge(kgrid_size, symmetries; kshift).kcoords
 
         @test length(irred_kcoords) == kirredsize
 
         # Try to reproduce all kcoords from irred_kcoords
         all_kcoords = Vector{Vec3{Rational{Int}}}()
         sym_preserving_grid = DFTK.symmetries_preserving_kgrid(symmetries, red_kcoords)
-        for (ik, k) in enumerate(irred_kcoords)
+        for k in irred_kcoords
             append!(all_kcoords, [symop.S * k for symop in sym_preserving_grid])
         end
 
diff --git a/test/bzmesh_symmetry.jl b/test/bzmesh_symmetry.jl
index ea38b4df0c..ac224e3428 100644
--- a/test/bzmesh_symmetry.jl
+++ b/test/bzmesh_symmetry.jl
@@ -4,11 +4,11 @@
     using LinearAlgebra
     testcase = TestCases.silicon
 
-    args = ((kgrid=[2, 2, 2], kshift=[1/2, 0, 0]),
-            (kgrid=[2, 2, 2], kshift=[1/2, 1/2, 0]),
-            (kgrid=[2, 2, 2], kshift=[0, 0, 0]),
-            (kgrid=[3, 2, 3], kshift=[0, 0, 0]),
-            (kgrid=[3, 2, 3], kshift=[0, 1/2, 1/2]))
+    args = ((; kgrid=[2, 2, 2], kshift=[1/2, 0, 0]),
+            (; kgrid=[2, 2, 2], kshift=[1/2, 1/2, 0]),
+            (; kgrid=[2, 2, 2], kshift=[0, 0, 0]),
+            (; kgrid=[3, 2, 3], kshift=[0, 0, 0]),
+            (; kgrid=[3, 2, 3], kshift=[0, 1/2, 1/2]))
     for case in args
         model_nosym = model_LDA(testcase.lattice, testcase.atoms, testcase.positions;
                                 symmetries=false)
diff --git a/test/chi0.jl b/test/chi0.jl
index 5284990a51..b558696757 100644
--- a/test/chi0.jl
+++ b/test/chi0.jl
@@ -38,7 +38,7 @@ function test_chi0(testcase; symmetries=false, temperature=0, spin_polarization=
         ham0  = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0).ham
         nbandsalg = is_εF_fixed ? FixedBands(; n_bands_converge=6) : AdaptiveBands(model)
         res = DFTK.next_density(ham0, nbandsalg; tol, eigensolver)
-        scfres = (ham=ham0, res...)
+        scfres = (; ham=ham0, res...)
 
         # create external small perturbation εδV
         n_spin = model.n_spin_components
diff --git a/test/compute_bands.jl b/test/compute_bands.jl
index 25cfff2fdc..8309e1595f 100644
--- a/test/compute_bands.jl
+++ b/test/compute_bands.jl
@@ -135,7 +135,7 @@ end
 
     eigres = diagonalize_all_kblocks(lobpcg_hyper, ham, n_bands + 3,
                                      n_conv_check=n_bands, tol=1e-5)
-    for ik in 1:length(basis.kpoints)
+    for ik = 1:length(basis.kpoints)
         @test eigres.λ[ik][1:n_bands] ≈ band_data.λ[ik] atol=1e-5
     end
 end
@@ -162,14 +162,14 @@ end
     @test ret.n_kcoord == 8
     @test ret.n_bands  == 4
 
-    for iband in 1:4
-        @test ret.λ[:, iband, 1] == [10ik .+ iband for ik in 1:8]
+    for iband = 1:4
+        @test ret.λ[:, iband, 1] == [10ik .+ iband for ik = 1:8]
         @test ret.λerror[:, iband, 1] == ret.λ[:, iband, 1] ./ 100
     end
 
     B = model.recip_lattice
     ref_kdist = [0.0]
-    for ik in 2:8
+    for ik = 2:8
         if ik != 4
             push!(ref_kdist, ref_kdist[end] + norm(B * (kinter[ik-1] - kinter[ik])))
         else
diff --git a/test/compute_density.jl b/test/compute_density.jl
index 10929dfc59..0b01dcb26c 100644
--- a/test/compute_density.jl
+++ b/test/compute_density.jl
@@ -14,7 +14,7 @@
         kwargs = ()
         n_bands = div(testcase.n_electrons, 2, RoundUp)
         if testcase.temperature !== nothing
-            kwargs = (temperature=testcase.temperature, smearing=DFTK.Smearing.FermiDirac())
+            kwargs = (; testcase.temperature, smearing=DFTK.Smearing.FermiDirac())
             n_bands = div(testcase.n_electrons, 2, RoundUp) + 4
         end
 
@@ -27,7 +27,7 @@
         occ, εF = DFTK.compute_occupation(basis, res.λ, FermiBisection())
         ρnew = compute_density(basis, res.X, occ)
 
-        for it in 1:n_rounds
+        for it = 1:n_rounds
             ham = Hamiltonian(basis; ρ=ρnew)
             res = diagonalize_all_kblocks(lobpcg_hyper, ham, n_bands; tol, ψguess=res.X)
 
diff --git a/test/compute_jacobian_eigen.jl b/test/compute_jacobian_eigen.jl
index 97fa821c8b..3823a6a140 100644
--- a/test/compute_jacobian_eigen.jl
+++ b/test/compute_jacobian_eigen.jl
@@ -22,7 +22,7 @@ function eigen_ΩplusK(basis::PlaneWaveBasis{T}, ψ, occupation, numval) where {
         precondprep!(Pks[ik], ψk)
     end
     function f_ldiv!(x, y)
-        for n in 1:size(y, 2)
+        for n = 1:size(y, 2)
             δψ = unpack(y[:, n])
             proj_tangent!(δψ, ψ)
             Pδψ = [ Pks[ik] \ δψk for (ik, δψk) in enumerate(δψ)]
diff --git a/test/energies_guess_density.jl b/test/energies_guess_density.jl
index 97224ed3a3..c471a45e33 100644
--- a/test/energies_guess_density.jl
+++ b/test/energies_guess_density.jl
@@ -25,7 +25,7 @@
     # Run one diagonalization and compute energies
     res = diagonalize_all_kblocks(lobpcg_hyper, H, n_bands, tol=1e-9)
     occupation = [[2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0]
-                  for i in 1:length(basis.kpoints)]
+                  for i = 1:length(basis.kpoints)]
     ρ = compute_density(H.basis, res.X, occupation)
     E, H = energy_hamiltonian(basis, res.X, occupation; ρ)
 
diff --git a/test/energy_cutoff_smearing.jl b/test/energy_cutoff_smearing.jl
index 7b43d540ec..a03d0023d7 100644
--- a/test/energy_cutoff_smearing.jl
+++ b/test/energy_cutoff_smearing.jl
@@ -22,7 +22,7 @@
     # Test irregularity of the standard band through its second finite diff derivative
     basis_std = PlaneWaveBasis(model, 5, silicon.kcoords, silicon.kweights; basis.fft_size)
     λ_std = vcat(compute_bands(basis_std, kcoords, n_bands=1; scfres.ρ).λ...)
-    ∂2λ_std = [(λ_std[i+1] - 2*λ_std[i] + λ_std[i-1])/δk^2 for i in 2:length(kcoords)-1]
+    ∂2λ_std = [(λ_std[i+1] - 2*λ_std[i] + λ_std[i-1])/δk^2 for i = 2:length(kcoords)-1]
 
     # Compute band for given blow-up and test regularity
     function test_blowup(kinetic_blowup)
@@ -30,7 +30,7 @@
         basis_mod = PlaneWaveBasis(model, 5, silicon.kcoords, silicon.kweights; basis.fft_size)
 
         λ_mod = vcat(compute_bands(basis_mod, kcoords, n_bands=1, ρ=scfres.ρ).λ...)
-        ∂2λ_mod = [(λ_mod[i+1] - 2*λ_mod[i] + λ_mod[i-1])/δk^2 for i in 2:length(kcoords)-1]
+        ∂2λ_mod = [(λ_mod[i+1] - 2*λ_mod[i] + λ_mod[i-1])/δk^2 for i = 2:length(kcoords)-1]
         @test norm(∂2λ_std) / norm(∂2λ_mod) > 1e4
         nothing
     end
diff --git a/test/external/atomsbase.jl b/test/external/atomsbase.jl
index 9e926414d6..b5aa05ef13 100644
--- a/test/external/atomsbase.jl
+++ b/test/external/atomsbase.jl
@@ -10,7 +10,7 @@
 
     lattice   = randn(3, 3)
     atoms     = [Si, C, H, C]
-    positions = [rand(3) for _ in 1:4]
+    positions = [rand(3) for _ = 1:4]
     magnetic_moments = rand(4)
 
     system = atomic_system(lattice, atoms, positions, magnetic_moments)
@@ -26,7 +26,7 @@
     parsed = DFTK.parse_system(system)
     @test parsed.lattice   ≈ lattice   atol=5e-13
     @test parsed.positions ≈ positions atol=5e-13
-    for i in 1:4
+    for i = 1:4
         @test iszero(parsed.magnetic_moments[i][1:2])
         @test parsed.magnetic_moments[i][3] == magnetic_moments[i]
     end
@@ -65,7 +65,7 @@ end
 
     lattice   = randn(3, 3)
     atoms     = [ElementCoulomb(:Si), ElementCoulomb(:C)]
-    positions = [rand(3) for _ in 1:2]
+    positions = [rand(3) for _ = 1:2]
     magnetic_moments = [rand(3), rand(3)]
     system = atomic_system(lattice, atoms, positions, magnetic_moments)
     @test system[:, :magnetic_moment] == magnetic_moments
@@ -78,14 +78,14 @@ end
     using AtomsBase
 
     @testset "Charged system" begin
-        lattice = [12u"bohr" * rand(3) for _ in 1:3]
+        lattice = [12u"bohr" * rand(3) for _ = 1:3]
         atoms   = [:C => rand(3), :Si => rand(3), :H => rand(3), :C => rand(3)]
         system  = periodic_system(atoms, lattice; fractional=true, charge=1.0u"e_au")
         @test_throws ErrorException Model(system)
     end
 
     @testset "Charged atoms, but neutral" begin
-        lattice = [12u"bohr" * rand(3) for _ in 1:3]
+        lattice = [12u"bohr" * rand(3) for _ = 1:3]
         atoms   = [Atom(:C,  rand(3) * 12u"bohr", charge=1.0u"e_au"),
                    Atom(:Si, rand(3) * 12u"bohr", charge=-1.0u"e_au")]
         system  = periodic_system(atoms, lattice)
@@ -94,7 +94,7 @@ end
     end
 
     @testset "Charged atoms and not neutral" begin
-        lattice = [12u"bohr" * rand(3) for _ in 1:3]
+        lattice = [12u"bohr" * rand(3) for _ = 1:3]
         atoms   = [Atom(:C,  rand(3) * 12u"bohr", charge=1.0u"e_au"),
                    Atom(:Si, rand(3) * 12u"bohr", charge=-2.0u"e_au")]
         system  = periodic_system(atoms, lattice)
@@ -108,7 +108,7 @@ end
     using UnitfulAtomic
     using AtomsBase
 
-    lattice     = [12u"bohr" * rand(3) for _ in 1:3]
+    lattice     = [12u"bohr" * rand(3) for _ = 1:3]
     weirdatom   = Atom(6, randn(3)u"Å"; atomic_symsymbol=:C1)
     atoms       = [:C => rand(3), :Si => rand(3), :H => rand(3), :C => rand(3)]
     pos_units   = last.(atoms)
@@ -179,7 +179,7 @@ end
     H  = ElementPsp(:H; psp=load_psp("hgh/lda/h-q1.hgh"))
     lattice   = randn(3, 3)
     atoms     = [Si, C, H, C]
-    positions = [rand(3) for _ in 1:4]
+    positions = [rand(3) for _ = 1:4]
     system    = atomic_system(lattice, atoms, positions)
 
     @test_throws ErrorException attach_psp(system; Si="hgh/pbe/si-q4.hgh")
@@ -197,7 +197,7 @@ end
     using UnitfulAtomic
     using AtomsBase
 
-    lattice     = [12u"bohr" * rand(3) for _ in 1:3]
+    lattice     = [12u"bohr" * rand(3) for _ = 1:3]
     atoms       = [Atom(6, randn(3)u"Å"; atomic_symbol=:C1),
                    Atom(6, randn(3)u"Å"; atomic_symbol=:C2)]
     system      = periodic_system(atoms, lattice)
diff --git a/test/forwarddiff.jl b/test/forwarddiff.jl
index 4bcdb22c8e..95fb1a7887 100644
--- a/test/forwarddiff.jl
+++ b/test/forwarddiff.jl
@@ -17,7 +17,7 @@
         end
         basis = PlaneWaveBasis(model; Ecut=5, kgrid=[2, 2, 2], kshift=[0, 0, 0])
 
-        response     = ResponseOptions(verbose=true)
+        response     = ResponseOptions(; verbose=true)
         is_converged = DFTK.ScfConvergenceForce(tol)
         scfres = self_consistent_field(basis; is_converged, response)
         compute_forces_cart(scfres)
@@ -78,7 +78,7 @@ end
         is_converged = DFTK.ScfConvergenceDensity(1e-10)
         scfres = self_consistent_field(basis; is_converged, mixing=KerkerMixing(),
                                        nbandsalg=FixedBands(; n_bands_converge=10),
-                                       damping=0.6, response=ResponseOptions(verbose=true))
+                                       damping=0.6, response=ResponseOptions(; verbose=true))
 
         ComponentArray(
            eigenvalues=hcat([ev[1:10] for ev in scfres.eigenvalues]...),
@@ -116,7 +116,7 @@ end
 
         is_converged = DFTK.ScfConvergenceDensity(1e-10)
         scfres = self_consistent_field(basis; is_converged,
-                                       response=ResponseOptions(verbose=true))
+                                       response=ResponseOptions(; verbose=true))
         compute_forces_cart(scfres)
     end
 
@@ -161,7 +161,7 @@ end
         ρ = zeros(Float64, basis.fft_size..., 1)
         is_converged = DFTK.ScfConvergenceDensity(1e-10)
         scfres = self_consistent_field(basis; ρ, is_converged,
-                                       response=ResponseOptions(verbose=true))
+                                       response=ResponseOptions(; verbose=true))
         compute_forces_cart(scfres)
     end
     derivative_ε = let ε = 1e-5
diff --git a/test/hamiltonian_consistency.jl b/test/hamiltonian_consistency.jl
index b64502cba6..11128503d8 100644
--- a/test/hamiltonian_consistency.jl
+++ b/test/hamiltonian_consistency.jl
@@ -40,7 +40,7 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2,
 
         ψ = [Matrix(qr(randn(ComplexF64, length(G_vectors(basis, kpt)), n_bands)).Q)
              for kpt in basis.kpoints]
-        occupation = [filled_occ * rand(n_bands) for _ in 1:length(basis.kpoints)]
+        occupation = [filled_occ * rand(n_bands) for _ = 1:length(basis.kpoints)]
         occ_scaling = n_electrons / sum(sum(occupation))
         occupation = [occ * occ_scaling for occ in occupation]
         ρ = with_logger(NullLogger()) do
@@ -63,7 +63,7 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2,
         diff = (compute_E(ε) - compute_E(-ε)) / (2ε)
 
         diff_predicted = 0.0
-        for ik in 1:length(basis.kpoints)
+        for ik = 1:length(basis.kpoints)
             Hψk = ham.blocks[ik]*ψ[ik]
             test_matrix_repr_operator(ham.blocks[ik], ψ[ik]; atol)
             δψkHψk = sum(occupation[ik][iband] * real(dot(δψ[ik][:, iband], Hψk[:, iband]))
diff --git a/test/helium_all_electron.jl b/test/helium_all_electron.jl
index fa4a9c01f3..ff85aa9602 100644
--- a/test/helium_all_electron.jl
+++ b/test/helium_all_electron.jl
@@ -13,7 +13,7 @@
         scfres.energies.total, DFTK.compute_forces(scfres)
     end
 
-    E, forces = energy_forces(Ecut=5, tol=1e-10)
+    E, forces = energy_forces(; Ecut=5, tol=1e-10)
     @test E ≈ -1.5869009433016852 atol=1e-12
     @test norm(forces) < 1e-7
 end
diff --git a/test/hessian.jl b/test/hessian.jl
index f992fbbdf3..342bfec0b2 100644
--- a/test/hessian.jl
+++ b/test/hessian.jl
@@ -92,9 +92,9 @@ end
     @testset "solve_ΩplusK_split agrees with solve_ΩplusK" begin
         scfres = self_consistent_field(basis; tol=1e-10)
         δψ1 = solve_ΩplusK_split(scfres, rhs).δψ
-        δψ1, _ = select_occupied_orbitals(basis, δψ1, scfres.occupation)
-        ψ, occupation = select_occupied_orbitals(basis, scfres.ψ, scfres.occupation)
-        rhs_trunc, _ = select_occupied_orbitals(basis, rhs, occupation)
+        δψ1 = select_occupied_orbitals(basis, δψ1, scfres.occupation).ψ
+        (; ψ, occupation) = select_occupied_orbitals(basis, scfres.ψ, scfres.occupation)
+        rhs_trunc = select_occupied_orbitals(basis, rhs, occupation).ψ
         δψ2 = solve_ΩplusK(basis, ψ, rhs_trunc, occupation).δψ
         @test norm(δψ1 - δψ2) < 1e-7
     end
diff --git a/test/kernel.jl b/test/kernel.jl
index cb90450b60..e153339db8 100644
--- a/test/kernel.jl
+++ b/test/kernel.jl
@@ -32,7 +32,7 @@
             ops_plus  = DFTK.ene_ops(term, basis, nothing, nothing; ρ=ρ_plus).ops
             δV = zero(ρ0)
 
-            for iσ in 1:model.n_spin_components
+            for iσ = 1:model.n_spin_components
                 # Index of the first spin-up or spin-down k-point
                 ifirst = first(DFTK.krange_spin(basis, iσ))
                 δV[:, :, :, iσ] = (ops_plus[ifirst].potential - ops_minus[ifirst].potential) / (2ε)
diff --git a/test/lobpcg.jl b/test/lobpcg.jl
index 7bfa6d1089..4de46634e8 100644
--- a/test/lobpcg.jl
+++ b/test/lobpcg.jl
@@ -28,7 +28,7 @@
                                       prec_type=nothing, interpolate_kpoints=false)
 
         @test res.converged
-        for ik in 1:length(basis.kpoints)
+        for ik = 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.n_iter[ik] < 200
@@ -40,7 +40,7 @@
                                       prec_type=PreconditionerTPA, interpolate_kpoints=false)
 
         @test res.converged
-        for ik in 1:length(basis.kpoints)
+        for ik = 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.n_iter[ik] < 50
@@ -73,7 +73,7 @@ end
         [-4.085991608422304, -4.085039856878318, -0.517299903754010,
          -0.513805498246478, -0.497036479690380]
     ]
-    for ik in 1:length(basis.kpoints)
+    for ik = 1:length(basis.kpoints)
         @test res.λ[ik][1:5] ≈ ref[basis.krange_thisproc[ik]] atol=5e-7
     end
 end
@@ -103,7 +103,7 @@ end
         [0.168662148987539, 0.238552367551507, 0.370743978236562,
          0.418387442903058, 0.619797227001203],
     ]
-    for ik in 1:length(basis.kpoints)
+    for ik = 1:length(basis.kpoints)
         @test res.λ[ik] ≈ ref[basis.krange_thisproc[ik]] atol=0.02
     end
 end
@@ -121,7 +121,7 @@ end
 
     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)
+    for ik = 1:length(basis.kpoints)
         @test res1.λ[ik] ≈ res2.λ[ik] atol=1e-6
     end
 end
diff --git a/test/occupation.jl b/test/occupation.jl
index 59b9f7f5c5..32ebc40997 100644
--- a/test/occupation.jl
+++ b/test/occupation.jl
@@ -46,12 +46,12 @@ end
     eigenvalues = [zeros(n_bands) for k in testcase.kcoords]
     n_occ = div(testcase.n_electrons, 2, RoundUp)
     n_k = length(testcase.kcoords)
-    for ik in 1:n_k
+    for ik = 1:n_k
         eigenvalues[ik] = sort(rand(n_bands))
         eigenvalues[ik][n_occ+1:end] .+= 2
     end
-    εHOMO = maximum(eigenvalues[ik][n_occ] for ik in 1:n_k)
-    εLUMO = minimum(eigenvalues[ik][n_occ + 1] for ik in 1:n_k)
+    εHOMO = maximum(eigenvalues[ik][n_occ] for ik = 1:n_k)
+    εLUMO = minimum(eigenvalues[ik][n_occ + 1] for ik = 1:n_k)
 
     # Occupation for zero temperature
     occupation0 = let
@@ -70,8 +70,8 @@ end
         model = Model(testcase.lattice, testcase.atoms, testcase.positions;
                       temperature, smearing, terms=[Kinetic()])
         basis = PlaneWaveBasis(model, Ecut, testcase.kcoords, testcase.kweights; fft_size)
-        occs, _ = with_logger(NullLogger()) do
-            DFTK.compute_occupation(basis, eigenvalues, alg; tol_n_elec=1e-12)
+        occs = with_logger(NullLogger()) do
+            DFTK.compute_occupation(basis, eigenvalues, alg; tol_n_elec=1e-12).occupation
         end
         @test sum(basis.kweights .* sum.(occs)) ≈ model.n_electrons
     end
@@ -82,9 +82,9 @@ end
         model = Model(testcase.lattice, testcase.atoms, testcase.positions;
                       temperature, smearing, terms=[Kinetic()])
         basis = PlaneWaveBasis(model, Ecut, testcase.kcoords, testcase.kweights; fft_size)
-        occupation, _ = DFTK.compute_occupation(basis, eigenvalues, alg; tol_n_elec=1e-6)
+        (; occupation) = DFTK.compute_occupation(basis, eigenvalues, alg; tol_n_elec=1e-6)
 
-        for ik in 1:n_k
+        for ik = 1:n_k
             @test all(isapprox.(occupation[ik], occupation0[ik], atol=1e-2))
         end
     end
diff --git a/test/phonon/helpers.jl b/test/phonon/helpers.jl
index 948d17872c..1787c06228 100644
--- a/test/phonon/helpers.jl
+++ b/test/phonon/helpers.jl
@@ -21,7 +21,7 @@ function ph_compute_reference(basis_supercell)
     n_dim = model_supercell.n_dim
     T = eltype(model_supercell.lattice)
     dynmat_ad = zeros(T, 3, n_atoms, 3, n_atoms)
-    for τ in 1:n_atoms, γ in 1:n_dim
+    for τ = 1:n_atoms, γ = 1:n_dim
         displacement = zero.(model_supercell.positions)
         displacement[τ] = setindex(displacement[τ], one(T), γ)
         dynmat_ad[:, :, γ, τ] = -ForwardDiff.derivative(zero(T)) do ε
@@ -51,7 +51,7 @@ end
 function generate_supercell_qpoints(; supercell_size=generate_random_supercell())
     qpoints_list = Iterators.product([1:n_sc for n_sc in supercell_size]...)
     qpoints = map(qpoints_list) do n_sc
-        normalize_kpoint_coordinate.([n_sc[i] / supercell_size[i] for i in 1:3])
+        normalize_kpoint_coordinate.([n_sc[i] / supercell_size[i] for i = 1:3])
     end |> vec
 
     (; supercell_size, qpoints)
diff --git a/test/phonon/lowlevel.jl b/test/phonon/lowlevel.jl
index 28bcb2a1bd..9c0181e0a6 100644
--- a/test/phonon/lowlevel.jl
+++ b/test/phonon/lowlevel.jl
@@ -8,7 +8,7 @@
     function transfer_blochwave_kpt_real(ψk_in, basis::PlaneWaveBasis, kpt_in, kpt_out, ΔG)
         ψk_out = zeros(eltype(ψk_in), length(kpt_out.G_vectors), size(ψk_in, 2))
         exp_ΔGr = cis2pi.(-dot.(Ref(ΔG), r_vectors(basis)))
-        for n in 1:size(ψk_in, 2)
+        for n = 1:size(ψk_in, 2)
             ψk_out[:, n] = fft(basis, kpt_out, exp_ΔGr .* ifft(basis, kpt_in, ψk_in[:, n]))
         end
         ψk_out
@@ -19,7 +19,7 @@
 
     positions = [[0.0, 0.0, 0.0]]
     n_scell = 2
-    for i in 1:n_scell-1
+    for i = 1:n_scell-1
         push!(positions, i * ones(3) / n_scell)
     end
     n_atoms = length(positions)
diff --git a/test/scf_compare.jl b/test/scf_compare.jl
index bae623368c..feb268551b 100644
--- a/test/scf_compare.jl
+++ b/test/scf_compare.jl
@@ -28,8 +28,7 @@
         @testset "Newton" begin
             scfres_start = self_consistent_field(basis, maxiter=1)
             # remove virtual orbitals
-            ψ0, _ = select_occupied_orbitals(basis, scfres_start.ψ,
-                                             scfres_start.occupation)
+            ψ0 = select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation).ψ
             ρ_newton = newton(basis, ψ0; tol).ρ
             @test maximum(abs, ρ_newton - ρ_def) < 5tol
         end
@@ -81,8 +80,7 @@ end
     ρ0    = guess_density(basis, magnetic_moments)
     ρ_def = self_consistent_field(basis; tol, ρ=ρ0).ρ
     scfres_start = self_consistent_field(basis, maxiter=1, ρ=ρ0)
-    ψ0, _ = select_occupied_orbitals(basis, scfres_start.ψ,
-                                     scfres_start.occupation)
+    ψ0 = select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation).ψ
 
     # Run DM
     if mpi_nprocs() == 1  # Distributed implementation not yet available
@@ -118,8 +116,8 @@ end
     ρ0    = guess_density(basis)
     ρ_ref = self_consistent_field(basis; ρ=ρ0, tol).ρ
 
-    for mixing_str in ("KerkerDosMixing()", "HybridMixing(RPA=true)",
-                       "LdosMixing(RPA=false)", "HybridMixing(εr=10, RPA=true)")
+    for mixing_str in ("KerkerDosMixing()", "HybridMixing(; RPA=true)",
+                       "LdosMixing(; RPA=false)", "HybridMixing(; εr=10, RPA=true)")
         @testset "Testing $mixing_str" begin
             mixing = eval(Meta.parse(mixing_str))
             ρ_mix = self_consistent_field(basis; ρ=ρ0, mixing, tol, damping=0.8).ρ
@@ -148,8 +146,8 @@ end
     scfres = self_consistent_field(basis; ρ=ρ0, tol)
     ρ_ref  = scfres.ρ
 
-    for mixing_str in ("KerkerMixing()", "KerkerDosMixing()", "DielectricMixing(εr=10)",
-                       "HybridMixing(εr=10)", "χ0Mixing(; χ0terms=[Applyχ0Model()], RPA=false)")
+    for mixing_str in ("KerkerMixing()", "KerkerDosMixing()", "DielectricMixing(; εr=10)",
+                       "HybridMixing(; εr=10)", "χ0Mixing(; χ0terms=[Applyχ0Model()], RPA=false)")
         @testset "Testing $mixing_str" begin
             mixing = eval(Meta.parse(mixing_str))
             scfres = self_consistent_field(basis; ρ=ρ0, mixing, tol, damping=0.8)
diff --git a/test/stresses.jl b/test/stresses.jl
index 988e5a12a0..7ac7678700 100644
--- a/test/stresses.jl
+++ b/test/stresses.jl
@@ -17,14 +17,14 @@
     function recompute_energy(lattice, symmetries, element)
         basis = make_basis(lattice, symmetries, element)
         scfres = self_consistent_field(basis; is_converged=DFTK.ScfConvergenceDensity(1e-13))
-        energies, H = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)
+        (; energies) = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)
         energies.total
     end
 
     function hellmann_feynman_energy(scfres, lattice, symmetries, element)
         basis = make_basis(lattice, symmetries, element)
         ρ = DFTK.compute_density(basis, scfres.ψ, scfres.occupation)
-        energies, H = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ)
+        (; energies) = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ)
         energies.total
     end