Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 : !!****m* ABINIT/m_FFTHyb
6 : !! NAME
7 : !! m_FFTHyb
8 : !!
9 : !! FUNCTION
10 : !! Almost useless. Just uses for FFT time evolution
11 : !! of number of electrons
12 : !!
13 : !! COPYRIGHT
14 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
15 : !! This file is distributed under the terms of the
16 : !! GNU General Public License, see ~abinit/COPYING
17 : !! or http://www.gnu.org/copyleft/gpl.txt .
18 : !!
19 : !! NOTES
20 : !!
21 : !! SOURCE
22 :
23 : #define FFTHyb_FORWARD 1
24 : #define FFTHyb_BACKWARD -1
25 : #include "defs.h"
26 : MODULE m_FFTHyb
27 : USE m_global
28 :
29 : IMPLICIT NONE
30 :
31 : !!***
32 :
33 : PRIVATE
34 :
35 : !!****t* m_FFTHyb/FFTHyb
36 : !! NAME
37 : !! FFTHyb
38 : !!
39 : !! FUNCTION
40 : !! This structured datatype contains the necessary data
41 : !!
42 : !! COPYRIGHT
43 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
44 : !! This file is distributed under the terms of the
45 : !! GNU General Public License, see ~abinit/COPYING
46 : !! or http://www.gnu.org/copyleft/gpl.txt .
47 : !!
48 : !! SOURCE
49 :
50 : TYPE, PUBLIC :: FFTHyb
51 : LOGICAL _PRIVATE :: set = .FALSE.
52 : INTEGER :: size
53 : DOUBLE PRECISION _PRIVATE :: Ts
54 : DOUBLE PRECISION _PRIVATE :: fs
55 : INTEGER , ALLOCATABLE, DIMENSION(:) _PRIVATE :: bit_rev
56 : COMPLEX(KIND=8) , ALLOCATABLE, DIMENSION(:) _PRIVATE :: data_inout
57 : END TYPE FFTHyb
58 : !!***
59 :
60 : PUBLIC :: FFTHyb_init
61 : PRIVATE :: FFTHyb_mirror
62 : PUBLIC :: FFTHyb_setData
63 : PUBLIC :: FFTHyb_run
64 : PUBLIC :: FFTHyb_getData
65 : PUBLIC :: FFTHyb_destroy
66 :
67 :
68 : CONTAINS
69 : !!***
70 :
71 : !!****f* ABINIT/m_FFTHyb/FFTHyb_init
72 : !! NAME
73 : !! FFTHyb_init
74 : !!
75 : !! FUNCTION
76 : !! Initialize ...
77 : !!
78 : !! COPYRIGHT
79 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
80 : !! This file is distributed under the terms of the
81 : !! GNU General Public License, see ~abinit/COPYING
82 : !! or http://www.gnu.org/copyleft/gpl.txt .
83 : !!
84 : !! INPUTS
85 : !! this=FFT
86 : !! n=Number of point (should be power of 2)
87 : !! samples_sec=number of samples per sec
88 : !!
89 : !! OUTPUT
90 : !!
91 : !! SIDE EFFECTS
92 : !!
93 : !! NOTES
94 : !!
95 : !! SOURCE
96 :
97 0 : SUBROUTINE FFTHyb_init(this,n,samples_sec)
98 :
99 : !Arguments ------------------------------------
100 : TYPE(FFTHyb), INTENT(INOUT) :: this
101 : INTEGER , INTENT(IN ) :: n
102 : DOUBLE PRECISION, INTENT(IN ) :: samples_sec
103 : !Local variables ------------------------------
104 : INTEGER :: i
105 : INTEGER :: inv_bit
106 : INTEGER :: total_size
107 :
108 : ! Check n is a power of 2
109 0 : IF ( n .LT. 2 .OR. IAND(n,n-1) .NE. 0 ) THEN
110 0 : CALL WARNALL("FFTHyb_init : array size is not a power of 2 -> auto fix")
111 0 : i = 1
112 0 : DO WHILE ( i .LT. n )
113 0 : i = ISHFT( i, 1 )
114 : END DO
115 0 : total_size = ISHFT(i, -1)
116 : ELSE
117 : total_size = n
118 : END IF
119 :
120 0 : this%size = total_size
121 0 : this%Ts = DBLE(total_size) / samples_sec
122 0 : this%fs = 1.d0 / DBLE(total_size)
123 :
124 0 : FREEIF(this%data_inout)
125 0 : MALLOC(this%data_inout,(0:total_size-1))
126 0 : this%data_inout(0:total_size-1) = CMPLX(0.d0,0.d0,8)
127 0 : FREEIF(this%bit_rev)
128 0 : MALLOC(this%bit_rev,(0:total_size-1))
129 0 : this%bit_rev = 0
130 :
131 0 : DO i = 1, total_size-1
132 0 : inv_bit = FFTHyb_mirror(i,total_size)
133 0 : this%bit_rev(inv_bit) = i
134 : END DO
135 :
136 0 : this%set = .TRUE.
137 :
138 0 : END SUBROUTINE FFTHyb_init
139 : !!***
140 :
141 : !!****f* ABINIT/m_FFTHyb/FFTHyb_mirror
142 : !! NAME
143 : !! FFTHyb_mirror
144 : !!
145 : !! FUNCTION
146 : !! mirror bits of an integer
147 : !!
148 : !! COPYRIGHT
149 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
150 : !! This file is distributed under the terms of the
151 : !! GNU General Public License, see ~abinit/COPYING
152 : !! or http://www.gnu.org/copyleft/gpl.txt .
153 : !!
154 : !! INPUTS
155 : !! i=bits
156 : !! n=integer
157 : !!
158 : !! OUTPUT
159 : !!
160 : !! SIDE EFFECTS
161 : !!
162 : !! NOTES
163 : !!
164 : !! SOURCE
165 :
166 : INTEGER FUNCTION FFTHyb_mirror(i,n)
167 :
168 : !Arguments ------------------------------------
169 : INTEGER, INTENT(IN ) :: i
170 : INTEGER, INTENT(IN ) :: n
171 : !Local variable -------------------------------
172 : INTEGER :: icp
173 : INTEGER :: ncp
174 :
175 : icp = i
176 : ncp = n
177 :
178 : FFTHyb_mirror = 0
179 :
180 0 : DO WHILE( ncp .GT. 1 )
181 0 : FFTHyb_mirror = IOR(ISHFT(FFTHyb_mirror,1) , IAND(icp,1))
182 0 : icp = ISHFT(icp,-1)
183 0 : ncp = ISHFT(ncp,-1)
184 : END DO
185 :
186 : END FUNCTION FFTHyb_mirror
187 : !!***
188 :
189 : !!****f* ABINIT/m_FFTHyb/FFTHyb_setData
190 : !! NAME
191 : !! FFTHyb_setData
192 : !!
193 : !! FUNCTION
194 : !! set input data (in time)
195 : !!
196 : !! COPYRIGHT
197 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
198 : !! This file is distributed under the terms of the
199 : !! GNU General Public License, see ~abinit/COPYING
200 : !! or http://www.gnu.org/copyleft/gpl.txt .
201 : !!
202 : !! INPUTS
203 : !! this=FFThyb
204 : !! array_in=data in
205 : !!
206 : !! OUTPUT
207 : !!
208 : !! SIDE EFFECTS
209 : !!
210 : !! NOTES
211 : !!
212 : !! SOURCE
213 :
214 0 : SUBROUTINE FFTHyb_setData(this, array_in)
215 : !Arguments ------------------------------------
216 : TYPE(FFTHyb), INTENT(INOUT) :: this
217 : DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: array_in
218 : INTEGER :: size_in
219 : INTEGER :: i
220 :
221 0 : size_in = SIZE(array_in)
222 0 : this%data_inout = CMPLX(0.d0,0.d0,KIND=4)
223 : ! IF ( size_in .NE. this%size ) &
224 : ! CALL WARNALL("FFTHyb_setData : size_in != size")
225 :
226 0 : DO i = 0, MIN(this%size,size_in)-1
227 0 : this%data_inout(i) = CMPLX(array_in(i+1), 0.d0,8)
228 : END DO
229 :
230 0 : END SUBROUTINE FFTHyb_setData
231 : !!***
232 :
233 : !!****f* ABINIT/m_FFTHyb/FFTHyb_run
234 : !! NAME
235 : !! FFTHyb_run
236 : !!
237 : !! FUNCTION
238 : !! perform FFT
239 : !!
240 : !! COPYRIGHT
241 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
242 : !! This file is distributed under the terms of the
243 : !! GNU General Public License, see ~abinit/COPYING
244 : !! or http://www.gnu.org/copyleft/gpl.txt .
245 : !!
246 : !! INPUTS
247 : !! this=FFT
248 : !! dir=direction FFTHyb_FORWARD or FFTHyb_BACKWARD
249 : !!
250 : !! OUTPUT
251 : !!
252 : !! SIDE EFFECTS
253 : !!
254 : !! NOTES
255 : !!
256 : !! SOURCE
257 :
258 0 : SUBROUTINE FFTHyb_run(this, dir)
259 :
260 : !Arguments ------------------------------------
261 : TYPE(FFTHyb), INTENT(INOUT) :: this
262 : INTEGER , INTENT(IN ) :: dir
263 : !Local variables ------------------------------
264 : INTEGER :: imax
265 : INTEGER :: istep
266 : INTEGER :: i
267 : INTEGER :: m
268 : INTEGER :: j
269 : DOUBLE PRECISION :: wtmp
270 : DOUBLE PRECISION :: wr
271 : DOUBLE PRECISION :: wpr
272 : DOUBLE PRECISION :: wpi
273 : DOUBLE PRECISION :: wi
274 : DOUBLE PRECISION :: theta
275 : DOUBLE PRECISION :: twoPi
276 : COMPLEX(KIND=8) :: tc
277 :
278 0 : imax = 1;
279 0 : istep = 2;
280 :
281 0 : twoPi = DBLE(dir)*2.d0*ACOS(-1.d0)
282 :
283 0 : DO WHILE ( imax .LT. this%size )
284 0 : istep = ISHFT(imax,1)
285 0 : theta = twoPi/DBLE(istep)
286 0 : Wtmp = SIN(0.5d0*theta)
287 0 : wpr = -2.d0*wtmp*wtmp
288 0 : wpi = SIN(theta)
289 0 : wr = 1.0d0
290 0 : wi = 0.d0
291 0 : DO m = 0, imax-1
292 0 : DO i = m, this%size-1, istep
293 0 : j= i+imax
294 : tc = CMPLX( wr* REAL (this%data_inout(this%bit_rev(j))) &
295 : -wi* AIMAG(this%data_inout(this%bit_rev(j))), &
296 : wr* AIMAG(this%data_inout(this%bit_rev(j))) &
297 0 : +wi* REAL (this%data_inout(this%bit_rev(j))), 8 )
298 0 : this%data_inout(this%bit_rev(j)) = this%data_inout(this%bit_rev(i)) - tc
299 0 : this%data_inout(this%bit_rev(i)) = this%data_inout(this%bit_rev(i)) + tc
300 : END DO
301 0 : wtmp = wr
302 0 : wr = wr*wpr - wi*wpi +wr
303 0 : wi = wi*wpr + wtmp*wpi+wi
304 : END DO
305 0 : imax = istep
306 : END DO
307 0 : IF ( dir .EQ. FFTHyb_FORWARD ) &
308 0 : this%data_inout = this%data_inout*this%fs
309 :
310 0 : END SUBROUTINE FFTHyb_run
311 : !!***
312 :
313 : !!****f* ABINIT/m_FFTHyb/FFTHyb_getData
314 : !! NAME
315 : !! FFTHyb_getData
316 : !!
317 : !! FUNCTION
318 : !! get result
319 : !!
320 : !! COPYRIGHT
321 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
322 : !! This file is distributed under the terms of the
323 : !! GNU General Public License, see ~abinit/COPYING
324 : !! or http://www.gnu.org/copyleft/gpl.txt .
325 : !!
326 : !! INPUTS
327 : !! this=FFT
328 : !!
329 : !! OUTPUT
330 : !! bound=number of frequencies
331 : !! array_out=output data
332 : !! freqs=output frequencies
333 : !!
334 : !! SIDE EFFECTS
335 : !!
336 : !! NOTES
337 : !!
338 : !! SOURCE
339 :
340 0 : SUBROUTINE FFTHyb_getData(this, bound, array_out, freqs)
341 :
342 : !Arguments ------------------------------------
343 : TYPE(FFTHyb), INTENT(IN) :: this
344 : INTEGER , INTENT(OUT) :: bound
345 : DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: array_out
346 : DOUBLE PRECISION, DIMENSION(:), OPTIONAL, INTENT(OUT) :: freqs
347 : !Local variables ------------------------------
348 : INTEGER :: i
349 : INTEGER :: size_out
350 : DOUBLE PRECISION :: re
351 : DOUBLE PRECISION :: im
352 :
353 0 : size_out = SIZE(array_out)
354 : ! IF ( SIZE(array_out) .NE. this%size/2 ) &
355 : ! CALL WARNALL("FFTHyb_getData : size_in != size")
356 0 : array_out = 0.d0
357 0 : bound = MIN(this%size/2, size_out)
358 0 : DO i=1, bound
359 0 : re = REAL(this%data_inout(this%bit_rev(i-1)))
360 0 : im = AIMAG(this%data_inout(this%bit_rev(i-1)))
361 0 : array_out(i) = SQRT(re*re + im*im)
362 : END DO
363 0 : IF ( PRESENT( freqs ) .AND. bound .LE. SIZE(freqs)) THEN
364 0 : DO i=1, bound
365 0 : freqs(i) = DBLE(i-1)/this%Ts
366 : END DO
367 : ! ELSE IF ( PRESENT( freqs ) .AND. bound .GT. SIZE(freqs) ) THEN
368 : ! CALL WARNALL("FFHyb_getData : freqs does is too small")
369 : END IF
370 :
371 0 : END SUBROUTINE FFTHyb_getData
372 : !!***
373 :
374 : !!****f* ABINIT/m_FFTHyb/FFTHyb_destroy
375 : !! NAME
376 : !! FFTHyb_destroy
377 : !!
378 : !! FUNCTION
379 : !! destroy every thing
380 : !!
381 : !! COPYRIGHT
382 : !! Copyright (C) 2013-2026 ABINIT group (J. Bieder)
383 : !! This file is distributed under the terms of the
384 : !! GNU General Public License, see ~abinit/COPYING
385 : !! or http://www.gnu.org/copyleft/gpl.txt .
386 : !!
387 : !! INPUTS
388 : !! this=FFT
389 : !!
390 : !! OUTPUT
391 : !!
392 : !! SIDE EFFECTS
393 : !!
394 : !! NOTES
395 : !!
396 : !! SOURCE
397 :
398 0 : SUBROUTINE FFTHyb_destroy(this)
399 :
400 : !Arguments ------------------------------------
401 : TYPE(FFTHyb), INTENT(INOUT) :: this
402 :
403 0 : FREEIF(this%bit_rev)
404 0 : FREEIF(this%data_inout)
405 0 : this%set = .FALSE.
406 0 : END SUBROUTINE FFTHyb_destroy
407 : !!***
408 :
409 0 : END MODULE m_FFTHyb
410 : !!***
411 :
|
