Addressed feedback to the pull request

src/screening.f90

@@ -127,8 +127,9 @@ subroutine calculate_Hilbert_weights(w_disc, w_cont, zeroplus, Hilbert_weights)
     end do
   end subroutine calculate_Hilbert_weights
 
-   subroutine hilbert_transform(Reeps_T, Omegas, spec_eps_T)
-    !! EQ (4), R. Balito et. al. 
+   subroutine Hilbert_transform(Reeps_T, Omegas, spec_eps_T)
+    !! Does Hilbert tranform of spectral head of bare polarizability 
+    !! Ref - EQ (4), R. Balito et. al. 
     !!
     !! Reeps_T Real part of bare polarizability
     !! Omega Energy of excitation in the electron gas
@@ -140,51 +141,54 @@ subroutine hilbert_transform(Reeps_T, Omegas, spec_eps_T)
     real(r64), intent(in) :: spec_eps_T(:)
 
     ! Local variables
-    real(r64) :: in1, in2, fin, bl !, domegas
-    integer(i64) :: nx, kh, np
-
-    nx = size(Omegas)
-    allocate(Reeps_T(nx))
-
-    !domegas = Omegas(2)-Omegas(1)
+    real(r64) :: sum1, sum2, diff, bh
+    integer(i64) :: nOmegas, kh, nh
+    ! Local constant
+    real(r64), parameter :: smalln = 1e-20  ! a number close to zero
 
-    do kh = 1, nx
-        in1=0.0_r64
-        in2=0.0_r64
-        do np=1,nx-1-kh
-            bl = log(real(np+1.0_r64)/real(np))
-            if ( (kh+np).lt.nx) then
-              in1 = in1 - (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh+np) + &
-                          (1.0_r64-np*bl)*spec_eps_T(kh+np+1)
-            else 
-              !print *,"Hilbert transform: Boundary breached",kh+np,nx
-              in1 = in1 - (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh+np) + &
-                      (1.0_r64-np*bl)*1e-12_r64
-            end if
+    nOmegas = size(Omegas)
+    allocate(Reeps_T(nOmegas))
+
+    !
+    do kh = 1, nOmegas
+        sum1 = 0.0_r64
+        sum2 = 0.0_r64
+        do nh = 1, nOmegas - 1 - kh
+           bh = log((nh + 1.0_r64)/nh)
+           ! To avoid edge error in the upper limit
+           if ((kh + nh) < nOmegas) then
+              sum1 = sum1 - (1.0_r64 - (nh + 1.0_r64)*bh)*spec_eps_T(kh + nh) + &
+                          (1.0_r64 - nh*bh)*spec_eps_T(kh + nh + 1)
+           else 
+              ! For outside the upper boundary, assume the function = small number -> 0
+              sum1 = sum1 - (1.0_r64 - (nh + 1.0_r64)*bh)*spec_eps_T(kh + nh) + &
+                      (1.0_r64 - nh*bh)*smalln
+           end if
         end do
 
-        do np=1,kh-1
-           bl = log(real(np+1.0_r64)/real(np))
-           if ((kh-np).gt.1) then
-             in2 = in2 + (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh-np) - &
-                         (1.0_r64-np*bl)*spec_eps_T(kh-np-1)
-            else 
-              !print *,"Hilbert transform: Boundary breached",kh+np,nx
-              in2 = in2 + (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh-np) - &
-                     (1.0_r64-np*bl)*1e-12_r64
-            end if
+        do nh = 1, kh - 1
+           bh = log((nh + 1.0_r64)/nh)
+           ! To avoid edge error in the lower limit
+           if ((kh - nh) > 1) then
+              sum2 = sum2 + (1.0_r64 - (nh + 1.0_r64)*bh)*spec_eps_T(kh - nh) - &
+                         (1.0_r64 - nh*bh)*spec_eps_T(kh - nh - 1)
+           else 
+              ! For outside the lower boundary, assume the function = small number -> 0
+              sum2 = sum2 + (1.0_r64 - (nh + 1.0_r64)*bh)*spec_eps_T(kh - nh) - &
+                     (1.0_r64 - nh*bh)*smalln
+           end if
         end do
 
-        if (kh.ge.2 .and. kh.le.(nx-1)) then
-           fin = spec_eps_T(kh+1)-spec_eps_T(kh-1)
-         else
-            fin = 1e-12_r64
-         end if
-        Reeps_T(kh) = (-1.0_r64/pi)*(fin+in1+in2)
-
+        ! To avoid edge error
+        if (kh >= 2 .and. kh <= (nOmegas-1)) then
+           diff = spec_eps_T(kh + 1) - spec_eps_T(kh - 1)
+        else
+        ! For outside the boundary
+            diff = smalln
+        end if
+        Reeps_T(kh) = (-1.0_r64/pi)*(diff + sum1 + sum2)
     end do
-
-  end subroutine hilbert_transform 
+  end subroutine Hilbert_transform 
 
   subroutine head_polarizability_real_3d_T(Reeps_T, Omegas, spec_eps_T, Hilbert_weights_T)
     !! Head of the bare real polarizability of the 3d Kohn-Sham system using