1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8!!@LICENSE 9! 10C ********************************************************************* 11C MODULE m_radfft 12C 13C Public procedures provided: 14C subroutine radfft ! Radial fast Fourier transform 15C 16C Public parameters, variables, and arrays: 17C none 18C 19C Used module procedures: 20C use m_bessph, only: bessph ! Spherical Bessel functions 21C use m_fft_gpfa,only: fft_gpfa_ez ! Fourier transform 22C use alloc, only: de_alloc ! Deallocation routines 23C use alloc, only: re_alloc ! (Re)allocation routines 24C 25C Used module parameters: 26C use precision, only: dp ! Double precision real kind 27C 28C ********************************************************************* 29C SUBROUTINE RADFFT( L, NR, RMAX, F, G ) 30C ********************************************************************* 31C Makes a fast Fourier transform of a radial function. 32C If function f is of the form 33C f(r_vec) = F(r_mod) * Ylm(theta,phi) 34C where Ylm is a spherical harmonic with l = argument L, and 35C argument F contains on input the real function F(r_mod), in a uniform 36C radial grid: 37C r_mod = ir * RMAX / NR, ir = 0,1,...,NR, 38C and if g is the 3-dimensional Fourier transform of f: 39C g(k_vec) = 1/(2*pi)**(3/2) * 40C Integral( d3_r_vec * exp(-i * k_vec * r_vec) * f(r_vec) ) 41C then g has the form 42C g(k_vec) = (-i)**L * G(k_mod) * Ylm(theta,phi) 43C where argument G contains on output the real function G(k_mod) in 44C a uniform radial grid: 45C k_mod = ik * k_max / NR, ik = 0,1,...,NR, k_max = NR*pi/RMAX 46C Ref: J.M.Soler notes of 16/08/95. 47C *************** INPUT *********************************************** 48C INTEGER L : Angular momentum quantum number 49C INTEGER NR : Number of radial intervals. 50C 2*NR must be an acceptable number of points for the 51C FFT routine used. 52C REAL*8 RMAX : Maximum radius 53C REAL*8 F(0:NR) : Function to be tranformed, in a radial mesh 54C *************** OUTPUT ********************************************** 55C REAL*8 G(0:NR) : Fourier transform of F (but see point 5 below) 56C *************** UNITS *********************************************** 57C Units of RMAX and F are arbitrary. 58C Units of k_max and G are related with those of RMAX and F in the 59C obvious way (see above). 60C *************** BEHAVIOUR ******************************************* 61C 1) F and G may be the same physical array, i.e. it is allowed: 62C CALL RADFFT( L, NR, RMAX, F, F ) 63C 2) It also works in the opposite direction, but then the factor 64C multiplying the output is (+i)**L. Thus, the following two calls 65C CALL RADFFT( L, NR, RMAX, F, G ) 66C CALL RADFFT( L, NR, NR*PI/RMAX, G, H ) 67C make H = F 68C 3) If you will divide the output by q**l, truncation errors may be 69C quite large for small k's if L and NR are large. Therefore, these 70C components are calculated by direct integration rather than FFT. 71C Parameter ERRFFT is the typical truncation error in the FFT, and 72C controls which k's are integrated directly. A good value is 1e-8. 73C If you will not divide by k**l, make ERRFFT=1.e-30. 74C 4) The function F is assumed to be zero at and beyond RMAX. The last 75C point F(NR) is not used to find G, except G(0) for L=0 (see 5) 76C 5) Because of the 'uncertainty principle', if f(r) is strictly zero 77C for r>RMAX, then g(k) cannot be strictly zero for k>kmax. 78C Therefore G(NR), which should be exactly zero, is used (for L=0) 79C as a 'reminder' term for the integral of G beyond kmax, to ensure 80C that F(0)=Sum[4*pi*r**2*dr*G(IR)] (this allows to recover F(0) 81C when called again in the inverse direction). Thus, the last value 82C G(NR) should be replaced by zero for any other use. NOTICE: this 83C is commented out in this version! 84C ********************************************************************* 85C Written by J.M.Soler. August 1996. 86! Work arrays handling by Rogeli Grima, ca 2009 87C ********************************************************************* 88 89 MODULE m_radfft 90 91 USE precision, only: dp ! Double precision real kind 92 USE m_bessph, only: bessph ! Spherical Bessel functions 93 use m_fft_gpfa,only: fft_gpfa_ez ! 1D fast Fourier transform 94 USE alloc, only: re_alloc, de_alloc 95! USE m_timer, only: timer_start ! Start counting CPU time 96! USE m_timer, only: timer_stop ! Stop counting CPU time 97 98 implicit none 99 100 PUBLIC :: radfft ! Radial fast Fourier transform 101 PUBLIC :: reset_radfft ! Deallocates work arrays 102 103 PRIVATE 104 105 ! Work arrays held in module to minimize reallocations 106 ! Note that we avoid "automatic" arrays, which may cause stack problems 107 real(dp), pointer :: GG(:) 108 real(dp), pointer :: FN(:,:) 109 real(dp), pointer :: P(:,:,:) 110 integer :: MAXL = -1 111 integer :: MAXNR = -1 112 113 CONTAINS 114 115 SUBROUTINE RADFFT( L, NR, RMAX, F, G ) 116 117 IMPLICIT NONE 118 119C Declare argument types and dimensions ----------------------------- 120 INTEGER, intent(in) :: L ! Angular momentum of function 121 INTEGER, intent(in) :: NR ! Number of radial points 122 real(dp),intent(in) :: RMAX ! Radius of last point 123 real(dp),intent(in) :: F(0:NR) ! Function to Fourier-transform 124 real(dp),intent(out):: G(0:NR) ! Fourier transform of F(r) 125C ------------------------------------------------------------------- 126 127C ERRFFT is the typical truncation error in the FFT routine --------- 128 real(dp), PARAMETER :: ERRFFT = 1.0E-8_dp 129C ------------------------------------------------------------------- 130 131C Internal variable types and dimensions ---------------------------- 132 INTEGER :: I, IQ, IR, JR, M, MQ, N, NQ 133 real(dp) :: C, DQ, DR, FR, PI, R, RN, Q, QMAX 134!! real(dp) :: GG(0:2*NR), FN(2,0:2*NR), P(2,0:L,0:L) 135C ------------------------------------------------------------------- 136 137C Start time counter ------------------------------------------------ 138* CALL TIMER_START( 'RADFFT' ) 139C 140C Allocate local memory 141 if (MAXL.eq.-1) nullify(P) 142 if (L.GT.MAXL) then 143 call re_alloc( P, 1, 2, 0, L, 0, L, 'P', 'RADFFT' ) 144 MAXL=L 145 endif 146 if (MAXNR.eq.-1) nullify(FN,GG) 147 if (NR.GT.MAXNR) then 148 call re_alloc( FN, 1, 2, 0, 2*NR, 'FN', 'RADFFT' ) 149 call re_alloc( GG, 0, 2*NR, 'GG', 'RADFFT' ) 150 MAXNR=NR 151 endif 152 153C Find some constants ----------------------------------------------- 154 PI = 4.D0 * ATAN( 1.D0 ) 155 NQ = NR 156 DR = RMAX / NR 157 DQ = PI / RMAX 158 QMAX = NQ * DQ 159 C = DR / SQRT( 2.D0*PI ) 160C ------------------------------------------------------------------- 161 162C Set up a complex polynomial such that the spherical Bessel function: 163C j_l(x) = Real( Sum_n( P(n,l) * x**n ) * exp(i*x) ) / x**(l+1) 164 P(1,0,0) = 0.D0 165 P(2,0,0) = -1.D0 166 if (l.gt.0) then 167 P(1,0,1) = 0.D0 168 P(2,0,1) = -1.D0 169 P(1,1,1) = -1.D0 170 P(2,1,1) = 0.D0 171 if (l.gt.1) then 172 DO M = 2,L 173 DO N = 0,M 174 DO I = 1,2 175 P(I,N,M) = 0.D0 176 IF (N .LT. M) P(I,N,M) = P(I,N,M) + (2*M-1) * P(I,N,M-1) 177 IF (N .GE. 2) P(I,N,M) = P(I,N,M) - P(I,N-2,M-2) 178 ENDDO 179 ENDDO 180 ENDDO 181 endif 182 endif 183C ------------------------------------------------------------------- 184 185C Initialize accumulation array ------------------------------------- 186 DO IQ = 0,NQ 187 GG(IQ) = 0.D0 188 ENDDO 189C ------------------------------------------------------------------- 190 191C Iterate on terms of the j_l(q*r) polynomial ----------------------- 192 DO N = 0,L 193 194C Set up function to be fast fourier transformed 195 FN(1,0) = 0.D0 196 FN(2,0) = 0.D0 197 DO JR = 1, 2*NR-1 198 199 IF (JR .LT. NR) THEN 200 IR = JR 201 R = IR * DR 202 FR = F(IR) 203 ELSEIF (JR .EQ. NR) THEN 204 IR = JR 205 R = IR * DR 206 FR = 0.D0 207 ELSE 208 IR = 2*NR - JR 209 R = - (IR * DR) 210 FR = F(IR) * (-1.D0)**L 211 ENDIF 212 213C Find r**2 * r**n / r**(l+1) 214 RN = R**(N-L+1) 215 216 FN(1,JR) = C * FR * RN * P(1,N,L) 217 FN(2,JR) = C * FR * RN * P(2,N,L) 218 ENDDO 219 220C Perform one-dimensional complex FFT 221! 222! Only the elements from 0 to 2*NR-1 of FN are used. 223! (a total of 2*NR). The fft routine will receive a one-dimensional 224! array of size 2*NR. 225! 226 CALL fft_gpfa_ez( FN, 2*NR, +1 ) 227 228C Accumulate contribution 229 DO IQ = 1,NQ 230 Q = IQ * DQ 231 GG(IQ) = ( GG(IQ) + FN(1,IQ) ) / Q 232 ENDDO 233 234 ENDDO 235C ------------------------------------------------------------------- 236 237C Special case for Q=0 --------------------------------------------- 238 GG(0) = 0.D0 239 IF ( L .EQ. 0 ) THEN 240 DO IR = 1,NR 241 R = IR * DR 242 GG(0) = GG(0) + R*R * F(IR) 243 ENDDO 244 GG(0) = GG(0) * 2.D0 * C 245 ENDIF 246C ------------------------------------------------------------------- 247 248C Direct integration for the smallest Q's --------------------------- 249 IF (L.EQ.0) THEN 250 MQ = 0 251 ELSE 252 MQ = NQ * ERRFFT**(1.D0/L) 253 ENDIF 254 DO IQ = 1,MQ 255 Q = IQ * DQ 256 GG(IQ) = 0.D0 257 DO IR = 1,NR 258 R = IR * DR 259 GG(IQ) = GG(IQ) + R*R * F(IR) * BESSPH(L,Q*R) 260 ENDDO 261 GG(IQ) = GG(IQ) * 2.D0 * C 262 ENDDO 263C ------------------------------------------------------------------- 264 265C Special case for Q=QMAX ------------------------------------------- 266* IF (L.EQ.0) THEN 267* GSUM = 0.D0 268* DO IQ = 1,NQ-1 269* Q = IQ * DQ 270* GSUM = GSUM + Q*Q * GG(IQ) 271* ENDDO 272* GSUM = GSUM * 4.D0 * PI * DQ 273* GG(NQ) = (2.D0*PI)**1.5D0 * F(0) - GSUM 274* GG(NQ) = GG(NQ) / (4.D0 * PI * DQ * QMAX**2) 275* ENDIF 276C ------------------------------------------------------------------- 277 278C Copy from local to output array ----------------------------------- 279 DO IQ = 0,NQ 280 G(IQ) = GG(IQ) 281 ENDDO 282C ------------------------------------------------------------------- 283 284C Stop time counter ------------------------------------------------ 285* CALL TIMER_STOP( 'RADFFT' ) 286C ------------------------------------------------------------------- 287 288 END SUBROUTINE radfft 289 290 SUBROUTINE RESET_RADFFT( ) 291 implicit none 292 call de_alloc( P, 'P', 'RADFFT' ) 293 call de_alloc( FN, 'FN', 'RADFFT' ) 294 call de_alloc( GG, 'GG', 'RADFFT' ) 295 MAXL = -1 296 MAXNR = -1 297 END SUBROUTINE RESET_RADFFT 298 299 END MODULE m_radfft 300 301