Merge branch 'master' into mpi-aware-omegaplusK

.github/workflows/ci.yaml

@@ -25,7 +25,7 @@ jobs:
           - {mode: stable,  os: macOS-latest,   payload: noslow         }
           - {mode: stable,  os: windows-latest, payload: noslow         }
           - {mode: stable,  os: ubuntu-latest,  payload: noslow-mpi     }
-          - {mode: nightly, os: ubuntu-latest,  payload: noslow         }
+          - {mode: latest,  os: ubuntu-latest,  payload: noslow         }
     env:
       GKS_ENCODING: utf8
       GKSwstype: 100       # Needed for Plots-related tests
@@ -41,12 +41,13 @@ jobs:
           version: '1.10'
           arch: x64
         if: ${{ matrix.mode == 'stable' }}
-      - name: Setup Julia nightly
+      - name: Setup Julia latest (pre-releases included)
         uses: julia-actions/setup-julia@v1
         with:
-          version: nightly
+          version: '1'
+          include-all-prereleases: true
           arch: x64
-        if: ${{ matrix.mode == 'nightly' }}
+        if: ${{ matrix.mode == 'latest' }}
 
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1

docs/src/guide/introductory_resources.md

@@ -34,7 +34,6 @@ see [Publications](@ref).
   in particular the [summary of DFT theory](https://michael-herbst.com/teaching/2022-mit-workshop-dftk/2022-mit-workshop-dftk/DFT_Theory.pdf).
 
 ## Textbooks and reviews
-
 - [Density Functional Theory](https://doi.org/10.1007/978-3-031-22340-2)
   edited by Eric Cancès and Gero Friesecke (Springer, 2023):
   Up to date textbook accessible to an interdisciplinary audience with contributions
@@ -58,14 +57,19 @@ see [Publications](@ref).
   Covers topics such as DFT, pseudos, SCF, response, ...
 
 ## Recordings
+- [DFTK.jl: 5 years of a multidisciplinary electronic-structure code](https://www.youtube.com/watch?v=ox_j2zKOuIk) by M. F. Herbst:
+  30-min talk at JuliaCon 2024 providing the state of DFTK 5 years after the project was started.
+  [Slides](https://michael-herbst.com/talks/2024.07.12_5years_DFTK.pdf)
+  [Pluto notebook](https://michael-herbst.com/talks/2024.07.12_5years_DFTK.html)
+
 - [Julia for Materials Modelling](https://www.youtube.com/watch?v=dujepKxxxkg) by M. F. Herbst:
   One-hour talk providing an overview of materials modelling tools for Julia.
   Key features of DFTK are highlighted as part of the talk.
   [Pluto notebooks](https://mfherbst.github.io/julia-for-materials/)
 
 - [DFTK: A Julian approach for simulating electrons in solids](https://www.youtube.com/watch?v=-RomkxjlIcQ) by M. F. Herbst:
   Pre-recorded talk for JuliaCon 2020.
-  Assumes no knowledge about DFT and gives the broad picture of DFTK.
+  Assumes no knowledge about DFT and gives the broad picture of DFTK. Starts to become a little outdated.
   [Slides](https://michael-herbst.com/talks/2020.07.29_juliacon_dftk.pdf).
 
 - [Juliacon 2021 DFT workshop](https://www.youtube.com/watch?v=HvpPMWVm8aw) by M. F. Herbst:

examples/dos.jl

@@ -0,0 +1,31 @@
+# # Densities of states (DOS)
+# In this example, we'll plot the DOS, local DOS, and projected DOS of Silicon.
+# DOS computation only supports finite temperature.
+
+using DFTK
+using Unitful
+using Plots
+using LazyArtifacts
+
+## Define the geometry and pseudopotential
+a = 10.26  # Silicon lattice constant in Bohr
+lattice = a / 2 * [[0 1 1.0];
+    [1 0 1.0];
+    [1 1 0.0]]
+Si = ElementPsp(:Si; psp=load_psp(artifact"pd_nc_sr_lda_standard_0.4.1_upf", "Si.upf"))
+atoms = [Si, Si]
+positions = [ones(3) / 8, -ones(3) / 8]
+
+## Run SCF
+model = model_LDA(lattice, atoms, positions; temperature=5e-3)
+basis = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4], symmetries_respect_rgrid=true)
+scfres = self_consistent_field(basis, tol=1e-8)
+
+## Plot the DOS
+plot_dos(scfres)
+
+## Plot the local DOS along one direction
+plot_ldos(scfres; n_points=100, ldos_xyz=[:, 10, 10])
+
+## Plot the projected DOS
+plot_pdos(scfres; εrange=(-0.3, 0.5))

ext/DFTKPlotsExt.jl

@@ -2,7 +2,7 @@ module DFTKPlotsExt
 using Brillouin: KPath
 using DFTK
 using DFTK: is_metal, data_for_plotting, spin_components, default_band_εrange
-import DFTK: plot_dos, plot_bandstructure
+import DFTK: plot_dos, plot_bandstructure, plot_ldos, plot_pdos
 using Plots
 using Unitful
 using UnitfulAtomic
@@ -108,4 +108,88 @@ function plot_dos(basis, eigenvalues; εF=nothing, unit=u"hartree",
 end
 plot_dos(scfres; kwargs...) = plot_dos(scfres.basis, scfres.eigenvalues; scfres.εF, kwargs...)
 
+function plot_ldos(basis, eigenvalues, ψ; εF=nothing, unit=u"hartree",
+                   temperature=basis.model.temperature,
+                   smearing=basis.model.smearing,
+                   εrange=default_band_εrange(eigenvalues; εF), 
+                   n_points=1000, ldos_xyz=[:, 1, 1], kwargs...)
+    eshift = something(εF, 0.0)
+    εs = range(austrip.(εrange)..., length=n_points)
+
+    # Constant to convert from AU to the desired unit
+    to_unit = ustrip(auconvert(unit, 1.0))
+
+    # LDε has three dimensions (x, y, z)
+    # map on a single axis to plot the variation with εs
+    LDεs = dropdims.(compute_ldos.(εs, Ref(basis), Ref(eigenvalues), Ref(ψ); smearing, temperature); dims=4)
+    LDεs_slice = similar(LDεs[1], n_points, length(LDεs[1][ldos_xyz...]))
+    for (i, LDε) in enumerate(LDεs)
+        LDεs_slice[i, :] = LDε[ldos_xyz...]
+    end
+    p = heatmap(1:size(LDεs_slice, 2), (εs .- eshift) .* to_unit, LDεs_slice; kwargs...)
+    if !isnothing(εF)
+        Plots.hline!(p, [0.0], label="εF", color=:green, lw=1.5)
+    end
+
+    if isnothing(εF)
+        Plots.ylabel!(p, "eigenvalues  ($(string(unit)))")
+    else
+        Plots.ylabel!(p, "eigenvalues -ε_F  ($(string(unit)))")
+    end
+    p
+end
+plot_ldos(scfres; kwargs...) = plot_ldos(scfres.basis, scfres.eigenvalues, scfres.ψ; scfres.εF, kwargs...)
+
+function plot_pdos(basis, eigenvalues, ψ, i, l, 
+                   psp, position, el::Symbol;
+                   εF=nothing, unit=u"hartree",
+                   temperature=basis.model.temperature,
+                   smearing=basis.model.smearing,
+                   εrange=default_band_εrange(eigenvalues; εF), 
+                   n_points=1000, p=nothing, kwargs...)
+    eshift = something(εF, 0.0)
+    εs = range(austrip.(εrange)..., length=n_points)
+
+    # Constant to convert from AU to the desired unit
+    to_unit = ustrip(auconvert(unit, 1.0))
+
+    # Calculate the projections of the atom with given i and l,
+    # and sum all angular momentums m=-l:l 
+    pdos = dropdims(sum(compute_pdos(εs, basis, eigenvalues, ψ, i, l, 
+                    psp, position; temperature, smearing), dims=2); dims=2)
+    label = String(el) * "-" * psp.pswfc_labels[l+1][i]
+
+    # Plot pdos 
+    p = something(p, Plots.plot(; kwargs...))
+    Plots.plot!(p, (εs .- eshift) .* to_unit, pdos; label)
+
+    p
+end
+
+function plot_pdos(scfres; kwargs...)
+    # Plot DOS
+    p = plot_dos(scfres; scfres.εF, kwargs...)
+
+    # TODO do the symmetrization instead of unfolding    
+    scfres_unfold = DFTK.unfold_bz(scfres)
+    basis = scfres_unfold.basis
+    psp_groups = [group for group in basis.model.atom_groups
+                  if basis.model.atoms[first(group)] isa ElementPsp]
+
+    # Plot PDOS for the first atom of each atom group
+    for group in psp_groups
+        psp = basis.model.atoms[first(group)].psp
+        position = basis.model.positions[first(group)]
+        el = basis.model.atoms[first(group)].symbol
+        for l = 0:psp.lmax
+            for i = 1:DFTK.count_n_pswfc_radial(psp, l)
+                plot_pdos(basis, scfres_unfold.eigenvalues, scfres_unfold.ψ, 
+                          i, l, psp, position, el; scfres.εF, p, kwargs...)
+            end
+        end                  
+    end
+
+    p
+end
+
 end

src/DFTK.jl

@@ -212,7 +212,10 @@ export compute_stresses_cart
 include("postprocess/stresses.jl")
 export compute_dos
 export compute_ldos
+export compute_pdos
 export plot_dos
+export plot_ldos
+export plot_pdos
 include("postprocess/dos.jl")
 export compute_χ0
 export apply_χ0

src/postprocess/dos.jl

@@ -7,6 +7,10 @@
 #
 # LDOS (local density of states)
 # LD(ε) = sum_n f_n' |ψn|^2 = sum_n δ(ε - ε_n) |ψn|^2
+#
+# PD(ε) = sum_n f_n' |<ψn,ϕ>|^2
+# ϕ = ∑_R ϕilm(r-pos-R) is the periodized atomic wavefunction, obtained from the pseudopotential
+# This is computed for a given (i,l) (eg i=2,l=2 for the 3p) and summed over all m
 
 """
 Total density of states at energy ε
@@ -63,7 +67,73 @@ function compute_ldos(scfres::NamedTuple; ε=scfres.εF, kwargs...)
     compute_ldos(ε, scfres.basis, scfres.eigenvalues, scfres.ψ; kwargs...)
 end
 
+"""
+Projected density of states at energy ε for an atom with given i and l.
+"""
+function compute_pdos(ε, basis, eigenvalues, ψ, i, l, psp, position;
+                      smearing=basis.model.smearing,
+                      temperature=basis.model.temperature)
+    if (temperature == 0) || smearing isa Smearing.None
+        error("compute_pdos only supports finite temperature")
+    end
+    filled_occ = filled_occupation(basis.model)
+
+    projs = compute_pdos_projs(basis, eigenvalues, ψ, i, l, psp, position)
+
+    D = zeros(typeof(ε[1]), length(ε), 2l+1)
+    for (i, iε) in enumerate(ε)
+        for (ik, projk) in enumerate(projs)
+            @views for (iband, εnk) in enumerate(eigenvalues[ik])
+                enred = (εnk - iε) / temperature
+                D[i, :] .-= (filled_occ .* basis.kweights[ik] .* projk[iband, :]
+                             ./ temperature
+                             .* Smearing.occupation_derivative(smearing, enred))       
+            end
+        end
+    end
+    D = mpi_sum(D, basis.comm_kpts)
+end
+
+function compute_pdos(scfres::NamedTuple, iatom, i, l; ε=scfres.εF, kwargs...)
+    psp = scfres.basis.model.atoms[iatom].psp
+    position = scfres.basis.model.positions[iatom]
+    # TODO do the symmetrization instead of unfolding    
+    scfres = unfold_bz(scfres)
+    compute_pdos(ε, scfres.basis, scfres.eigenvalues, scfres.ψ, i, l, psp, position; kwargs...)
+end
+
+# Build atomic orbitals projections projs[ik][iband,m] = |<ψnk, ϕ>|^2 for a single atom.
+# TODO optimization ? accept a range of positions, in case we want to compute the PDOS
+# for all atoms of the same type (and reuse the computation of the atomic orbitals)
+function compute_pdos_projs(basis, eigenvalues, ψ, i, l, psp::NormConservingPsp, position)
+    # Precompute the form factors on all kpoints at once so we can better use the cache (memory-compute tradeoff).
+    # Revisit this (pass the cache around explicitly) if RAM becomes an issue.
+    G_plus_k_all = [Gplusk_vectors(basis, basis.kpoints[ik])
+                    for ik = 1:length(basis.kpoints)]
+    G_plus_k_all_cart = [map(recip_vector_red_to_cart(basis.model), gpk) 
+                         for gpk in G_plus_k_all]
+
+    # Build form factors of pseudo-wavefunctions centered at 0.
+    fun(p) = eval_psp_pswfc_fourier(psp, i, l, p)
+    # form_factors_all[ik][p,m]
+    form_factors_all = build_form_factors(fun, l, G_plus_k_all_cart)
+
+    projs = Vector{Matrix}(undef, length(basis.kpoints))
+    for (ik, ψk) in enumerate(ψ)
+        structure_factor = [cis2pi(-dot(position, p)) for p in G_plus_k_all[ik]]
+        # TODO orthogonalize pseudo-atomic wave functions?
+        proj_vectors = structure_factor .* form_factors_all[ik] ./ sqrt(basis.model.unit_cell_volume)
+        projs[ik] = abs2.(ψk' * proj_vectors) # contract on p -> projs[ik][iband,m]
+    end
+
+    projs
+end
+
 """
 Plot the density of states over a reasonable range. Requires to load `Plots.jl` beforehand.
 """
 function plot_dos end
+
+function plot_ldos end
+
+function plot_pdos end

src/pseudo/NormConservingPsp.jl

@@ -11,20 +11,24 @@ abstract type NormConservingPsp end
 #    description::String         # Descriptive string
 
 #### Methods:
-# charge_ionic(psp::NormConservingPsp)
-# has_valence_density(psp:NormConservingPsp)
-# has_core_density(psp:NormConservingPsp)
-# eval_psp_projector_real(psp::NormConservingPsp, i, l, r::Real)
-# eval_psp_projector_fourier(psp::NormConservingPsp, i, l, p::Real)
-# eval_psp_local_real(psp::NormConservingPsp, r::Real)
-# eval_psp_local_fourier(psp::NormConservingPsp, p::Real)
-# eval_psp_energy_correction(T::Type, psp::NormConservingPsp, n_electrons::Integer)
+# charge_ionic(psp)
+# has_valence_density(psp)
+# has_core_density(psp)
+# eval_psp_projector_real(psp, i, l, r::Real)
+# eval_psp_projector_fourier(psp, i, l, p::Real)
+# eval_psp_local_real(psp, r::Real)
+# eval_psp_local_fourier(psp, p::Real)
+# eval_psp_energy_correction(T::Type, psp, n_electrons::Integer)
 
 #### Optional methods:
-# eval_psp_density_valence_real(psp::NormConservingPsp, r::Real)
-# eval_psp_density_valence_fourier(psp::NormConservingPsp, p::Real)
-# eval_psp_density_core_real(psp::NormConservingPsp, r::Real)
-# eval_psp_density_core_fourier(psp::NormConservingPsp, p::Real)
+# eval_psp_density_valence_real(psp, r::Real)
+# eval_psp_density_valence_fourier(psp, p::Real)
+# eval_psp_density_core_real(psp, r::Real)
+# eval_psp_density_core_fourier(psp, p::Real)
+# eval_psp_pswfc_real(psp, i::Int, l::Int, p::Real)
+# eval_psp_pswfc_fourier(psp, i::Int, l::Int, p::Real)
+# count_n_pswfc(psp, l::Integer)
+# count_n_pswfc_radial(psp, l::Integer)
 
 """
     eval_psp_projector_real(psp, i, l, r)
@@ -203,3 +207,14 @@ function count_n_proj(psps, psp_positions)
     sum(count_n_proj(psp) * length(positions)
         for (psp, positions) in zip(psps, psp_positions))
 end
+
+count_n_pswfc_radial(psp::NormConservingPsp, l) = error("Pseudopotential $psp does not implement atomic wavefunctions.")
+
+function count_n_pswfc_radial(psp::NormConservingPsp)
+    sum(l -> count_n_pswfc_radial(psp, l), 0:psp.lmax; init=0)::Int 
+end
+
+count_n_pswfc(psp::NormConservingPsp, l) = count_n_pswfc_radial(psp, l) * (2l + 1)
+function count_n_pswfc(psp::NormConservingPsp)
+    sum(l -> count_n_pswfc(psp, l), 0:psp.lmax; init=0)::Int
+end

src/pseudo/PspUpf.jl

@@ -171,6 +171,8 @@ function eval_psp_projector_fourier(psp::PspUpf, i, l, p::T)::T where {T<:Real}
     hankel(rgrid, r2_proj, l, p)
 end
 
+count_n_pswfc_radial(psp::PspUpf, l) = length(psp.r2_pswfcs[l+1])
+
 function eval_psp_pswfc_real(psp::PspUpf, i, l, r::T)::T where {T<:Real}
     psp.r2_pswfcs_interp[l+1][i](r) / r^2
 end