[DRAFT] Fourthorder lineshape

src/libolle/mpi_wrappers.f90

@@ -54,8 +54,8 @@ module mpi_wrappers
         !> mpi allreduce
         generic :: allreduce=>allreduce_int,allreduce_real,allreduce_1d_real,allreduce_2d_real,allreduce_3d_real,&
                               allreduce_4d_real,allreduce_5d_real,allreduce_6d_real,allreduce_1d_int,allreduce_2d_int,allreduce_3d_int,&
-                              allreduce_2d_complex,allreduce_3d_complex,allreduce_4d_complex,allreduce_5d_complex,allreduce_1d_logical,&
-                              allreduce_int8
+                              allreduce_2d_complex,allreduce_3d_complex,allreduce_4d_complex,allreduce_5d_complex,allreduce_6d_complex,&
+                              allreduce_1d_logical,allreduce_int8
         procedure, private :: allreduce_int
         procedure, private :: allreduce_1d_int
         procedure, private :: allreduce_2d_int
@@ -71,6 +71,7 @@ module mpi_wrappers
         procedure, private :: allreduce_3d_complex
         procedure, private :: allreduce_4d_complex
         procedure, private :: allreduce_5d_complex
+        procedure, private :: allreduce_6d_complex
         procedure, private :: allreduce_1d_logical
         procedure, private :: allreduce_int8
         !> mpi broadcast
@@ -1442,6 +1443,42 @@ subroutine allreduce_5d_complex(mw,operation,x,y,filename,linenumber)
         endif
     endif
 end subroutine
+subroutine allreduce_6d_complex(mw,operation,x,y,filename,linenumber)
+    !> MPI helper
+    class(lo_mpi_helper), intent(inout) :: mw
+    !> what operation to do
+    character(len=*), intent(in) :: operation
+    !> array to allreduce
+    complex(r8), dimension(:,:,:,:,:,:), intent(inout) :: x
+    !> array to allreduce to, if omitted default to in-place
+    complex(r8), dimension(:,:,:,:,:,:), intent(out), optional :: y
+    !> filename we call from
+    character(len=*), intent(in), optional :: filename
+    !> line number we call from
+    integer, intent(in), optional :: linenumber
+
+    integer :: mpiop
+    ! Fetch proper operation code
+    mpiop=mpi_operation_code(operation)
+    ! Do the actual communication
+    if ( present(y) ) then
+        if ( size(x) .ne. size(y) ) then
+            call lo_stop_gracefully(['mpi_allreduce inconsistent array sizes'],lo_exitcode_mpi,filename,linenumber,mw%comm)
+        endif
+        y=0.0_r8
+        call mpi_allreduce(x,y,size(x),MPI_DOUBLE_COMPLEX,mpiop,mw%comm,mw%error)
+    else
+        call mpi_allreduce(MPI_IN_PLACE,x,size(x),MPI_DOUBLE_COMPLEX,mpiop,mw%comm,mw%error)
+    endif
+    ! Check that things went ok
+    if ( mw%error .ne. 0 ) then
+        if ( present(filename) .and. present(linenumber) ) then
+            call lo_stop_gracefully(['mpi_allreduce exit code '//tochar(mw%error)],lo_exitcode_mpi,filename,linenumber,mw%comm)
+        else
+            call lo_stop_gracefully(['mpi_allreduce exit code '//tochar(mw%error)],lo_exitcode_mpi,__FILE__,__LINE__,mw%comm)
+        endif
+    endif
+end subroutine
 subroutine allreduce_1d_logical(mw,operation,i,j,filename,linenumber)
     !> MPI helper
     class(lo_mpi_helper), intent(inout) :: mw

src/lineshape/dielscatter.f90

@@ -248,7 +248,7 @@ subroutine generate(dir, wp, di, p, qp, dr, fc, fct, se, isf, opts, tmr, mw, mem
                 ! fetch probing energy and real and imaginary part of the self-energy
                 Z = se%energy_axis(ie)
                 dl = se%re_3ph(ie, imode) + se%re_4ph(ie, imode)
-                gm = se%im_3ph(ie, imode) + se%im_iso(ie, imode)
+                gm = se%im_3ph(ie, imode) + se%im_4ph(ie, imode) + se%im_iso(ie, imode)
 
                 ! Zeroth order IR, selfenergy=0
                 c0 = diagonal_greensfunction(om, Z, 0.0_r8, smallsmearing)
@@ -346,7 +346,7 @@ subroutine generate(dir, wp, di, p, qp, dr, fc, fct, se, isf, opts, tmr, mw, mem
                 ! fetch probing energy and real and imaginary part of the self-energy
                 Z = se%energy_axis(ie)
                 dl = se%re_3ph(ie, b1) + se%re_4ph(ie, b1)
-                gm = se%im_3ph(ie, b1) + se%im_iso(ie, b1)
+                gm = se%im_3ph(ie, b1) + se%im_4ph(ie, b1) + se%im_iso(ie, b1)
                 grfun(ie, b1) = aimag(diagonal_greensfunction(om1, Z, dl, gm))
             end do
         end do

src/lineshape/io.f90

@@ -161,6 +161,7 @@ subroutine write_lineshape_to_hdf5(se, p, qpt, wp, dir, di, qp, enhet, temperatu
         call h5%store_data(se%energy_axis*unitfactor, h5%group_id, 'frequency', enhet=trim(unitname))
         call h5%store_data(se%im_iso*unitfactor, h5%group_id, 'imaginary_isotope_selfenergy', enhet=trim(unitname))
         call h5%store_data(se%im_3ph*unitfactor, h5%group_id, 'imaginary_threephonon_selfenergy', enhet=trim(unitname))
+        call h5%store_data(se%im_4ph*unitfactor, h5%group_id, 'imaginary_fourphonon_selfenergy', enhet=trim(unitname))
         call h5%store_data(se%re_3ph*unitfactor, h5%group_id, 'real_threephonon_selfenergy', enhet=trim(unitname))
         call h5%store_data(se%re_4ph*unitfactor, h5%group_id, 'real_fourphonon_selfenergy', enhet=trim(unitname))
 
@@ -179,7 +180,7 @@ subroutine write_lineshape_to_hdf5(se, p, qpt, wp, dir, di, qp, enhet, temperatu
         ! Imaginary part of Green's function
         do i = 1, se%n_mode
             if (wp%omega(i) .gt. lo_freqtol*10) then
-                call evaluate_spectral_function(se%energy_axis, se%im_3ph(:, i) + se%im_iso(:, i), se%re_3ph(:, i) + se%re_4ph(:, i), wp%omega(i), dr0(:, i))
+                call evaluate_spectral_function(se%energy_axis, se%im_3ph(:, i) + se%im_4ph(:, i) + se%im_iso(:, i), se%re_3ph(:, i) + se%re_4ph(:, i), wp%omega(i), dr0(:, i))
             else
                 dr0(:, i) = 0.0_r8
             end if

src/lineshape/main.f90

@@ -72,7 +72,7 @@ program lineshape
         call fct%readfromfile(uc, 'infile.forceconstant_thirdorder')
         if (mw%talk) write (*, *) '... read third order forceconstant'
     end if
-    if (opts%fourthorder) then
+    if (opts%fourthorder_real .or. opts%fourthorder) then
         call fcf%readfromfile(uc, 'infile.forceconstant_fourthorder')
         if (mw%talk) write (*, *) '... read fourth order forceconstant'
     end if

src/lineshape/options.f90

@@ -20,6 +20,7 @@ module options
     logical :: isotopescattering = .false.
     logical :: thirdorder = .false.
     logical :: fourthorder = .false.
+    logical :: fourthorder_real = .false.
     logical :: mct = .false.
     logical :: slightsmearing = .false.
     integer :: integrationtype = -lo_hugeint
@@ -126,9 +127,12 @@ subroutine parse(opts)
                  help='Switch of three-phonon scattering', &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
-    call cli%add(switch='--fourthorder', &
+    call cli%add(switch='--fourthorder_real', &
                  help='Consider four-phonon contributions to the real part of the self-energy.', hidden=.true., &
                  required=.false., act='store_true', def='.false.', error=lo_status)
+    call cli%add(switch='--fourthorder', &
+                 help='Consider four-phonon contributions to the dynamic part of the self-energy.', hidden=.true., &
+                 required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
     call cli%add(switch='--remove_static_selfenergy', &
                  help='Remove the static contribution to the self-energy as in the mode-coupling theory.', &
@@ -237,6 +241,7 @@ subroutine parse(opts)
     opts%isotopescattering = .not. dumlog
     call cli%get(switch='--no_thirdorder_scattering', val=dumlog)
     opts%thirdorder = .not. dumlog
+    call cli%get(switch='--fourthorder_real', val=opts%fourthorder_real)
     call cli%get(switch='--fourthorder', val=opts%fourthorder)
     call cli%get(switch='--remove_static_selfenergy', val=opts%mct)
     call cli%get(switch='--nondiagonal', val=dumlog)
@@ -336,7 +341,7 @@ subroutine parse(opts)
     end if
 
     ! If we are in the mode-coupling approach, the real part four phonon makes no sense
-    if (opts%mct .and. opts%fourthorder) then
+    if (opts%mct .and. opts%fourthorder_real) then
         write(*, *) 'There is no real part from the fourth order in the mode-coupling theory'
         stop
     end if

src/lineshape/phonondamping.f90

@@ -41,7 +41,7 @@ module phonondamping
     ! energy axis for the self-energy
     real(r8), dimension(:), allocatable :: energy_axis
     ! Imaginary part of self energy (energy,mode)
-    real(r8), dimension(:, :), allocatable :: im_3ph, im_iso
+    real(r8), dimension(:, :), allocatable :: im_3ph, im_iso, im_4ph
     ! Real part of self energy (energy,mode)
     real(r8), dimension(:, :), allocatable :: re_3ph, re_4ph
     !> Direction of probe
@@ -54,6 +54,7 @@ module phonondamping
     logical :: isotope_scattering = .false.
     logical :: thirdorder_scattering = .false.
     logical :: fourthorder_scattering = .false.
+    logical :: fourthorder_real = .false.
     logical :: diagonal = .false.
     logical :: skipsym = .false.
     ! auxiliary information about lineshapes:
@@ -75,7 +76,8 @@ module phonondamping
 
 ! Prefactors, to make sure they are the same all the time.
 real(r8), parameter :: threephonon_prefactor = lo_pi/16.0_r8
-real(r8), parameter :: fourphonon_prefactor = 1.0_r8/8.0_r8
+real(r8), parameter :: fourphonon_real_prefactor = 1.0_r8/8.0_r8
+real(r8), parameter :: fourphonon_imag_prefactor = lo_pi/96.0_r8
 real(r8), parameter :: isotope_prefactor = lo_pi/4.0_r8
 
 ! dos interfaces
@@ -199,6 +201,7 @@ function se_size_in_mem(se) result(mem)
     if (allocated(se%im_3ph)) mem = mem + storage_size(se%im_3ph)*size(se%im_3ph)
     if (allocated(se%im_iso)) mem = mem + storage_size(se%im_iso)*size(se%im_iso)
     if (allocated(se%re_3ph)) mem = mem + storage_size(se%re_3ph)*size(se%re_3ph)
+    if (allocated(se%im_4ph)) mem = mem + storage_size(se%im_4ph)*size(se%im_4ph)
     if (allocated(se%re_4ph)) mem = mem + storage_size(se%re_4ph)*size(se%re_4ph)
     if (allocated(se%xmid)) mem = mem + storage_size(se%xmid)*size(se%xmid)
     if (allocated(se%xlo)) mem = mem + storage_size(se%xlo)*size(se%xlo)

src/lineshape/phonondamping_dos.f90

@@ -146,7 +146,7 @@ module subroutine get_intensity_as_dos(pd, qpd, drd, uc, fc, fct, fcf, ise, sf,
                 if (drd%iq(iq)%omega(imode) .lt. lo_freqtol) cycle
                 ! Get the spectral function
                 buf1 = 0.0_r8
-                call evaluate_spectral_function(se%energy_axis, se%im_3ph(:, imode) + se%im_iso(:, imode), se%re_3ph(:, imode) + se%re_4ph(:, imode), drd%iq(iq)%omega(imode), buf1)
+                call evaluate_spectral_function(se%energy_axis, se%im_3ph(:, imode) + se%im_4ph(:, imode) + se%im_iso(:, imode), se%re_3ph(:, imode) + se%re_4ph(:, imode), drd%iq(iq)%omega(imode), buf1)
                 buf1 = buf1/lo_trapezoid_integration(se%energy_axis, buf1)
 
                 ! Stash the unsmeared spectral function
@@ -173,7 +173,9 @@ module subroutine get_intensity_as_dos(pd, qpd, drd, uc, fc, fct, fcf, ise, sf,
                 ! Store shifts and linewidth
                 buf_shift_3rd(imode, iq) = lo_linear_interpolation(se%energy_axis, se%re_3ph(:, imode), drd%iq(iq)%omega(imode))
                 buf_shift_4th(imode, iq) = lo_linear_interpolation(se%energy_axis, se%re_4ph(:, imode), drd%iq(iq)%omega(imode))
-                buf_linewidth(imode, iq) = lo_linear_interpolation(se%energy_axis, se%im_3ph(:, imode), drd%iq(iq)%omega(imode)) + lo_linear_interpolation(se%energy_axis, se%im_iso(:, imode), drd%iq(iq)%omega(imode))
+                buf_linewidth(imode, iq) = lo_linear_interpolation(se%energy_axis, se%im_3ph(:, imode), drd%iq(iq)%omega(imode)) + &
+                                           lo_linear_interpolation(se%energy_axis, se%im_iso(:, imode), drd%iq(iq)%omega(imode)) + &
+                                           lo_linear_interpolation(se%energy_axis, se%im_4ph(:, imode), drd%iq(iq)%omega(imode))
             end do
 
             ! Store spectral function

src/lineshape/phonondamping_fourthorder.f90

@@ -13,7 +13,7 @@ subroutine fourphonon_selfenergy(qpoint, wp, gp, qp, uc, temperature, dr, fcf, d
     type(lo_crystalstructure), intent(in) :: uc
     !> harmonic properties for q',q'',q''
     type(lo_phonon_dispersions), intent(in) :: dr
-    !> third order force constant
+    !> fourth order force constant
     type(lo_forceconstant_fourthorder), intent(in) :: fcf
     !> temperature
     real(r8), intent(in) :: temperature
@@ -238,7 +238,7 @@ subroutine fourphonon_selfenergy(qpoint, wp, gp, qp, uc, temperature, dr, fcf, d
                 ctr = ctr + 1
                 if (mod(ctr, mw%n) .ne. mw%r) cycle
 
-                prefactor = qp%ip(iq)%integration_weight*fourphonon_prefactor
+                prefactor = qp%ip(iq)%integration_weight*fourphonon_real_prefactor
                 ! pre-transform the matrix element
                 call pretransform_phi4(fcf, qpoint%r, qp%ip(iq)%r, ptf)
 
@@ -278,7 +278,7 @@ subroutine fourphonon_selfenergy(qpoint, wp, gp, qp, uc, temperature, dr, fcf, d
                 ctr = ctr + 1
                 if (mod(ctr, mw%n) .ne. mw%r) cycle
 
-                prefactor = qp%ap(iq)%integration_weight*fourphonon_prefactor
+                prefactor = qp%ap(iq)%integration_weight*fourphonon_real_prefactor
                 ! pre-transform the matrix element
                 call pretransform_phi4(fcf, qpoint%r, qp%ap(iq)%r, ptf)
 
@@ -432,7 +432,7 @@ subroutine fourphonon_selfenergy(qpoint, wp, gp, qp, uc, temperature, dr, fcf, d
     !             if ( minval(omega) .gt. omegathres ) then
     !                 f0=real(fcf%scatteringamplitude(omega,egv,-qv2,-qv3,-qv4),r8)
     !                 f0=f0*(1.0_r8+lo_planck(temperature,omega(3))+lo_planck(temperature,omega(4)) )
-    !                 f0=f0*qp%ap(q)%integration_weight*fourphonon_prefactor
+    !                 f0=f0*qp%ap(q)%integration_weight*fourphonon_real_prefactor
     !             else
     !                 f0=0.0_r8
     !             endif

src/lineshape/phonondamping_gaussian.f90

@@ -372,3 +372,143 @@ subroutine threephonon_imaginary_selfenergy_gaussian(wp, se, sr, qp, dr, tempera
     call mem%deallocate(buf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
     if (verbosity .gt. 0) call lo_progressbar(' ... threephonon imaginary selfenergy', dr%n_mode*qp%n_full_point, dr%n_mode*qp%n_full_point, walltime() - t0)
 end subroutine
+
+
+!> fourphonon self energy gaussian
+subroutine fourphonon_imaginary_selfenergy_gaussian(wp, se, sr, qp, dr, temperature, mw, mem, verbosity)
+    !> harmonic properties at this q-point
+    type(lo_phonon_dispersions_qpoint), intent(in) :: wp
+    !> self-energy
+    type(lo_phonon_selfenergy), intent(inout) :: se
+    !> scattering rates
+    type(lo_listofscatteringrates), intent(inout) :: sr
+    !> qpoint mesh
+    class(lo_qpoint_mesh), intent(in) :: qp
+    !> harmonic dispersions
+    type(lo_phonon_dispersions), intent(in) :: dr
+    !> temperature
+    real(r8), intent(in) :: temperature
+    !> MPI helper
+    type(lo_mpi_helper), intent(inout) :: mw
+    !> memory tracker
+    type(lo_mem_helper), intent(inout) :: mem
+    !> talk a lot?
+    integer, intent(in) :: verbosity
+
+    real(r8), dimension(:), allocatable :: buf
+    real(r8) :: psisquare
+    real(r8) :: sigma, dom, om2, om3, om4, invf, s2, s3, s4
+    real(r8) :: plf1, plf2, t0, pref
+    real(r8) :: n2, n3, n4
+    integer :: q2, q3, ctr, i, ii, jj, b1, b2, b3, b4, ilo, ihi
+
+    ! set the unit
+    t0 = walltime()
+
+    invf = se%n_energy/se%energy_axis(se%n_energy)
+
+    ctr = 0
+    se%im_4ph = 0.0_r8
+    call mem%allocate(buf, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+    buf = 0.0_r8
+
+    qloopfull: do q2 = 1, qp%n_full_point
+    do q3=1, qp%n_full_point
+        ! Prefactor for this q-point
+        pref = fourphonon_imag_prefactor * qp%ap(q2)%integration_weight * qp%ap(q3)%integration_weight
+        do b1 = 1, dr%n_mode
+            do b2 = 1, dr%n_mode
+                ! Make it parallel
+                ctr = ctr + 1
+                if (mod(ctr, mw%n) .ne. mw%r) cycle
+
+                do b3 = 1, dr%n_mode
+                    do b4 = 1, dr%n_mode
+                        ! reset buffer
+                        buf = 0.0_r8
+                        ilo = lo_hugeint
+                        ihi = -lo_hugeint
+                        ! Get the smearing parameter, in case it is adaptive
+                        select case (se%integrationtype)
+                        case (1)
+                            s2 = dr%default_smearing(b2)
+                            s3 = dr%default_smearing(b3)
+                            s4 = dr%default_smearing(b4)
+                            sigma = sqrt(s2**2 + s3**2 + s4**2)
+                        case (2)
+                            s2 = qp%adaptive_sigma(qp%ap(q2)%radius, dr%aq(q2)%vel(:, b2), dr%default_smearing(b2), se%smearing_prefactor)
+                            s3 = qp%adaptive_sigma(qp%ap(q3)%radius, dr%aq(q3)%vel(:, b3), dr%default_smearing(b3), se%smearing_prefactor)
+                            s4 = qp%adaptive_sigma(qp%ap(q3)%radius, sr%vel4(:, b4, q2, q3), dr%default_smearing(b4), se%smearing_prefactor)
+                            sigma = sqrt(s2**2 + s3**2 + s4**2)
+                        end select
+                        ! Fetch frequencies
+                        om2 = dr%aq(q2)%omega(b2)
+                        om3 = dr%aq(q3)%omega(b3)
+                        om4 = sr%omega4(b4, q2, q3)
+
+                        ! Define occupation numbers
+                        n2 = lo_planck(temperature, om2)
+                        n3 = lo_planck(temperature, om3)
+                        n4 = lo_planck(temperature, om4)
+
+                        plf1 = (n2 + 1) * (n3 + 1) * (n4 + 1)
+                        plf2 = n2 * n3 * n4
+                        dom = om2 + om3 + om4
+                        ! Delta(BigOmega - om2 - om3 - om4)
+                        ii = max(floor((dom - 4 * sigma)*invf), 1)
+                        jj = min(ceiling((dom + 4 * sigma) * invf), se%n_energy)
+                        ilo = min(ilo, ii)
+                        ihi = max(ihi, jj)
+                        do i = ii, jj
+                            buf(i) = buf(i) + (plf1 - plf2) * lo_gauss(se%energy_axis(i), dom, sigma)
+                        end do
+
+                        plf1 = n2 * (n3 + 1) * (n4 + 1)
+                        plf2 = (n2 + 1) * n3 * n4
+                        dom = -om2 + om3 + om4
+                        ! Delta(BigOmega - om2 + om3 + om4)
+                        ii = max(floor((dom - 4 * sigma)*invf), 1)
+                        jj = min(ceiling((dom + 4 * sigma)*invf), se%n_energy)
+                        ilo = min(ilo, ii)
+                        ihi = max(ihi, jj)
+                        do i = ii, jj
+                            buf(i) = buf(i) + 3 * (plf1 - plf2) * lo_gauss(se%energy_axis(i), dom, sigma)
+                        end do
+
+                        ! Increment the self-energy
+                        if (ilo .lt. ihi) then
+                            psisquare = abs(sr%psi_4ph(b1, b2, b3, b4, q2, q3) * conjg(sr%psi_4ph(b1, b2, b3, b4, q2, q3))) * pref
+                            se%im_4ph(ilo:ihi, b1) = se%im_4ph(ilo:ihi, b1) + buf(ilo:ihi)*psisquare
+                        end if
+                    end do
+                end do
+            end do
+        end do
+        if (verbosity .gt. 0) then
+            if (lo_trueNtimes(ctr, 127, qp%n_full_point*qp%n_full_point*dr%n_mode*dr%n_mode)) then
+                call lo_progressbar(' ... fourphonon imaginary selfenergy', ctr, dr%n_mode*dr%n_mode*qp%n_full_point*qp%n_full_point)
+            end if
+        end if
+    end do
+    end do qloopfull
+
+    ! Sum it up
+    call mw%allreduce('sum', se%im_4ph)
+
+    ! Fix degeneracies ?
+    do b1 = 1, dr%n_mode
+        buf = 0.0_r8
+        do i = 1, wp%degeneracy(b1)
+            b2 = wp%degenmode(i, b1)
+            buf = buf + se%im_4ph(:, b2)
+        end do
+        buf = buf / real(wp%degeneracy(b1), r8)
+        do i = 1, wp%degeneracy(b1)
+            b2 = wp%degenmode(i, b1)
+            se%im_4ph(:, b2) = buf
+        end do
+    end do
+    call mem%deallocate(buf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+    if (verbosity .gt. 0) call lo_progressbar(' ... fourphonon imaginary selfenergy', dr%n_mode*qp%n_full_point,&
+    dr%n_mode*qp%n_full_point, walltime() - t0)
+end subroutine