Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 12 : 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 12 : POINTER :: ft_samples, samples, samples_input
112 : INTEGER :: handle, i, i0, j, nsamples, nseries, stat
113 : LOGICAL :: subtract_initial
114 : REAL(kind=dp) :: damping, 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 12 : t0 = 0.0_dp
122 12 : IF (PRESENT(t0_opt)) t0 = t0_opt
123 : ! Determine zero index
124 12 : i0 = 1
125 12 : DO i = 1, SIZE(time_series)
126 12 : IF (time_series(i) >= t0) THEN
127 : i0 = i
128 : EXIT
129 : END IF
130 : END DO
131 : ! Determine nsamples
132 12 : nsamples = SIZE(time_series) - i0 + 1
133 : ! Determine total time
134 12 : t_total = time_series(SIZE(time_series)) - time_series(i0)
135 : ! Now can determine default damping
136 12 : damping = 4.0_dp/(t_total)
137 : ! Damping option supplied in au units of time
138 12 : IF (PRESENT(damping_opt)) THEN
139 12 : IF (damping_opt > 0.0_dp) THEN
140 0 : damping = 1.0_dp/damping_opt
141 12 : ELSE IF (damping_opt == 0.0_dp) THEN
142 : ! Special case - zero damping
143 0 : damping = 0.0_dp
144 : END IF
145 : END IF
146 : ! subtract initial
147 12 : subtract_initial = .TRUE.
148 12 : subtract_value = 0.0_dp
149 12 : IF (PRESENT(subtract_initial_opt)) subtract_initial = subtract_initial_opt
150 :
151 : ! Determine nseries
152 12 : nseries = SIZE(value_series, 1)
153 : ! Reallocate results if nsamples lower than current size
154 12 : IF (nsamples /= SIZE(result_series, 2)) THEN
155 0 : DEALLOCATE (result_series)
156 0 : ALLOCATE (result_series(nseries, nsamples), source=CMPLX(0.0, 0.0, kind=dp))
157 : END IF
158 :
159 : ! Calculate the omega series values, ordered from negative to positive
160 12 : IF (PRESENT(omega_series)) THEN
161 10 : CALL fft_freqs(nsamples, t_total, omega_series, fft_ordering_opt=.FALSE.)
162 : END IF
163 :
164 : ! Use FFTW3 library
165 : ! Allocate the in-out arrays (on every rank)
166 12 : CALL timeset(routineN, handle)
167 12 : NULLIFY (samples)
168 12 : NULLIFY (samples_input)
169 12 : NULLIFY (ft_samples)
170 24 : CALL fft_alloc(samples, [nsamples*nseries])
171 24 : CALL fft_alloc(samples_input, [nsamples*nseries])
172 24 : CALL fft_alloc(ft_samples, [nsamples*nseries])
173 : ! Fill the samples with data
174 48 : DO i = 1, nseries
175 4884 : DO j = 1, nsamples
176 : ! Subtract initial value if required
177 4836 : IF (subtract_initial) THEN
178 4836 : subtract_value = value_series(i, 1)
179 : END IF
180 4836 : samples_input(j + (i - 1)*nsamples) = value_series(i, i0 + j - 1) - subtract_value
181 : ! Apply damping
182 : samples_input(j + (i - 1)*nsamples) = samples_input(j + (i - 1)*nsamples)* &
183 4872 : EXP(-damping*(time_series(i0 + j - 1) - time_series(i0)))
184 : END DO
185 : END DO
186 : ! Create the plan (this overwrites samples and ft_samples with planning data)
187 12 : CALL fft_create_plan_1dm(fft_plan, fft_library("FFTW3"), -1, .FALSE., nsamples, nseries, samples, ft_samples, 3)
188 : ! Carry out the transform
189 : ! Scale by dt - to transform to an integral
190 12 : CALL fft_1dm(fft_plan, samples_input, ft_samples, time_series(2) - time_series(1), stat)
191 12 : IF (stat /= 0) THEN
192 : ! Failed fftw3 - go to backup
193 : ! Uses value_series and result_series - no need to reassign data
194 : ! TODO : OMP parallel for different series?
195 0 : DO i = 1, nseries
196 0 : IF (subtract_initial) THEN
197 0 : subtract_value = value_series(i, 1)
198 : END IF
199 : CALL ft_simple(time_series(i0:SIZE(time_series)), &
200 : value_series(i, i0:SIZE(value_series, 2)), result_series(i, 1:nsamples), &
201 0 : damping, subtract_value)
202 : END DO
203 : ELSE
204 : ! Successful FT requires shift
205 48 : DO i = 1, nseries
206 36 : CALL fft_shift(ft_samples((i - 1)*nsamples + 1:i*nsamples))
207 4884 : result_series(i, :) = ft_samples((i - 1)*nsamples + 1:i*nsamples)
208 : END DO
209 : END IF
210 : ! Deallocate
211 12 : CALL fft_dealloc(samples)
212 12 : CALL fft_dealloc(ft_samples)
213 12 : CALL fft_dealloc(samples_input)
214 12 : CALL fft_destroy_plan(fft_plan)
215 12 : CALL timestop(handle)
216 48 : END SUBROUTINE multi_fft
217 : ! **************************************************************************************************
218 : !> \brief Switches the order in result of FT, so that negative frequencies go first
219 : !> \param source Array containing the FT - buffer is used to reorder it
220 : !> \date 10.2025
221 : !> \author Stepan Marek
222 : ! **************************************************************************************************
223 36 : SUBROUTINE fft_shift(source)
224 : COMPLEX(kind=dp), DIMENSION(:) :: source
225 :
226 36 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer
227 : INTEGER :: n
228 : INTEGER, DIMENSION(2) :: neg_lower, neg_upper, pos_lower, &
229 : pos_upper
230 :
231 : ! Boundary indices for positive/negative part of the spectrum
232 : ! Index 1 : 1 = transformed order, 2 = FFT order
233 :
234 36 : n = SIZE(source)
235 36 : IF (MOD(n, 2) == 0) THEN
236 : ! Even case
237 0 : pos_lower(1) = n/2 + 1
238 0 : pos_upper(2) = n/2
239 0 : neg_lower(2) = n/2 + 1
240 0 : neg_upper(1) = n/2
241 : ELSE
242 36 : pos_lower(1) = (n + 1)/2
243 36 : pos_upper(2) = (n + 1)/2
244 36 : neg_lower(2) = (n + 1)/2 + 1
245 36 : neg_upper(1) = (n - 1)/2
246 : END IF
247 : ! Parity independent positions
248 36 : pos_lower(2) = 1
249 36 : pos_upper(1) = n
250 36 : neg_lower(1) = 1
251 36 : neg_upper(2) = n
252 :
253 108 : ALLOCATE (buffer(n))
254 2436 : buffer(neg_lower(1):neg_upper(1)) = source(neg_lower(2):neg_upper(2))
255 2472 : buffer(pos_lower(1):pos_upper(1)) = source(pos_lower(2):pos_upper(2))
256 4872 : source(:) = buffer(:)
257 36 : DEALLOCATE (buffer)
258 :
259 36 : END SUBROUTINE fft_shift
260 : ! **************************************************************************************************
261 : !> \brief Switches the order in result of FT, so that negative frequencies go first
262 : !> \param n Number of frequencies
263 : !> \param t_total Total corresponding propagation time
264 : !> \param omegas Array of frequencies
265 : !> \param fft_ordering_opt Whether to switch to FFT ordering
266 : !> \date 10.2025
267 : !> \author Stepan Marek
268 : ! **************************************************************************************************
269 10 : SUBROUTINE fft_freqs(n, t_total, omegas, fft_ordering_opt)
270 : ! Number of FT samples
271 : INTEGER :: n
272 : REAL(kind=dp) :: t_total
273 : REAL(kind=dp), DIMENSION(:) :: omegas
274 : LOGICAL, OPTIONAL :: fft_ordering_opt
275 :
276 : INTEGER :: finish, i, start
277 : LOGICAL :: fft_ordering
278 :
279 : ! Total window time, dt = nsamples / t_total
280 :
281 : ! Determine the order, by default, use physics order,
282 : ! i.e. negative frequencies before positive ones
283 10 : fft_ordering = .FALSE.
284 10 : IF (PRESENT(fft_ordering_opt)) fft_ordering = fft_ordering_opt
285 :
286 10 : IF (.NOT. fft_ordering) THEN
287 : ! Physics order case
288 : ! Unit frequencies at
289 : ! - for even n : -n/2, -n/2 + 1, -n/2 + 2, ..., -1, 0, 1, ..., n/2 - 1
290 : ! - for odd n : -(n-1)/2, -(n-1)/2 + 1, ..., -1, 0, 1, ..., (n-1)/2
291 10 : IF (MOD(n, 2) == 0) THEN
292 0 : start = -n/2
293 : ELSE
294 10 : start = -(n - 1)/2
295 : END IF
296 1520 : DO i = 1, n
297 1520 : omegas(i) = start + i - 1
298 : END DO
299 : ELSE
300 : ! FFT order case
301 : ! Unit frequencies at
302 : ! - for even n : 0, 1, ..., n/2 - 1, -n/2, -n/2 + 1, -n/2 + 2, ..., -1
303 : ! - for odd n : 0, 1, ..., (n-1)/2, -(n-1)/2, -(n-1)/2 + 1, ..., -1
304 0 : IF (MOD(n, 2) == 0) THEN
305 0 : finish = n/2 - 1
306 0 : start = -n/2
307 : ELSE
308 0 : finish = (n - 1)/2
309 0 : start = -(n - 1)/2
310 : END IF
311 : ! Positive frequencies
312 0 : DO i = 1, finish + 1
313 0 : omegas(i) = (i - 1)
314 : END DO
315 : ! Negative frequencies
316 0 : DO i = finish + 2, n
317 0 : omegas(i) = start + i - finish - 2
318 : END DO
319 : END IF
320 :
321 : ! Finally, multiply by the factor to translate to angular frequency
322 1520 : omegas(:) = omegas(:)*twopi/t_total
323 10 : END SUBROUTINE fft_freqs
324 :
325 : END MODULE rt_propagation_ft
|
