Fix broadening for thermal conductivity

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
@@ -256,13 +256,17 @@ 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
                     ! Skip gamma for acoustic branches
                     om2 = dr%iq(iq)%omega(kmode)
                     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
                     ! This comes from the fact that we don't work with creation/annihilation but
@@ -411,7 +415,12 @@ 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)
+                ! 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
+                    Fnb(:, b1, q1) = -buf(:, il)/dr%iq(q1)%qs(b1)
+                end if
             end do
             call mw%allreduce('sum', Fnb)
         end block applyXi

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
@@ -39,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, thresh_sigma
 
 contains
     !> Generate the scattering amplitudes
@@ -78,13 +83,18 @@ 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
         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
 
         ! grid dimensions
         select type (qp)
@@ -116,20 +126,47 @@ 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__)
+
+        sigavg = 0.0_r8
+
         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
+                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), opts%sigma)**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
+                sigavg(b1) = sigavg(b1) + sigma * qp%ip(q1)%integration_weight
             end do
         end do
 
+        ! 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 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__)
+
         if (mw%talk) write (*, *) '... distributing q-point/modes on MPI ranks'
         ctr = 0
         my_nqpoints = 0
@@ -244,7 +281,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)

src/thermal_conductivity/scattering_fourphonon.f90

@@ -32,30 +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
+    !> 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
 
-    real(r8), dimension(:, :), allocatable :: od_terms
+    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__)
@@ -82,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)
@@ -146,18 +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)
 
-                    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 +165,10 @@ 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
+                    f3 = mult3*psisq*plf3*d3
 
                     fall = f0 + f1 + f2 + f3
 
@@ -268,6 +261,70 @@ 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)
+    !> 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, j
+
+    select case (integrationtype)
+    ! Gaussian smearing with fixed width
+    case (1)
+        sigma = sr%sigma
+    ! 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
+        sigma = sigma + sr%thresh_sigma
+    end select
+
+    om1 = dr%iq(q1)%omega(b1)
+    om2 = dr%aq(q2)%omega(b2)
+    om3 = dr%aq(q3)%omega(b3)
+    om4 = dr%aq(q4)%omega(b4)
+
+    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
 
 subroutine quartet_is_irreducible(qp, uc, q1, q2, q3, q4, isred, red_quartet, mw, mem)
@@ -349,6 +406,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]

src/thermal_conductivity/scattering_isotope.f90

@@ -26,10 +26,12 @@ 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
-    integer :: q1, b1, q2, b2, i, niso
+    integer :: q1, b1, q2, b2, i2, niso
 
     q1 = sr%my_qpoints(il)
     b1 = sr%my_modes(il)
@@ -49,19 +51,20 @@ 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 = sigma + sr%thresh_sigma
             end select
 
-            i = (q2 - 1)*dr%n_mode + b2
+            ! 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)
             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