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
+  use params, only: kB, twopi, pi
 
   implicit none
 
@@ -1653,4 +1653,45 @@ 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"
+     !!
+
+     real(r64), intent(in) :: fx(:)
+     real(r64), allocatable, intent(out) :: Hfx(:)
+
+     ! 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
+     ! Note: In the reference, we have N + 1 points and indexing is 0-based
+     ! whereas here we have N points and indexing is 1-based
+     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) = -(fx(k + 2) - fx(k) + term2 + term3)/pi
+     end do
+  end subroutine Hilbert_transform
 end module misc

src/screening.f90

@@ -7,7 +7,8 @@ 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
+       distribute_points, sort, qdist, operator(.umklapp.), Bose, &
+       Hilbert_transform
   use wannier_module, only: wannier
   use delta, only: delta_fn, get_delta_fn_pointer
 
@@ -58,53 +59,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.

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
+       interpolate_using_precomputed, operator(.umklapp.), shrink, Hilbert_transform
 
   implicit none
 
   integer :: itest
-  integer, parameter :: num_tests = 28
+  integer, parameter :: num_tests = 32
   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,6 +23,9 @@ 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>>'
 
@@ -334,8 +337,82 @@ program test_misc
   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: f(x) = 1/(1 + x^2), even points")
+  allocate(x_even(n_even), hfx1_even(n_even))
+  call linspace(x_even, xmin, xmax, n_even)
+  call Hilbert_transform(fx1_el(x_even), hfx1_even)
+  call test_array(itest)%assert(hfx1_even(ind_even), hfx1_el(x_even(ind_even)), &
+       tol = 2e-4_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), even points")
+  allocate(hfx2_even(n_even))
+  call Hilbert_transform(fx2_el(x_even), hfx2_even)
+  call test_array(itest)%assert(hfx2_even(ind_even), hfx2_el(x_even(ind_even)), &
+       tol = 1e-4_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("Hilbert transform: f(x) = 1/(1 + x^2), odd points")
+  allocate(x_odd(n_odd), hfx1_odd(n_odd))
+  call linspace(x_odd, xmin, xmax, n_odd)
+  call Hilbert_transform(fx1_el(x_odd), hfx1_odd)
+  call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_el(x_odd(ind_odd)), &
+       tol = 4e-4_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), odd points")
+  allocate(hfx2_odd(n_odd))
+  call Hilbert_transform(fx2_el(x_odd), hfx2_odd)
+  call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_el(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
+     elemental function fx1_el(x) result(fx1)
+        real(r64), intent(in) :: x
+        real(r64) :: fx1
+        fx1 = 1/(1.0_r64 + x**2)
+     end function fx1_el
+
+     ! Hfx1 is actual hilbert transform of fx1, is an odd function
+     elemental function hfx1_el(x) result(hfx1)
+        real(r64), intent(in) :: x
+        real(r64) :: hfx1
+        hfx1 = x/(1.0_r64 + x**2)
+     end function hfx1_el
+
+     ! fx2 is an odd function
+     elemental function fx2_el(x) result(fx2)
+        real(r64), intent(in) :: x
+        real(r64) :: fx2
+        fx2 = sin(x)/(1.0_r64 + x**2)
+     end function fx2_el
+
+     ! Hfx2 is actual hilbert transform of fx2, is an even function
+     elemental function hfx2_el(x) result(hfx2)
+        real(r64), intent(in) :: x
+        real(r64) :: hfx2
+        hfx2 = (1/exp(1.0_r64) - cos(x))/(1.0_r64 + x**2)
+     end function hfx2_el
 end program test_misc