diff --git a/src/misc.f90 b/src/misc.f90 index 156d22e..30254ed 100644 --- a/src/misc.f90 +++ b/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 diff --git a/src/screening.f90 b/src/screening.f90 index 9d80cea..4734fcb 100644 --- a/src/screening.f90 +++ b/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. diff --git a/test/test_misc.f90 b/test/test_misc.f90 index e363557..5cab0e4 100644 --- a/test/test_misc.f90 +++ b/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*, '<>' @@ -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