Band structure for non-standard lattice, cleanup (#555)

Co-authored-by: Michael F. Herbst <[email protected]>

Project.toml

@@ -48,7 +48,7 @@ spglib_jll = "ac4a9f1e-bdb2-5204-990c-47c8b2f70d4e"
 [compat]
 AbstractFFTs = "1"
 AtomsBase = "0.2.2"
-Brillouin = "0.5"
+Brillouin = "0.5.4"
 ChainRulesCore = "1.15"
 Conda = "1"
 DftFunctionals = "0.2"

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
+Brillouin = "23470ee3-d0df-4052-8b1a-8cbd6363e7f0"
 DFTK = "acf6eb54-70d9-11e9-0013-234b7a5f5337"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

examples/graphene.jl

@@ -32,27 +32,7 @@ model = model_PBE(lattice, atoms, positions; temperature)
 basis = PlaneWaveBasis(model; Ecut, kgrid)
 scfres = self_consistent_field(basis)
 
-## Choose the points of the band diagram, in reduced coordinates (in the (b1,b2) basis)
-Γ  = [0, 0, 0]
-K  = [ 1, 1, 0]/3
-Kp = [-1, 2, 0]/3
-M  = (K + Kp)/2
-kpath_coords = [Γ, K, M, Γ]
-kpath_names  = ["Γ", "K", "M", "Γ"]
-
-## Build the path manually for now
-kline_density = 20
-function build_path(k1, k2)
-    target_Δk = 1/kline_density  # the actual Δk is |k2-k1|/npt
-    npt = ceil(Int, norm(model.recip_lattice * (k2-k1)) / target_Δk)
-    [k1 + t * (k2-k1) for t in range(0, 1, length=npt)]
-end
-kcoords = []
-for i = 1:length(kpath_coords)-1
-    append!(kcoords, build_path(kpath_coords[i], kpath_coords[i+1]))
-end
-klabels = Dict(zip(kpath_names, kpath_coords))
-
-## Plot the bands
-band_data = compute_bands(basis, kcoords; scfres.ρ)
-DFTK.plot_band_data(band_data; scfres.εF, klabels)
+## Construct 2D path through Brillouin zone
+sgnum = 13  # Graphene space group number
+kpath = irrfbz_path(model; dim=2, sgnum)
+plot_bandstructure(scfres, kpath; kline_density=20)

examples/publications/2022_cazalis.jl

@@ -59,30 +59,9 @@ basis = PlaneWaveBasis(model; Ecut, kgrid)
 ## Run SCF
 scfres = self_consistent_field(basis, tol=1e-10)
 
-## Choose the points of the band diagram, in reduced coordinates (in the (b1,b2) basis)
-Γ  = [0, 0, 0]
-K  = [ 1, 1, 0]/3
-Kp = [-1, 2, 0]/3
-M  = (K + Kp)/2
-kpath_coords = [Γ, K, M, Γ]
-kpath_names  = ["Γ", "K", "M", "Γ"]
-
-## Build the path
-kline_density = 20
-function build_path(k1, k2)
-    target_Δk = 1/kline_density  # the actual Δk is |k2-k1|/npt
-    npt = ceil(Int, norm(model.recip_lattice * (k2-k1)) / target_Δk)
-    [k1 + t * (k2-k1) for t in range(0, 1, length=npt)]
-end
-kcoords = []
-for i = 1:length(kpath_coords)-1
-    append!(kcoords, build_path(kpath_coords[i], kpath_coords[i+1]))
-end
-klabels = Dict(kpath_names[i] => kpath_coords[i] for i=1:length(kpath_coords))
-
-## Plot the bands
-band_data = compute_bands(basis, kcoords; n_bands=5, scfres.ρ)
-p = DFTK.plot_band_data(band_data; klabels, markersize=nothing)
+## Plot bands
+sgnum = 13  # Graphene space group number
+p = plot_bandstructure(scfres, irrfbz_path(model; sgnum); n_bands=5)
 Plots.hline!(p, [scfres.εF], label="", color="black")
-Plots.ylims!(p, -Inf,Inf)
+Plots.ylims!(p, (-Inf, Inf))
 p

src/DFTK.jl

@@ -177,8 +177,8 @@ include("external/pymatgen.jl")
 include("external/stubs.jl")  # Function stubs for conditionally defined methods
 
 export compute_bands
-export high_symmetry_kpath
 export plot_bandstructure
+export irrfbz_path
 include("postprocess/band_structure.jl")
 
 export compute_forces

src/common/zeros_like.jl

@@ -7,3 +7,4 @@ function zeros_like(X::AbstractArray, T::Type=eltype(X), dims::Integer...=size(X
 end
 zeros_like(X::AbstractArray, dims::Integer...) = zeros_like(X, eltype(X), dims...)
 zeros_like(X::Array, T::Type=eltype(X), dims::Integer...=size(X)...) = zeros(T, dims...)
+zeros_like(X::StaticArray, T::Type=eltype(X), dims::Integer...=size(X)...) = @SArray zeros(T, dims...)

src/external/spglib.jl

@@ -56,6 +56,9 @@ function spglib_cell(lattice, atom_groups, positions, magnetic_moments)
     spg = spglib_atoms(atom_groups, positions, magnetic_moments)
     (; cell=Spglib.Cell(lattice, spg.positions, spg.numbers, spg.spins), spg.collinear)
 end
+function spglib_cell(model::Model, magnetic_moments)
+    spglib_cell(model.lattice, model.atom_groups, model.positions, magnetic_moments)
+end
 
 
 @timing function spglib_get_symmetry(lattice::AbstractMatrix{<:AbstractFloat}, atom_groups,
@@ -180,6 +183,6 @@ end
 function spglib_spacegroup_number(model, magnetic_moments=[]; tol_symmetry=SYMMETRY_TOLERANCE)
     # Get spacegroup number according to International Tables for Crystallography (ITA)
     # TODO Time-reversal symmetry disabled? (not yet available in DFTK)
-    cell, _ = spglib_cell(model.lattice, model.atom_groups, model.positions, magnetic_moments)
+    cell, _ = spglib_cell(model, magnetic_moments)
     Spglib.get_spacegroup_number(cell, tol_symmetry)
 end

src/external/wannier90.jl

@@ -42,16 +42,16 @@ function write_w90_win(fileprefix::String, basis::PlaneWaveBasis;
         println(fp, "!"^20 * " k_points\n")
 
         if bands_plot
-            kpathdata = high_symmetry_kpath(basis.model)
-            length(kpathdata.kpath) > 1 || @warn(
+            kpath  = irrfbz_path(basis.model)
+            length(kpath.paths) > 1 || @warn(
                 "Only first kpath branch considered in write_w90_win")
-            path = kpathdata.kpath[1]
+            path = kpath.paths[1]
 
             println(fp, "begin kpoint_path")
             for i in 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.(kpathdata.klabels[A], digits=5)...)
-                @printf(fp, "%s %10.6f %10.6f %10.6f\n", B, round.(kpathdata.klabels[B], digits=5)...)
+                @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)...)
             end
             println(fp, "end kpoint_path")
             println(fp, "bands_plot = true\n")

src/plotting.jl

@@ -21,38 +21,34 @@ function ScfPlotTrace(plt=Plots.plot(yaxis=:log); kwargs...)
 end
 
 
-function plot_band_data(band_data; εF=nothing,
-                        klabels=Dict{String, Vector{Float64}}(), unit=u"hartree", kwargs...)
+function plot_band_data(kpath::KPathInterpolant, band_data;
+                        εF=nothing, unit=u"hartree", kwargs...)
     eshift = isnothing(εF) ? 0.0 : εF
-    data = prepare_band_data(band_data, klabels=klabels)
+    data = data_for_plotting(kpath, band_data)
 
     # Constant to convert from AU to the desired unit
     to_unit = ustrip(auconvert(unit, 1.0))
 
-    markerargs = ()
-    if !(:markersize in keys(kwargs)) && !(:markershape in keys(kwargs))
-        if length(krange_spin(band_data.basis, 1)) < 70
-            markerargs = (markersize=2, markershape=:circle)
+    # Plot all bands, spins and errors
+    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
+        yerror = nothing
+        if hasproperty(data, :λerror)
+            yerror = data.λerror[:, iband, σ][branch] .* to_unit
         end
+        energies = (data.λ[:, iband, σ][branch] .- eshift) .* to_unit
+        Plots.plot!(p, data.kdistances[branch], energies;
+                    label="", yerror, color=(:blue, :red)[σ], margs..., kwargs...)
     end
 
-    # For each branch, plot all bands, spins and errors
-    p = Plots.plot(xlabel="wave vector")
-    for branch in data.branches
-        for σ in 1:data.n_spin, iband = 1:data.n_bands
-            yerror = nothing
-            if hasproperty(branch, :λerror)
-                yerror = branch.λerror[:, iband, σ] .* to_unit
-            end
-            energies = (branch.λ[:, iband, σ] .- eshift) .* to_unit
-            color = (:blue, :red)[σ]
-            Plots.plot!(p, branch.kdistances, energies; color, label="",
-                        yerror, markerargs..., kwargs...)
-        end
+    # Delimiter for branches
+    for branch in data.kbranches[1:end-1]
+        Plots.vline!(p, [data.kdistances[last(branch)]], color=:black, label="")
     end
 
     # X-range: 0 to last kdistance value
-    Plots.xlims!(p, (0, data.branches[end].kdistances[end]))
+    Plots.xlims!(p, (0, data.kdistances[end]))
     Plots.xticks!(p, data.ticks.distances, data.ticks.labels)
 
     ylims = [-0.147, 0.147]