From 761629db57b9ad98c4f47c785136493c1ec5ba0c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Thu, 13 Mar 2025 10:15:13 +0100
Subject: [PATCH 01/14] WIP broadening + viscosity

---
 src/thermal_conductivity/kappa.f90            | 59 +++++++++++++++-
 src/thermal_conductivity/scattering.f90       | 69 +++++++++++++++++--
 .../scattering_threephonon.f90                | 16 ++++-
 3 files changed, 137 insertions(+), 7 deletions(-)

diff --git a/src/thermal_conductivity/kappa.f90 b/src/thermal_conductivity/kappa.f90
index bc324c4a..d3f44324 100644
--- a/src/thermal_conductivity/kappa.f90
+++ b/src/thermal_conductivity/kappa.f90
@@ -1,7 +1,7 @@
 #include "precompilerdefinitions"
 module kappa
 use konstanter, only: r8, lo_sqtol, lo_kb_hartree, lo_freqtol, lo_kappa_au_to_SI, &
-                      lo_groupvel_Hartreebohr_to_ms
+                      lo_groupvel_Hartreebohr_to_ms, lo_twopi
 use gottochblandat, only: lo_sqnorm, lo_planck, lo_outerproduct, lo_chop, lo_harmonic_oscillator_cv
 use mpi_wrappers, only: lo_mpi_helper
 use lo_memtracker, only: lo_mem_helper
@@ -22,6 +22,7 @@ module kappa
 public :: get_kappa_offdiag
 public :: iterative_solution
 public :: symmetrize_kappa
+public :: get_viscosity
 contains
 
 !> Calculate the thermal conductivity
@@ -484,4 +485,60 @@ subroutine iterative_solution(sr, dr, qp, uc, temperature, niter, tol, classical
     call mem%deallocate(Fbb, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(buf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 end subroutine
+
+!> Calculate the thermal conductivity
+subroutine get_viscosity(dr, qp, uc, temperature, classical, mu)
+    !> dispersions
+    type(lo_phonon_dispersions), intent(inout) :: dr
+    !> q-mesh
+    class(lo_qpoint_mesh), intent(in) :: qp
+    !> structure
+    type(lo_crystalstructure), intent(in) :: uc
+    !> temperature
+    real(r8), intent(in) :: temperature
+    !> Are we in the classical limit ?
+    logical, intent(in) :: classical
+    !> thermal conductivity tensor
+    real(r8), dimension(3, 3, 3, 3), intent(out) :: mu
+
+    real(r8), dimension(3) :: v0, v1, kr
+    real(r8) :: om1, cv, pref, n
+    integer :: j
+
+    integer :: q1, b1, iq, a, b, c, d, iop
+    real(r8), dimension(3, 3) :: v2, buf
+
+    mu = 0.0_r8
+    do q1=1, qp%n_full_point
+        iq = qp%ap(q1)%irreducible_index
+        kr = qp%ap(q1)%r * lo_twopi
+        iop = qp%ap(q1)%operation_from_irreducible
+        do b1=1, dr%n_mode
+            om1 = dr%aq(q1)%omega(b1)
+            if (om1 .lt. lo_freqtol) cycle
+            n = lo_planck(temperature, om1)
+            pref = n * (n + 1.0_r8) / uc%volume / lo_kb_hartree / temperature
+            v2 = 0.0_r8
+            if (iop .gt. 0) then
+                v0 = lo_operate_on_vector(uc%sym%op(iop), dr%iq(iq)%Fn(:, b1))
+                v1 = lo_operate_on_vector(uc%sym%op(iop), dr%iq(iq)%vel(:, b1))
+            else
+                v0 = -lo_operate_on_vector(uc%sym%op(abs(iop)), dr%iq(iq)%Fn(:, b1))
+                v1 = -lo_operate_on_vector(uc%sym%op(abs(iop)), dr%iq(iq)%vel(:, b1))
+            end if
+            v2 = lo_outerproduct(v0, v1)
+            do a=1, 3
+            do b=1, 3
+            do c=1, 3
+            do d=1, 3
+                mu(a, b, c, d) = mu(a, b, c, d) + pref * kr(a) * v0(b) * kr(c) * v1(d) / qp%n_full_point
+!               mu(a, b, c, d) = mu(a, b, c, d) + pref * v2(a, b) * kr(c) * kr(d) / qp%n_full_point
+            end do
+            end do
+            end do
+            end do
+        end do
+    end do
+    mu = lo_chop(mu, sum(abs(mu))*1e-6_r8)
+end subroutine
 end module
diff --git a/src/thermal_conductivity/scattering.f90 b/src/thermal_conductivity/scattering.f90
index 45f1ad1b..7cc4aa4c 100644
--- a/src/thermal_conductivity/scattering.f90
+++ b/src/thermal_conductivity/scattering.f90
@@ -1,7 +1,8 @@
 #include "precompilerdefinitions"
 module scattering
 use konstanter, only: r8, i8, lo_freqtol, lo_twopi, lo_exitcode_param, lo_hugeint, lo_pi, lo_tol, &
-                      lo_phonongroupveltol, lo_tol, lo_frequency_THz_to_Hartree, lo_kb_hartree, lo_huge
+                      lo_phonongroupveltol, lo_tol, lo_frequency_THz_to_Hartree, lo_kb_hartree, lo_huge, &
+                      lo_frequency_Hartree_to_meV
 use gottochblandat, only: walltime, lo_trueNtimes, lo_progressbar_init, lo_progressbar, lo_gauss, lo_planck, lo_return_unique
 use mpi_wrappers, only: lo_mpi_helper, lo_stop_gracefully
 use lo_memtracker, only: lo_mem_helper
@@ -83,8 +84,13 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
         !> The q-point grid dimension
         integer, dimension(3) :: dims
         !> Some integers for the do loop/indices
-        integer :: q1, b1, il, j, my_nqpoints, ctr
-        !> The seed for the random number generator for the Monte-Carlo integration
+        integer :: q1, q2, b1, il, j, my_nqpoints, ctr
+
+        real(r8), dimension(:), allocatable :: sigavg, bla
+        real(r8), dimension(3, 3) :: reclat
+        real(r8), dimension(6) :: w
+        real(r8) :: sigma
+        integer, dimension(3) :: fq, fqp
 
         ! grid dimensions
         select type (qp)
@@ -116,8 +122,15 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
         ! We can start some precomputation
         allocate (sr%be(qp%n_irr_point, dr%n_mode))
         allocate (sr%sigsq(qp%n_irr_point, dr%n_mode))
+        call mem%allocate(sigavg, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+        call mem%allocate(bla, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
         sr%be = 0.0_r8
         sr%sigsq = 0.0_r8
+        bla = 0.0_r8
+        sigavg = 0.0_r8
+        do j=1, 3
+            reclat(:, j) = uc%reciprocal_latticevectors(:, j) / dims(j)
+        end do
         do q1 = 1, qp%n_irr_point
             do b1 = 1, dr%n_mode
                 if (opts%classical) then
@@ -126,9 +139,57 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
                     sr%be(q1, b1) = lo_planck(opts%temperature, dr%iq(q1)%omega(b1))
                 end if
 
-                sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), dr%default_smearing(b1), opts%sigma)**2
+                do j=1, 3
+                    w(j) = dot_product(dr%iq(q1)%vel(:, b1), reclat(:, j))**2
+                end do
+!               sigma = sum(w(1:3))
+                sigma = maxval(w(1:3))
+                sr%sigsq(q1, b1) = sigma
+                if (sigma .gt. 0) then
+                    sigavg(b1) = sigavg(b1) + sigma * qp%ip(q1)%integration_weight
+                    bla(b1) = bla(b1) + 1.0_r8
+                end if
+
+                fq = singlet_to_triplet(qp%ip(q1)%full_index, dims(2), dims(3))
+                do j=1, 3
+                    fqp = fq
+                    fqp(j) = mod(fqp(j) + 1, dims(j))
+                    if (fqp(j) .eq. 0) fqp(j) = dims(j)
+                    q2 = triplet_to_singlet(fqp, dims(2), dims(3))
+!                   w(j) = dot_product(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b1), reclat(:, j))**2
+!                   w(j) = sum((dr%iq(q1)%vel(j, b1) * reclat(:, j))**2)
+                    w(j) = 0.5_r8 * (dr%iq(q1)%omega(b1) - dr%aq(q2)%omega(b1))
+
+                    fqp = fq
+                    fqp(j) = fqp(j) - 1
+                    if (fqp(j) .eq. 0) fqp(j) = dims(j)
+                    q2 = triplet_to_singlet(fqp, dims(2), dims(3))
+!                   w(j+3) = dot_product(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b1), reclat(:, j))**2
+!                   w(j+3) = sum((dr%iq(q1)%vel(j, b1) * reclat(:, j))**2)
+                    w(j+3) = 0.5_r8 * (dr%iq(q1)%omega(b1) - dr%aq(q2)%omega(b1))
+
+                end do
+!               sigma = (sum(w * lo_twopi / sqrt(2.0_r8)) / 6.0_r8)**2
+                sigma = sum(w)**2 / 6.0_r8
+!               sigma = norm2(dr%iq(q1)%vel(:, b1)) * qp%ip(q1)%radius
+!               sigma = sigma**2
+!               sigma = maxval(w(1:3))
+                sigavg(b1) = sigavg(b1) + sigma * qp%ip(q1)%integration_weight
+                sr%sigsq(q1, b1) = sigma
+
+!               sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), minval(dr%default_smearing), opts%sigma)**2
+!               sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), dr%default_smearing(b1), opts%sigma)**2
+            end do
+        end do
+        sigavg = sigavg / bla
+        do q1=1, qp%n_irr_point
+            do b1=1, dr%n_mode
+                if (sqrt(sr%sigsq(q1, b1)) .lt. lo_freqtol) sr%sigsq(q1, b1) = sigavg(b1)
             end do
         end do
+        if (mw%talk) write(*, *) dr%default_smearing * lo_frequency_Hartree_to_meV
+        if (mw%talk) write(*, *) sqrt(sigavg) * lo_frequency_Hartree_to_meV
+        call mem%deallocate(sigavg, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 
         if (mw%talk) write (*, *) '... distributing q-point/modes on MPI ranks'
         ctr = 0
diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90
index f164acdc..e0850542 100644
--- a/src/thermal_conductivity/scattering_threephonon.f90
+++ b/src/thermal_conductivity/scattering_threephonon.f90
@@ -49,6 +49,13 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     !> buff to keep the off diagonal scattering matrix elements
     real(r8), dimension(:, :), allocatable :: od_terms
 
+    real(r8), dimension(3, 3) :: reclat
+    real(r8), dimension(3) :: w
+
+    do i=1, 3
+        reclat(:, i) = uc%reciprocal_latticevectors(:, i) / mcg%full_dims(i)
+    end do
+
     ! We start by allocating everything
     call mem%allocate(ptf, dr%n_mode**3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%allocate(evp1, dr%n_mode**2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
@@ -107,8 +114,13 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
                                  sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
                                  sr%sigsq(qp%ap(q3)%irreducible_index, b3))
                 case (6)
-                    sigma = qp%smearingparameter(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), &
-                                                 dr%default_smearing(b3), smearing)
+!                   sigma = qp%smearingparameter(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), &
+!                                                dr%default_smearing(b3), smearing)
+
+                    do i=1, 3
+                        w(i) = dot_product(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), reclat(:, i))**2
+                    end do
+                    sigma = sqrt(maxval(w) * 0.5_r8)
                 end select
 
                 ! This is the multiplication of eigv of phonons 1 and 2 and now 3

From 1b04be47672c3815fa6183b171e692d87a350c18 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Sun, 27 Apr 2025 23:13:56 +0200
Subject: [PATCH 02/14] Still a work in progress

---
 src/thermal_conductivity/scattering.f90       | 129 ++++++++++++------
 .../scattering_fourphonon.f90                 | 121 +++++++++++++---
 .../scattering_isotope.f90                    |   9 +-
 .../scattering_threephonon.f90                | 110 ++++++++++++---
 4 files changed, 291 insertions(+), 78 deletions(-)

diff --git a/src/thermal_conductivity/scattering.f90 b/src/thermal_conductivity/scattering.f90
index 7cc4aa4c..42728fc9 100644
--- a/src/thermal_conductivity/scattering.f90
+++ b/src/thermal_conductivity/scattering.f90
@@ -40,6 +40,10 @@ module scattering
     real(r8), dimension(:, :), allocatable :: be, sigsq
     !> The scattering matrix
     real(r8), dimension(:, :), allocatable :: Xi
+    !> The reciprocal lattice vectors scaled by the qgrid, for adaptive broadening
+    real(r8), dimension(3, 3) :: reclat = -lo_huge
+    !> The smearing parameter
+    real(r8) :: sigma
 
 contains
     !> Generate the scattering amplitudes
@@ -92,6 +96,12 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
         real(r8) :: sigma
         integer, dimension(3) :: fq, fqp
 
+        real(r8), dimension(3) :: v0
+        real(r8), dimension(:, :), allocatable :: allsig
+        real(r8), dimension(:), allocatable :: sigsort, sigmin, sigmax
+        real(r8) :: f0
+        integer :: i25, i75
+
         ! grid dimensions
         select type (qp)
         class is (lo_fft_mesh)
@@ -123,72 +133,70 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
         allocate (sr%be(qp%n_irr_point, dr%n_mode))
         allocate (sr%sigsq(qp%n_irr_point, dr%n_mode))
         call mem%allocate(sigavg, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
-        call mem%allocate(bla, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+        call mem%allocate(allsig, [qp%n_irr_point, dr%n_mode], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+        call mem%allocate(sigsort, qp%n_irr_point, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+        call mem%allocate(sigmin, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+        call mem%allocate(sigmax, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+
+        sigavg = 0.0_r8
+        allsig = 0.0_r8
+        sigsort = 0.0_r8
+        sigmin = lo_huge
+        sigmax = -lo_huge
+
         sr%be = 0.0_r8
         sr%sigsq = 0.0_r8
-        bla = 0.0_r8
-        sigavg = 0.0_r8
+        sr%sigma = opts%sigma * lo_frequency_THz_to_Hartree
         do j=1, 3
-            reclat(:, j) = uc%reciprocal_latticevectors(:, j) / dims(j)
+            sr%reclat(:, j) = uc%reciprocal_latticevectors(:, j) / dims(j)
         end do
         do q1 = 1, qp%n_irr_point
             do b1 = 1, dr%n_mode
+                ! Skip gamma point
+                if (dr%iq(q1)%omega(b1) .lt. lo_freqtol) cycle
+
                 if (opts%classical) then
                     sr%be(q1, b1) = lo_kb_hartree*opts%temperature/dr%iq(q1)%omega(b1)
                 else
                     sr%be(q1, b1) = lo_planck(opts%temperature, dr%iq(q1)%omega(b1))
                 end if
+                sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), dr%default_smearing(b1), 1.0_r8)**2
 
-                do j=1, 3
-                    w(j) = dot_product(dr%iq(q1)%vel(:, b1), reclat(:, j))**2
-                end do
-!               sigma = sum(w(1:3))
-                sigma = maxval(w(1:3))
-                sr%sigsq(q1, b1) = sigma
-                if (sigma .gt. 0) then
-                    sigavg(b1) = sigavg(b1) + sigma * qp%ip(q1)%integration_weight
-                    bla(b1) = bla(b1) + 1.0_r8
-                end if
-
-                fq = singlet_to_triplet(qp%ip(q1)%full_index, dims(2), dims(3))
-                do j=1, 3
-                    fqp = fq
-                    fqp(j) = mod(fqp(j) + 1, dims(j))
-                    if (fqp(j) .eq. 0) fqp(j) = dims(j)
-                    q2 = triplet_to_singlet(fqp, dims(2), dims(3))
-!                   w(j) = dot_product(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b1), reclat(:, j))**2
-!                   w(j) = sum((dr%iq(q1)%vel(j, b1) * reclat(:, j))**2)
-                    w(j) = 0.5_r8 * (dr%iq(q1)%omega(b1) - dr%aq(q2)%omega(b1))
-
-                    fqp = fq
-                    fqp(j) = fqp(j) - 1
-                    if (fqp(j) .eq. 0) fqp(j) = dims(j)
-                    q2 = triplet_to_singlet(fqp, dims(2), dims(3))
-!                   w(j+3) = dot_product(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b1), reclat(:, j))**2
-!                   w(j+3) = sum((dr%iq(q1)%vel(j, b1) * reclat(:, j))**2)
-                    w(j+3) = 0.5_r8 * (dr%iq(q1)%omega(b1) - dr%aq(q2)%omega(b1))
+!               sigma = norm2(dr%iq(q1)%vel(:, b1)) * qp%ip(q1)%radius * lo_twopi / sqrt(2.0_r8)
+!               sr%sigsq(q1, b1) = sigma**2
 
-                end do
-!               sigma = (sum(w * lo_twopi / sqrt(2.0_r8)) / 6.0_r8)**2
-                sigma = sum(w)**2 / 6.0_r8
-!               sigma = norm2(dr%iq(q1)%vel(:, b1)) * qp%ip(q1)%radius
-!               sigma = sigma**2
-!               sigma = maxval(w(1:3))
-                sigavg(b1) = sigavg(b1) + sigma * qp%ip(q1)%integration_weight
+                v0 = matmul(abs(dr%iq(q1)%vel(:, b1)), sr%reclat)**2
+                sigma = maxval(v0*0.5_r8)
                 sr%sigsq(q1, b1) = sigma
 
-!               sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), minval(dr%default_smearing), opts%sigma)**2
-!               sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), dr%default_smearing(b1), opts%sigma)**2
+                f0 = norm2(dr%iq(q1)%vel(:, b1)) * qp%ip(q1)%radius * lo_twopi / sqrt(2.0_r8)
+                if (f0 .gt. 0.0_r8) then
+                    sigmin(b1) = min(sigmin(b1), f0**2)
+                    sigmax(b1) = max(sigmax(b1), f0**2)
+                    sigavg(b1) = sigavg(b1) + f0**2 * qp%ip(q1)%integration_weight
+                end if
             end do
         end do
-        sigavg = sigavg / bla
+
+        do b1=1, dr%n_mode
+!           sigavg(b1) = max(sigavg(b1) - 1.5_r8 * (sigmax(b1) - sigmin(b1)), sigmin(b1))
+            sigavg(b1) = sigavg(b1)
+        end do
+
         do q1=1, qp%n_irr_point
             do b1=1, dr%n_mode
-                if (sqrt(sr%sigsq(q1, b1)) .lt. lo_freqtol) sr%sigsq(q1, b1) = sigavg(b1)
+                if (sr%sigsq(q1, b1) .lt. 1e-20_r8) then
+                sr%sigsq(q1, b1) = max(sr%sigsq(q1, b1), (0.25_r8 * dr%default_smearing(b1))**2)
+!               sr%sigsq(q1, b1) = min(sr%sigsq(q1, b1), 4.0_r8 * dr%default_smearing(b1))
+!               sr%sigsq(q1, b1) = max(sr%sigsq(q1, b1), sigavg(b1))
+!               sr%sigsq(q1, b1) = max(sr%sigsq(q1, b1), sigmin(b1))
+                end if
             end do
         end do
+
         if (mw%talk) write(*, *) dr%default_smearing * lo_frequency_Hartree_to_meV
         if (mw%talk) write(*, *) sqrt(sigavg) * lo_frequency_Hartree_to_meV
+        if (mw%talk) write(*, *) sqrt(sigmin) * lo_frequency_Hartree_to_meV
         call mem%deallocate(sigavg, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 
         if (mw%talk) write (*, *) '... distributing q-point/modes on MPI ranks'
@@ -413,4 +421,39 @@ function sr_size_in_mem(sr) result(mem)
     if (allocated(sr%sigsq)) mem = mem + storage_size(sr%sigsq)*size(sr%sigsq)
     mem = mem/8
 end function
+
+subroutine approximate_hessian(q1, b1, qp, dims, hess)
+    !> The qpoint-mesh
+    type(lo_qpoint_mesh), intent(in) :: qp
+    !> The qpoint
+    integer, intent(in) :: q1
+    !> The mode
+    integer, intent(in) :: b1
+    !> The dimension of the grid
+    integer, dimension(3), intent(in) :: dims
+    !> The Hessian
+    real(r8), dimension(3, 3), intent(out) :: hess
+
+!   fq = qp%ip(q1)%full_index
+!   om1 = dr%iq(q1)%omega(b1)
+!   vel1 = dr%iq(q1)%vel(:, b1)
+!   qtriplet = singlet_to_triplet(q1, dims(2), dims(3))
+
+!   hess = 0.0_r8
+!   do i=1, 3
+!       qt2 = qtriplet
+!       qt2(i) = qt2(i) + 1
+!       if (qt2(i) .gt. dims(i)) qt2(i) = 1
+!       q2 = triplet_to_singlet(qt2, dims(2), dims(3))
+!       vel2 = dr%aq(q2)%vel(:, b2)
+
+!       dq = norm2(qp%ap(q2)%r - qp%ap(fq)%r)
+
+!       do j=1, 3
+!           hess(i, j) = hess(i, j) + (vel2(i) - vel1(i)) / dq
+!       end do
+!   end do
+!   hess = hess / 3.0_r8
+
+end subroutine
 end module
diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90
index 32ba10e9..28421071 100644
--- a/src/thermal_conductivity/scattering_fourphonon.f90
+++ b/src/thermal_conductivity/scattering_fourphonon.f90
@@ -54,6 +54,8 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     real(r8) :: f0, f1, f2, f3, f4, f5, f6, f7, fall
     !> The reducible triplet corresponding to the currently computed quartet
     integer, dimension(:, :), allocatable :: red_quartet
+    !> The approximated dirac
+    real(r8) :: d0, d1, d2, d3, d4
 
     real(r8), dimension(:, :), allocatable :: od_terms
 
@@ -146,18 +148,20 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
                     n4p = n4 + 1.0_r8
                     egv4 = dr%aq(q4)%egv(:, b4)/sqrt(om4)
 
-                    select case (integrationtype)
-                    case (1)
-                        sigma = smearing*lo_frequency_THz_to_Hartree
-                    case (2)
-                        sigma = sqrt(sr%sigsq(q1, b1) + &
-                                     sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
-                                     sr%sigsq(qp%ap(q3)%irreducible_index, b3) + &
-                                     sr%sigsq(qp%ap(q4)%irreducible_index, b4))
-                    case (6)
-                        sigma = qp%smearingparameter(dr%aq(q3)%vel(:, b3) - dr%aq(q4)%vel(:, b4), &
-                                                     dr%default_smearing(b3), smearing)
-                    end select
+!                   select case (integrationtype)
+!                   case (1)
+!                       sigma = smearing*lo_frequency_THz_to_Hartree
+!                   case (2)
+!                       sigma = sqrt(sr%sigsq(q1, b1) + &
+!                                    sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
+!                                    sr%sigsq(qp%ap(q3)%irreducible_index, b3) + &
+!                                    sr%sigsq(qp%ap(q4)%irreducible_index, b4))
+!                   case (6)
+!                       sigma = qp%smearingparameter(dr%aq(q3)%vel(:, b3) - dr%aq(q4)%vel(:, b4), &
+!                                                    dr%default_smearing(b3), smearing)
+!                   end select
+
+                    call get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3)
 
                     evp3 = 0.0_r8
                     call zgeru(dr%n_mode, dr%n_mode**3, (1.0_r8, 0.0_r8), egv4, 1, evp2, 1, evp3, dr%n_mode)
@@ -172,10 +176,14 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
                     plf3 = 3.0_r8*n4*n3p*n2p - n4p*n3*n2
 
                     ! Prefactors, including the matrix elements and dirac
-                    f0 = mult0*psisq*plf0*(lo_gauss(om1, om2 + om3 + om4, sigma) - lo_gauss(om1, -om2 - om3 - om4, sigma))
-                    f1 = mult1*psisq*plf1*(lo_gauss(om1, -om2 + om3 + om4, sigma) - lo_gauss(om1, om2 - om3 - om4, sigma))
-                    f2 = mult2*psisq*plf2*(lo_gauss(om1, -om3 + om2 + om4, sigma) - lo_gauss(om1, om3 - om2 - om4, sigma))
-                    f3 = mult3*psisq*plf3*(lo_gauss(om1, -om4 + om3 + om2, sigma) - lo_gauss(om1, om4 - om3 - om2, sigma))
+!                   f0 = mult0*psisq*plf0*(lo_gauss(om1, om2 + om3 + om4, sigma) - lo_gauss(om1, -om2 - om3 - om4, sigma))
+!                   f1 = mult1*psisq*plf1*(lo_gauss(om1, -om2 + om3 + om4, sigma) - lo_gauss(om1, om2 - om3 - om4, sigma))
+!                   f2 = mult2*psisq*plf2*(lo_gauss(om1, -om3 + om2 + om4, sigma) - lo_gauss(om1, om3 - om2 - om4, sigma))
+!                   f3 = mult3*psisq*plf3*(lo_gauss(om1, -om4 + om3 + om2, sigma) - lo_gauss(om1, om4 - om3 - om2, sigma))
+                    f0 = mult0*psisq*plf0*d0
+                    f1 = mult1*psisq*plf1*d1
+                    f2 = mult2*psisq*plf2*d2
+                    f3 = mult3*psisq*plf3*d3
 
                     fall = f0 + f1 + f2 + f3
 
@@ -268,6 +276,87 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     call mem%deallocate(qgridfull1, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(qgridfull2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(od_terms, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+
+    contains
+subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3)
+    !> The scattering amplitudes
+    type(lo_scattering_rates), intent(in) :: sr
+    !> The qpoint mesh
+    class(lo_qpoint_mesh), intent(in) :: qp
+    !> Harmonic dispersions
+    type(lo_phonon_dispersions), intent(in) :: dr
+    !> The qpoints
+    integer, intent(in) :: q1, q2, q3, q4
+    !> The modes
+    integer, intent(in) :: b1, b2, b3, b4
+    !> THe integration type
+    integer, intent(in) :: integrationtype
+    !> The approximated dirac
+    real(r8), intent(out) :: d0, d1, d2, d3
+
+    !> The frequencies
+    real(r8) :: om1, om2, om3, om4
+    ! Buffer to contains all the sigmas
+    real(r8), dimension(3, 6) :: allsig
+    !> An integer for the do loop
+    integer :: i
+
+    select case (integrationtype)
+    ! Gaussian smearing with fixed width
+    case (1)
+        sigma = smearing * lo_frequency_THz_to_Hartree
+    ! Adaptive Gaussian smearing with smeared frequency approach
+    case (2)
+        sigma = sqrt(sr%sigsq(q1, b1) + &
+                     sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
+                     sr%sigsq(qp%ap(q3)%irreducible_index, b3) + &
+                     sr%sigsq(qp%ap(q4)%irreducible_index, b4))
+    ! Adaptive Gaussian smearing with velocity difference approach
+    case (6)
+        ! We take an average to respect the symmetry of the scattering matrix
+        allsig = 0.0_r8
+
+        allsig(:, 1) = matmul(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b2), sr%reclat)**2
+        allsig(:, 2) = matmul(dr%iq(q1)%vel(:, b1) - dr%aq(q3)%vel(:, b3), sr%reclat)**2
+        allsig(:, 3) = matmul(dr%iq(q1)%vel(:, b1) - dr%aq(q4)%vel(:, b4), sr%reclat)**2
+        allsig(:, 4) = matmul(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), sr%reclat)**2
+        allsig(:, 5) = matmul(dr%aq(q2)%vel(:, b2) - dr%aq(q4)%vel(:, b4), sr%reclat)**2
+        allsig(:, 6) = matmul(dr%aq(q3)%vel(:, b3) - dr%aq(q4)%vel(:, b4), sr%reclat)**2
+        sigma = 0.0_r8
+        do i=1, 6
+            sigma = sigma + sqrt(maxval(allsig(:, i)) * 0.5_r8) / 6.0_r8
+        end do
+    end select
+    ! To catch edge cases, where the sum rule is actually working
+!   if (sigma .lt. lo_freqtol) sigma = 1.0_r8 / sqrt(lo_twopi)
+!   if (sigma .lt. lo_freqtol) sigma = minval(dr%default_smearing)
+
+    d0 = 0.0_r8
+    d1 = 0.0_r8
+    d2 = 0.0_r8
+    d3 = 0.0_r8
+
+    om1 = dr%iq(q1)%omega(b1)
+    om2 = dr%aq(q2)%omega(b2)
+    om3 = dr%aq(q3)%omega(b3)
+    om4 = dr%aq(q4)%omega(b4)
+    ! We can have some cases actually respecting strictly the sum rules
+    if (abs(om1 - om2) .lt. lo_freqtol .and. abs(om3 - om4) .lt. lo_freqtol) then
+        d2 = 1.0_r8
+        d3 = 1.0_r8
+    elseif (abs(om1 - om3) .lt. lo_freqtol .and. abs(om2 - om4) .lt. lo_freqtol) then
+        d1 = 1.0_r8
+        d3 = 1.0_r8
+    elseif (abs(om1 - om4) .lt. lo_freqtol .and. abs(om2 - om3) .lt. lo_freqtol) then
+        d1 = 1.0_r8
+        d2 = 1.0_r8
+    else
+        d0 = lo_gauss(om1, om2 + om3 + om4, sigma)  - lo_gauss(om1, -om2 - om3 - om4, sigma)
+        d1 = lo_gauss(om1, -om2 + om3 + om4, sigma) - lo_gauss(om1, om2 - om3 - om4, sigma)
+        d2 = lo_gauss(om1, -om3 + om2 + om4, sigma) - lo_gauss(om1, om3 - om2 - om4, sigma)
+        d3 = lo_gauss(om1, -om4 + om3 + om2, sigma) - lo_gauss(om1, om4 - om3 - om2, sigma)
+    end if
+end subroutine
 end subroutine
 
 subroutine quartet_is_irreducible(qp, uc, q1, q2, q3, q4, isred, red_quartet, mw, mem)
diff --git a/src/thermal_conductivity/scattering_isotope.f90 b/src/thermal_conductivity/scattering_isotope.f90
index f0855cbc..2a67b68f 100644
--- a/src/thermal_conductivity/scattering_isotope.f90
+++ b/src/thermal_conductivity/scattering_isotope.f90
@@ -26,6 +26,8 @@ subroutine compute_isotope_scattering(il, sr, qp, dr, uc, temperature, &
 
     ! Eigenvectors
     complex(r8), dimension(uc%na*3, 2) :: egviso
+    !> For the broadening calculation
+    real(r8), dimension(3) :: allsig
     ! prefactor and phonon buffers
     real(r8) :: om1, om2, sigma, psisq, prefactor, f0
     ! Integers for do loops
@@ -49,8 +51,11 @@ subroutine compute_isotope_scattering(il, sr, qp, dr, uc, temperature, &
                 sigma = sqrt(sr%sigsq(q1, b1) + &
                              sr%sigsq(qp%ap(q2)%irreducible_index, b2))
             case (6)
-                sigma = qp%smearingparameter(dr%aq(q2)%vel(:, b2), &
-                                             dr%default_smearing(b2), smearing)
+                allsig = matmul(dr%aq(q2)%vel(:, b2), sr%reclat)**2
+                sigma = sqrt(maxval(allsig) * 0.5_r8)
+                sigma = max(sigma, dr%default_smearing(b2) * 0.25_r8)
+!               sigma = qp%smearingparameter(dr%aq(q2)%vel(:, b2), &
+!                                            dr%default_smearing(b2), smearing)
             end select
 
             i = (q2 - 1)*dr%n_mode + b2
diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90
index e0850542..64d6f1eb 100644
--- a/src/thermal_conductivity/scattering_threephonon.f90
+++ b/src/thermal_conductivity/scattering_threephonon.f90
@@ -48,9 +48,11 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     integer, dimension(:, :), allocatable :: red_triplet
     !> buff to keep the off diagonal scattering matrix elements
     real(r8), dimension(:, :), allocatable :: od_terms
+    !> The approximated dirac
+    real(r8) :: d0, d1
 
-    real(r8), dimension(3, 3) :: reclat
-    real(r8), dimension(3) :: w
+    real(r8), dimension(3, 3) :: reclat, wig
+    real(r8), dimension(3) :: w, sigsig
 
     do i=1, 3
         reclat(:, i) = uc%reciprocal_latticevectors(:, i) / mcg%full_dims(i)
@@ -106,22 +108,33 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
                 if (om3 .lt. lo_freqtol) cycle
                 egv3 = dr%aq(q3)%egv(:, b3)/sqrt(om3)
 
-                select case (integrationtype)
-                case (1)
-                    sigma = smearing*lo_frequency_THz_to_Hartree
-                case (2)
-                    sigma = sqrt(sr%sigsq(q1, b1) + &
-                                 sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
-                                 sr%sigsq(qp%ap(q3)%irreducible_index, b3))
-                case (6)
+!               select case (integrationtype)
+!               case (1)
+!                   sigma = smearing*lo_frequency_THz_to_Hartree
+!               case (2)
+!                   sigma = sqrt(sr%sigsq(q1, b1) + &
+!                                sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
+!                                sr%sigsq(qp%ap(q3)%irreducible_index, b3))
+!               case (6)
 !                   sigma = qp%smearingparameter(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), &
 !                                                dr%default_smearing(b3), smearing)
 
-                    do i=1, 3
-                        w(i) = dot_product(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), reclat(:, i))**2
-                    end do
-                    sigma = sqrt(maxval(w) * 0.5_r8)
-                end select
+!                   wig = 0.0_r8
+!                   sigsig = 0.0_r8
+!                   do i=1, 3
+!                       w(i) = dot_product(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), reclat(:, i))**2
+!                       wig(i, 1) = dot_product(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b2), reclat(:, i))**2
+!                       wig(i, 2) = dot_product(dr%iq(q1)%vel(:, b1) - dr%aq(q3)%vel(:, b3), reclat(:, i))**2
+!                       wig(i, 3) = dot_product(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), reclat(:, i))**2
+!                   end do
+!                   do i=1, 3
+!                       sigsig(i) = sqrt(maxval(wig(:, i)) * 0.5_r8)
+!                   end do
+!                   sigma = sum(sigsig) / 3.0_r8
+!                   ! sigma = sqrt(maxval(w) * 0.5_r8)
+!               end select
+
+                call get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1)
 
                 ! This is the multiplication of eigv of phonons 1 and 2 and now 3
                 evp2 = 0.0_r8
@@ -138,8 +151,10 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
                 ! The prefactor for the scattering
                 plf0 = n2 - n3
                 plf1 = n2 + n3 + 1.0_r8
-                f0 = perm*psisq*plf0*(lo_gauss(om1, -om2 + om3, sigma) - lo_gauss(om1, om2 - om3, sigma))
-                f1 = perm*psisq*plf1*(lo_gauss(om1, om2 + om3, sigma) - lo_gauss(om1, -om2 - om3, sigma))
+!               f0 = perm*psisq*plf0*(lo_gauss(om1, -om2 + om3, sigma) - lo_gauss(om1, om2 - om3, sigma))
+!               f1 = perm*psisq*plf1*(lo_gauss(om1, om2 + om3, sigma) - lo_gauss(om1, -om2 - om3, sigma))
+                f0 = perm*psisq*plf0*d0
+                f1 = perm*psisq*plf1*d1
 
                 ! And we add everything for each triplet equivalent to the one we are actually computing
                 do i = 1, size(red_triplet, 2)
@@ -234,6 +249,67 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     call mem%deallocate(qgridfull, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(od_terms, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     deallocate (red_triplet)
+
+    contains
+subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1)
+    !> The scattering amplitudes
+    type(lo_scattering_rates), intent(in) :: sr
+    !> The qpoint mesh
+    class(lo_qpoint_mesh), intent(in) :: qp
+    !> Harmonic dispersions
+    type(lo_phonon_dispersions), intent(in) :: dr
+    !> The qpoints
+    integer, intent(in) :: q1, q2, q3
+    !> The modes
+    integer, intent(in) :: b1, b2, b3
+    !> THe integration type
+    integer, intent(in) :: integrationtype
+    !> The approximated dirac
+    real(r8), intent(out) :: d0, d1
+
+    !> The frequencies
+    real(r8) :: om1, om2, om3
+    ! Buffer to contains all the sigmas
+    real(r8), dimension(3, 3) :: allsig
+    !> An integer for the do loop
+    integer :: i
+
+    select case (integrationtype)
+    ! Gaussian smearing with fixed width
+    case (1)
+        sigma = smearing * lo_frequency_THz_to_Hartree
+    ! Adaptive Gaussian smearing with smeared frequency approach
+    case (2)
+        sigma = sqrt(sr%sigsq(q1, b1) + &
+                     sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
+                     sr%sigsq(qp%ap(q3)%irreducible_index, b3))
+    ! Adaptive Gaussian smearing with velocity difference approach
+    case (6)
+        allsig = 0.0_r8
+        allsig(:, 1) = matmul(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b2), sr%reclat)**2
+        allsig(:, 2) = matmul(dr%iq(q1)%vel(:, b1) - dr%aq(q3)%vel(:, b3), sr%reclat)**2
+        allsig(:, 3) = matmul(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), sr%reclat)**2
+        sigma = 0.0_r8
+        do i=1, 3
+            sigma = sigma + sqrt(maxval(allsig(:, i)) * 0.5_r8) / 3.0_r8
+        end do
+
+        sigma = sqrt(maxval(matmul(dr%aq(q2)%vel(:, b1) - dr%aq(q3)%vel(:, b3), sr%reclat)**2) * 0.5_r8)
+    end select
+
+    om1 = dr%iq(q1)%omega(b1)
+    om2 = dr%aq(q2)%omega(b2)
+    om3 = dr%aq(q3)%omega(b3)
+    d0 = 0.0_r8
+    d1 = 0.0_r8
+    ! We cannot have interaction with same mode (or degenerate)
+    if (abs(om1 - om2) .gt. lo_freqtol .and. &
+        abs(om1 - om3) .gt. lo_freqtol .and. &
+        abs(om2 - om3) .gt. lo_freqtol) then
+        d0 = lo_gauss(om1, -om2 + om3, sigma) - lo_gauss(om1, om2 - om3, sigma)
+        d1 = lo_gauss(om1, om2 + om3, sigma) - lo_gauss(om1, -om2 - om3, sigma)
+    end if
+end subroutine
 end subroutine
 
 !> Get the Fourier transform of the third order matrix element

From e112b448fb11dce843435d53a368249d6de60766 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Thu, 8 May 2025 14:22:11 +0200
Subject: [PATCH 03/14] Working broadening parameters

---
 src/thermal_conductivity/scattering.f90 | 60 +++++++------------------
 1 file changed, 16 insertions(+), 44 deletions(-)

diff --git a/src/thermal_conductivity/scattering.f90 b/src/thermal_conductivity/scattering.f90
index 42728fc9..f6914bc1 100644
--- a/src/thermal_conductivity/scattering.f90
+++ b/src/thermal_conductivity/scattering.f90
@@ -83,6 +83,12 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
     type(lo_montecarlo_grid) :: mcg3, mcg4
 
     init: block
+        !> The average broadening factor for each branch
+        real(r8), dimension(:), allocatable :: sigavg
+        !> A buffer for a vector
+        real(r8), dimension(3) :: v0
+        !> The calculated sigma for a mode
+        real(r8) :: sigma
         !> To initialize the random number generator and timing
         real(r8) :: rseed
         !> The q-point grid dimension
@@ -90,18 +96,6 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
         !> Some integers for the do loop/indices
         integer :: q1, q2, b1, il, j, my_nqpoints, ctr
 
-        real(r8), dimension(:), allocatable :: sigavg, bla
-        real(r8), dimension(3, 3) :: reclat
-        real(r8), dimension(6) :: w
-        real(r8) :: sigma
-        integer, dimension(3) :: fq, fqp
-
-        real(r8), dimension(3) :: v0
-        real(r8), dimension(:, :), allocatable :: allsig
-        real(r8), dimension(:), allocatable :: sigsort, sigmin, sigmax
-        real(r8) :: f0
-        integer :: i25, i75
-
         ! grid dimensions
         select type (qp)
         class is (lo_fft_mesh)
@@ -133,23 +127,18 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
         allocate (sr%be(qp%n_irr_point, dr%n_mode))
         allocate (sr%sigsq(qp%n_irr_point, dr%n_mode))
         call mem%allocate(sigavg, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
-        call mem%allocate(allsig, [qp%n_irr_point, dr%n_mode], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
-        call mem%allocate(sigsort, qp%n_irr_point, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
-        call mem%allocate(sigmin, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
-        call mem%allocate(sigmax, dr%n_mode, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 
         sigavg = 0.0_r8
-        allsig = 0.0_r8
-        sigsort = 0.0_r8
-        sigmin = lo_huge
-        sigmax = -lo_huge
 
         sr%be = 0.0_r8
         sr%sigsq = 0.0_r8
+        ! This could be useful if integrationtype=1 -> fixed gaussian smearing
         sr%sigma = opts%sigma * lo_frequency_THz_to_Hartree
+        ! This could be useful if integrationtype=6 -> Adapative smearing from groupvel diff
         do j=1, 3
             sr%reclat(:, j) = uc%reciprocal_latticevectors(:, j) / dims(j)
         end do
+        ! And we get the values in the array
         do q1 = 1, qp%n_irr_point
             do b1 = 1, dr%n_mode
                 ! Skip gamma point
@@ -160,43 +149,26 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
                 else
                     sr%be(q1, b1) = lo_planck(opts%temperature, dr%iq(q1)%omega(b1))
                 end if
-                sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), dr%default_smearing(b1), 1.0_r8)**2
-
-!               sigma = norm2(dr%iq(q1)%vel(:, b1)) * qp%ip(q1)%radius * lo_twopi / sqrt(2.0_r8)
-!               sr%sigsq(q1, b1) = sigma**2
+!               sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), dr%default_smearing(b1), 1.0_r8)**2
 
                 v0 = matmul(abs(dr%iq(q1)%vel(:, b1)), sr%reclat)**2
+                ! This allows to work around problem with 2D materials
                 sigma = maxval(v0*0.5_r8)
                 sr%sigsq(q1, b1) = sigma
 
-                f0 = norm2(dr%iq(q1)%vel(:, b1)) * qp%ip(q1)%radius * lo_twopi / sqrt(2.0_r8)
-                if (f0 .gt. 0.0_r8) then
-                    sigmin(b1) = min(sigmin(b1), f0**2)
-                    sigmax(b1) = max(sigmax(b1), f0**2)
-                    sigavg(b1) = sigavg(b1) + f0**2 * qp%ip(q1)%integration_weight
-                end if
+                ! Let's accumulate the average broadening factor for each mode, for a sanity check
+                sigavg(b1) = sigavg(b1) + sigma * qp%ip(q1)%integration_weight
             end do
         end do
 
-        do b1=1, dr%n_mode
-!           sigavg(b1) = max(sigavg(b1) - 1.5_r8 * (sigmax(b1) - sigmin(b1)), sigmin(b1))
-            sigavg(b1) = sigavg(b1)
-        end do
-
+        ! We add a little sanity check, in case a broadening factor is zero
         do q1=1, qp%n_irr_point
             do b1=1, dr%n_mode
-                if (sr%sigsq(q1, b1) .lt. 1e-20_r8) then
-                sr%sigsq(q1, b1) = max(sr%sigsq(q1, b1), (0.25_r8 * dr%default_smearing(b1))**2)
-!               sr%sigsq(q1, b1) = min(sr%sigsq(q1, b1), 4.0_r8 * dr%default_smearing(b1))
-!               sr%sigsq(q1, b1) = max(sr%sigsq(q1, b1), sigavg(b1))
-!               sr%sigsq(q1, b1) = max(sr%sigsq(q1, b1), sigmin(b1))
+                if (sr%sigsq(q1, b1) .lt. sigavg(b1) / 10.0_r8) then
+                    sr%sigsq(q1, b1) = sigavg(b1) / 10.0_r8
                 end if
             end do
         end do
-
-        if (mw%talk) write(*, *) dr%default_smearing * lo_frequency_Hartree_to_meV
-        if (mw%talk) write(*, *) sqrt(sigavg) * lo_frequency_Hartree_to_meV
-        if (mw%talk) write(*, *) sqrt(sigmin) * lo_frequency_Hartree_to_meV
         call mem%deallocate(sigavg, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 
         if (mw%talk) write (*, *) '... distributing q-point/modes on MPI ranks'

From 27e5cdf7e3f1c9ee0e8fd35a2a1bbbe1d6bccb5f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Fri, 9 May 2025 09:03:05 +0200
Subject: [PATCH 04/14] Little speed-up by skipping last IFC projection if
 omega diff larger than 4 sigma

---
 src/thermal_conductivity/scattering.f90       |  2 +-
 .../scattering_fourphonon.f90                 | 73 +++++++++++++------
 .../scattering_threephonon.f90                | 65 +++++++++--------
 3 files changed, 84 insertions(+), 56 deletions(-)

diff --git a/src/thermal_conductivity/scattering.f90 b/src/thermal_conductivity/scattering.f90
index f6914bc1..f0deca9e 100644
--- a/src/thermal_conductivity/scattering.f90
+++ b/src/thermal_conductivity/scattering.f90
@@ -152,7 +152,7 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
 !               sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), dr%default_smearing(b1), 1.0_r8)**2
 
                 v0 = matmul(abs(dr%iq(q1)%vel(:, b1)), sr%reclat)**2
-                ! This allows to work around problem with 2D materials
+                ! This allows to work around problems with 2D materials
                 sigma = maxval(v0*0.5_r8)
                 sr%sigsq(q1, b1) = sigma
 
diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90
index 28421071..d299d77e 100644
--- a/src/thermal_conductivity/scattering_fourphonon.f90
+++ b/src/thermal_conductivity/scattering_fourphonon.f90
@@ -56,6 +56,8 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     integer, dimension(:, :), allocatable :: red_quartet
     !> The approximated dirac
     real(r8) :: d0, d1, d2, d3, d4
+    !> Can we skip the process ?
+    logical :: isok
 
     real(r8), dimension(:, :), allocatable :: od_terms
 
@@ -148,20 +150,8 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
                     n4p = n4 + 1.0_r8
                     egv4 = dr%aq(q4)%egv(:, b4)/sqrt(om4)
 
-!                   select case (integrationtype)
-!                   case (1)
-!                       sigma = smearing*lo_frequency_THz_to_Hartree
-!                   case (2)
-!                       sigma = sqrt(sr%sigsq(q1, b1) + &
-!                                    sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
-!                                    sr%sigsq(qp%ap(q3)%irreducible_index, b3) + &
-!                                    sr%sigsq(qp%ap(q4)%irreducible_index, b4))
-!                   case (6)
-!                       sigma = qp%smearingparameter(dr%aq(q3)%vel(:, b3) - dr%aq(q4)%vel(:, b4), &
-!                                                    dr%default_smearing(b3), smearing)
-!                   end select
-
-                    call get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3)
+                    call get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, isok)
+                    if (.not. isok) cycle
 
                     evp3 = 0.0_r8
                     call zgeru(dr%n_mode, dr%n_mode**3, (1.0_r8, 0.0_r8), egv4, 1, evp2, 1, evp3, dr%n_mode)
@@ -176,10 +166,6 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
                     plf3 = 3.0_r8*n4*n3p*n2p - n4p*n3*n2
 
                     ! Prefactors, including the matrix elements and dirac
-!                   f0 = mult0*psisq*plf0*(lo_gauss(om1, om2 + om3 + om4, sigma) - lo_gauss(om1, -om2 - om3 - om4, sigma))
-!                   f1 = mult1*psisq*plf1*(lo_gauss(om1, -om2 + om3 + om4, sigma) - lo_gauss(om1, om2 - om3 - om4, sigma))
-!                   f2 = mult2*psisq*plf2*(lo_gauss(om1, -om3 + om2 + om4, sigma) - lo_gauss(om1, om3 - om2 - om4, sigma))
-!                   f3 = mult3*psisq*plf3*(lo_gauss(om1, -om4 + om3 + om2, sigma) - lo_gauss(om1, om4 - om3 - om2, sigma))
                     f0 = mult0*psisq*plf0*d0
                     f1 = mult1*psisq*plf1*d1
                     f2 = mult2*psisq*plf2*d2
@@ -278,7 +264,7 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     call mem%deallocate(od_terms, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 
     contains
-subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3)
+subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, isok)
     !> The scattering amplitudes
     type(lo_scattering_rates), intent(in) :: sr
     !> The qpoint mesh
@@ -293,13 +279,15 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
     integer, intent(in) :: integrationtype
     !> The approximated dirac
     real(r8), intent(out) :: d0, d1, d2, d3
+    !> Can we skip the rest of the calculation ?
+    logical, intent(out) :: isok
 
     !> The frequencies
     real(r8) :: om1, om2, om3, om4
     ! Buffer to contains all the sigmas
     real(r8), dimension(3, 6) :: allsig
     !> An integer for the do loop
-    integer :: i
+    integer :: i, j
 
     select case (integrationtype)
     ! Gaussian smearing with fixed width
@@ -344,18 +332,55 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
     if (abs(om1 - om2) .lt. lo_freqtol .and. abs(om3 - om4) .lt. lo_freqtol) then
         d2 = 1.0_r8
         d3 = 1.0_r8
+        j = j + 1
     elseif (abs(om1 - om3) .lt. lo_freqtol .and. abs(om2 - om4) .lt. lo_freqtol) then
         d1 = 1.0_r8
         d3 = 1.0_r8
+        j = j + 1
     elseif (abs(om1 - om4) .lt. lo_freqtol .and. abs(om2 - om3) .lt. lo_freqtol) then
         d1 = 1.0_r8
         d2 = 1.0_r8
+        j = j + 1
     else
-        d0 = lo_gauss(om1, om2 + om3 + om4, sigma)  - lo_gauss(om1, -om2 - om3 - om4, sigma)
-        d1 = lo_gauss(om1, -om2 + om3 + om4, sigma) - lo_gauss(om1, om2 - om3 - om4, sigma)
-        d2 = lo_gauss(om1, -om3 + om2 + om4, sigma) - lo_gauss(om1, om3 - om2 - om4, sigma)
-        d3 = lo_gauss(om1, -om4 + om3 + om2, sigma) - lo_gauss(om1, om4 - om3 - om2, sigma)
+!       d0 = lo_gauss(om1, om2 + om3 + om4, sigma)  - lo_gauss(om1, -om2 - om3 - om4, sigma)
+!       d1 = lo_gauss(om1, -om2 + om3 + om4, sigma) - lo_gauss(om1, om2 - om3 - om4, sigma)
+!       d2 = lo_gauss(om1, -om3 + om2 + om4, sigma) - lo_gauss(om1, om3 - om2 - om4, sigma)
+!       d3 = lo_gauss(om1, -om4 + om3 + om2, sigma) - lo_gauss(om1, om4 - om3 - om2, sigma)
+
+        if (abs(om1 - om2 - om3 - om4) .lt. 4.0_r8 * sigma) then
+            d0 = d0 + lo_gauss(om1, om2 + om3 + om4, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 + om2 + om3 + om4) .lt. 4.0_r8 * sigma) then
+            d0 = d0 - lo_gauss(om1, -om2 - om3 - om4, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 + om2 - om3 - om4) .lt. 4.0_r8 * sigma) then
+            d1 = d1 + lo_gauss(om1, -om2 + om3 + om4, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 - om2 + om3 + om4) .lt. 4.0_r8 * sigma) then
+            d1 = d1 - lo_gauss(om1, om2 - om3 - om4, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 -om2 + om3 - om4) .lt. 4.0_r8 * sigma) then
+            d2 = d2 + lo_gauss(om1, om2 - om3 + om4, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 + om2 - om3 + om4) .lt. 4.0_r8 * sigma) then
+            d2 = d2 - lo_gauss(om1, -om2 + om3 - om4, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 - om2 - om3 + om4) .lt. 4.0_r8 * sigma) then
+            d3 = d3 + lo_gauss(om1, om2 + om3 - om4, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 + om2 + om3 - om4) .lt. 4.0_r8 * sigma) then
+            d3 = d3 - lo_gauss(om1, -om2 - om3 + om4, sigma)
+            j = j + 1
+        end if
     end if
+    if (j .gt. 0) isok = .true.
 end subroutine
 end subroutine
 
diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90
index 64d6f1eb..57b78f57 100644
--- a/src/thermal_conductivity/scattering_threephonon.f90
+++ b/src/thermal_conductivity/scattering_threephonon.f90
@@ -54,6 +54,8 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     real(r8), dimension(3, 3) :: reclat, wig
     real(r8), dimension(3) :: w, sigsig
 
+    logical :: isok
+
     do i=1, 3
         reclat(:, i) = uc%reciprocal_latticevectors(:, i) / mcg%full_dims(i)
     end do
@@ -108,33 +110,9 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
                 if (om3 .lt. lo_freqtol) cycle
                 egv3 = dr%aq(q3)%egv(:, b3)/sqrt(om3)
 
-!               select case (integrationtype)
-!               case (1)
-!                   sigma = smearing*lo_frequency_THz_to_Hartree
-!               case (2)
-!                   sigma = sqrt(sr%sigsq(q1, b1) + &
-!                                sr%sigsq(qp%ap(q2)%irreducible_index, b2) + &
-!                                sr%sigsq(qp%ap(q3)%irreducible_index, b3))
-!               case (6)
-!                   sigma = qp%smearingparameter(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), &
-!                                                dr%default_smearing(b3), smearing)
-
-!                   wig = 0.0_r8
-!                   sigsig = 0.0_r8
-!                   do i=1, 3
-!                       w(i) = dot_product(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), reclat(:, i))**2
-!                       wig(i, 1) = dot_product(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b2), reclat(:, i))**2
-!                       wig(i, 2) = dot_product(dr%iq(q1)%vel(:, b1) - dr%aq(q3)%vel(:, b3), reclat(:, i))**2
-!                       wig(i, 3) = dot_product(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), reclat(:, i))**2
-!                   end do
-!                   do i=1, 3
-!                       sigsig(i) = sqrt(maxval(wig(:, i)) * 0.5_r8)
-!                   end do
-!                   sigma = sum(sigsig) / 3.0_r8
-!                   ! sigma = sqrt(maxval(w) * 0.5_r8)
-!               end select
-
-                call get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1)
+                call get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, isok)
+
+                if (.not. isok) cycle
 
                 ! This is the multiplication of eigv of phonons 1 and 2 and now 3
                 evp2 = 0.0_r8
@@ -251,7 +229,7 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     deallocate (red_triplet)
 
     contains
-subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1)
+subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, isok)
     !> The scattering amplitudes
     type(lo_scattering_rates), intent(in) :: sr
     !> The qpoint mesh
@@ -267,6 +245,8 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     !> The approximated dirac
     real(r8), intent(out) :: d0, d1
 
+    logical, intent(out) :: isok
+
     !> The frequencies
     real(r8) :: om1, om2, om3
     ! Buffer to contains all the sigmas
@@ -274,6 +254,8 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     !> An integer for the do loop
     integer :: i
 
+    integer :: j
+
     select case (integrationtype)
     ! Gaussian smearing with fixed width
     case (1)
@@ -289,12 +271,13 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
         allsig(:, 1) = matmul(dr%iq(q1)%vel(:, b1) - dr%aq(q2)%vel(:, b2), sr%reclat)**2
         allsig(:, 2) = matmul(dr%iq(q1)%vel(:, b1) - dr%aq(q3)%vel(:, b3), sr%reclat)**2
         allsig(:, 3) = matmul(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), sr%reclat)**2
+
         sigma = 0.0_r8
         do i=1, 3
             sigma = sigma + sqrt(maxval(allsig(:, i)) * 0.5_r8) / 3.0_r8
         end do
 
-        sigma = sqrt(maxval(matmul(dr%aq(q2)%vel(:, b1) - dr%aq(q3)%vel(:, b3), sr%reclat)**2) * 0.5_r8)
+!       sigma = sqrt(maxval(matmul(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), sr%reclat)**2) * 0.5_r8)
     end select
 
     om1 = dr%iq(q1)%omega(b1)
@@ -302,13 +285,33 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     om3 = dr%aq(q3)%omega(b3)
     d0 = 0.0_r8
     d1 = 0.0_r8
+    j = 0
+    isok = .false.
     ! We cannot have interaction with same mode (or degenerate)
     if (abs(om1 - om2) .gt. lo_freqtol .and. &
         abs(om1 - om3) .gt. lo_freqtol .and. &
         abs(om2 - om3) .gt. lo_freqtol) then
-        d0 = lo_gauss(om1, -om2 + om3, sigma) - lo_gauss(om1, om2 - om3, sigma)
-        d1 = lo_gauss(om1, om2 + om3, sigma) - lo_gauss(om1, -om2 - om3, sigma)
+!       d0 = lo_gauss(om1, -om2 + om3, sigma) - lo_gauss(om1, om2 - om3, sigma)
+!       d1 = lo_gauss(om1, om2 + om3, sigma) - lo_gauss(om1, -om2 - om3, sigma)
+
+        if (abs(om1 + om2 - om3) .lt. 4.0_r8 * sigma) then
+            d0 = d0 + lo_gauss(om1, -om2 + om3, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 - om2 + om3) .lt. 4.0_r8 * sigma) then
+            d0 = d0 - lo_gauss(om1, om2 - om3, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 - om2 - om3) .lt. 4.0_r8 * sigma) then
+            d1 = d1 + lo_gauss(om1, om2 + om3, sigma)
+            j = j + 1
+        end if
+        if (abs(om1 + om2 + om3) .lt. 4.0_r8 * sigma) then
+            d1 = d1 - lo_gauss(om1, -om2 - om3, sigma)
+            j = j + 1
+        end if
     end if
+    if (j .gt. 0) isok = .true.
 end subroutine
 end subroutine
 

From 1b89eda9f222b3f17194faade09bb30376ab008e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Fri, 9 May 2025 11:59:06 +0200
Subject: [PATCH 05/14] Little fix for the broadening parameter

---
 src/thermal_conductivity/scattering.f90       | 18 +++++++-------
 .../scattering_fourphonon.f90                 | 10 ++------
 .../scattering_isotope.f90                    |  4 +---
 .../scattering_threephonon.f90                | 24 ++++++-------------
 4 files changed, 18 insertions(+), 38 deletions(-)

diff --git a/src/thermal_conductivity/scattering.f90 b/src/thermal_conductivity/scattering.f90
index f0deca9e..178c8398 100644
--- a/src/thermal_conductivity/scattering.f90
+++ b/src/thermal_conductivity/scattering.f90
@@ -43,7 +43,7 @@ module scattering
     !> The reciprocal lattice vectors scaled by the qgrid, for adaptive broadening
     real(r8), dimension(3, 3) :: reclat = -lo_huge
     !> The smearing parameter
-    real(r8) :: sigma
+    real(r8) :: sigma, thresh_sigma
 
 contains
     !> Generate the scattering amplitudes
@@ -138,6 +138,7 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
         do j=1, 3
             sr%reclat(:, j) = uc%reciprocal_latticevectors(:, j) / dims(j)
         end do
+
         ! And we get the values in the array
         do q1 = 1, qp%n_irr_point
             do b1 = 1, dr%n_mode
@@ -161,14 +162,10 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
             end do
         end do
 
-        ! We add a little sanity check, in case a broadening factor is zero
-        do q1=1, qp%n_irr_point
-            do b1=1, dr%n_mode
-                if (sr%sigsq(q1, b1) .lt. sigavg(b1) / 10.0_r8) then
-                    sr%sigsq(q1, b1) = sigavg(b1) / 10.0_r8
-                end if
-            end do
-        end do
+        ! We add a little baseline as a sanity check, to avoid pathological case when group velocities are near zero
+        sr%sigsq = sr%sigsq + maxval(sigavg) * 1e-4_r8
+        ! This is to have the sanity check in memory, for integrationtype 6
+        sr%thresh_sigma = maxval(sigavg) * 1e-4_r8
         call mem%deallocate(sigavg, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 
         if (mw%talk) write (*, *) '... distributing q-point/modes on MPI ranks'
@@ -285,7 +282,8 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
                 ! While we are at it, we can set other things
                 dr%iq(q1)%qs(b1) = 2.0_r8*dr%iq(q1)%linewidth(b1)
                 velnorm = norm2(dr%iq(q1)%vel(:, b1))
-                if (velnorm .gt. lo_phonongroupveltol) then
+!               if (velnorm .gt. lo_phonongroupveltol) then
+                if (dr%iq(q1)%linewidth(b1) .gt. lo_freqtol) then
                     dr%iq(q1)%mfp(:, b1) = dr%iq(q1)%vel(:, b1)/dr%iq(q1)%qs(b1)
                     dr%iq(q1)%scalar_mfp(b1) = velnorm/dr%iq(q1)%qs(b1)
                     dr%iq(q1)%F0(:, b1) = dr%iq(q1)%mfp(:, b1)
diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90
index d299d77e..9ba96ac2 100644
--- a/src/thermal_conductivity/scattering_fourphonon.f90
+++ b/src/thermal_conductivity/scattering_fourphonon.f90
@@ -292,7 +292,7 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
     select case (integrationtype)
     ! Gaussian smearing with fixed width
     case (1)
-        sigma = smearing * lo_frequency_THz_to_Hartree
+        sigma = sr%sigma
     ! Adaptive Gaussian smearing with smeared frequency approach
     case (2)
         sigma = sqrt(sr%sigsq(q1, b1) + &
@@ -314,10 +314,8 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
         do i=1, 6
             sigma = sigma + sqrt(maxval(allsig(:, i)) * 0.5_r8) / 6.0_r8
         end do
+        sigma = sigma + sr%thresh_sigma
     end select
-    ! To catch edge cases, where the sum rule is actually working
-!   if (sigma .lt. lo_freqtol) sigma = 1.0_r8 / sqrt(lo_twopi)
-!   if (sigma .lt. lo_freqtol) sigma = minval(dr%default_smearing)
 
     d0 = 0.0_r8
     d1 = 0.0_r8
@@ -342,10 +340,6 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
         d2 = 1.0_r8
         j = j + 1
     else
-!       d0 = lo_gauss(om1, om2 + om3 + om4, sigma)  - lo_gauss(om1, -om2 - om3 - om4, sigma)
-!       d1 = lo_gauss(om1, -om2 + om3 + om4, sigma) - lo_gauss(om1, om2 - om3 - om4, sigma)
-!       d2 = lo_gauss(om1, -om3 + om2 + om4, sigma) - lo_gauss(om1, om3 - om2 - om4, sigma)
-!       d3 = lo_gauss(om1, -om4 + om3 + om2, sigma) - lo_gauss(om1, om4 - om3 - om2, sigma)
 
         if (abs(om1 - om2 - om3 - om4) .lt. 4.0_r8 * sigma) then
             d0 = d0 + lo_gauss(om1, om2 + om3 + om4, sigma)
diff --git a/src/thermal_conductivity/scattering_isotope.f90 b/src/thermal_conductivity/scattering_isotope.f90
index 2a67b68f..2aeb7342 100644
--- a/src/thermal_conductivity/scattering_isotope.f90
+++ b/src/thermal_conductivity/scattering_isotope.f90
@@ -53,9 +53,7 @@ subroutine compute_isotope_scattering(il, sr, qp, dr, uc, temperature, &
             case (6)
                 allsig = matmul(dr%aq(q2)%vel(:, b2), sr%reclat)**2
                 sigma = sqrt(maxval(allsig) * 0.5_r8)
-                sigma = max(sigma, dr%default_smearing(b2) * 0.25_r8)
-!               sigma = qp%smearingparameter(dr%aq(q2)%vel(:, b2), &
-!                                            dr%default_smearing(b2), smearing)
+                sigma = sigma + sr%thresh_sigma
             end select
 
             i = (q2 - 1)*dr%n_mode + b2
diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90
index 57b78f57..ecd577d2 100644
--- a/src/thermal_conductivity/scattering_threephonon.f90
+++ b/src/thermal_conductivity/scattering_threephonon.f90
@@ -50,16 +50,9 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     real(r8), dimension(:, :), allocatable :: od_terms
     !> The approximated dirac
     real(r8) :: d0, d1
-
-    real(r8), dimension(3, 3) :: reclat, wig
-    real(r8), dimension(3) :: w, sigsig
-
+    !> Do we have to compute the scattering process ?
     logical :: isok
 
-    do i=1, 3
-        reclat(:, i) = uc%reciprocal_latticevectors(:, i) / mcg%full_dims(i)
-    end do
-
     ! We start by allocating everything
     call mem%allocate(ptf, dr%n_mode**3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%allocate(evp1, dr%n_mode**2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
@@ -110,8 +103,9 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
                 if (om3 .lt. lo_freqtol) cycle
                 egv3 = dr%aq(q3)%egv(:, b3)/sqrt(om3)
 
+                ! Get the weight for each processes
                 call get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, isok)
-
+                ! If all the weights are zero, skip the last part
                 if (.not. isok) cycle
 
                 ! This is the multiplication of eigv of phonons 1 and 2 and now 3
@@ -252,14 +246,12 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     ! Buffer to contains all the sigmas
     real(r8), dimension(3, 3) :: allsig
     !> An integer for the do loop
-    integer :: i
-
-    integer :: j
+    integer :: i, j
 
     select case (integrationtype)
     ! Gaussian smearing with fixed width
     case (1)
-        sigma = smearing * lo_frequency_THz_to_Hartree
+        sigma = sr%sigma
     ! Adaptive Gaussian smearing with smeared frequency approach
     case (2)
         sigma = sqrt(sr%sigsq(q1, b1) + &
@@ -276,8 +268,6 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
         do i=1, 3
             sigma = sigma + sqrt(maxval(allsig(:, i)) * 0.5_r8) / 3.0_r8
         end do
-
-!       sigma = sqrt(maxval(matmul(dr%aq(q2)%vel(:, b2) - dr%aq(q3)%vel(:, b3), sr%reclat)**2) * 0.5_r8)
     end select
 
     om1 = dr%iq(q1)%omega(b1)
@@ -291,8 +281,6 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     if (abs(om1 - om2) .gt. lo_freqtol .and. &
         abs(om1 - om3) .gt. lo_freqtol .and. &
         abs(om2 - om3) .gt. lo_freqtol) then
-!       d0 = lo_gauss(om1, -om2 + om3, sigma) - lo_gauss(om1, om2 - om3, sigma)
-!       d1 = lo_gauss(om1, om2 + om3, sigma) - lo_gauss(om1, -om2 - om3, sigma)
 
         if (abs(om1 + om2 - om3) .lt. 4.0_r8 * sigma) then
             d0 = d0 + lo_gauss(om1, -om2 + om3, sigma)
@@ -441,6 +429,8 @@ subroutine triplet_is_irreducible(qp, uc, q1, q2, q3, isred, red_triplet, mw, me
             if (qpp(1) .gt. q2) isred = .true.
         end do
 
+        if (isred) return
+
         ! For stability, replace points not on the grid with starting q-point
         ! Should not be needed if the grid respect the symmetries of the lattice though
         if (minval(newqp) .lt. 0) then

From 00eed06d07f1b06ccd1285a3dba29ed542b00fe6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Fri, 9 May 2025 12:32:59 +0200
Subject: [PATCH 06/14] Correcting deallocation bug

---
 src/thermal_conductivity/scattering_fourphonon.f90  | 3 +++
 src/thermal_conductivity/scattering_threephonon.f90 | 4 +---
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90
index 9ba96ac2..87f4afd1 100644
--- a/src/thermal_conductivity/scattering_fourphonon.f90
+++ b/src/thermal_conductivity/scattering_fourphonon.f90
@@ -262,6 +262,7 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     call mem%deallocate(qgridfull1, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(qgridfull2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(od_terms, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+    if (allocated(red_quartet)) deallocate(red_quartet)
 
     contains
 subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, isok)
@@ -457,6 +458,8 @@ subroutine quartet_is_irreducible(qp, uc, q1, q2, q3, q4, isred, red_quartet, mw
             if (qpp(1) .gt. q2 .or. qpp(2) .gt. q3) isred = .true.
         end do
 
+        if (isred) return
+
         if (minval(newqp) .lt. 0) then
             do j = 1, size(newqp, 2)
                 if (any(newqp(:, j) .lt. 0)) newqp(:, j) = [q2, q3, q4]
diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90
index ecd577d2..e4cd1762 100644
--- a/src/thermal_conductivity/scattering_threephonon.f90
+++ b/src/thermal_conductivity/scattering_threephonon.f90
@@ -123,8 +123,6 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
                 ! The prefactor for the scattering
                 plf0 = n2 - n3
                 plf1 = n2 + n3 + 1.0_r8
-!               f0 = perm*psisq*plf0*(lo_gauss(om1, -om2 + om3, sigma) - lo_gauss(om1, om2 - om3, sigma))
-!               f1 = perm*psisq*plf1*(lo_gauss(om1, om2 + om3, sigma) - lo_gauss(om1, -om2 - om3, sigma))
                 f0 = perm*psisq*plf0*d0
                 f1 = perm*psisq*plf1*d1
 
@@ -220,7 +218,7 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     call mem%deallocate(egv3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(qgridfull, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(od_terms, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
-    deallocate (red_triplet)
+    if (allocated(red_triplet)) deallocate (red_triplet)
 
     contains
 subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, isok)

From 87013ca7d3d2c3d0dfd76498aee84ed892bf0343 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Fri, 9 May 2025 14:57:09 +0200
Subject: [PATCH 07/14] Fix broadening for some pathological cases

---
 src/thermal_conductivity/kappa.f90                  | 8 +++++++-
 src/thermal_conductivity/scattering_fourphonon.f90  | 1 -
 src/thermal_conductivity/scattering_isotope.f90     | 6 +++---
 src/thermal_conductivity/scattering_threephonon.f90 | 1 -
 4 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/thermal_conductivity/kappa.f90 b/src/thermal_conductivity/kappa.f90
index d3f44324..cf96d20d 100644
--- a/src/thermal_conductivity/kappa.f90
+++ b/src/thermal_conductivity/kappa.f90
@@ -257,6 +257,7 @@ subroutine get_kappa_offdiag(dr, qp, uc, fc, temperature, classical, mem, mw, ka
                 if (dr%iq(iq)%omega(jmode) .lt. lo_freqtol) cycle
                 om1 = dr%iq(iq)%omega(jmode)
                 tau1 = dr%iq(iq)%linewidth(jmode)
+                if (tau1 .lt. lo_freqtol) cycle
                 do kmode = 1, dr%n_mode
                     if (jmode .eq. kmode) cycle  ! We only want the off diagonal contribution
                     ! Skip gamma for acoustic branches
@@ -264,6 +265,7 @@ subroutine get_kappa_offdiag(dr, qp, uc, fc, temperature, classical, mem, mw, ka
                     if (om2 .lt. lo_freqtol) cycle
 
                     tau2 = dr%iq(iq)%linewidth(kmode)
+                    if (tau2 .lt. lo_freqtol) cycle
 
                     ! This is consistent with the paper, but a bit different from QHGK
                     ! This comes from the fact that we don't work with creation/annihilation but
@@ -412,7 +414,11 @@ subroutine iterative_solution(sr, dr, qp, uc, temperature, niter, tol, classical
             do il = 1, sr%my_nqpoints
                 q1 = sr%my_qpoints(il)
                 b1 = sr%my_modes(il)
-                Fnb(:, b1, q1) = -buf(:, il)/dr%iq(q1)%qs(b1)
+                if (dr%iq(q1)%qs(b1) .lt. lo_freqtol) then
+                    Fnb(:, b1, q1) = 0.0_r8
+                else
+                    Fnb(:, b1, q1) = -buf(:, il)/dr%iq(q1)%qs(b1)
+                end if
             end do
             call mw%allreduce('sum', Fnb)
         end block applyXi
diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90
index 87f4afd1..8433a67d 100644
--- a/src/thermal_conductivity/scattering_fourphonon.f90
+++ b/src/thermal_conductivity/scattering_fourphonon.f90
@@ -341,7 +341,6 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
         d2 = 1.0_r8
         j = j + 1
     else
-
         if (abs(om1 - om2 - om3 - om4) .lt. 4.0_r8 * sigma) then
             d0 = d0 + lo_gauss(om1, om2 + om3 + om4, sigma)
             j = j + 1
diff --git a/src/thermal_conductivity/scattering_isotope.f90 b/src/thermal_conductivity/scattering_isotope.f90
index 2aeb7342..dc692288 100644
--- a/src/thermal_conductivity/scattering_isotope.f90
+++ b/src/thermal_conductivity/scattering_isotope.f90
@@ -31,7 +31,7 @@ subroutine compute_isotope_scattering(il, sr, qp, dr, uc, temperature, &
     ! prefactor and phonon buffers
     real(r8) :: om1, om2, sigma, psisq, prefactor, f0
     ! Integers for do loops
-    integer :: q1, b1, q2, b2, i, niso
+    integer :: q1, b1, q2, b2, i2, niso
 
     q1 = sr%my_qpoints(il)
     b1 = sr%my_modes(il)
@@ -56,7 +56,7 @@ subroutine compute_isotope_scattering(il, sr, qp, dr, uc, temperature, &
                 sigma = sigma + sr%thresh_sigma
             end select
 
-            i = (q2 - 1)*dr%n_mode + b2
+            i2 = (q2 - 1)*dr%n_mode + b2
 
             egviso(:, 2) = dr%aq(q2)%egv(:, b2)
 
@@ -64,7 +64,7 @@ subroutine compute_isotope_scattering(il, sr, qp, dr, uc, temperature, &
 
             f0 = psisq*om1*om2*lo_gauss(om1, om2, sigma)
             g0 = g0 + f0
-            sr%Xi(il, i) = sr%Xi(il, i) + f0*om2/om1
+            sr%Xi(il, i2) = sr%Xi(il, i2) + f0*om2/om1
         end do
     end do
 end subroutine
diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90
index e4cd1762..f71e3433 100644
--- a/src/thermal_conductivity/scattering_threephonon.f90
+++ b/src/thermal_conductivity/scattering_threephonon.f90
@@ -200,7 +200,6 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
         end do allq2
 
         ! And now we add things, with the normalization
-        od_terms = od_terms
         do q2 = 1, qp%n_full_point
             do b2 = 1, dr%n_mode
                 i2 = (q2 - 1)*dr%n_mode + b2

From fa7ca7acc17814671203370d37b230690a1f944b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Fri, 9 May 2025 15:32:07 +0200
Subject: [PATCH 08/14] Some cleaning

---
 src/thermal_conductivity/kappa.f90            |  3 ++
 src/thermal_conductivity/scattering.f90       |  7 ++--
 .../scattering_fourphonon.f90                 | 38 +++++++++---------
 .../scattering_isotope.f90                    |  2 +-
 .../scattering_threephonon.f90                | 39 +++++++++----------
 5 files changed, 44 insertions(+), 45 deletions(-)

diff --git a/src/thermal_conductivity/kappa.f90 b/src/thermal_conductivity/kappa.f90
index cf96d20d..67daa70d 100644
--- a/src/thermal_conductivity/kappa.f90
+++ b/src/thermal_conductivity/kappa.f90
@@ -257,6 +257,7 @@ subroutine get_kappa_offdiag(dr, qp, uc, fc, temperature, classical, mem, mw, ka
                 if (dr%iq(iq)%omega(jmode) .lt. lo_freqtol) cycle
                 om1 = dr%iq(iq)%omega(jmode)
                 tau1 = dr%iq(iq)%linewidth(jmode)
+                ! Not needed almost always, but needed for some extremely rare pathological cases
                 if (tau1 .lt. lo_freqtol) cycle
                 do kmode = 1, dr%n_mode
                     if (jmode .eq. kmode) cycle  ! We only want the off diagonal contribution
@@ -265,6 +266,7 @@ subroutine get_kappa_offdiag(dr, qp, uc, fc, temperature, classical, mem, mw, ka
                     if (om2 .lt. lo_freqtol) cycle
 
                     tau2 = dr%iq(iq)%linewidth(kmode)
+                    ! Not needed almost always, but needed for some extremely rare pathological cases
                     if (tau2 .lt. lo_freqtol) cycle
 
                     ! This is consistent with the paper, but a bit different from QHGK
@@ -414,6 +416,7 @@ subroutine iterative_solution(sr, dr, qp, uc, temperature, niter, tol, classical
             do il = 1, sr%my_nqpoints
                 q1 = sr%my_qpoints(il)
                 b1 = sr%my_modes(il)
+                ! The check is needed for some extremely rare pathological cases
                 if (dr%iq(q1)%qs(b1) .lt. lo_freqtol) then
                     Fnb(:, b1, q1) = 0.0_r8
                 else
diff --git a/src/thermal_conductivity/scattering.f90 b/src/thermal_conductivity/scattering.f90
index 178c8398..7f2c9bc2 100644
--- a/src/thermal_conductivity/scattering.f90
+++ b/src/thermal_conductivity/scattering.f90
@@ -150,21 +150,20 @@ subroutine generate(sr, qp, dr, uc, fct, fcf, opts, tmr, mw, mem)
                 else
                     sr%be(q1, b1) = lo_planck(opts%temperature, dr%iq(q1)%omega(b1))
                 end if
-!               sr%sigsq(q1, b1) = qp%smearingparameter(dr%iq(q1)%vel(:, b1), dr%default_smearing(b1), 1.0_r8)**2
 
                 v0 = matmul(abs(dr%iq(q1)%vel(:, b1)), sr%reclat)**2
                 ! This allows to work around problems with 2D materials
                 sigma = maxval(v0*0.5_r8)
                 sr%sigsq(q1, b1) = sigma
 
-                ! Let's accumulate the average broadening factor for each mode, for a sanity check
+                ! Let's accumulate the average broadening factor for each mode
                 sigavg(b1) = sigavg(b1) + sigma * qp%ip(q1)%integration_weight
             end do
         end do
 
-        ! We add a little baseline as a sanity check, to avoid pathological case when group velocities are near zero
+        ! We add a little baseline, to avoid pathological case when group velocities are near zero
         sr%sigsq = sr%sigsq + maxval(sigavg) * 1e-4_r8
-        ! This is to have the sanity check in memory, for integrationtype 6
+        ! This is to have a sanity check in memory, for integrationtype 6
         sr%thresh_sigma = maxval(sigavg) * 1e-4_r8
         call mem%deallocate(sigavg, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 
diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90
index 8433a67d..310d218b 100644
--- a/src/thermal_conductivity/scattering_fourphonon.f90
+++ b/src/thermal_conductivity/scattering_fourphonon.f90
@@ -32,34 +32,32 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     complex(r8), dimension(:), allocatable :: egv1, egv2, egv3, egv4
     !> Helper for Fourier transform of psi3
     complex(r8), dimension(:), allocatable :: ptf, evp1, evp2, evp3
-    !> The full qpoint grid, to be shuffled
-    integer, dimension(:), allocatable :: qgridfull1, qgridfull2
+    !> Buffer for the off-diagonal terms in the scattering matrix
+    real(r8), dimension(:, :), allocatable :: od_terms
     !> The qpoints in cartesian coordinates
     real(r8), dimension(3) :: qv2, qv3, qv4
+    !> The reducible triplet corresponding to the currently computed quartet
+    integer, dimension(:, :), allocatable :: red_quartet
+    !> The full qpoint grid, to be shuffled
+    integer, dimension(:), allocatable :: qgridfull1, qgridfull2
     !> The complex scattering amplitude
     complex(r8) :: c0
     !> Frequencies, bose-einstein occupation and scattering strength
-    real(r8) :: om1, om2, om3, om4, psisq
-    ! The gaussian integration width
-    real(r8) :: sigma
+    real(r8) :: om1, om2, om3, om4, psisq, sigma
     !> Stuff for the linewidths
     real(r8) :: n2, n3, n4, n2p, n3p, n4p, plf0, plf1, plf2, plf3
-    !> Integers for do loops
-    integer :: q1, q2, q3, q4, q2p, q3p, q4p, b1, b2, b3, b4, qi, qj, i
-    !> Is the quartet irreducible ?
-    logical :: isred
     !> If so, what is its multiplicity
     real(r8) :: mult0, mult1, mult2, mult3
     !> All the prefactors for the scattering
     real(r8) :: f0, f1, f2, f3, f4, f5, f6, f7, fall
-    !> The reducible triplet corresponding to the currently computed quartet
-    integer, dimension(:, :), allocatable :: red_quartet
     !> The approximated dirac
     real(r8) :: d0, d1, d2, d3, d4
-    !> Can we skip the process ?
-    logical :: isok
-
-    real(r8), dimension(:, :), allocatable :: od_terms
+    !> Integers for do loops
+    integer :: q1, q2, q3, q4, q2p, q3p, q4p, b1, b2, b3, b4, qi, qj, i
+    !> Is the quartet irreducible ?
+    logical :: isred
+    !> Is there any scattering happening ?
+    logical :: is_scatter
 
     ! We start by allocating everything
     call mem%allocate(ptf, dr%n_mode**4, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
@@ -150,8 +148,8 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
                     n4p = n4 + 1.0_r8
                     egv4 = dr%aq(q4)%egv(:, b4)/sqrt(om4)
 
-                    call get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, isok)
-                    if (.not. isok) cycle
+                    call get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, is_scatter)
+                    if (.not. is_scatter) cycle
 
                     evp3 = 0.0_r8
                     call zgeru(dr%n_mode, dr%n_mode**3, (1.0_r8, 0.0_r8), egv4, 1, evp2, 1, evp3, dr%n_mode)
@@ -265,7 +263,7 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     if (allocated(red_quartet)) deallocate(red_quartet)
 
     contains
-subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, isok)
+subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, is_scatter)
     !> The scattering amplitudes
     type(lo_scattering_rates), intent(in) :: sr
     !> The qpoint mesh
@@ -281,7 +279,7 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
     !> The approximated dirac
     real(r8), intent(out) :: d0, d1, d2, d3
     !> Can we skip the rest of the calculation ?
-    logical, intent(out) :: isok
+    logical, intent(out) :: is_scatter
 
     !> The frequencies
     real(r8) :: om1, om2, om3, om4
@@ -374,7 +372,7 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
             j = j + 1
         end if
     end if
-    if (j .gt. 0) isok = .true.
+    if (j .gt. 0) is_scatter = .true.
 end subroutine
 end subroutine
 
diff --git a/src/thermal_conductivity/scattering_isotope.f90 b/src/thermal_conductivity/scattering_isotope.f90
index dc692288..2fcd3963 100644
--- a/src/thermal_conductivity/scattering_isotope.f90
+++ b/src/thermal_conductivity/scattering_isotope.f90
@@ -56,10 +56,10 @@ subroutine compute_isotope_scattering(il, sr, qp, dr, uc, temperature, &
                 sigma = sigma + sr%thresh_sigma
             end select
 
+            ! The off-diagonal index in the scattering matrix
             i2 = (q2 - 1)*dr%n_mode + b2
 
             egviso(:, 2) = dr%aq(q2)%egv(:, b2)
-
             psisq = isotope_scattering_strength(uc, egviso)*prefactor
 
             f0 = psisq*om1*om2*lo_gauss(om1, om2, sigma)
diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90
index f71e3433..6c2e063d 100644
--- a/src/thermal_conductivity/scattering_threephonon.f90
+++ b/src/thermal_conductivity/scattering_threephonon.f90
@@ -32,26 +32,24 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     complex(r8), dimension(:), allocatable :: egv1, egv2, egv3
     !> Helper for Fourier transform of psi3
     complex(r8), dimension(:), allocatable :: ptf, evp1, evp2
+    !> buff to keep the off diagonal scattering matrix elements
+    real(r8), dimension(:, :), allocatable :: od_terms
+    !> The reducible triplet corresponding to the currently computed triplet
+    integer, dimension(:, :), allocatable :: red_triplet
     !> The full qpoint grid, to be shuffled
     integer, dimension(:), allocatable :: qgridfull
-    !> The qpoints and the dimension of the qgrid
-    real(r8), dimension(3) :: qv2, qv3
-    !> Frequencies, bose-einstein occupation and scattering strength and some other buffer
-    real(r8) :: sigma, om1, om2, om3, n2, n3, psisq, f0, f1, plf0, plf1, perm
     !> The complex threephonon matrix element
     complex(r8) :: c0
+    !> Frequencies, bose-einstein occupation and scattering strength and some other buffer
+    real(r8) :: sigma, om1, om2, om3, n2, n3, psisq, f0, f1, plf0, plf1, perm
+    !> The approximated dirac
+    real(r8) :: d0, d1
     !> Integers for do loops and counting
     integer :: qi, q1, q2, q3, q2p, q3p, b1, b2, b3, i
     !> Is the triplet irreducible ?
     logical :: isred
-    !> The reducible triplet corresponding to the currently computed triplet
-    integer, dimension(:, :), allocatable :: red_triplet
-    !> buff to keep the off diagonal scattering matrix elements
-    real(r8), dimension(:, :), allocatable :: od_terms
-    !> The approximated dirac
-    real(r8) :: d0, d1
     !> Do we have to compute the scattering process ?
-    logical :: isok
+    logical :: is_scatter
 
     ! We start by allocating everything
     call mem%allocate(ptf, dr%n_mode**3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
@@ -104,9 +102,9 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
                 egv3 = dr%aq(q3)%egv(:, b3)/sqrt(om3)
 
                 ! Get the weight for each processes
-                call get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, isok)
+                call get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, is_scatter)
                 ! If all the weights are zero, skip the last part
-                if (.not. isok) cycle
+                if (.not. is_scatter) cycle
 
                 ! This is the multiplication of eigv of phonons 1 and 2 and now 3
                 evp2 = 0.0_r8
@@ -148,7 +146,7 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
         !> To keep track of the number of mode we actually add
         integer, dimension(dr%n_mode) :: nn
         !> To hold the q-point in reduced coordinates
-        real(r8), dimension(3) :: qv2p
+        real(r8), dimension(3) :: qv2, qv2p
         !> To get the index of the new triplet on the fft_grid
         integer, dimension(3) :: gi
         !> Some buffer
@@ -220,7 +218,7 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     if (allocated(red_triplet)) deallocate (red_triplet)
 
     contains
-subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, isok)
+subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, is_scatter)
     !> The scattering amplitudes
     type(lo_scattering_rates), intent(in) :: sr
     !> The qpoint mesh
@@ -235,8 +233,8 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     integer, intent(in) :: integrationtype
     !> The approximated dirac
     real(r8), intent(out) :: d0, d1
-
-    logical, intent(out) :: isok
+    !> Is there any scattering happening ?
+    logical, intent(out) :: is_scatter
 
     !> The frequencies
     real(r8) :: om1, om2, om3
@@ -265,6 +263,7 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
         do i=1, 3
             sigma = sigma + sqrt(maxval(allsig(:, i)) * 0.5_r8) / 3.0_r8
         end do
+        sigma = sigma + sr%thresh_sigma
     end select
 
     om1 = dr%iq(q1)%omega(b1)
@@ -273,8 +272,8 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     d0 = 0.0_r8
     d1 = 0.0_r8
     j = 0
-    isok = .false.
-    ! We cannot have interaction with same mode (or degenerate)
+    is_scatter = .false.
+    ! We cannot have interactions if two modes are degenerate (or the same)
     if (abs(om1 - om2) .gt. lo_freqtol .and. &
         abs(om1 - om3) .gt. lo_freqtol .and. &
         abs(om2 - om3) .gt. lo_freqtol) then
@@ -296,7 +295,7 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
             j = j + 1
         end if
     end if
-    if (j .gt. 0) isok = .true.
+    if (j .gt. 0) is_scatter = .true.
 end subroutine
 end subroutine
 

From 0689e9effad25ad9727abac25ab77ea7186d4323 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Fri, 9 May 2025 15:44:45 +0200
Subject: [PATCH 09/14] Update test for thermal conductivity

---
 .../reference/outfile.thermal_conductivity                | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/thermal_conductivity/reference/outfile.thermal_conductivity b/tests/thermal_conductivity/reference/outfile.thermal_conductivity
index 64f54c76..14dc96e7 100644
--- a/tests/thermal_conductivity/reference/outfile.thermal_conductivity
+++ b/tests/thermal_conductivity/reference/outfile.thermal_conductivity
@@ -2,13 +2,13 @@
 # Temperature:          0.300000000000E+03
 # Single mode approximation
 #                      kxx                      kyy                      kzz                      kxy                      kxz                      kyz
-        0.487382103035E+02       0.487382103035E+02       0.487382103035E+02       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
+        0.557829741287E+02       0.557829741287E+02       0.557829741287E+02       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
 # Collective contribution
 #                      kxx                      kyy                      kzz                      kxy                      kxz                      kyz
-        0.255832668695E+01       0.255832668695E+01       0.255832668695E+01       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
+        0.241886547641E+01       0.241886547641E+01       0.241886547641E+01       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
 # Off diagonal (coherence) contribution
 #                      kxx                      kyy                      kzz                      kxy                      kxz                      kyz
-        0.215161320755E-01       0.215161320755E-01       0.215161320755E-01       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
+        0.223142292275E-01       0.223142292275E-01       0.223142292275E-01       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
 # Total thermal conductivity
 #                      kxx                      kyy                      kzz                      kxy                      kxz                      kyz
-        0.513180531225E+02       0.513180531225E+02       0.513180531225E+02       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
+        0.582241538344E+02       0.582241538344E+02       0.582241538344E+02       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00

From 868c47250840e5aeb436fcafcb71d566dd5ce9ff Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Fri, 9 May 2025 18:08:13 +0200
Subject: [PATCH 10/14] little fix

---
 src/thermal_conductivity/scattering_fourphonon.f90 | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90
index 310d218b..ae6c2860 100644
--- a/src/thermal_conductivity/scattering_fourphonon.f90
+++ b/src/thermal_conductivity/scattering_fourphonon.f90
@@ -327,14 +327,17 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
     om4 = dr%aq(q4)%omega(b4)
     ! We can have some cases actually respecting strictly the sum rules
     if (abs(om1 - om2) .lt. lo_freqtol .and. abs(om3 - om4) .lt. lo_freqtol) then
+        d0 = 1.0_r8
         d2 = 1.0_r8
         d3 = 1.0_r8
         j = j + 1
     elseif (abs(om1 - om3) .lt. lo_freqtol .and. abs(om2 - om4) .lt. lo_freqtol) then
+        d0 = 1.0_r8
         d1 = 1.0_r8
         d3 = 1.0_r8
         j = j + 1
     elseif (abs(om1 - om4) .lt. lo_freqtol .and. abs(om2 - om3) .lt. lo_freqtol) then
+        d0 = 1.0_r8
         d1 = 1.0_r8
         d2 = 1.0_r8
         j = j + 1

From 4dd5e9e8d0b247b944db158a995c2f6856f18627 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Sat, 10 May 2025 11:38:12 +0200
Subject: [PATCH 11/14] remove useless things

---
 src/thermal_conductivity/kappa.f90 | 57 ------------------------------
 1 file changed, 57 deletions(-)

diff --git a/src/thermal_conductivity/kappa.f90 b/src/thermal_conductivity/kappa.f90
index 67daa70d..71747e38 100644
--- a/src/thermal_conductivity/kappa.f90
+++ b/src/thermal_conductivity/kappa.f90
@@ -22,7 +22,6 @@ module kappa
 public :: get_kappa_offdiag
 public :: iterative_solution
 public :: symmetrize_kappa
-public :: get_viscosity
 contains
 
 !> Calculate the thermal conductivity
@@ -494,60 +493,4 @@ subroutine iterative_solution(sr, dr, qp, uc, temperature, niter, tol, classical
     call mem%deallocate(Fbb, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     call mem%deallocate(buf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
 end subroutine
-
-!> Calculate the thermal conductivity
-subroutine get_viscosity(dr, qp, uc, temperature, classical, mu)
-    !> dispersions
-    type(lo_phonon_dispersions), intent(inout) :: dr
-    !> q-mesh
-    class(lo_qpoint_mesh), intent(in) :: qp
-    !> structure
-    type(lo_crystalstructure), intent(in) :: uc
-    !> temperature
-    real(r8), intent(in) :: temperature
-    !> Are we in the classical limit ?
-    logical, intent(in) :: classical
-    !> thermal conductivity tensor
-    real(r8), dimension(3, 3, 3, 3), intent(out) :: mu
-
-    real(r8), dimension(3) :: v0, v1, kr
-    real(r8) :: om1, cv, pref, n
-    integer :: j
-
-    integer :: q1, b1, iq, a, b, c, d, iop
-    real(r8), dimension(3, 3) :: v2, buf
-
-    mu = 0.0_r8
-    do q1=1, qp%n_full_point
-        iq = qp%ap(q1)%irreducible_index
-        kr = qp%ap(q1)%r * lo_twopi
-        iop = qp%ap(q1)%operation_from_irreducible
-        do b1=1, dr%n_mode
-            om1 = dr%aq(q1)%omega(b1)
-            if (om1 .lt. lo_freqtol) cycle
-            n = lo_planck(temperature, om1)
-            pref = n * (n + 1.0_r8) / uc%volume / lo_kb_hartree / temperature
-            v2 = 0.0_r8
-            if (iop .gt. 0) then
-                v0 = lo_operate_on_vector(uc%sym%op(iop), dr%iq(iq)%Fn(:, b1))
-                v1 = lo_operate_on_vector(uc%sym%op(iop), dr%iq(iq)%vel(:, b1))
-            else
-                v0 = -lo_operate_on_vector(uc%sym%op(abs(iop)), dr%iq(iq)%Fn(:, b1))
-                v1 = -lo_operate_on_vector(uc%sym%op(abs(iop)), dr%iq(iq)%vel(:, b1))
-            end if
-            v2 = lo_outerproduct(v0, v1)
-            do a=1, 3
-            do b=1, 3
-            do c=1, 3
-            do d=1, 3
-                mu(a, b, c, d) = mu(a, b, c, d) + pref * kr(a) * v0(b) * kr(c) * v1(d) / qp%n_full_point
-!               mu(a, b, c, d) = mu(a, b, c, d) + pref * v2(a, b) * kr(c) * kr(d) / qp%n_full_point
-            end do
-            end do
-            end do
-            end do
-        end do
-    end do
-    mu = lo_chop(mu, sum(abs(mu))*1e-6_r8)
-end subroutine
 end module

From f0745e8f9daef22d2768da57a1376ef2767c81e0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Sat, 10 May 2025 11:46:03 +0200
Subject: [PATCH 12/14] remove useless things

---
 src/thermal_conductivity/scattering.f90 | 35 -------------------------
 1 file changed, 35 deletions(-)

diff --git a/src/thermal_conductivity/scattering.f90 b/src/thermal_conductivity/scattering.f90
index 7f2c9bc2..22c567ff 100644
--- a/src/thermal_conductivity/scattering.f90
+++ b/src/thermal_conductivity/scattering.f90
@@ -390,39 +390,4 @@ function sr_size_in_mem(sr) result(mem)
     if (allocated(sr%sigsq)) mem = mem + storage_size(sr%sigsq)*size(sr%sigsq)
     mem = mem/8
 end function
-
-subroutine approximate_hessian(q1, b1, qp, dims, hess)
-    !> The qpoint-mesh
-    type(lo_qpoint_mesh), intent(in) :: qp
-    !> The qpoint
-    integer, intent(in) :: q1
-    !> The mode
-    integer, intent(in) :: b1
-    !> The dimension of the grid
-    integer, dimension(3), intent(in) :: dims
-    !> The Hessian
-    real(r8), dimension(3, 3), intent(out) :: hess
-
-!   fq = qp%ip(q1)%full_index
-!   om1 = dr%iq(q1)%omega(b1)
-!   vel1 = dr%iq(q1)%vel(:, b1)
-!   qtriplet = singlet_to_triplet(q1, dims(2), dims(3))
-
-!   hess = 0.0_r8
-!   do i=1, 3
-!       qt2 = qtriplet
-!       qt2(i) = qt2(i) + 1
-!       if (qt2(i) .gt. dims(i)) qt2(i) = 1
-!       q2 = triplet_to_singlet(qt2, dims(2), dims(3))
-!       vel2 = dr%aq(q2)%vel(:, b2)
-
-!       dq = norm2(qp%ap(q2)%r - qp%ap(fq)%r)
-
-!       do j=1, 3
-!           hess(i, j) = hess(i, j) + (vel2(i) - vel1(i)) / dq
-!       end do
-!   end do
-!   hess = hess / 3.0_r8
-
-end subroutine
 end module

From 93aa84f55dd165a925a0257479988436afc54078 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Sat, 17 May 2025 12:47:47 +0200
Subject: [PATCH 13/14] Final fix for the broadening

---
 .../scattering_fourphonon.f90                 | 74 +++----------------
 .../scattering_threephonon.f90                | 42 ++---------
 2 files changed, 18 insertions(+), 98 deletions(-)

diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90
index ae6c2860..05eedb6b 100644
--- a/src/thermal_conductivity/scattering_fourphonon.f90
+++ b/src/thermal_conductivity/scattering_fourphonon.f90
@@ -56,8 +56,8 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     integer :: q1, q2, q3, q4, q2p, q3p, q4p, b1, b2, b3, b4, qi, qj, i
     !> Is the quartet irreducible ?
     logical :: isred
-    !> Is there any scattering happening ?
-    logical :: is_scatter
+
+    integer, dimension(3) :: qq
 
     ! We start by allocating everything
     call mem%allocate(ptf, dr%n_mode**4, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
@@ -84,6 +84,8 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     call mcg%generate_grid(qgridfull1, rng)
     call mcg%generate_grid(qgridfull2, rng)
 
+    allocate(red_quartet(3, 1))
+
     compute_loop: do qi = 1, mcg%npoints
     do qj = 1, mcg%npoints
         q2 = qgridfull1(qi)
@@ -148,8 +150,7 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
                     n4p = n4 + 1.0_r8
                     egv4 = dr%aq(q4)%egv(:, b4)/sqrt(om4)
 
-                    call get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, is_scatter)
-                    if (.not. is_scatter) cycle
+                    call get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3)
 
                     evp3 = 0.0_r8
                     call zgeru(dr%n_mode, dr%n_mode**3, (1.0_r8, 0.0_r8), egv4, 1, evp2, 1, evp3, dr%n_mode)
@@ -263,7 +264,7 @@ subroutine compute_fourphonon_scattering(il, sr, qp, dr, uc, fcf, mcg, rng, &
     if (allocated(red_quartet)) deallocate(red_quartet)
 
     contains
-subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3, is_scatter)
+subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype, d0, d1, d2, d3)
     !> The scattering amplitudes
     type(lo_scattering_rates), intent(in) :: sr
     !> The qpoint mesh
@@ -278,8 +279,6 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
     integer, intent(in) :: integrationtype
     !> The approximated dirac
     real(r8), intent(out) :: d0, d1, d2, d3
-    !> Can we skip the rest of the calculation ?
-    logical, intent(out) :: is_scatter
 
     !> The frequencies
     real(r8) :: om1, om2, om3, om4
@@ -316,66 +315,15 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, q4, b1, b2, b3, b4, integrationtype
         sigma = sigma + sr%thresh_sigma
     end select
 
-    d0 = 0.0_r8
-    d1 = 0.0_r8
-    d2 = 0.0_r8
-    d3 = 0.0_r8
-
     om1 = dr%iq(q1)%omega(b1)
     om2 = dr%aq(q2)%omega(b2)
     om3 = dr%aq(q3)%omega(b3)
     om4 = dr%aq(q4)%omega(b4)
-    ! We can have some cases actually respecting strictly the sum rules
-    if (abs(om1 - om2) .lt. lo_freqtol .and. abs(om3 - om4) .lt. lo_freqtol) then
-        d0 = 1.0_r8
-        d2 = 1.0_r8
-        d3 = 1.0_r8
-        j = j + 1
-    elseif (abs(om1 - om3) .lt. lo_freqtol .and. abs(om2 - om4) .lt. lo_freqtol) then
-        d0 = 1.0_r8
-        d1 = 1.0_r8
-        d3 = 1.0_r8
-        j = j + 1
-    elseif (abs(om1 - om4) .lt. lo_freqtol .and. abs(om2 - om3) .lt. lo_freqtol) then
-        d0 = 1.0_r8
-        d1 = 1.0_r8
-        d2 = 1.0_r8
-        j = j + 1
-    else
-        if (abs(om1 - om2 - om3 - om4) .lt. 4.0_r8 * sigma) then
-            d0 = d0 + lo_gauss(om1, om2 + om3 + om4, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 + om2 + om3 + om4) .lt. 4.0_r8 * sigma) then
-            d0 = d0 - lo_gauss(om1, -om2 - om3 - om4, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 + om2 - om3 - om4) .lt. 4.0_r8 * sigma) then
-            d1 = d1 + lo_gauss(om1, -om2 + om3 + om4, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 - om2 + om3 + om4) .lt. 4.0_r8 * sigma) then
-            d1 = d1 - lo_gauss(om1, om2 - om3 - om4, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 -om2 + om3 - om4) .lt. 4.0_r8 * sigma) then
-            d2 = d2 + lo_gauss(om1, om2 - om3 + om4, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 + om2 - om3 + om4) .lt. 4.0_r8 * sigma) then
-            d2 = d2 - lo_gauss(om1, -om2 + om3 - om4, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 - om2 - om3 + om4) .lt. 4.0_r8 * sigma) then
-            d3 = d3 + lo_gauss(om1, om2 + om3 - om4, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 + om2 + om3 - om4) .lt. 4.0_r8 * sigma) then
-            d3 = d3 - lo_gauss(om1, -om2 - om3 + om4, sigma)
-            j = j + 1
-        end if
-    end if
-    if (j .gt. 0) is_scatter = .true.
+
+    d0 = lo_gauss(om1, om2 + om3 + om4, sigma) - lo_gauss(om1, -om2 - om3 - om4, sigma)
+    d1 = lo_gauss(om1, -om2 + om3 + om4, sigma) - lo_gauss(om1, om2 - om3 - om4, sigma)
+    d2 = lo_gauss(om1, om2 - om3 + om4, sigma) - lo_gauss(om1, -om2 + om3 - om4, sigma)
+    d3 = lo_gauss(om1, om2 + om3 - om4, sigma) - lo_gauss(om1, -om2 - om3 + om4, sigma)
 end subroutine
 end subroutine
 
diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90
index 6c2e063d..b147b06a 100644
--- a/src/thermal_conductivity/scattering_threephonon.f90
+++ b/src/thermal_conductivity/scattering_threephonon.f90
@@ -48,8 +48,6 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     integer :: qi, q1, q2, q3, q2p, q3p, b1, b2, b3, i
     !> Is the triplet irreducible ?
     logical :: isred
-    !> Do we have to compute the scattering process ?
-    logical :: is_scatter
 
     ! We start by allocating everything
     call mem%allocate(ptf, dr%n_mode**3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
@@ -90,6 +88,7 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
         do b2 = 1, dr%n_mode
             om2 = dr%aq(q2)%omega(b2)
             if (om2 .lt. lo_freqtol) cycle
+            if (abs(om1 - om2) .lt. lo_freqtol) cycle
 
             egv2 = dr%aq(q2)%egv(:, b2)/sqrt(om2)
 
@@ -99,12 +98,11 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
             do b3 = 1, dr%n_mode
                 om3 = dr%aq(q3)%omega(b3)
                 if (om3 .lt. lo_freqtol) cycle
+                if (abs(om1 - om3) .lt. lo_freqtol .or. abs(om2 - om3) .lt. lo_freqtol) cycle
                 egv3 = dr%aq(q3)%egv(:, b3)/sqrt(om3)
 
                 ! Get the weight for each processes
-                call get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, is_scatter)
-                ! If all the weights are zero, skip the last part
-                if (.not. is_scatter) cycle
+                call get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1)
 
                 ! This is the multiplication of eigv of phonons 1 and 2 and now 3
                 evp2 = 0.0_r8
@@ -218,7 +216,7 @@ subroutine compute_threephonon_scattering(il, sr, qp, dr, uc, fct, mcg, rng, &
     if (allocated(red_triplet)) deallocate (red_triplet)
 
     contains
-subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1, is_scatter)
+subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1)
     !> The scattering amplitudes
     type(lo_scattering_rates), intent(in) :: sr
     !> The qpoint mesh
@@ -233,8 +231,6 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     integer, intent(in) :: integrationtype
     !> The approximated dirac
     real(r8), intent(out) :: d0, d1
-    !> Is there any scattering happening ?
-    logical, intent(out) :: is_scatter
 
     !> The frequencies
     real(r8) :: om1, om2, om3
@@ -269,33 +265,9 @@ subroutine get_dirac(sr, qp, dr, q1, q2, q3, b1, b2, b3, integrationtype, d0, d1
     om1 = dr%iq(q1)%omega(b1)
     om2 = dr%aq(q2)%omega(b2)
     om3 = dr%aq(q3)%omega(b3)
-    d0 = 0.0_r8
-    d1 = 0.0_r8
-    j = 0
-    is_scatter = .false.
-    ! We cannot have interactions if two modes are degenerate (or the same)
-    if (abs(om1 - om2) .gt. lo_freqtol .and. &
-        abs(om1 - om3) .gt. lo_freqtol .and. &
-        abs(om2 - om3) .gt. lo_freqtol) then
-
-        if (abs(om1 + om2 - om3) .lt. 4.0_r8 * sigma) then
-            d0 = d0 + lo_gauss(om1, -om2 + om3, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 - om2 + om3) .lt. 4.0_r8 * sigma) then
-            d0 = d0 - lo_gauss(om1, om2 - om3, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 - om2 - om3) .lt. 4.0_r8 * sigma) then
-            d1 = d1 + lo_gauss(om1, om2 + om3, sigma)
-            j = j + 1
-        end if
-        if (abs(om1 + om2 + om3) .lt. 4.0_r8 * sigma) then
-            d1 = d1 - lo_gauss(om1, -om2 - om3, sigma)
-            j = j + 1
-        end if
-    end if
-    if (j .gt. 0) is_scatter = .true.
+
+    d0 = lo_gauss(om1, -om2 + om3, sigma) - lo_gauss(om1, om2 - om3, sigma)
+    d1 = lo_gauss(om1, om2 + om3, sigma) - lo_gauss(om1, -om2 - om3, sigma)
 end subroutine
 end subroutine
 

From 330c9b78d1458cd268dea11e29cce5c34fa49ed1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Alo=C3=AFs=20Castellano?= <alois.castel@gmail.com>
Date: Sat, 17 May 2025 13:12:59 +0200
Subject: [PATCH 14/14] tests for thermal conductivity

---
 .../reference/outfile.thermal_conductivity                | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/thermal_conductivity/reference/outfile.thermal_conductivity b/tests/thermal_conductivity/reference/outfile.thermal_conductivity
index 14dc96e7..db1eac5a 100644
--- a/tests/thermal_conductivity/reference/outfile.thermal_conductivity
+++ b/tests/thermal_conductivity/reference/outfile.thermal_conductivity
@@ -2,13 +2,13 @@
 # Temperature:          0.300000000000E+03
 # Single mode approximation
 #                      kxx                      kyy                      kzz                      kxy                      kxz                      kyz
-        0.557829741287E+02       0.557829741287E+02       0.557829741287E+02       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
+        0.545593928549E+02       0.545593928549E+02       0.545593928549E+02       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
 # Collective contribution
 #                      kxx                      kyy                      kzz                      kxy                      kxz                      kyz
-        0.241886547641E+01       0.241886547641E+01       0.241886547641E+01       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
+        0.332468484893E+01       0.332468484893E+01       0.332468484893E+01       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
 # Off diagonal (coherence) contribution
 #                      kxx                      kyy                      kzz                      kxy                      kxz                      kyz
-        0.223142292275E-01       0.223142292275E-01       0.223142292275E-01       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
+        0.225646351121E-01       0.225646351121E-01       0.225646351121E-01       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
 # Total thermal conductivity
 #                      kxx                      kyy                      kzz                      kxy                      kxz                      kyz
-        0.582241538344E+02       0.582241538344E+02       0.582241538344E+02       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00
+        0.579066423389E+02       0.579066423389E+02       0.579066423389E+02       0.000000000000E+00       0.000000000000E+00       0.000000000000E+00