diff --git a/GX-NAOIntegrals/CMakeLists.txt b/GX-NAOIntegrals/CMakeLists.txt index 2165e65..f8481b0 100644 --- a/GX-NAOIntegrals/CMakeLists.txt +++ b/GX-NAOIntegrals/CMakeLists.txt @@ -24,6 +24,7 @@ target_sources(LibGXNAO PRIVATE src/gaunt.f90 src/spline.f90 src/gauss_quadrature.f90 + src/legendre_polynomial.f90 ) target_link_libraries(LibGXNAO GXCommon) diff --git a/GX-NAOIntegrals/src/legendre_polynomial.f90 b/GX-NAOIntegrals/src/legendre_polynomial.f90 new file mode 100644 index 0000000..155d663 --- /dev/null +++ b/GX-NAOIntegrals/src/legendre_polynomial.f90 @@ -0,0 +1,66 @@ +!> @brief evaluation of standard legendre polinomials +!! +!! using the recurrence relation +!! P(0,x) = 1 +!! P(1,x) = x +!! P(n,x) = (2*n-1)/n * x * P(n-1,x) - (n-1)/n * P(n-2,x) +!! +!! see https://people.sc.fsu.edu/~jburkardt/f_src/legendre_polynomial/legendre_polynomial.html +module legendre_polynomial + + use kinds, only: dp + + implicit none + + private + public :: evaluate_legendre_polinomial, evaluate_legendre_polinomial_batch + + contains + + !> @brief evaluate all legegendre pol. up to dregree n at point x + !! + !! @param[in] n -- degree of legendre polynomial + !! @param[in] x -- point + !! @param[out] p -- all legendre polynomials P_0, ..., P_n + subroutine evaluate_legendre_polinomial(n, x, p) + integer, intent(in) :: n + real(kind=dp), intent(in) :: x + real(kind=dp), dimension(n+1), intent(out) :: p + + ! internal variables + integer :: i + + p(1) = 1.0_dp + p(2) = x + do i = 2, n + p(i+1) = (2.0_dp*n-1.0_dp)/dble(n) * x * p(i) - (n-1.0_dp)/dble(n) * p(i-1) + end do + + end subroutine + + + + !> @brief evaluate all legegendre pol. up to dregree n at all points x + !! + !! @param[in] n -- degree of legendre polynomial + !! @param[in] n_points -- number of points in batch + !! @param[in] x -- batch of points + !! @param[out] p -- all legendre polynomials P_0, ..., P_n at pints x + subroutine evaluate_legendre_polinomial_batch(n, n_points, x, p) + integer, intent(in) :: n + integer, intent(in) :: n_points + real(kind=dp), dimension(n_points), intent(in) :: x + real(kind=dp), dimension(n_points, n+1), intent(out) :: p + + ! internal variables + integer :: i + + p(:, 1) = 1.0_dp + p(:, 2) = x(:) + do i = 2, n + p(:, i+1) = (2.0_dp*n-1.0_dp)/dble(n) * x(:) * p(:, i) - (n-1.0_dp)/dble(n) * p(:, i-1) + end do + + end subroutine + +end module \ No newline at end of file diff --git a/GX-NAOIntegrals/test/test_nao.f90 b/GX-NAOIntegrals/test/test_nao.f90 index e043916..0337cd9 100644 --- a/GX-NAOIntegrals/test/test_nao.f90 +++ b/GX-NAOIntegrals/test/test_nao.f90 @@ -4,6 +4,7 @@ program test use gaunt, only: r_gaunt use spline, only: cubic_spline use gauss_quadrature, only: get_gauss_legendre_grid, gauss_legendre_integrator + use legendre_polynomial, only: evaluate_legendre_polinomial_batch implicit none @@ -20,24 +21,32 @@ program test !a = r_gaunt(10, 10, 10, 0, 0, 0) - call get_grid_function(0.0001d0, 10.0d0, 200, r_grid, slater) - call my_spline%create(r_grid, slater) - !call get_grid_function(0.01d0, 9.0d0, 122, r_grid_out, slater_out) - !call my_spline%evaluate_batch(r_grid_out, y_out) - !do i = 1, 122 - ! !a = my_spline%evaluate(r_grid_out(i)) - ! print *, i, r_grid_out(i), slater_out(i), y_out(i), slater_out(i) - y_out(i) - !end do - - - call get_gauss_legendre_grid(n, 0.001d0, 10.0d0, gauleg_grid, gauleg_weight) - call my_spline%evaluate_batch(gauleg_grid, gauleg_func) - do i = 1, n - gauleg_func(i) = (gauleg_func(i)*gauleg_grid(i))**2 - print *, i, gauleg_grid(i), gauleg_weight(i), gauleg_func(i) - end do - a = gauss_legendre_integrator(n, gauleg_weight, gauleg_func) - print *, a + ! test splines + ! call get_grid_function(0.0001d0, 10.0d0, 200, r_grid, slater) + ! call my_spline%create(r_grid, slater) + ! call get_grid_function(0.01d0, 9.0d0, 122, r_grid_out, slater_out) + ! call my_spline%evaluate_batch(r_grid_out, y_out) + ! do i = 1, 122 + ! !a = my_spline%evaluate(r_grid_out(i)) + ! print *, i, r_grid_out(i), slater_out(i), y_out(i), slater_out(i) - y_out(i) + ! end do + + + ! test gauss integration + ! call get_grid_function(0.0001d0, 10.0d0, 200, r_grid, slater) + ! call my_spline%create(r_grid, slater) + ! call get_gauss_legendre_grid(n, 0.001d0, 10.0d0, gauleg_grid, gauleg_weight) + ! call my_spline%evaluate_batch(gauleg_grid, gauleg_func) + ! do i = 1, n + ! gauleg_func(i) = (gauleg_func(i)*gauleg_grid(i))**2 + ! print *, i, gauleg_grid(i), gauleg_weight(i), gauleg_func(i) + ! end do + ! a = gauss_legendre_integrator(n, gauleg_weight, gauleg_func) + ! print *, a + + + ! test the legendre polinomials + call test_legendre_polinomials() contains @@ -69,5 +78,35 @@ subroutine get_grid_function(r_start, r_end, n, grid, func) end subroutine + subroutine test_legendre_polinomials() + + integer, parameter :: n = 7 + integer, parameter :: n_points = 300 + real(kind=8) :: p(n_points, n+1), x(n_points) + + real(kind=8), parameter :: r_start = -1.0d0 + real(kind=8), parameter :: r_end = 1.0d0 + + real(kind=8) :: step + integer :: i, j + + step = (r_end - r_start)/(n_points-1) + do i = 1, n_points + x(i) = r_start + (i-1) * step + end do + + call evaluate_legendre_polinomial_batch(n, n_points, x, p) + + do i = 1, n_points + do j = 1, n+1 + write(*, "(E20.8)", advance="no") p(i, j) + write(*, "(A)", advance="no") " " + end do + print * + end do + + + end subroutine test_legendre_polinomials + end program test \ No newline at end of file