Moved Hilbert_transform to misc.f90

src/misc.f90

@@ -1653,4 +1653,49 @@ subroutine subtitle(text)
     string2print(75 - length + 1 : 75) = text
     if(this_image() == 1) write(*,'(A75)') string2print
   end subroutine subtitle
+
+  subroutine Hilbert_transform(fx, Hfx)
+     !! Does Hilbert tranform of spectral head of bare polarizability 
+     !! Ref - EQ (4)3, R. Balito et. al. 
+     !!
+     !! fx is the function 
+     !! Hfx is the Hilbert transform of the function
+     !! 
+
+     real(r64), allocatable, intent(out) :: Hfx(:)
+     real(r64), intent(in) :: fx(:)
+
+     ! Local variables
+     integer :: nfx
+
+     integer :: n, k
+     real(r64) :: term2, term3, b
+
+     nfx = size(fx)
+     allocate(Hfx(nfx))
+
+     ! Hilbert function is zero at the edges
+     Hfx(1) = 0.0_r64
+     Hfx(nfx) = 0.0_r64
+
+     ! Both Hfx and fx follow 1-indexing system unlike the paper
+     do k = 1, nfx-1      ! Run over the internal points
+        term2 = 0.0_r64  ! 2nd term in Bilato Eq. 4
+        term3 = 0.0_r64  ! 3rd term in Bilato Eq. 4
+
+        do n = 1, nfx-1-k  ! Partial sum over internal points
+           b = log((n + 1.0_r64)/n)
+           term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*fx(k + n + 1) + &
+                          (1.0_r64 - n*b)*fx(k + n + 2)
+        end do
+
+        do n = 1, k-1      ! Partial sum over internal points
+           b = log((n + 1.0_r64)/n)
+           term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*fx(k - n + 1) - &
+                         (1.0_r64 - n*b)*fx(k - n)
+        end do
+
+        Hfx(k + 1) = (-1.0_r64/pi)*(fx(k + 2) - fx(k) + term2 + term3)
+     end do
+  end subroutine Hilbert_transform
 end module misc

src/screening.f90

@@ -58,53 +58,6 @@ subroutine calculate_qTF(crys, el)
     end if
   end subroutine calculate_qTF
 
-  subroutine Hilbert_transform(f, Hf)
-    !! Does Hilbert tranform of spectral head of bare polarizability 
-    !! Ref - EQ (4), R. Balito et. al. 
-    !!
-    !! Hf H.f(x), the Hilbert transform
-    !! f(x) the function
-
-    real(r64), intent(in) :: f(:)
-    real(r64), allocatable, intent(out) :: Hf(:)
-
-    ! Locals
-    real(r64) :: term2, term3, term1, b
-    integer(i64) :: nxs, k, n
-
-    !Number points on domain grid
-    nxs = size(f)
-
-    allocate(Hf(nxs))
-
-    !Assume that f vanishes at the edges, and Hf also
-    Hf(1) = 0.0_r64
-    Hf(nxs) = 0.0_r64
-
-    do k = 2, nxs - 1 !Run over the internal points
-       term2 = 0.0_r64 !2nd term in Bilato Eq. 4
-       term3 = 0.0_r64 !3rd term in Bilato Eq. 4
-
-       do n = 2, nxs - 1 - k !Partial sum over internal points
-          b = log((n + 1.0_r64)/n)
-
-          term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*f(k + n) + &
-               (1.0_r64 - n*b)*f(k + n + 1)
-       end do
-
-       do n = 2, k - 2 !Partial sum over internal points
-          b = log((n + 1.0_r64)/n)
-
-          term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*f(k - n) - &
-               (1.0_r64 - n*b)*f(k - n - 1)
-       end do
-
-       term1 = f(k + 1) - f(k - 1)
-
-       Hf(k) = (-1.0_r64/pi)*(term1 + term2 + term3)
-    end do
-  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
 !!$    !! Hilbert transform for a given set of temperature-dependent quantities.