From 225622a4a61489d7828ee72ef78e1085921af26e Mon Sep 17 00:00:00 2001 From: Nakib Haider Date: Wed, 27 Sep 2023 16:24:15 +0200 Subject: [PATCH] Implemented the electron density states postproc tool. --- postproc/dos.jl | 128 +++++++++++++++++++++++++++++++---------------- src/electron.f90 | 6 ++- src/numerics.f90 | 13 +++++ src/phonon.f90 | 20 ++------ 4 files changed, 106 insertions(+), 61 deletions(-) diff --git a/postproc/dos.jl b/postproc/dos.jl index 7c4d1fea..5f0a1d2d 100644 --- a/postproc/dos.jl +++ b/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 diff --git a/src/electron.f90 b/src/electron.f90 index a4fe06b0..474943c5 100644 --- a/src/electron.f90 +++ b/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 diff --git a/src/numerics.f90 b/src/numerics.f90 index 1d6be472..83291df9 100644 --- a/src/numerics.f90 +++ b/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 diff --git a/src/phonon.f90 b/src/phonon.f90 index 8f47b538..fac303df 100644 --- a/src/phonon.f90 +++ b/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)