diff --git a/plugins/local/spher_harm/README.rst b/plugins/local/spher_harm/README.rst index bf897f73..9c9b12a6 100644 --- a/plugins/local/spher_harm/README.rst +++ b/plugins/local/spher_harm/README.rst @@ -2,3 +2,6 @@ spher_harm ========== +Routines for spherical Harmonics evaluation in real space. +The main routine is "spher_harm_func_r3(r,l,m,re_ylm, im_ylm)". +The test routine is "test_spher_harm" where everything is explained in details. diff --git a/plugins/local/spher_harm/routines_test.irp.f b/plugins/local/spher_harm/routines_test.irp.f index 6f7cbc1c..fe8fc422 100644 --- a/plugins/local/spher_harm/routines_test.irp.f +++ b/plugins/local/spher_harm/routines_test.irp.f @@ -1,10 +1,93 @@ +subroutine test_spher_harm + implicit none + BEGIN_DOC + ! routine to test the generic spherical harmonics routine "spher_harm_func_r3" from R^3 --> C + ! + ! We test = delta_m1,m2 delta_l1,l2 + ! + ! The test is done through the integration on a sphere with the Lebedev grid. + END_DOC + include 'constants.include.F' + integer :: l1,m1,i,l2,m2,lmax + double precision :: r(3),weight,accu_re, accu_im,accu + double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 + double precision :: theta,phi,r_abs + lmax = 5 ! Maximum angular momentum until which we are going to test orthogonality conditions + do l1 = 0,lmax + do m1 = -l1 ,l1 + do l2 = 0,lmax + do m2 = -l2 ,l2 + accu_re = 0.d0 ! accumulator for the REAL part of + accu_im = 0.d0 ! accumulator for the IMAGINARY part of + accu = 0.d0 ! accumulator for the weights ==> should be \int dOmega == 4 pi + ! = \int dOmega Y_l1,m1^* Y_l2,m2 + ! \approx \sum_i W_i Y_l1,m1^*(r_i) Y_l2,m2(r_i) WITH r_i being on the spher of radius 1 + do i = 1, n_points_integration_angular + r(1:3) = angular_quadrature_points(i,1:3) ! ith Lebedev point (x,y,z) on the sphere of radius 1 + weight = weights_angular_points(i) ! associated Lebdev weight not necessarily positive + +!!!!!!!!!!! Test of the Cartesian --> Spherical coordinates + ! theta MUST belong to [0,pi] and phi to [0,2pi] + ! gets the cartesian to spherical change of coordinates + call cartesian_to_spherical(r,theta,phi,r_abs) + if(theta.gt.pi.or.theta.lt.0.d0)then + print*,'pb with theta, it should be in [0,pi]',theta + print*,r + endif + if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then + print*,'pb with phi, it should be in [0,2 pi]',phi/pi + print*,r + endif + +!!!!!!!!!!! Routines returning the Spherical harmonics on the grid point + call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) + call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) + +!!!!!!!!!!! Integration of Y_l1,m1^*(r) Y_l2,m2(r) + ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) + ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) + accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu += weight + enddo + ! Test that the sum of the weights is 4 pi + if(dabs(accu - dfour_pi).gt.1.d-6)then + print*,'Problem !! The sum of the Lebedev weight is not 4 pi ..' + print*,accu + stop + endif + ! Test for the delta l1,l2 and delta m1,m2 + ! + ! Test for the off-diagonal part of the Kronecker delta + if(l1.ne.l2.or.m1.ne.m2)then + if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then + print*,'pb OFF DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + ! Test for the diagonal part of the Kronecker delta + if(l1==l2.and.m1==m2)then + if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then + print*,'pb DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + enddo + enddo + enddo + enddo +end subroutine test_cart implicit none BEGIN_DOC ! test for the cartesian --> spherical change of coordinates ! - ! simple test such that the polar angle theta ranges in [0,pi] + ! test the routine "cartesian_to_spherical" such that the polar angle theta ranges in [0,pi] ! ! and the asymuthal angle phi ranges in [0,2pi] END_DOC @@ -40,97 +123,18 @@ subroutine test_cart print*,phi/pi end -subroutine test_spher_harm - implicit none - BEGIN_DOC - ! routine to test the spherical harmonics integration on a sphere with the grid. - ! - ! We test = delta_m1,m2 delta_l1,l2 - END_DOC - include 'constants.include.F' - integer :: l1,m1,i,l2,m2,lmax - double precision :: r(3),weight,accu_re, accu_im,accu - double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 - l1 = 0 - m1 = 0 - l2 = 0 - m2 = 0 - lmax = 5 - do l1 = 0,lmax - do m1 = -l1 ,l1 - do l2 = 0,lmax - do m2 = -l2 ,l2 - accu_re = 0.d0 - accu_im = 0.d0 - ! = \int dOmega Y_l1,m1^* Y_l2,m2 - ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) - ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu = 0.d0 - do i = 1, n_points_integration_angular - double precision :: theta,phi,r_abs - r(1:3) = angular_quadrature_points(i,1:3) - weight = weights_angular_points(i) - call cartesian_to_spherical(r,theta,phi,r_abs) - if(theta.gt.pi.or.theta.lt.0.d0)then - print*,'pb with theta',theta - print*,r - endif - if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then - print*,'pb with phi',phi/pi - print*,r - endif - call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) - call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) - accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) - accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu += weight - write(33,'(10(F16.10,X))')phi/pi - enddo - ! Test for the delta l1,l2 and delta m1,m2 - if(l1.ne.l2.or.m1.ne.m2)then - if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then - print*,'pb OFF DIAG !!!!! ' - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - endif - endif - if(l1==l2.and.m1==m2)then - if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then - print*,'pb DIAG !!!!! ' - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - endif - endif - enddo - enddo - enddo - enddo - double precision :: x,dx,xmax,xmin - integer:: nx - nx = 10000 - xmin = -5.d0 - xmax = 5.d0 - dx = (xmax - xmin)/dble(nx) - x = xmin - do i = 1, nx - write(34,*)x,datan(x),dacos(x) - x += dx - enddo -end subroutine test_brutal_spheric implicit none include 'constants.include.F' BEGIN_DOC - ! test for the = delta_m1,m2 delta_l1,l2 using a two dimentional integration + ! Test for the = delta_m1,m2 delta_l1,l2 using the following two dimentional integration ! ! \int_0^2pi d Phi \int_-1^+1 d(cos(Theta)) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi) ! != \int_0^2pi d Phi \int_0^pi dTheta sin(Theta) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi) ! - ! Allows to test for the general functions spher_harm_func_m_pos with spher_harm_func_expl + ! Allows to test for the general functions "spher_harm_func_m_pos" with "spher_harm_func_expl" END_DOC integer :: itheta, iphi,ntheta,nphi double precision :: theta_min, theta_max, dtheta,theta @@ -147,7 +151,7 @@ subroutine test_brutal_spheric dphi = (phi_max - phi_min)/dble(nphi) dtheta = (theta_max - theta_min)/dble(ntheta) - lmax = 3 + lmax = 2 do l1 = 0,lmax do m1 = 0 ,l1 do l2 = 0,lmax @@ -196,7 +200,7 @@ end subroutine test_assoc_leg_pol implicit none BEGIN_DOC -! TODO : Put the documentation of the program here +! Test for the associated Legendre Polynoms. The test is done through the orthogonality condition. END_DOC print *, 'Hello world' integer :: l1,m1,ngrid,i,l2,m2 diff --git a/plugins/local/spher_harm/spher_harm.irp.f b/plugins/local/spher_harm/spher_harm.irp.f index e8deafb9..7a2eea06 100644 --- a/plugins/local/spher_harm/spher_harm.irp.f +++ b/plugins/local/spher_harm/spher_harm.irp.f @@ -1,7 +1,7 @@ program spher_harm implicit none - call test_spher_harm +! call test_spher_harm ! call test_cart -! call test_brutal_spheric + call test_brutal_spheric end