add legendre polinomial module

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)
 

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 

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