From 075b1e5761960ab581976b3c10ec1e1e41a955b0 Mon Sep 17 00:00:00 2001
From: Nakib Haider <nakib.haider.protik@gmail.com>
Date: Wed, 27 Mar 2024 16:48:42 +0100
Subject: [PATCH] Worked a bit on the spectral transport coeff julia script.

---
 postproc/spectral_trans_coeffs.jl | 72 ++++++++++++++++++++++++++-----
 1 file changed, 62 insertions(+), 10 deletions(-)

diff --git a/postproc/spectral_trans_coeffs.jl b/postproc/spectral_trans_coeffs.jl
index 2c26576..00602a6 100644
--- a/postproc/spectral_trans_coeffs.jl
+++ b/postproc/spectral_trans_coeffs.jl
@@ -1,9 +1,14 @@
 include("parameters.jl")
 include("statistics.jl")
 
+using ThreadsX
 using ArgParse
 using DelimitedFiles
 
+function Heaviside(x::Float64)
+    (1.0 + sign(x))/2.0
+end
+
 function calculate_mfp_cumulative_kappa(rundir, outdir, T)
     # TODO ...
 
@@ -15,27 +20,74 @@ function calculate_mfp_cumulative_kappa(rundir, outdir, T)
     
     #TODO Generate T-dependent directory from T
     Tdir = "T0.300E+03/"
+
+    #Read map from full Brillouin zone to the irreducible irreducible one.
+    println("ϟ Reading FBZ->IBZ map...")
+    fbz2ibz = readdlm(rundir*species*".fbz2ibz_map", Int32)
     
-    #Read full-Brillouin zone energies and velocities. Reshape the latter appropriately.
-    println("ϟ Reading full-Brillouin zone energies...")
+    #Read full Brillouin zone energies and velocities. Reshape the latter appropriately.
+    println("ϟ Reading FBZ 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))
+
     εs_shape = size(εs)
     nq_fbz, nb = εs_shape[1], εs_shape[2] #number of FBZ points, bands
     
-    #Read full-Brillouin zone energies and velocities. Reshape the latter appropriately.
+    #Read full Brillouin zone energies.
     println("ϟ Reading full-Brillouin zone mean-free-paths (mfps)...")
     λs = readdlm(rundir*Tdir*"nodrag_iterated_ph_mfps_ibz") #nm
 
-    #Create mfp sampling grid (log scale)
-    low = -6 #1e-6 nm, a really tiny number I'd say.
-    high = log10(1.5*maximum(λs)) #50% higher than largest value in λs
+    #Read full Brillouin zone reponse functions.
 
-    N = 25 # TODO should read this from user input
-    #λ' = exp10.(range(low, high, length = N))
-    λp = exp10.(range(low, high, length = N))
+    response = zeros(nb, nq_fbz, 3)
+    for ib in 1:nb
+        fname = "nodrag_F0_"*string(ib)
+        response[ib, :, :] = readdlm(rundir*Tdir*fname)
+    end
     
+    #Create mfp sampling grid (log scale)
+    #low = -6 
+    #high = log10(1.5*maximum(λs)) 
+    nsamp = 25 # TODO should read this from user input
+    λsamp = exp10.(range(-6, #1e-6 nm, a really tiny number I'd say.
+                         log10(1.5*maximum(λs)) , #50% higher than largest value in λs
+                         nsamp))
+
+    #Calculate band resolved DOS
+    println("ϟ Calculating phonon thermal conductivity accumulation wrt mean-free-path...")
+    
+    #kappa_accum = reshape(ThreadsX.map(l -> mfp_cumulative_kappa(l, εs), λsamp), [nsamp, nb])
+    
+    kappa_accum = zeros(nsamp, nb, 3, 3)
+    for isamp in 1:nsamp
+        for ib in 1:nb
+            for iq_fbz in 1:nq_fbz
+                iq_ibz = fbz2ibz[iq_fbz]
+                if(λs[iq_ibz, ib] <= λsamp[isamp])    
+                    ε = εs[iq_fbz, ib]
+
+                    population_fac = Bose(ε, T)
+                    population_fac *= 1.0 + population_fac
+
+                    for icart in 1:3
+                        kappa_accum[isamp, ib, icart, :] +=
+                            ε*population_fac*vs[iq_fbz, ib, icart]*response[ib, iq_fbz, :]
+                    end
+                end
+            end
+        end
+    end
     #Common multiplicative factor
-    fac = 1.0e21/kB/T/volume/nq_fbz
+    fac = qe*1.0e21/kB/T/volume/nq_fbz
+    kappa_accum *= fac
+    
+    #Write to file
+    #println("ϟ Writing results out to file...")
+    #open(outdir*species*".kappa_accum_wrt_mfp", "w") do w
+    #    writedlm(w, kappa_accum)
+    #end
 end
 
 function parse_commandline()