diff --git a/src/thermal_conductivity/kappa.f90 b/src/thermal_conductivity/kappa.f90 index bc324c4a..71747e38 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 @@ -256,6 +256,8 @@ 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 @@ -263,6 +265,8 @@ 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 ! 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 diff --git a/src/thermal_conductivity/scattering.f90 b/src/thermal_conductivity/scattering.f90 index 45f1ad1b..22c567ff 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 @@ -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) diff --git a/src/thermal_conductivity/scattering_fourphonon.f90 b/src/thermal_conductivity/scattering_fourphonon.f90 index 32ba10e9..05eedb6b 100644 --- a/src/thermal_conductivity/scattering_fourphonon.f90 +++ b/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] diff --git a/src/thermal_conductivity/scattering_isotope.f90 b/src/thermal_conductivity/scattering_isotope.f90 index f0855cbc..2fcd3963 100644 --- a/src/thermal_conductivity/scattering_isotope.f90 +++ b/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 diff --git a/src/thermal_conductivity/scattering_threephonon.f90 b/src/thermal_conductivity/scattering_threephonon.f90 index f164acdc..b147b06a 100644 --- a/src/thermal_conductivity/scattering_threephonon.f90 +++ b/src/thermal_conductivity/scattering_threephonon.f90 @@ -32,22 +32,22 @@ 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 ! We start by allocating everything call mem%allocate(ptf, dr%n_mode**3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) @@ -88,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) @@ -97,19 +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) - 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) - end select + ! Get the weight for each processes + 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 @@ -126,8 +119,8 @@ 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 ! And we add everything for each triplet equivalent to the one we are actually computing do i = 1, size(red_triplet, 2) @@ -151,7 +144,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 @@ -203,7 +196,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 @@ -221,7 +213,62 @@ 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) + !> 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, 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)) + ! 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 = sigma + sr%thresh_sigma + end select + + om1 = dr%iq(q1)%omega(b1) + om2 = dr%aq(q2)%omega(b2) + om3 = dr%aq(q3)%omega(b3) + + 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 !> Get the Fourier transform of the third order matrix element @@ -350,6 +397,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 diff --git a/tests/thermal_conductivity/reference/outfile.thermal_conductivity b/tests/thermal_conductivity/reference/outfile.thermal_conductivity index 64f54c76..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.487382103035E+02 0.487382103035E+02 0.487382103035E+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.255832668695E+01 0.255832668695E+01 0.255832668695E+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.215161320755E-01 0.215161320755E-01 0.215161320755E-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.513180531225E+02 0.513180531225E+02 0.513180531225E+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