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_1d,&
14 : fft_alloc,&
15 : fft_create_plan_1d,&
16 : fft_dealloc,&
17 : fft_destroy_plan
18 : USE fft_plan, ONLY: fft_plan_type
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: twopi
21 : #include "../base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 :
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_ft'
28 :
29 : PUBLIC :: multi_fft, &
30 : fft_freqs, &
31 : fft_shift
32 :
33 : CONTAINS
34 : ! **************************************************************************************************
35 : !> \brief Naively calculates the Fourier transform - it is not the bottleneck of this calculation
36 : !> \param time_series Timestamps in atomic units of time
37 : !> \param value_series Values to be Fourier transformed - moments, field etc.
38 : !> \param result_series FT of the value series
39 : !> \param damping Applied exponential damping
40 : !> \param subtract_value Value to be subtracted from the value_series (for example initial value)
41 : !> \par History
42 : !> 10.2025 Refactored for use with multi_fft routine, moved to separate file [Stepan Marek]
43 : !> 09.2024 Initial version [Stepan Marek]
44 : !> \author Stepan Marek
45 : !> \note Uses physics ordering in frequencies, those can be constructed by fft_freq
46 : ! **************************************************************************************************
47 0 : SUBROUTINE ft_simple(time_series, value_series, result_series, damping, subtract_value)
48 : REAL(kind=dp), DIMENSION(:) :: time_series
49 : COMPLEX(kind=dp), DIMENSION(:) :: value_series, result_series
50 : REAL(kind=dp) :: damping
51 : COMPLEX(kind=dp) :: subtract_value
52 :
53 : CHARACTER(len=*), PARAMETER :: routineN = 'ft_simple'
54 :
55 : INTEGER :: handle, i, j, N, start
56 : REAL(kind=dp) :: delta_t
57 :
58 0 : CALL timeset(routineN, handle)
59 :
60 0 : N = SIZE(time_series)
61 :
62 0 : delta_t = time_series(2) - time_series(1)
63 :
64 0 : IF (MOD(N, 2) == 0) THEN
65 0 : start = -N/2
66 : ELSE
67 0 : start = -(N - 1)/2
68 : END IF
69 :
70 : ! TODO : At least OMP, but ideally even MPI parallelize, or handle this on higher level?
71 0 : DO i = 1, N
72 0 : result_series(i) = CMPLX(0.0, 0.0, kind=dp)
73 0 : DO j = 1, N
74 : result_series(i) = result_series(i) + EXP(CMPLX(0.0, twopi*(start + i - 1)*(j - 1)/N, kind=dp))* &
75 0 : EXP(-damping*delta_t*(j - 1))*(value_series(j) - subtract_value)
76 : END DO
77 : END DO
78 0 : result_series(:) = delta_t*result_series(:)
79 :
80 0 : CALL timestop(handle)
81 :
82 0 : END SUBROUTINE ft_simple
83 : ! **************************************************************************************************
84 : !> \brief Calculates the Fourier transform - couples to FFT libraries in CP2K, if available
85 : !> \param time_series Timestamps in atomic units of time
86 : !> \param value_series Values to be Fourier transformed - moments, field etc. Real only. Many series can be provided.
87 : !> \param result_series FT of the value series - complex numbers
88 : !> \param omega_series ...
89 : !> \param damping_opt Supply custom exponential damping - default is 4.0/totalTime, i.e. ratio
90 : !> of last and first element in windowed value series is reduced by e^(-4)
91 : !> \param t0_opt Carry the FT only starting from certain time - allows for exclusion of trace before
92 : !> the pulse application etc.
93 : !> \param subtract_initial_opt Subtract the value at the start of the array
94 : !> \date 10.2025
95 : !> \author Stepan Marek
96 : ! **************************************************************************************************
97 72 : SUBROUTINE multi_fft(time_series, value_series, result_series, omega_series, &
98 : damping_opt, t0_opt, subtract_initial_opt)
99 : REAL(kind=dp), DIMENSION(:) :: time_series
100 : COMPLEX(kind=dp), DIMENSION(:, :) :: value_series
101 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: result_series
102 : REAL(kind=dp), DIMENSION(:), OPTIONAL :: omega_series
103 : REAL(kind=dp), OPTIONAL :: damping_opt, t0_opt
104 : LOGICAL, OPTIONAL :: subtract_initial_opt
105 :
106 : CHARACTER(len=*), PARAMETER :: routineN = 'multi_fft'
107 :
108 : COMPLEX(kind=dp) :: subtract_value
109 : COMPLEX(kind=dp), CONTIGUOUS, DIMENSION(:), &
110 72 : POINTER :: ft_samples, samples, samples_input
111 : INTEGER :: handle, i, i0, j, nsamples, nseries, stat
112 : LOGICAL :: subtract_initial
113 : REAL(kind=dp) :: damping, delta_t, t0, t_total
114 : TYPE(fft_plan_type) :: fft_plan
115 :
116 : ! For value and result series: Index 1 - different series, Index 2 - single series entry
117 :
118 : ! Evaluate optional arguments
119 : ! Start with t0
120 72 : t0 = 0.0_dp
121 72 : IF (PRESENT(t0_opt)) t0 = t0_opt
122 72 : IF (SIZE(time_series) < 2) THEN
123 0 : CPABORT("multi_fft requires at least two time samples.")
124 : END IF
125 : ! Determine zero index
126 72 : i0 = 1
127 72 : DO i = 1, SIZE(time_series)
128 72 : IF (time_series(i) >= t0) THEN
129 : i0 = i
130 : EXIT
131 : END IF
132 : END DO
133 : ! Determine nsamples
134 72 : nsamples = SIZE(time_series) - i0 + 1
135 72 : IF (nsamples < 2) THEN
136 0 : CPABORT("multi_fft requires at least two samples in the selected time window.")
137 : END IF
138 : ! Determine total time
139 72 : t_total = time_series(SIZE(time_series)) - time_series(i0)
140 72 : delta_t = time_series(i0 + 1) - time_series(i0)
141 72 : IF (t_total /= t_total .OR. ABS(t_total) >= HUGE(t_total) .OR. t_total <= 0.0_dp) THEN
142 0 : CPABORT("multi_fft detected an abnormal total time window (NaN/Inf/non-positive).")
143 : END IF
144 72 : IF (delta_t /= delta_t .OR. ABS(delta_t) >= HUGE(delta_t) .OR. delta_t <= 0.0_dp) THEN
145 0 : CPABORT("multi_fft detected an abnormal timestep (NaN/Inf/non-positive).")
146 : END IF
147 : ! Now can determine default damping
148 72 : damping = 4.0_dp/(t_total)
149 : ! Damping option supplied in au units of time
150 72 : IF (PRESENT(damping_opt)) THEN
151 72 : IF (damping_opt > 0.0_dp) THEN
152 58 : damping = 1.0_dp/damping_opt
153 14 : ELSE IF (damping_opt == 0.0_dp) THEN
154 : ! Special case - zero damping
155 0 : damping = 0.0_dp
156 : END IF
157 : END IF
158 72 : IF (damping /= damping .OR. ABS(damping) >= HUGE(damping)) THEN
159 0 : CPABORT("multi_fft detected an abnormal damping factor (NaN/Inf).")
160 : END IF
161 : ! subtract initial
162 72 : subtract_initial = .TRUE.
163 72 : subtract_value = 0.0_dp
164 72 : IF (PRESENT(subtract_initial_opt)) subtract_initial = subtract_initial_opt
165 :
166 : ! Determine nseries
167 72 : nseries = SIZE(value_series, 1)
168 : ! Reallocate results if nsamples lower than current size
169 72 : IF (nsamples /= SIZE(result_series, 2)) THEN
170 0 : DEALLOCATE (result_series)
171 0 : ALLOCATE (result_series(nseries, nsamples), source=CMPLX(0.0, 0.0, kind=dp))
172 : END IF
173 :
174 : ! Calculate the omega series values, ordered from negative to positive
175 72 : IF (PRESENT(omega_series)) THEN
176 68 : CALL fft_freqs(nsamples, t_total, omega_series, fft_ordering_opt=.FALSE.)
177 3160 : IF (ANY(omega_series /= omega_series) .OR. &
178 : ANY(ABS(omega_series) >= HUGE(omega_series))) THEN
179 0 : CPABORT("multi_fft produced abnormal frequencies (NaN/Inf).")
180 : END IF
181 : END IF
182 :
183 : ! Use FFTW3 library
184 : ! Allocate the in-out arrays (on every rank)
185 72 : CALL timeset(routineN, handle)
186 72 : NULLIFY (samples)
187 72 : NULLIFY (samples_input)
188 72 : NULLIFY (ft_samples)
189 144 : CALL fft_alloc(samples, [nsamples*nseries])
190 144 : CALL fft_alloc(samples_input, [nsamples*nseries])
191 144 : CALL fft_alloc(ft_samples, [nsamples*nseries])
192 : ! Fill the samples with data
193 312 : DO i = 1, nseries
194 5904 : DO j = 1, nsamples
195 : ! Subtract initial value if required
196 5592 : IF (subtract_initial) THEN
197 5592 : subtract_value = value_series(i, 1)
198 : END IF
199 5592 : samples_input(j + (i - 1)*nsamples) = value_series(i, i0 + j - 1) - subtract_value
200 : ! Apply damping
201 : samples_input(j + (i - 1)*nsamples) = samples_input(j + (i - 1)*nsamples)* &
202 5832 : EXP(-damping*(time_series(i0 + j - 1) - time_series(i0)))
203 : END DO
204 : END DO
205 : ! Create the plan (this overwrites samples and ft_samples with planning data)
206 : CALL fft_create_plan_1d(fft_plan, -1, .FALSE., .FALSE., &
207 72 : nsamples, nsamples, nsamples, nseries, samples, ft_samples)
208 : ! Carry out the transform
209 : ! Scale by dt - to transform to an integral
210 72 : CALL fft_1d(fft_plan, samples_input, ft_samples, delta_t, stat)
211 72 : 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 312 : DO i = 1, nseries
226 240 : CALL fft_shift(ft_samples((i - 1)*nsamples + 1:i*nsamples))
227 5904 : 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 29520 : 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 72 : CALL fft_dealloc(samples)
238 72 : CALL fft_dealloc(ft_samples)
239 72 : CALL fft_dealloc(samples_input)
240 72 : CALL fft_destroy_plan(fft_plan)
241 72 : CALL timestop(handle)
242 288 : 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 240 : SUBROUTINE fft_shift(source)
250 : COMPLEX(kind=dp), DIMENSION(:) :: source
251 :
252 240 : 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 240 : n = SIZE(source)
261 240 : 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 228 : pos_lower(1) = (n + 1)/2
269 228 : pos_upper(2) = (n + 1)/2
270 228 : neg_lower(2) = (n + 1)/2 + 1
271 228 : neg_upper(1) = (n - 1)/2
272 : END IF
273 : ! Parity independent positions
274 240 : pos_lower(2) = 1
275 240 : pos_upper(1) = n
276 240 : neg_lower(1) = 1
277 240 : neg_upper(2) = n
278 :
279 720 : ALLOCATE (buffer(n))
280 2922 : buffer(neg_lower(1):neg_upper(1)) = source(neg_lower(2):neg_upper(2))
281 3150 : buffer(pos_lower(1):pos_upper(1)) = source(pos_lower(2):pos_upper(2))
282 5832 : source(:) = buffer(:)
283 240 : DEALLOCATE (buffer)
284 :
285 240 : 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 68 : 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 68 : fft_ordering = .FALSE.
310 68 : IF (PRESENT(fft_ordering_opt)) fft_ordering = fft_ordering_opt
311 :
312 68 : 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 68 : IF (MOD(n, 2) == 0) THEN
318 4 : start = -n/2
319 : ELSE
320 64 : start = -(n - 1)/2
321 : END IF
322 1580 : DO i = 1, n
323 1580 : 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 1580 : omegas(:) = omegas(:)*twopi/t_total
349 68 : END SUBROUTINE fft_freqs
350 :
351 : END MODULE rt_propagation_ft
|
