Implemented the electron density states postproc tool.

postproc/dos.jl

@@ -22,7 +22,7 @@ function dos(ω, εs, vs, Bs, wvmesh)
 
     #Number of wave vectors and bands
     nq, nb = size(εs)
-
+    
     #Hard-code a minimum smearing value.
     σ_min = 1.0e-6 # 1 μeV
 
@@ -44,16 +44,93 @@ function dos(ω, εs, vs, Bs, wvmesh)
             dos += δ_gaussian(ω, εs[iq′, ib′], σ)
         end
     end
-    return dos/nq
+    return dos/prod(wvmesh)
+end
+
+function calculate_dos(species, rundir, outdir, chempot)
+    #Calculate the density of states of a given 
+    #species ∈ ["ph", "el"].
+
+    #Read full-Brillouin zone energies and velocities. Reshape the latter appropriately.
+    println("ϟ Reading full-Brillouin zone energies and velocities...")
+    εs = readdlm(rundir*species*".ens_fbz") #eV
+    vs = readdlm(rundir*species*".vels_fbz") #Km/s
+    vs_shape = size(vs)
+    vs = reshape(vs, (vs_shape[1], vs_shape[2]÷3, 3))
+
+    #Read irreducible Brillouin zone energies and flatten.
+    println("ϟ Reading irreducible Brillouin zone energies...")
+    εs_irred = readdlm(rundir*species*".ens_ibz") #eV
+    εs_irred_shape = size(εs_irred)
+    ωs = reshape(εs_irred, (prod(εs_irred_shape)))
+
+    #Read reciprocal lattice vectors and their discretization grid
+    println("ϟ Reading reciprocal lattice vector data...")
+    reclattvecs_data = readdlm(rundir*species*".reclattvecs") #nm^-1
+    Bs = reclattvecs_data[1:3, :]
+    wvmesh = convert.(Int32, reclattvecs_data[4, :])
+
+    #Calculate band resolved DOS
+    println("ϟ Calculating density of states...")
+    species_dos = reshape(ThreadsX.map(ω -> dos(ω, εs, vs, Bs, wvmesh), ωs), εs_irred_shape)
+    if(species == "ph")
+        species_dos[1, 1:3] .= 0.0 #Take care of the 3 Γ-point acoustic phonons
+    end
+
+    #Write DOS to file
+    println("ϟ Writing results out to file...")
+    open(outdir*species*".dos_branches", "w") do w
+        writedlm(w, species_dos)
+    end
+
+    #Plot in meV^-1-mev units
+    println("ϟ Plotting results...")
+    if(species == "ph")
+        dos_plot = plot(1.0e3*εs_irred, 1.0e-3*species_dos,
+                        legend = false, seriestype = :scatter, mc = :green, ma = 0.6,
+                        xlabel = "phonon energy [meV]",
+                        ylabel = "DOS [meV]"*L"$^{-1}$")
+    else
+        #Set hald of Fermi shell width
+        half_Fermi_shell = 0.2 #eV
+
+        #Find (closest) Cartesian index to the chemical potential
+        chempot_index = argmin(abs.(εs_irred .- chempot))
+
+        #Set vertical range
+        padding = 1.25 #extra vertical space
+        ymin = (2.0 - padding)*species_dos[chempot_index]
+        ymax = padding*species_dos[chempot_index]
+
+        dos_plot = plot(εs_irred .- chempot, species_dos,
+                        xlims = [-half_Fermi_shell, half_Fermi_shell],
+                        ylims = [ymin, ymax],
+                        legend = false, seriestype = :scatter, mc = :green, ma = 0.6,
+                        xlabel = "electron energy - "*L"\mu"*" [eV]",
+                        ylabel = "DOS [eV.spin]"*L"$^{-1}$")
+    end
+    savefig(dos_plot, outdir*species*"dos_bands.pdf")
+
+    #display(dos_plot)
 end
 
 function parse_commandline()
     args = ArgParseSettings()
     @add_arg_table args begin
         "--rundir"
-            help = "elphbolt run directory"
-            arg_type = String
-            default = "./"
+        help = "elphbolt run directory"
+        arg_type = String
+        default = "./"
+
+        "--dosof"
+        help = "(ph)onon or (el)ectron density of states"
+        arg_type = String
+        default = "ph"
+
+        "--chempot"
+        help = "Chemical potential"
+        arg_type = Float64
+        default = 0.0
     end
     return parse_args(args)
 end
@@ -73,45 +150,10 @@ function main()
     outdir = rundir*"postproc_results/"
     mkpath(outdir)
 
-    #Read full-Brillouin zone energies and velocities. Reshape the latter appropriately.
-    println("ϟ Reading full-Brillouin zone energies and velocities...")
-    εs = readdlm(rundir*"ph.ens_fbz") #eV
-    vs = readdlm(rundir*"ph.vels_fbz") #Km/s
-    vs_shape = size(vs)
-    vs = reshape(vs, (vs_shape[1], vs_shape[2]÷3, 3))
-
-    #Read irreducible Brillouin zone energies and flatten.
-    println("ϟ Reading irreducible Brillouin zone energies...")
-    εs_irred = readdlm(rundir*"ph.ens_ibz") #eV
-    εs_irred_shape = size(εs_irred)
-    ωs = reshape(εs_irred, (prod(size(εs_irred))))
+    #Set chemical potential
+    chempot = parsed_args["chempot"]
 
-    #Read reciprocal lattice vectors and their discretization grid
-    println("ϟ Reading reciprocal lattice vector data...")
-    reclattvecs_data = readdlm(rundir*"ph.reclattvecs") #nm^-1
-    Bs = reclattvecs_data[1:3, :]
-    wvmesh = reclattvecs_data[4, :]
-
-    #Calculate band resolved phonon DOS
-    println("ϟ Calculating phonon density of states...")
-    ph_dos = reshape(ThreadsX.map(ω -> dos(ω, εs, vs, Bs, wvmesh), ωs), εs_irred_shape)
-    ph_dos[1, 1:3] .= 0.0 #Take care of the 3 Γ-point acoustic phonons
-
-    #Write DOS to file
-    println("ϟ Writing results out to file...")
-    open(outdir*"ph.dos_branches", "w") do w
-        writedlm(w, ph_dos)
-    end
-
-    #Plot in meV^-1-mev units
-    println("ϟ Plotting results...")
-    dos_plot = plot(1.0e3*εs_irred, 1.0e-3*ph_dos,
-                    legend = false, seriestype = :scatter, mc = :green, ma = 0.6,
-                    xlabel = "phonon energy [meV]",
-                    ylabel = "DOS [meV]"*L"$^{-1}$")
-    savefig(dos_plot, outdir*"phdos_bands.pdf")
-
-    #display(dos_plot)
+    calculate_dos(parsed_args["dosof"], rundir, outdir, chempot)
 
     println("ϟ All done!")
 end

src/electron.f90

@@ -281,7 +281,7 @@ subroutine read_input_and_setup(self, wann, crys, sym, num)
        if (scissor .ne. 0.0_r64) then
           write(*, "(A, 1E16.8, A)") "Scissor operator = ", &
             self%scissor(self%indlowconduction) , " eV"
-       end if     
+       end if
     end if
 
     !Calculate electrons
@@ -573,6 +573,10 @@ subroutine calculate_electrons(self, wann, crys, sym, num)
     !Print out irreducible electron energies and velocities
     call write2file_rank2_real("el.ens_ibz", self%ens_irred)
     call write2file_rank3_real("el.vels_ibz", self%vels_irred)
+
+    !Print out full BZ electron energies and velocities
+    call write2file_rank2_real("el.ens_fbz", self%ens)
+    call write2file_rank3_real("el.vels_fbz", self%vels)
 
     !Calculate electron tetrahedra
     if(num%tetrahedra) then

src/numerics.f90

@@ -362,6 +362,19 @@ subroutine read_input_and_setup(self, crys)
        end do
        write(1, *) self%qmesh
        close(1)
+
+       open(1, file = "el.reclattvecs", status = "replace")
+       do i = 1, 3
+          write(1, "(" // trim(adjustl(numcols)) // "E20.10)") &
+               crys%reclattvecs(:, i)
+       end do
+       if(crys%twod) then
+          write(1, *) self%mesh_ref*self%qmesh(1), self%mesh_ref*self%qmesh(2), 1
+       else
+          write(1, *) self%mesh_ref*self%qmesh(1), self%mesh_ref*self%qmesh(2), &
+               self%mesh_ref*self%qmesh(3)
+       end if
+       close(1)
 
        write(*, "(A, (3I5,x))") "q-mesh = ", self%qmesh
        if(crys%twod) then

src/phonon.f90

@@ -23,7 +23,7 @@ module phonon_module
   use particle_module, only: particle
   use misc, only: print_message, subtitle, expi, distribute_points, &
        write2file_rank2_real, exit_with_message, create_set, coarse_grain, &
-       mux_state
+       mux_state, write2file_rank3_real
   use numerics_module, only: numerics
   use crystal_module, only: crystal, calculate_wavevectors_full
   use symmetry_module, only: symmetry, find_irred_wedge, create_fbz2ibz_map
@@ -326,23 +326,9 @@ subroutine calculate_phonons(self, crys, sym, num, wann)
     end if
 
     !Print out full BZ phonon energies and velocities
+    call write2file_rank2_real("ph.ens_fbz", self%ens)
+    call write2file_rank3_real("ph.vels_fbz", self%vels)
     if(this_image() == 1) then
-       write(numcols, "(I0)") self%numbands
-       open(1, file = "ph.ens_fbz", status = "replace")
-       do iq = 1, self%nwv
-          write(1, "(" // trim(adjustl(numcols)) // "E20.10)") &
-               self%ens(iq, :)
-       end do
-       close(1)
-
-       write(numcols, "(I0)") 3*self%numbands
-       open(1, file = "ph.vels_fbz", status = "replace")
-       do iq = 1, self%nwv
-          write(1, "(" // trim(adjustl(numcols)) // "E20.10)") &
-               self%vels(iq, :, :)
-       end do
-       close(1)
-
        open(1, file = "ph.fbz2ibz_map", status = "replace")
        do iq = 1, self%nwv
           write(1, "(I10)") fbz2ibz_map(iq)