Consistency and documentation fixes (#920)

.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

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

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",

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)

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)]));
 

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;

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

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:

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

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)

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

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

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)

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)

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
 

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"Å"),
             )
 

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,