Add a flag in lineshape to remove the omega=0 real part of the self-energy

docs/program/lineshape.md

@@ -86,6 +86,10 @@ Optional switches:
     default value .false.  
     Read the q-point mesh from file. To generate a q-mesh file, see the genkpoints utility.
 
+* `--remove_static_selfenergy`
+    default value .false.
+    Remove the static contribution to the self-energy as in the mode-coupling theory.
+
 * `--help`, `-h`  
     Print this help message
 

src/lineshape/options.f90

@@ -20,6 +20,7 @@ module options
     logical :: isotopescattering = .false.
     logical :: thirdorder = .false.
     logical :: fourthorder = .false.
+    logical :: mct = .false.
     logical :: slightsmearing = .false.
     integer :: integrationtype = -lo_hugeint
     logical :: readiso = .false.
@@ -129,6 +130,10 @@ subroutine parse(opts)
                  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)
     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.', &
+                 required=.false., act='store_true', def='.false.', error=lo_status)
+    if (lo_status .ne. 0) stop
     call cli%add(switch='--nondiagonal', &
                  help='Consider non-diagonal contributions to the self-energy.', hidden=.true., &
                  required=.false., act='store_true', def='.false.', error=lo_status)
@@ -233,6 +238,7 @@ subroutine parse(opts)
     call cli%get(switch='--no_thirdorder_scattering', val=dumlog)
     opts%thirdorder = .not. dumlog
     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)
     opts%diagonal = .not. dumlog
     call cli%get(switch='--minsmear', val=opts%minsmear)
@@ -329,6 +335,12 @@ subroutine parse(opts)
         end if
     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
+        write(*, *) 'There is no real part from the fourth order in the mode-coupling theory'
+        stop
+    end if
+
     ! Not really options
     opts%timereversal = .true.
     opts%slightsmearing = .false.

src/lineshape/phonondamping_generation.f90

@@ -247,10 +247,77 @@ subroutine generate(se, qpoint, qdir, wp, uc, fc, fct, fcf, ise, isf, qp, dr, op
             call mem%deallocate(xs, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
             call mem%deallocate(z0, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
             call mem%deallocate(y0, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+
+            ! If we are in the mode-coupling approximation, we need to remove the static contribution
+            if (opts%mct) then
+                do imode=1, se%n_mode
+                    se%re_3ph(:, imode) = se%re_3ph(:, imode) - se%re_3ph(1, imode)
+                end do
+            end if
         end block kktransform
         call tmr%tock('Kramers-Kronig transformation')
     end if
 
+    ! Correct the memory matrix to get back to the orthogonal time one
+    if (opts%mct .and. se%thirdorder_scattering) then
+        kktransform2: block
+            real(r8) :: pref = 2.0_r8/lo_pi
+            complex(r8), dimension(:), allocatable :: z0
+            complex(r8) :: eta
+            real(r8), dimension(:), allocatable :: x, xs, y0
+            real(r8) :: xp, dlx
+            integer :: imode, ie, ctr
+
+            call mem%allocate(x, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+            call mem%allocate(xs, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+            call mem%allocate(z0, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+            call mem%allocate(y0, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+
+            x = se%energy_axis                  ! copy of x-axis
+            xs = x**2                           ! x^2, precalculated
+            dlx = x(2) - x(1)                     ! prefactor for riemann integral
+            eta = lo_imag*(x(2) - x(1))*1E-8_r8   ! small imaginary thing to avoid divergence
+
+            do imode=1, se%n_mode
+                se%im_3ph(2:, imode) = (x(2:)**2 * se%im_3ph(2:, imode)) / ((x(2:) - se%re_3ph(2:, imode))**2 + se%im_3ph(2:, imode)**2)
+            end do
+
+            ctr = 0
+            se%re_3ph = 0.0_r8
+            do imode = 1, se%n_mode
+            do ie = 1, se%n_energy
+                ctr = ctr + 1
+                if (mod(ctr, mw%n) .ne. mw%r) cycle
+                xp = x(ie)
+                y0 = se%im_3ph(:, imode)*x
+                ! To anyone reading this code:
+                ! The correct way to compute the principal value would be
+                ! z0 = (xp + eta)**2 - xs
+                ! However, it doesn't matter for the Im -> Re transform because we
+                ! don't really hit 0.
+                z0 = xp**2 - xs + eta
+                y0 = real(y0/z0, r8)
+                y0(1) = y0(1)*0.5_r8
+                y0(se%n_energy) = y0(se%n_energy)*0.5_r8
+                se%re_3ph(ie, imode) = sum(y0)*dlx
+            end do
+            end do
+            call mw%allreduce('sum', se%re_3ph)
+            se%re_3ph = se%re_3ph*pref
+
+            call mem%deallocate(x, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+            call mem%deallocate(xs, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+            call mem%deallocate(z0, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+            call mem%deallocate(y0, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
+
+            ! If we are in the mode-coupling approximation, we need to remove the static contribution
+            do imode=1, se%n_mode
+                se%re_3ph(:, imode) = se%re_3ph(:, imode) - se%re_3ph(1, imode)
+            end do
+        end block kktransform2
+        call tmr%tock('Kramers-Kronig transformation')
+    end if
+
     ! Normalize the spectral functions
     normalize: block
         real(r8), parameter :: integraltol = 1E-13_r8