Revert "Wrote unit test for Hilbert transform"

src/misc.f90

@@ -22,7 +22,7 @@ module misc
 #endif
 
   use precision, only: r128, r64, i64
-  use params, only: kB, twopi, pi
+  use params, only: kB, twopi
 
   implicit none
 
@@ -1653,49 +1653,4 @@ 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 for a given function 
-     !! Ref - EQ (4)3, R. Balito et. al.
-     !! - An algorithm for fast Hilbert transform of real
-     !!   functions
-     !!
-     !! 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 :: n, k, nfx
-     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 - 2      ! 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 - 2 - 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

@@ -7,8 +7,7 @@ module screening_module
   use numerics_module, only: numerics
   use misc, only: linspace, mux_vector, binsearch, Fermi, print_message, &
        compsimps, twonorm, write2file_rank2_real, write2file_rank1_real, &
-       distribute_points, sort, qdist, operator(.umklapp.), Bose, &
-       hilbert_transform
+       distribute_points, sort, qdist, operator(.umklapp.), Bose
   use wannier_module, only: wannier
   use delta, only: delta_fn, get_delta_fn_pointer
 
@@ -59,6 +58,53 @@ 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.

test/test_misc.f90

@@ -7,12 +7,12 @@ program test_misc
        twonorm, binsearch, mux_vector, demux_vector, interpolate, coarse_grained, &
        unique, linspace, compsimps, mux_state, demux_state, demux_mesh, expm1, &
        Fermi, Bose, Pade_continued, precompute_interpolation_corners_and_weights, &
-       interpolate_using_precomputed, operator(.umklapp.), shrink, hilbert_transform
+       interpolate_using_precomputed, operator(.umklapp.), shrink
 
   implicit none
 
   integer :: itest
-  integer, parameter :: num_tests = 32
+  integer, parameter :: num_tests = 28
   type(testify) :: test_array(num_tests), tests_all
   integer(i64) :: index, quotient, remainder, int_array(5), v1(3), v2(3), &
        v1_muxed, v2_muxed, ik, ik1, ik2, ik3, ib1, ib2, ib3, wvmesh(3), &
@@ -23,9 +23,6 @@ program test_misc
        real_array(5), result, q1(3, 4), q2(3, 4), q3(3, 4)
   real(r64), allocatable :: integrand(:), domain(:), im_axis(:), real_func(:), &
        widc(:, :), f_coarse(:), f_interp(:), array_of_reals(:)
-  real(r64), allocatable :: hfx1_even(:), hfx1_odd(:), hfx2_even(:), hfx2_odd(:), &
-       ind_even(:), ind_odd(:), x_even(:), x_odd(:), xmin, xmax
-  integer(i64) :: n_even, n_odd  
 
   print*, '<<module misc unit tests>>'
 
@@ -336,88 +333,9 @@ program test_misc
   array_of_reals = [1, 2, 3, 4, 5]*1.0_r64
   call shrink(array_of_reals, 2_i64)
   call test_array(itest)%assert(array_of_reals, [1, 2]*1.0_r64)
-
-  ! Hilbert transform tests (H)
-  ! fx1 -> function 1, fx2 -> function 2
-  ! - hfx1_even stores hilbert transform calculated for fx1, and for even number
-  ! - of points 
-  xmin = -30.0
-  xmax = 30.0
-  n_even = 4000
-  n_odd = 4001
-  ! ind_even are indices to compare in case of even number of points
-  ! ind_odd are indices to compare in case of odd number of points 
-  allocate(ind_even(6),ind_odd(5))
-  ind_even = [801, 1201, 1601, 2001, 2401, 2801]
-  ind_odd = [889, 1333, 1777, 2221, 2665]
-
-  itest = itest + 1
-  test_array(itest) = testify("Hilbert transform: even function, even points")
-  allocate(x_even(n_even), hfx1_even(n_even))
-  call linspace(x_even, xmin, xmax, n_even)
-  call Hilbert_transform(fx1_array(x_even), hfx1_even)
-  call test_array(itest)%assert(hfx1_even(ind_even), hfx1_array(x_even(ind_even)), &
-       tol = 2e-4_r64)
-
-  itest = itest + 1
-  test_array(itest) = testify("Hilbert transform: odd function, even points")
-  allocate(hfx2_even(n_even))
-  call Hilbert_transform(fx2_array(x_even), hfx2_even)
-  call test_array(itest)%assert(hfx2_even(ind_even), hfx2_array(x_even(ind_even)), &
-       tol = 1e-4_r64)
-
-  itest = itest + 1
-  test_array(itest) = testify("Hilbert transform: even function, odd points")
-  allocate(x_odd(n_odd), hfx1_odd(n_odd))
-  call linspace(x_odd, xmin, xmax, n_odd)
-  call Hilbert_transform(fx1_array(x_odd), hfx1_odd)
-  call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_array(x_odd(ind_odd)), &
-       tol = 4e-4_r64)
-
-  itest = itest + 1
-  test_array(itest) = testify("Hilbert transform: odd function, odd points")
-  allocate(hfx2_odd(n_odd))
-  call Hilbert_transform(fx2_array(x_odd), hfx2_odd)
-  call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_array(x_odd(ind_odd)), &
-       tol = 1e-5_r64)
 
   tests_all = testify(test_array)
   call tests_all%report
 
   if(tests_all%get_status() .eqv. .false.) error stop -1
-
-  contains
-     ! reference functions for the Hilbert transform test
-     ! fx1 is an even function
-     function fx1_array(x) result(fx1)
-        real(r64), intent(in) :: x(:)
-        real(r64), allocatable :: fx1(:)
-        allocate(fx1(size(x)))
-        fx1 = 1/(1 + x**2)
-     end function fx1_array
-
-     ! Hfx1 is actual hilbert transform of fx1, is an odd function
-     function hfx1_array(x) result(hfx1)
-        real(r64), intent(in) :: x(:)
-        real(r64), allocatable :: hfx1(:)
-        allocate(hfx1(size(x)))
-        hfx1 = x/(1 + x**2)
-     end function hfx1_array
-
-     ! fx2 is an odd function
-     function fx2_array(x) result(fx2)
-        real(r64), intent(in) :: x(:)
-        real(r64), allocatable :: fx2(:)
-        allocate(fx2(size(x)))
-        fx2 = sin(x)/(1 + x**2)
-     end function fx2_array
-
-     ! Hfx2 is actual hilbert transform of fx2, is an even function
-     function hfx2_array(x) result(hfx2)
-        real(r64), intent(in) :: x(:)
-        real(r64), allocatable :: hfx2(:)
-        real(r64), parameter :: e = 2.718281
-        allocate(hfx2(size(x)))
-        hfx2 = (1/e - cos(x))/(1 + x**2)
-     end function hfx2_array
 end program test_misc