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