From 36c6f48b0fa118adeadca214c9c595224d30584a Mon Sep 17 00:00:00 2001 From: sbraccia Date: Wed, 17 Nov 2004 12:05:53 +0000 Subject: [PATCH] Cleanup of spherical bessel functions. Added a routine to compute the derivativs of spherical bessel functions (used by calc_btq). C.S. git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1450 c92efa57-630b-4861-b058-cf58834340f0 --- flib/Makefile | 1 + flib/sph_bes.f90 | 398 +++++++++++++++++++++++++--------------------- flib/sph_dbes.f90 | 109 +++++++++++++ 3 files changed, 327 insertions(+), 181 deletions(-) create mode 100644 flib/sph_dbes.f90 diff --git a/flib/Makefile b/flib/Makefile index cb87b1482..ebf07c9d0 100644 --- a/flib/Makefile +++ b/flib/Makefile @@ -22,6 +22,7 @@ scnds.o \ simpsn.o \ sort.o \ sph_bes.o \ +sph_dbes.o \ transto.o \ date_and_tim.o \ sort_gvec.o \ diff --git a/flib/sph_bes.f90 b/flib/sph_bes.f90 index e2e68200f..25791e44c 100644 --- a/flib/sph_bes.f90 +++ b/flib/sph_bes.f90 @@ -1,247 +1,283 @@ ! -! Copyright (C) 2001 PWSCF group +! Copyright (C) 2001-2004 PWSCF-FPMD-CP90 group ! This file is distributed under the terms of the ! GNU General Public License. See the file `License' ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! -! -!-------------------------------------------------------------------- -subroutine sph_bes (msh, r, q, l, jl) - !-------------------------------------------------------------------- +!---------------------------------------------------------------------------- +SUBROUTINE sph_bes( msh, r, q, l, jl ) + !---------------------------------------------------------------------------- ! - ! input: - ! msh = number of grid points points - ! r(1:msh)= radial grid - ! q = q - ! l = angular momentum (-1 <= l <= 6) - ! output: - ! jl(1:msh) = j_l(q*r(i)) (j_l = spherical bessel function) + ! ... input: + ! ... msh = number of grid points points + ! ... r(1:msh)= radial grid + ! ... q = q + ! ... l = angular momentum (-1 <= l <= 6) + ! ... output: + ! ... jl(1:msh) = j_l(q*r(i)) (j_l = spherical bessel function) ! - use kinds, only: DP - implicit none + USE kinds, ONLY : DP + USE constants, ONLY : eps8 ! - integer :: msh, l - real(kind=DP) :: r (msh), q, jl (msh) + IMPLICIT NONE ! - real(kind=DP), parameter :: eps = 1.d-8 - integer :: ir, ir0 + INTEGER :: msh, l + REAL(kind=DP) :: r(msh), q, jl(msh) ! -#if defined __MASS - real(kind=DP) :: qr(msh), sin_qr(msh), cos_qr(msh) + INTEGER :: ir, ir0 + ! +#if defined (__MASS) + REAL(kind=DP) :: qr(msh), sin_qr(msh), cos_qr(msh) #endif ! - if (abs (q) < eps) then - if (l == -1) then - call errore ('sph_bes', 'j_{-1}(0) ?!?', 1) - elseif (l == 0) then - jl(:) = 1.d0 - else - jl(:) = 0.d0 - endif - else - if (abs (q * r (1) ) > eps) then + IF ( ABS( q ) < eps8 ) THEN + ! + IF ( l == -1 ) THEN + ! + CALL errore( 'sph_bes', 'j_{-1}(0) ?!?', 1 ) + ! + ELSE IF ( l == 0 ) THEN + ! + jl(:) = 1.D0 + ! + ELSE + ! + jl(:) = 0.D0 + ! + END IF + ! + ELSE + ! + IF ( ABS( q*r(1) ) > eps8 ) THEN + ! ir0 = 1 - else - if (l == -1) then - call errore ('sph_bes', 'j_{-1}(0) ?!?', 2) - elseif (l == 0) then - jl (1) = 1.d0 - else - jl (1) = 0.d0 - endif + ! + ELSE + ! + IF (l == -1) THEN + ! + CALL errore( 'sph_bes', 'j_{-1}(0) ?!?', 2 ) + ! + ELSE IF ( l == 0 ) THEN + ! + jl(1) = 1.D0 + ! + ELSE + ! + jl(1) = 0.D0 + ! + ENDIF + ! ir0 = 2 - endif - - if (l == - 1) then - -#if defined __MASS - + ! + END IF + ! + IF ( l == - 1 ) THEN + ! +#if defined (__MASS) + ! qr = q * r - call vcos( cos_qr, qr, msh) - do ir = ir0, msh + ! + CALL vcos( cos_qr, qr, msh) + ! + DO ir = ir0, msh jl (ir) = cos_qr(ir) / (q * r (ir) ) - enddo - + END DO + ! #else - - do ir = ir0, msh - jl (ir) = cos (q * r (ir) ) / (q * r (ir) ) - enddo - + ! + DO ir = ir0, msh + jl (ir) = COS (q * r (ir) ) / (q * r (ir) ) + END DO + ! #endif - - elseif (l == 0) then - -#if defined __MASS - + ! + ELSE IF ( l == 0 ) THEN + ! +#if defined (__MASS) + ! qr = q * r - call vsin( sin_qr, qr, msh) - do ir = ir0, msh + ! + CALL vsin( sin_qr, qr, msh) + ! + DO ir = ir0, msh jl (ir) = sin_qr(ir) / (q * r (ir) ) - enddo - + END DO + ! #else - - do ir = ir0, msh - jl (ir) = sin (q * r (ir) ) / (q * r (ir) ) - enddo - + ! + DO ir = ir0, msh + jl (ir) = SIN (q * r (ir) ) / (q * r (ir) ) + END DO + ! #endif - - elseif (l == 1) then - -#if defined __MASS - + ! + ELSE IF ( l == 1 ) THEN + ! +#if defined (__MASS) + ! qr = q * r - call vcos( cos_qr, qr, msh) - call vsin( sin_qr, qr, msh) - do ir = ir0, msh + ! + CALL vcos( cos_qr, qr, msh) + CALL vsin( sin_qr, qr, msh) + ! + DO ir = ir0, msh jl (ir) = ( sin_qr(ir) / (q * r (ir) ) - cos_qr(ir) ) / (q * r (ir) ) - enddo - + END DO + ! #else - - do ir = ir0, msh - jl (ir) = (sin (q * r (ir) ) / (q * r (ir) ) - cos (q * r ( & + ! + DO ir = ir0, msh + jl (ir) = (SIN (q * r (ir) ) / (q * r (ir) ) - COS (q * r ( & ir) ) ) / (q * r (ir) ) - enddo - + END DO + ! #endif - - elseif (l == 2) then - -#if defined __MASS - + ! + ELSE IF ( l == 2 ) THEN + ! +#if defined (__MASS) + ! qr = q * r - call vcos( cos_qr, qr, msh) - call vsin( sin_qr, qr, msh) - do ir = ir0, msh + ! + CALL vcos( cos_qr, qr, msh) + CALL vsin( sin_qr, qr, msh) + ! + DO ir = ir0, msh jl (ir) = ( (3.d0 / (q * r (ir) ) - (q * r (ir) ) ) * sin_qr( & ir ) - 3.d0 * cos_qr(ir) ) / (q * r (ir) ) ** 2 - enddo - + END DO + ! #else - - do ir = ir0, msh - jl (ir) = ( (3.d0 / (q * r (ir) ) - (q * r (ir) ) ) * sin ( & - q * r (ir) ) - 3.d0 * cos (q * r (ir) ) ) / (q * r (ir) ) ** & - 2 - enddo - + ! + DO ir = ir0, msh + jl (ir) = ( (3.d0 / (q * r (ir) ) - (q * r (ir) ) ) * SIN ( & + q * r (ir) ) - 3.d0 * COS (q * r (ir) ) ) / (q * r (ir) ) ** 2 + END DO + ! #endif - - elseif (l == 3) then - -#if defined __MASS - + ! + ELSE IF ( l == 3 ) THEN + ! +#if defined (__MASS) + ! qr = q * r - call vcos( cos_qr, qr, msh) - call vsin( sin_qr, qr, msh) - do ir = ir0, msh + ! + CALL vcos( cos_qr, qr, msh) + CALL vsin( sin_qr, qr, msh) + ! + DO ir = ir0, msh jl (ir) = (sin_qr (ir) * (15.d0 / (q * r (ir) ) & - 6.d0 * (q * r (ir) ) ) + cos_qr (ir) * ( (q * r (ir) & ) **2 - 15.d0) ) / (q * r (ir) ) **3 - enddo - + ENDDO + ! #else - - do ir = ir0, msh - jl (ir) = (sin (q * r (ir) ) * (15.d0 / (q * r (ir) ) & - - 6.d0 * (q * r (ir) ) ) + cos (q * r (ir) ) * ( (q * r (ir) & - ) **2 - 15.d0) ) / (q * r (ir) ) **3 - enddo - + ! + DO ir = ir0, msh + jl (ir) = (SIN (q * r (ir) ) * (15.d0 / (q * r (ir) ) & + - 6.d0 * (q * r (ir) ) ) + COS (q * r (ir) ) * ( (q * r (ir) & + ) **2 - 15.d0) ) / (q * r (ir) ) ** 3 + END DO + ! #endif - - elseif (l == 4) then - -#if defined __MASS - + ! + ELSE IF ( l == 4 ) THEN + ! +#if defined (__MASS) + ! qr = q * r - call vcos( cos_qr, qr, msh) - call vsin( sin_qr, qr, msh) - do ir = ir0, msh + ! + CALL vcos( cos_qr, qr, msh) + CALL vsin( sin_qr, qr, msh) + ! + DO ir = ir0, msh jl (ir) = (sin_qr (ir) * (105.d0 - 45.d0 * (q * r (ir) & ) **2 + (q * r (ir) ) **4) + cos_qr (ir) * (10.d0 * & (q * r (ir) ) **3 - 105.d0 * (q * r (ir) ) ) ) / (q * r (ir) & ) **5 - enddo - + END DO + ! #else - - do ir = ir0, msh - jl (ir) = (sin (q * r (ir) ) * (105.d0 - 45.d0 * (q * r (ir) & - ) **2 + (q * r (ir) ) **4) + cos (q * r (ir) ) * (10.d0 * & + ! + DO ir = ir0, msh + jl (ir) = (SIN (q * r (ir) ) * (105.d0 - 45.d0 * (q * r (ir) & + ) **2 + (q * r (ir) ) **4) + COS (q * r (ir) ) * (10.d0 * & (q * r (ir) ) **3 - 105.d0 * (q * r (ir) ) ) ) / (q * r (ir) & ) **5 - enddo - + END DO + ! #endif - - elseif (l == 5) then - -#if defined __MASS - + ! + ELSE IF ( l == 5 ) THEN + ! +#if defined (__MASS) + ! qr = q * r - call vcos( cos_qr, qr, msh) - call vsin( sin_qr, qr, msh) - do ir = ir0, msh + ! + CALL vcos( cos_qr, qr, msh ) + CALL vsin( sin_qr, qr, msh ) + ! + DO ir = ir0, msh jl (ir) = (-cos_qr(ir) - (945.d0*cos_qr(ir))/(q*r(ir)) ** 4 + & (105.d0*cos_qr(ir))/ (q*r(ir)) ** 2 + & (945.d0*sin_qr(ir))/ (q*r(ir)) ** 5 - & (420.d0*sin_qr(ir))/(q*r(ir)) ** 3 + & (15.d0*sin_qr(ir))/ (q*r(ir)))/(q*r(ir)) - end do - + END DO + ! #else - - do ir = ir0, msh - jl (ir) = (-cos(q*r(ir)) - (945.d0*cos(q*r(ir)))/(q*r(ir)) ** 4 + & - (105.d0*cos(q*r(ir)))/ (q*r(ir)) ** 2 + & - (945.d0*sin(q*r(ir)))/ (q*r(ir)) ** 5 - & - (420.d0*sin(q*r(ir)))/(q*r(ir)) ** 3 + & - (15.d0*sin(q*r(ir)))/ (q*r(ir)))/(q*r(ir)) - end do - + ! + DO ir = ir0, msh + jl (ir) = (-COS(q*r(ir)) - (945.d0*COS(q*r(ir)))/(q*r(ir)) ** 4 + & + (105.d0*COS(q*r(ir)))/ (q*r(ir)) ** 2 + & + (945.d0*SIN(q*r(ir)))/ (q*r(ir)) ** 5 - & + (420.d0*SIN(q*r(ir)))/(q*r(ir)) ** 3 + & + (15.d0*SIN(q*r(ir)))/ (q*r(ir)))/(q*r(ir)) + END DO + ! #endif - - elseif (l == 6) then - -#if defined __MASS - + ! + ELSE IF ( l == 6 ) THEN + ! +#if defined (__MASS) + ! qr = q * r - call vcos( cos_qr, qr, msh) - call vsin( sin_qr, qr, msh) - do ir = ir0, msh + ! + CALL vcos( cos_qr, qr, msh ) + CALL vsin( sin_qr, qr, msh ) + ! + DO ir = ir0, msh jl (ir) = ((-10395.d0*cos_qr(ir))/(q*r(ir)) ** 5 + & (1260.d0*cos_qr(ir))/(q*r(ir)) ** 3 - & (21.d0*cos_qr(ir))/ (q*r(ir)) - sin_qr(ir) + & (10395.d0*sin_qr(ir))/(q*r(ir)) ** 6 - & (4725.d0*sin_qr(ir))/ (q*r(ir)) ** 4 + & (210.d0*sin_qr(ir))/(q*r(ir)) ** 2)/(q*r(ir)) - end do - + END DO + ! #else - - do ir = ir0, msh - jl (ir) = ((-10395.d0*cos(q*r(ir)))/(q*r(ir)) ** 5 + & - (1260.d0*cos(q*r(ir)))/(q*r(ir)) ** 3 - & - (21.d0*cos(q*r(ir)))/ (q*r(ir)) - sin(q*r(ir)) + & - (10395.d0*sin(q*r(ir)))/(q*r(ir)) ** 6 - & - (4725.d0*sin(q*r(ir)))/ (q*r(ir)) ** 4 + & - (210.d0*sin(q*r(ir)))/(q*r(ir)) ** 2)/(q*r(ir)) - end do - + ! + DO ir = ir0, msh + jl (ir) = ((-10395.d0*COS(q*r(ir)))/(q*r(ir)) ** 5 + & + (1260.d0*COS(q*r(ir)))/(q*r(ir)) ** 3 - & + (21.d0*COS(q*r(ir)))/ (q*r(ir)) - SIN(q*r(ir)) + & + (10395.d0*SIN(q*r(ir)))/(q*r(ir)) ** 6 - & + (4725.d0*SIN(q*r(ir)))/ (q*r(ir)) ** 4 + & + (210.d0*SIN(q*r(ir)))/(q*r(ir)) ** 2)/(q*r(ir)) + END DO + ! #endif - - else - - call errore ('sph_bes', 'not implemented', l) - - endif - - endif + ! + ELSE + ! + CALL errore ('sph_bes', 'not implemented', l) + ! + END IF + ! + END IF ! - return -end subroutine sph_bes + RETURN + ! +END SUBROUTINE sph_bes diff --git a/flib/sph_dbes.f90 b/flib/sph_dbes.f90 new file mode 100644 index 000000000..178e72c50 --- /dev/null +++ b/flib/sph_dbes.f90 @@ -0,0 +1,109 @@ +! +! Copyright (C) 2001-2004 PWSCF-FPMD-CP90 group +! This file is distributed under the terms of the +! GNU General Public License. See the file `License' +! in the root directory of the present distribution, +! or http://www.gnu.org/copyleft/gpl.txt . +! +!---------------------------------------------------------------------------- +SUBROUTINE sph_dbes( MMAX, R, XG, L, DJL ) + !---------------------------------------------------------------------------- + ! + ! ... calculates derivatives of spherical bessel functions j_l(Gr) + ! ... with respect to h_alpha,beta (without the factor GAGK(KK,IG)*HTM1) + ! ... i.e. -x * D(jl(x))/dx + ! + USE kinds, ONLY : DP + USE constants, ONLY : eps8 + ! + IMPLICIT NONE + ! + INTEGER :: MMAX, L + REAL(KIND=DP) :: XG + REAL(KIND=DP) :: DJL(MMAX), R(MMAX) + ! + INTEGER :: IR + REAL(KIND=DP) :: XRG, XRG2 + ! + ! + IF ( L == 1 ) THEN ! S PART + IF( XG < eps8 ) THEN + DO IR=1,MMAX + DJL(IR) = 0.D0 + END DO + ELSE + DJL(1) = 0.D0 + DO IR=2,MMAX + XRG=R(IR)*XG + DJL(IR) = SIN(XRG)/XRG-COS(XRG) + END DO + ENDIF + ENDIF + ! + IF ( L == 2 ) THEN ! P PART + IF( XG < eps8 ) THEN + DO IR=1,MMAX + DJL(IR) = 0.D0 + END DO + ELSE + DJL(1) = 0.D0 + DO IR=2,MMAX + XRG=R(IR)*XG + DJL(IR) = 2.D0*(SIN(XRG)/XRG-COS(XRG))/XRG - SIN(XRG) + END DO + ENDIF + ENDIF + ! + IF ( L == 3 ) THEN ! D PART + IF ( XG < eps8 ) THEN + DO IR=1,MMAX + DJL(IR) = 0.D0 + END DO + ELSE + DJL(1) = 0.D0 + DO IR=2,MMAX + XRG=R(IR)*XG + DJL(IR) = ( SIN(XRG)*(9.D0/(XRG*XRG)-4.D0) - & + 9.D0*COS(XRG)/XRG ) /XRG + COS(XRG) + END DO + END IF + END IF + ! + IF ( L == 4 ) THEN ! F PART + IF ( XG < eps8 ) THEN + DO IR=1,MMAX + DJL(IR) = 0.D0 + END DO + ELSE + DJL(1) = 0.D0 + DO IR=2,MMAX + XRG=R(IR)*XG + XRG2=XRG*XRG + DJL(IR) = SIN(XRG)*(60.D0/(XRG2*XRG2)-27.D0/XRG2+1.d0) - & + COS(XRG)*(60.D0/XRG2-7.D0)/XRG + END DO + END IF + END IF + ! + IF ( L == 5 ) THEN ! G PART + IF ( XG < eps8 ) THEN + DO IR=1,MMAX + DJL(IR) = 0.D0 + END DO + ELSE + DJL(1) = 0.D0 + DO IR=2,MMAX + XRG=R(IR)*XG + XRG2=XRG*XRG + DJL(IR) = SIN(XRG)*(525.D0/(XRG2*XRG2)-240.D0/XRG2+11.D0)/XRG - & + COS(XRG)*(525.D0/(XRG2*XRG2)-65.D0/XRG2+1.D0) + END DO + END IF + END IF + ! + IF ( L <= 0 .OR. L >= 6 ) & + CALL errore( 'sph_dbes', ' L NOT PROGRAMMED, L= ',L ) + ! + RETURN + ! +END SUBROUTINE sph_dbes