Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Separation of Fourier transform utilities into separate file
10 : !> \author Stepan Marek (08.24)
11 : ! **************************************************************************************************
12 : MODULE rt_propagation_ft
13 : USE fft_lib, ONLY: fft_1dm,&
14 : fft_alloc,&
15 : fft_create_plan_1dm,&
16 : fft_dealloc,&
17 : fft_destroy_plan,&
18 : fft_library
19 : USE fft_plan, ONLY: fft_plan_type
20 : USE kinds, ONLY: dp
21 : USE mathconstants, ONLY: twopi
22 : #include "../base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 :
28 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_ft'
29 :
30 : PUBLIC :: multi_fft, &
31 : fft_freqs, &
32 : fft_shift
33 :
34 : CONTAINS
35 : ! **************************************************************************************************
36 : !> \brief Naively calculates the Fourier transform - it is not the bottleneck of this calculation
37 : !> \param time_series Timestamps in atomic units of time
38 : !> \param value_series Values to be Fourier transformed - moments, field etc.
39 : !> \param result_series FT of the value series
40 : !> \param damping Applied exponential damping
41 : !> \param subtract_value Value to be subtracted from the value_series (for example initial value)
42 : !> \par History
43 : !> 10.2025 Refactored for use with multi_fft routine, moved to separate file [Stepan Marek]
44 : !> 09.2024 Initial version [Stepan Marek]
45 : !> \author Stepan Marek
46 : !> \note Uses physics ordering in frequencies, those can be constructed by fft_freq
47 : ! **************************************************************************************************
48 0 : SUBROUTINE ft_simple(time_series, value_series, result_series, damping, subtract_value)
49 : REAL(kind=dp), DIMENSION(:) :: time_series
50 : COMPLEX(kind=dp), DIMENSION(:) :: value_series, result_series
51 : REAL(kind=dp) :: damping
52 : COMPLEX(kind=dp) :: subtract_value
53 :
54 : CHARACTER(len=*), PARAMETER :: routineN = 'ft_simple'
55 :
56 : INTEGER :: handle, i, j, N, start
57 : REAL(kind=dp) :: delta_t
58 :
59 0 : CALL timeset(routineN, handle)
60 :
61 0 : N = SIZE(time_series)
62 :
63 0 : delta_t = time_series(2) - time_series(1)
64 :
65 0 : IF (MOD(N, 2) == 0) THEN
66 0 : start = -N/2
67 : ELSE
68 0 : start = -(N - 1)/2
69 : END IF
70 :
71 : ! TODO : At least OMP, but ideally even MPI parallelize, or handle this on higher level?
72 0 : DO i = 1, N
73 0 : result_series(i) = CMPLX(0.0, 0.0, kind=dp)
74 0 : DO j = 1, N
75 : result_series(i) = result_series(i) + EXP(CMPLX(0.0, twopi*(start + i - 1)*(j - 1)/N, kind=dp))* &
76 0 : EXP(-damping*delta_t*(j - 1))*(value_series(j) - subtract_value)
77 : END DO
78 : END DO
79 0 : result_series(:) = delta_t*result_series(:)
80 :
81 0 : CALL timestop(handle)
82 :
83 0 : END SUBROUTINE ft_simple
84 : ! **************************************************************************************************
85 : !> \brief Calculates the Fourier transform - couples to FFT libraries in CP2K, if available
86 : !> \param time_series Timestamps in atomic units of time
87 : !> \param value_series Values to be Fourier transformed - moments, field etc. Real only. Many series can be provided.
88 : !> \param result_series FT of the value series - complex numbers
89 : !> \param omega_series ...
90 : !> \param damping_opt Supply custom exponential damping - default is 4.0/totalTime, i.e. ratio
91 : !> of last and first element in windowed value series is reduced by e^(-4)
92 : !> \param t0_opt Carry the FT only starting from certain time - allows for exclusion of trace before
93 : !> the pulse application etc.
94 : !> \param subtract_initial_opt Subtract the value at the start of the array
95 : !> \date 10.2025
96 : !> \author Stepan Marek
97 : ! **************************************************************************************************
98 66 : SUBROUTINE multi_fft(time_series, value_series, result_series, omega_series, &
99 : damping_opt, t0_opt, subtract_initial_opt)
100 : REAL(kind=dp), DIMENSION(:) :: time_series
101 : COMPLEX(kind=dp), DIMENSION(:, :) :: value_series
102 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: result_series
103 : REAL(kind=dp), DIMENSION(:), OPTIONAL :: omega_series
104 : REAL(kind=dp), OPTIONAL :: damping_opt, t0_opt
105 : LOGICAL, OPTIONAL :: subtract_initial_opt
106 :
107 : CHARACTER(len=*), PARAMETER :: routineN = 'multi_fft'
108 :
109 : COMPLEX(kind=dp) :: subtract_value
110 : COMPLEX(kind=dp), CONTIGUOUS, DIMENSION(:), &
111 66 : POINTER :: ft_samples, samples, samples_input
112 : INTEGER :: handle, i, i0, j, nsamples, nseries, stat
113 : LOGICAL :: subtract_initial
114 : REAL(kind=dp) :: damping, delta_t, t0, t_total
115 : TYPE(fft_plan_type) :: fft_plan
116 :
117 : ! For value and result series: Index 1 - different series, Index 2 - single series entry
118 :
119 : ! Evaluate optional arguments
120 : ! Start with t0
121 66 : t0 = 0.0_dp
122 66 : IF (PRESENT(t0_opt)) t0 = t0_opt
123 66 : IF (SIZE(time_series) < 2) THEN
124 0 : CPABORT("multi_fft requires at least two time samples.")
125 : END IF
126 : ! Determine zero index
127 66 : i0 = 1
128 66 : DO i = 1, SIZE(time_series)
129 66 : IF (time_series(i) >= t0) THEN
130 : i0 = i
131 : EXIT
132 : END IF
133 : END DO
134 : ! Determine nsamples
135 66 : nsamples = SIZE(time_series) - i0 + 1
136 66 : IF (nsamples < 2) THEN
137 0 : CPABORT("multi_fft requires at least two samples in the selected time window.")
138 : END IF
139 : ! Determine total time
140 66 : t_total = time_series(SIZE(time_series)) - time_series(i0)
141 66 : delta_t = time_series(i0 + 1) - time_series(i0)
142 66 : IF (t_total /= t_total .OR. ABS(t_total) >= HUGE(t_total) .OR. t_total <= 0.0_dp) THEN
143 0 : CPABORT("multi_fft detected an abnormal total time window (NaN/Inf/non-positive).")
144 : END IF
145 66 : IF (delta_t /= delta_t .OR. ABS(delta_t) >= HUGE(delta_t) .OR. delta_t <= 0.0_dp) THEN
146 0 : CPABORT("multi_fft detected an abnormal timestep (NaN/Inf/non-positive).")
147 : END IF
148 : ! Now can determine default damping
149 66 : damping = 4.0_dp/(t_total)
150 : ! Damping option supplied in au units of time
151 66 : IF (PRESENT(damping_opt)) THEN
152 66 : IF (damping_opt > 0.0_dp) THEN
153 52 : damping = 1.0_dp/damping_opt
154 14 : ELSE IF (damping_opt == 0.0_dp) THEN
155 : ! Special case - zero damping
156 0 : damping = 0.0_dp
157 : END IF
158 : END IF
159 66 : IF (damping /= damping .OR. ABS(damping) >= HUGE(damping)) THEN
160 0 : CPABORT("multi_fft detected an abnormal damping factor (NaN/Inf).")
161 : END IF
162 : ! subtract initial
163 66 : subtract_initial = .TRUE.
164 66 : subtract_value = 0.0_dp
165 66 : IF (PRESENT(subtract_initial_opt)) subtract_initial = subtract_initial_opt
166 :
167 : ! Determine nseries
168 66 : nseries = SIZE(value_series, 1)
169 : ! Reallocate results if nsamples lower than current size
170 66 : IF (nsamples /= SIZE(result_series, 2)) THEN
171 0 : DEALLOCATE (result_series)
172 0 : ALLOCATE (result_series(nseries, nsamples), source=CMPLX(0.0, 0.0, kind=dp))
173 : END IF
174 :
175 : ! Calculate the omega series values, ordered from negative to positive
176 66 : IF (PRESENT(omega_series)) THEN
177 62 : CALL fft_freqs(nsamples, t_total, omega_series, fft_ordering_opt=.FALSE.)
178 2896 : IF (ANY(omega_series /= omega_series) .OR. &
179 : ANY(ABS(omega_series) >= HUGE(omega_series))) THEN
180 0 : CPABORT("multi_fft produced abnormal frequencies (NaN/Inf).")
181 : END IF
182 : END IF
183 :
184 : ! Use FFTW3 library
185 : ! Allocate the in-out arrays (on every rank)
186 66 : CALL timeset(routineN, handle)
187 66 : NULLIFY (samples)
188 66 : NULLIFY (samples_input)
189 66 : NULLIFY (ft_samples)
190 132 : CALL fft_alloc(samples, [nsamples*nseries])
191 132 : CALL fft_alloc(samples_input, [nsamples*nseries])
192 132 : CALL fft_alloc(ft_samples, [nsamples*nseries])
193 : ! Fill the samples with data
194 288 : DO i = 1, nseries
195 5502 : DO j = 1, nsamples
196 : ! Subtract initial value if required
197 5214 : IF (subtract_initial) THEN
198 5214 : subtract_value = value_series(i, 1)
199 : END IF
200 5214 : samples_input(j + (i - 1)*nsamples) = value_series(i, i0 + j - 1) - subtract_value
201 : ! Apply damping
202 : samples_input(j + (i - 1)*nsamples) = samples_input(j + (i - 1)*nsamples)* &
203 5436 : EXP(-damping*(time_series(i0 + j - 1) - time_series(i0)))
204 : END DO
205 : END DO
206 : ! Create the plan (this overwrites samples and ft_samples with planning data)
207 66 : CALL fft_create_plan_1dm(fft_plan, fft_library("FFTW3"), -1, .FALSE., nsamples, nseries, samples, ft_samples, 3)
208 : ! Carry out the transform
209 : ! Scale by dt - to transform to an integral
210 66 : CALL fft_1dm(fft_plan, samples_input, ft_samples, delta_t, stat)
211 66 : IF (stat /= 0) THEN
212 : ! Failed fftw3 - go to backup
213 : ! Uses value_series and result_series - no need to reassign data
214 : ! TODO : OMP parallel for different series?
215 0 : DO i = 1, nseries
216 0 : IF (subtract_initial) THEN
217 0 : subtract_value = value_series(i, 1)
218 : END IF
219 : CALL ft_simple(time_series(i0:SIZE(time_series)), &
220 : value_series(i, i0:SIZE(value_series, 2)), result_series(i, 1:nsamples), &
221 0 : damping, subtract_value)
222 : END DO
223 : ELSE
224 : ! Successful FT requires shift
225 288 : DO i = 1, nseries
226 222 : CALL fft_shift(ft_samples((i - 1)*nsamples + 1:i*nsamples))
227 5502 : result_series(i, :) = ft_samples((i - 1)*nsamples + 1:i*nsamples)
228 : END DO
229 : END IF
230 : IF (ANY(REAL(result_series, kind=dp) /= REAL(result_series, kind=dp)) .OR. &
231 : ANY(AIMAG(result_series) /= AIMAG(result_series)) .OR. &
232 27480 : ANY(ABS(REAL(result_series, kind=dp)) >= HUGE(1.0_dp)) .OR. &
233 : ANY(ABS(AIMAG(result_series)) >= HUGE(1.0_dp))) THEN
234 0 : CPABORT("multi_fft produced abnormal Fourier amplitudes (NaN/Inf).")
235 : END IF
236 : ! Deallocate
237 66 : CALL fft_dealloc(samples)
238 66 : CALL fft_dealloc(ft_samples)
239 66 : CALL fft_dealloc(samples_input)
240 66 : CALL fft_destroy_plan(fft_plan)
241 66 : CALL timestop(handle)
242 264 : END SUBROUTINE multi_fft
243 : ! **************************************************************************************************
244 : !> \brief Switches the order in result of FT, so that negative frequencies go first
245 : !> \param source Array containing the FT - buffer is used to reorder it
246 : !> \date 10.2025
247 : !> \author Stepan Marek
248 : ! **************************************************************************************************
249 222 : SUBROUTINE fft_shift(source)
250 : COMPLEX(kind=dp), DIMENSION(:) :: source
251 :
252 222 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer
253 : INTEGER :: n
254 : INTEGER, DIMENSION(2) :: neg_lower, neg_upper, pos_lower, &
255 : pos_upper
256 :
257 : ! Boundary indices for positive/negative part of the spectrum
258 : ! Index 1 : 1 = transformed order, 2 = FFT order
259 :
260 222 : n = SIZE(source)
261 222 : IF (MOD(n, 2) == 0) THEN
262 : ! Even case
263 12 : pos_lower(1) = n/2 + 1
264 12 : pos_upper(2) = n/2
265 12 : neg_lower(2) = n/2 + 1
266 12 : neg_upper(1) = n/2
267 : ELSE
268 210 : pos_lower(1) = (n + 1)/2
269 210 : pos_upper(2) = (n + 1)/2
270 210 : neg_lower(2) = (n + 1)/2 + 1
271 210 : neg_upper(1) = (n - 1)/2
272 : END IF
273 : ! Parity independent positions
274 222 : pos_lower(2) = 1
275 222 : pos_upper(1) = n
276 222 : neg_lower(1) = 1
277 222 : neg_upper(2) = n
278 :
279 666 : ALLOCATE (buffer(n))
280 2724 : buffer(neg_lower(1):neg_upper(1)) = source(neg_lower(2):neg_upper(2))
281 2934 : buffer(pos_lower(1):pos_upper(1)) = source(pos_lower(2):pos_upper(2))
282 5436 : source(:) = buffer(:)
283 222 : DEALLOCATE (buffer)
284 :
285 222 : END SUBROUTINE fft_shift
286 : ! **************************************************************************************************
287 : !> \brief Switches the order in result of FT, so that negative frequencies go first
288 : !> \param n Number of frequencies
289 : !> \param t_total Total corresponding propagation time
290 : !> \param omegas Array of frequencies
291 : !> \param fft_ordering_opt Whether to switch to FFT ordering
292 : !> \date 10.2025
293 : !> \author Stepan Marek
294 : ! **************************************************************************************************
295 62 : SUBROUTINE fft_freqs(n, t_total, omegas, fft_ordering_opt)
296 : ! Number of FT samples
297 : INTEGER :: n
298 : REAL(kind=dp) :: t_total
299 : REAL(kind=dp), DIMENSION(:) :: omegas
300 : LOGICAL, OPTIONAL :: fft_ordering_opt
301 :
302 : INTEGER :: finish, i, start
303 : LOGICAL :: fft_ordering
304 :
305 : ! Total window time, dt = nsamples / t_total
306 :
307 : ! Determine the order, by default, use physics order,
308 : ! i.e. negative frequencies before positive ones
309 62 : fft_ordering = .FALSE.
310 62 : IF (PRESENT(fft_ordering_opt)) fft_ordering = fft_ordering_opt
311 :
312 62 : IF (.NOT. fft_ordering) THEN
313 : ! Physics order case
314 : ! Unit frequencies at
315 : ! - for even n : -n/2, -n/2 + 1, -n/2 + 2, ..., -1, 0, 1, ..., n/2 - 1
316 : ! - for odd n : -(n-1)/2, -(n-1)/2 + 1, ..., -1, 0, 1, ..., (n-1)/2
317 62 : IF (MOD(n, 2) == 0) THEN
318 4 : start = -n/2
319 : ELSE
320 58 : start = -(n - 1)/2
321 : END IF
322 1448 : DO i = 1, n
323 1448 : omegas(i) = start + i - 1
324 : END DO
325 : ELSE
326 : ! FFT order case
327 : ! Unit frequencies at
328 : ! - for even n : 0, 1, ..., n/2 - 1, -n/2, -n/2 + 1, -n/2 + 2, ..., -1
329 : ! - for odd n : 0, 1, ..., (n-1)/2, -(n-1)/2, -(n-1)/2 + 1, ..., -1
330 0 : IF (MOD(n, 2) == 0) THEN
331 0 : finish = n/2 - 1
332 0 : start = -n/2
333 : ELSE
334 0 : finish = (n - 1)/2
335 0 : start = -(n - 1)/2
336 : END IF
337 : ! Positive frequencies
338 0 : DO i = 1, finish + 1
339 0 : omegas(i) = (i - 1)
340 : END DO
341 : ! Negative frequencies
342 0 : DO i = finish + 2, n
343 0 : omegas(i) = start + i - finish - 2
344 : END DO
345 : END IF
346 :
347 : ! Finally, multiply by the factor to translate to angular frequency
348 1448 : omegas(:) = omegas(:)*twopi/t_total
349 62 : END SUBROUTINE fft_freqs
350 :
351 : END MODULE rt_propagation_ft
|
