Support band structures for 1D systems (JuliaMolSim#492)

.github/workflows/documentation.yaml

@@ -41,7 +41,7 @@ jobs:
             doctest(DFTK)'
       - name: Generate docs
         run: julia --project=docs docs/make.jl
-        timeout-minutes: 20
+        timeout-minutes: 30
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

examples/free_electron.jl

@@ -0,0 +1,14 @@
+using DFTK
+using Plots
+
+# This example plots the free-electron bands for a 1D system.
+# TODO Convert into an example in the docs
+#      show how one could add a potential to this and see the bands change.
+
+lattice = diagm([5., 0, 0])
+model   = Model(lattice; n_electrons=4, terms=[Kinetic()])
+basis   = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1));
+
+n_bands = 6
+ρ0 = guess_density(basis)  # Just dummy, has no meaning in this model
+p = plot_bandstructure(basis, ρ0, n_bands, kline_density=15)

src/postprocess/band_structure.jl

@@ -3,6 +3,17 @@ using PyCall
 # Functionality for computing band structures
 
 function high_symmetry_kpath(model; kline_density=20)
+    if model.n_dim == 1  # Return fast for 1D model
+        # kline_density = Number of kpoint per Angström
+        # Length of the kpath is lattice[1, 1] in 1D
+        n_points = ceil(Int, kline_density / austrip(1u"Å") * model.lattice[1, 1])
+        return (
+            kcoords = [[coord, 0, 0] for coord in range(-1//2, 1//2, length=1+n_points)],
+            klabels=Dict("Γ" => zeros(3), "-1/2" => [-0.5, 0.0, 0.0], "1/2" => [0.5, 0, 0]),
+            kpath=[[raw"-1/2", raw"1/2"]],
+        )
+    end
+
     # TODO This is the last function that hard-depends on pymatgen. The way to solve this
     # is to use the julia version implemented in
     # https://github.com/louisponet/DFControl.jl/blob/master/src/structure.jl

src/pseudo/PspHgh.jl

@@ -182,7 +182,7 @@ and `Q` is a polynomial. This function returns `Q`.
     @assert 0 <= l <= length(psp.rp) - 1
     @assert i > 0
     rp::T = psp.rp[l + 1]
-    common::T = 4T(π)^(5 / T(4)) * sqrt(T(2)^(l + 1) * rp^3)
+    common::T = 4T(π)^(5 / T(4)) * sqrt(T(2^(l + 1)) * rp^3)
 
     # Note: In the (l == 0 && i == 2) case the HGH paper has an error.
     #       The first 8 in equation (8) should not be under the sqrt-sign

test/compute_bands.jl

@@ -130,6 +130,18 @@ if mpi_nprocs() == 1  # not easy to distribute
     @test kpath[2] == ["U", "X"]
 end
 
+@testset "High-symmetry kpath construction for 1D system" begin
+    lattice = diagm([1.0, 0, 0])
+    model = Model(lattice, n_electrons=1, terms=[Kinetic()])
+    kcoords, klabels, kpath = high_symmetry_kpath(model; kline_density=10)
+
+    @test length(kcoords) == 7
+    @test kcoords[1] ≈ [-1/2, 0, 0]
+    @test kcoords[4] ≈ [   0, 0, 0]
+    @test kcoords[7] ≈ [ 1/2, 0, 0]
+    @test length(kpath) == 1
+end
+
 @testset "Compute bands for silicon" begin
     testcase = silicon
     Ecut = 7