From 9e1067867cbccd874a6ee02171c4d60d17bd43c2 Mon Sep 17 00:00:00 2001
From: Kun Chen <chenkun0228@gmail.com>
Date: Fri, 11 Nov 2022 23:02:51 -0500
Subject: [PATCH] map DFTK band to BZMesh

---
 example/sodium/sodium.jl | 59 +++++++++++++++-------------------------
 1 file changed, 22 insertions(+), 37 deletions(-)

diff --git a/example/sodium/sodium.jl b/example/sodium/sodium.jl
index 8138cd0..1fd82c1 100644
--- a/example/sodium/sodium.jl
+++ b/example/sodium/sodium.jl
@@ -4,6 +4,8 @@ using UnitfulAtomic
 using Plots
 using PlotlyJS
 using Printf
+using PyCall
+using BrillouinZoneMeshes
 
 # Setup a sodium lattice, see the following link for more information
 # http://lampx.tugraz.at/~hadley/ss1/crystalstructure/structures/bcc/bcc.php
@@ -38,46 +40,29 @@ _, data = compute_bands(basis, [[0.0, 0.0, 0.0],], ρ=scfres.ρ);
 bandwidth = uconvert.(u"eV", (data[1][5] - scfres.εF) * u"hartree")
 println(bandwidth)
 
-# analysis the wavefunction
-# the 1s band is the 2s^2 orbital
-# the 2-4 bands are the 2p^6 orbitals
-# the >5 bands are the valence electrons
+# # interpolation of bands
+# interp = pyimport("scipy.interpolate")
+# np = pyimport("numpy")
+# griddata = interp.griddata
 
-kidx = 1
+# build BZ mesh and maps
+cell = Cells.Cell(; lattice=austrip.(lattice), atoms=[1,], positions)
+bzmesh = BZMeshes.Monkhorst_Pack(cell=cell, size=kgrid, shift=[false, false, false])
+meshmap = MeshMaps.MeshMap(bzmesh)
 
-for (gi, g) in enumerate(scfres.basis.kpoints[kidx].G_vectors)
-    ψ=scfres.ψ[kidx]
-    @printf("%32s    %16.8f    %16.8f    %16.8f\n", "$g", abs(ψ[gi, 1]), abs(ψ[gi, 2]), abs(ψ[gi, 3])) 
+println(length(meshmap.irreducible_indices))
+for idx in meshmap.irreducible_indices
+    k = inv_lattice_vector(bzmesh) * bzmesh[idx]
+    println(k)
 end
 
-Gz = 0
-bandidx = 1 
-
-x, y, psi = [], [], []
-
-for (gi, g) in enumerate(scfres.basis.kpoints[kidx].G_vectors)
-    if g[3]==Gz
-        push!(x, g[1])
-        push!(y, g[2])
-        push!(psi, abs(scfres.ψ[kidx][gi, bandidx]))
-    end
-end
-Gxmax = maximum(abs.(x))
-Gymax = maximum(abs.(y))
-
-z = zeros(Gxmax*2+1, Gymax*2+1)
-for (gi, g) in enumerate(scfres.basis.kpoints[kidx].G_vectors)
-    if g[3]==Gz
-        z[g[1]+Gxmax+1, g[2]+Gymax+1] =abs(scfres.ψ[kidx][gi, bandidx])
-    end
+# band structure with BZ mesh and maps
+band = Dict{Int,Float64}()
+bandix = 5 # the fifth band is the valence band
+for (ki, kpoints) in enumerate(basis.kpoints)
+    idx = AbstractMeshes.locate(bzmesh, basis.kpoints[ki].coordinate)
+    band[idx] = scfres.eigenvalues[ki][bandix] - scfres.εF
 end
 
-data = PlotlyJS.contour(x=[i for i in -Gxmax:Gxmax], y=[i for i in -Gymax:Gymax], z=z)
-PlotlyJS.plot(data)
-# pyplot()
-# Plots.contour(x, y, psi, xlabel="Gx", ylabel="Gy", zlabel="|ψ|", legend=false)
-# using PlotlyJS
-# trace = PlotlyJS.surface(x=x, y=y, z=psi)
-# layout = Layout(title="Mt. Bruno Elevation", autosize=false, width=500,
-#                     height=500)
-# PlotlyJS.plot(trace, layout)
\ No newline at end of file
+# band fourier transform
+