(git:691081d)
Loading...
Searching...
No Matches
time_frequency_grids.F
Go to the documentation of this file.
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 Definition and construction of time/frequency grids for correlation methods.
10!> \par History
11!> 05.2019 Refactored from rpa_ri_gpw [Frederick Stein]
12! **************************************************************************************************
15 USE kinds, ONLY: dp
16 USE mathconstants, ONLY: pi
21#include "./base/base_uses.f90"
22
23 IMPLICIT NONE
24
25 PRIVATE
26
27 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'time_frequency_grids'
28
30 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: frequency
31 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: frequency_weights
32 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: imaginary_time
33 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: time_weights_at_zero_frequency
34 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: cosine_time_to_frequency_weights
35 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: cosine_frequency_to_time_weights
36 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: sine_time_to_frequency_weights
38
39 INTEGER, PARAMETER :: cosine_time_to_frequency = 1, &
40 sine_time_to_frequency = 2, &
41 cosine_frequency_to_time = 3
42
45
46CONTAINS
47
48! **************************************************************************************************
49!> \brief Release all data owned by a time_frequency_grid_type object.
50!> \param grid The grid object to release.
51! **************************************************************************************************
53
54 TYPE(time_frequency_grid_type), INTENT(INOUT) :: grid
55
56 IF (ALLOCATED(grid%frequency)) DEALLOCATE (grid%frequency)
57 IF (ALLOCATED(grid%frequency_weights)) DEALLOCATE (grid%frequency_weights)
58 IF (ALLOCATED(grid%imaginary_time)) DEALLOCATE (grid%imaginary_time)
59 IF (ALLOCATED(grid%time_weights_at_zero_frequency)) DEALLOCATE (grid%time_weights_at_zero_frequency)
60 IF (ALLOCATED(grid%cosine_time_to_frequency_weights)) DEALLOCATE (grid%cosine_time_to_frequency_weights)
61 IF (ALLOCATED(grid%cosine_frequency_to_time_weights)) DEALLOCATE (grid%cosine_frequency_to_time_weights)
62 IF (ALLOCATED(grid%sine_time_to_frequency_weights)) DEALLOCATE (grid%sine_time_to_frequency_weights)
63
64 END SUBROUTINE time_frequency_grid_release
65
66! **************************************************************************************************
67!> \brief Build a Clenshaw-Curtis frequency grid.
68!> \param num_points Number of integration points.
69!> \param grid Grid object to fill.
70! **************************************************************************************************
71 SUBROUTINE build_clenshaw_grid(num_points, grid)
72
73 INTEGER, INTENT(IN) :: num_points
74 TYPE(time_frequency_grid_type), INTENT(OUT) :: grid
75
76 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_clenshaw_grid'
77
78 INTEGER :: handle, jquad
79
80 CALL timeset(routinen, handle)
81
82 cpassert(num_points > 0)
83
84 ALLOCATE (grid%frequency(num_points), grid%frequency_weights(num_points))
85 grid%frequency = 0.0_dp
86 grid%frequency_weights = 0.0_dp
87
88 DO jquad = 1, num_points - 1
89 grid%frequency(jquad) = jquad*pi/(2.0_dp*num_points)
90 grid%frequency_weights(jquad) = pi/(num_points*sin(grid%frequency(jquad))**2)
91 END DO
92 grid%frequency(num_points) = pi/2.0_dp
93 grid%frequency_weights(num_points) = &
94 pi/(2.0_dp*num_points*sin(grid%frequency(num_points))**2)
95
96 CALL timestop(handle)
97
98 END SUBROUTINE build_clenshaw_grid
99
100! **************************************************************************************************
101!> \brief Build a minimax time/frequency grid through the common backend boundary.
102!> \param num_points Number of minimax points.
103!> \param energy_min Lower end of the physical energy interval.
104!> \param energy_max Upper end of the physical energy interval.
105!> \param regularization Regularization used for the fitted transform weights.
106!> \param num_points_per_magnitude Number of fitting points per decade.
107!> \param grid Grid object to fill.
108!> \param build_frequency Whether to construct the frequency components.
109!> \param build_time Whether to construct the imaginary-time components.
110!> \param build_transforms Whether to construct the time/frequency transform weights.
111!> \param build_sine Whether to construct the sine time-to-frequency weights.
112!> \param time_scaling Scaling applied to the imaginary-time abscissas.
113!> \param time_weight_scaling Scaling applied to the imaginary-time weights.
114!> \param max_fit_error Maximum fitting error across all requested transforms.
115!> \param print_warning Whether the minimax coefficient routine prints warnings.
116!> \param unit_nr Output unit used by an external backend.
117!> \param prefer_external_backend Whether to try an available external backend first.
118!> \param used_external_backend Whether the external backend supplied the grid.
119! **************************************************************************************************
120 SUBROUTINE build_minimax_time_frequency_grid(num_points, energy_min, energy_max, regularization, &
121 num_points_per_magnitude, grid, build_frequency, build_time, &
122 build_transforms, build_sine, time_scaling, time_weight_scaling, &
123 max_fit_error, print_warning, unit_nr, prefer_external_backend, &
124 used_external_backend)
125
126 INTEGER, INTENT(IN) :: num_points
127 REAL(kind=dp), INTENT(IN) :: energy_min, energy_max, regularization
128 INTEGER, INTENT(IN) :: num_points_per_magnitude
129 TYPE(time_frequency_grid_type), INTENT(OUT) :: grid
130 LOGICAL, INTENT(IN) :: build_frequency, build_time, &
131 build_transforms, build_sine
132 REAL(kind=dp), INTENT(IN) :: time_scaling, time_weight_scaling
133 REAL(kind=dp), INTENT(OUT) :: max_fit_error
134 LOGICAL, INTENT(IN) :: print_warning
135 INTEGER, INTENT(IN) :: unit_nr
136 LOGICAL, INTENT(IN) :: prefer_external_backend
137 LOGICAL, INTENT(OUT), OPTIONAL :: used_external_backend
138
139 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_minimax_time_frequency_grid'
140
141 INTEGER :: external_ierr, handle, ierr
142 REAL(kind=dp) :: e_range, max_error
143 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: points_and_weights
144
145 CALL timeset(routinen, handle)
146
147 IF (PRESENT(used_external_backend)) used_external_backend = .false.
148 IF (prefer_external_backend) THEN
149 CALL greenx_get_minimax_grid(unit_nr, num_points, energy_min, energy_max, regularization, &
150 grid%imaginary_time, grid%time_weights_at_zero_frequency, &
151 grid%frequency, grid%frequency_weights, &
152 grid%cosine_time_to_frequency_weights, &
153 grid%cosine_frequency_to_time_weights, &
154 grid%sine_time_to_frequency_weights, external_ierr)
155 IF (external_ierr == 0) THEN
156 max_fit_error = 0.0_dp
157 IF (PRESENT(used_external_backend)) used_external_backend = .true.
158 CALL timestop(handle)
159 RETURN
160 END IF
162 END IF
163
164 cpassert(num_points > 0)
165 cpassert(energy_min > 0.0_dp)
166 cpassert(energy_max >= energy_min)
167 cpassert(time_scaling > 0.0_dp)
168 cpassert(time_weight_scaling > 0.0_dp)
169 cpassert(.NOT. build_transforms .OR. (build_frequency .AND. build_time))
170 cpassert(.NOT. build_sine .OR. build_transforms)
171 max_error = 0.0_dp
172
173 e_range = energy_max/energy_min
174 ALLOCATE (points_and_weights(2*num_points))
175
176 IF (build_frequency) THEN
177 IF (num_points <= 20) THEN
178 CALL get_rpa_minimax_coeff(num_points, e_range, points_and_weights, ierr, print_warning)
179 ELSE
180 CALL get_rpa_minimax_coeff_larger_grid(num_points, e_range, points_and_weights)
181 END IF
182
183 ALLOCATE (grid%frequency(num_points))
184 ALLOCATE (grid%frequency_weights(num_points))
185 grid%frequency(:) = points_and_weights(1:num_points)*energy_min
186 grid%frequency_weights(:) = points_and_weights(num_points + 1:)*energy_min
187 IF (num_points >= 26) grid%frequency_weights(:) = grid%frequency_weights(:)*4.0_dp
188 END IF
189
190 IF (build_time) THEN
191 IF (num_points <= 20) THEN
192 CALL get_exp_minimax_coeff(num_points, e_range, points_and_weights)
193 ELSE
194 CALL get_exp_minimax_coeff_gw(num_points, e_range, points_and_weights)
195 END IF
196
197 ALLOCATE (grid%imaginary_time(num_points))
198 ALLOCATE (grid%time_weights_at_zero_frequency(num_points))
199 grid%imaginary_time(:) = points_and_weights(1:num_points)/time_scaling/energy_min
200 grid%time_weights_at_zero_frequency(:) = points_and_weights(num_points + 1:)/time_weight_scaling/energy_min
201
202 IF (build_transforms) THEN
203 ALLOCATE (grid%cosine_time_to_frequency_weights(num_points, num_points))
204 CALL fit_l_sq_weights(cosine_time_to_frequency, num_points, grid%imaginary_time, &
205 grid%cosine_time_to_frequency_weights, grid%frequency, energy_min, energy_max, &
206 max_error, num_points_per_magnitude, regularization)
207
208 ALLOCATE (grid%cosine_frequency_to_time_weights(num_points, num_points))
209 CALL fit_l_sq_weights(cosine_frequency_to_time, num_points, grid%imaginary_time, &
210 grid%cosine_frequency_to_time_weights, grid%frequency, energy_min, energy_max, &
211 max_error, num_points_per_magnitude, regularization)
212
213 IF (build_sine) THEN
214 ALLOCATE (grid%sine_time_to_frequency_weights(num_points, num_points))
215 CALL fit_l_sq_weights(sine_time_to_frequency, num_points, grid%imaginary_time, &
216 grid%sine_time_to_frequency_weights, grid%frequency, energy_min, energy_max, &
217 max_error, num_points_per_magnitude, regularization)
218 END IF
219 END IF
220 END IF
221
222 max_fit_error = max_error
223
224 DEALLOCATE (points_and_weights)
225
226 CALL timestop(handle)
227
229
230! **************************************************************************************************
231!> \brief Calculate least-squares weights for a time/frequency transform.
232!> \param transform_kind Type of transform to fit.
233!> \param num_integ_points Number of integration points.
234!> \param tau_tj Imaginary-time integration points.
235!> \param weights Transform weights to construct.
236!> \param omega_tj Frequency integration points.
237!> \param E_min Lower end of the fitting interval.
238!> \param E_max Upper end of the fitting interval.
239!> \param max_error Maximum fitting error.
240!> \param num_points_per_magnitude Number of fitting points per decade.
241!> \param regularization Regularization parameter for the pseudoinverse.
242! **************************************************************************************************
243 SUBROUTINE fit_l_sq_weights(transform_kind, num_integ_points, tau_tj, weights, omega_tj, E_min, E_max, &
244 max_error, num_points_per_magnitude, regularization)
245
246 INTEGER, INTENT(IN) :: transform_kind, num_integ_points
247 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: tau_tj
248 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: weights
249 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: omega_tj
250 REAL(kind=dp), INTENT(IN) :: e_min, e_max
251 REAL(kind=dp), INTENT(OUT) :: max_error
252 INTEGER, INTENT(IN) :: num_points_per_magnitude
253 REAL(kind=dp), INTENT(IN) :: regularization
254
255 CHARACTER(LEN=*), PARAMETER :: routinen = 'fit_l_sq_weights'
256
257 INTEGER :: handle, iii, info, jjj, jquad, lwork, &
258 num_x_nodes
259 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
260 LOGICAL :: sine_transform, time_to_frequency
261 REAL(kind=dp) :: func_val, max_error_tmp, multiplicator, &
262 omega, tau, x_value
263 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, vec_uty, weight_work, work, &
264 x_values, y_values
265 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_a, mat_sinvvsinvsigma, &
266 mat_sinvvsinvt, mat_u
267
268 CALL timeset(routinen, handle)
269
270 cpassert(transform_kind >= cosine_time_to_frequency)
271 cpassert(transform_kind <= cosine_frequency_to_time)
272 cpassert(SIZE(tau_tj) == num_integ_points)
273 cpassert(SIZE(omega_tj) == num_integ_points)
274 cpassert(SIZE(weights, 1) == num_integ_points)
275 cpassert(SIZE(weights, 2) == num_integ_points)
276
277 time_to_frequency = transform_kind /= cosine_frequency_to_time
278 sine_transform = transform_kind == sine_time_to_frequency
279
280 ! take num_points_per_magnitude points per 10-interval
281 num_x_nodes = (int(log10(e_max/e_min)) + 1)*num_points_per_magnitude
282
283 ! take at least as many x points as integration points to have clear
284 ! input for the singular value decomposition
285 num_x_nodes = max(num_x_nodes, num_integ_points)
286
287 ALLOCATE (x_values(num_x_nodes), y_values(num_x_nodes), mat_a(num_x_nodes, num_integ_points), &
288 weight_work(num_integ_points), sing_values(num_integ_points), mat_u(num_x_nodes, num_x_nodes), &
289 mat_sinvvsinvt(num_x_nodes, num_integ_points), work(8*num_integ_points*num_integ_points + &
290 12*num_integ_points + 2*num_x_nodes), iwork(8*num_integ_points), &
291 mat_sinvvsinvsigma(num_integ_points, num_x_nodes), vec_uty(num_x_nodes))
292 mat_sinvvsinvsigma = 0.0_dp
293
294 ! double the value nessary for 'A' to achieve good performance
295 lwork = SIZE(work)
296
297 multiplicator = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
298 DO iii = 1, num_x_nodes
299 x_values(iii) = e_min*multiplicator**(iii - 1)
300 END DO
301
302 max_error = 0.0_dp
303
304 DO jquad = 1, num_integ_points
305
306 IF (time_to_frequency) THEN
307 omega = omega_tj(jquad)
308
309 IF (sine_transform) THEN
310 ! y=2*omega/(x^2+omega^2)
311 DO iii = 1, num_x_nodes
312 y_values(iii) = 2.0_dp*omega/(x_values(iii)**2 + omega**2)
313 END DO
314 ELSE
315 ! y=2*x/(x^2+omega^2)
316 DO iii = 1, num_x_nodes
317 y_values(iii) = 2.0_dp*x_values(iii)/(x_values(iii)**2 + omega**2)
318 END DO
319 END IF
320
321 DO jjj = 1, num_integ_points
322 DO iii = 1, num_x_nodes
323 IF (sine_transform) THEN
324 mat_a(iii, jjj) = sin(omega*tau_tj(jjj))*exp(-x_values(iii)*tau_tj(jjj))
325 ELSE
326 mat_a(iii, jjj) = cos(omega*tau_tj(jjj))*exp(-x_values(iii)*tau_tj(jjj))
327 END IF
328 END DO
329 END DO
330 ELSE
331 tau = tau_tj(jquad)
332
333 ! y=exp(-x*|tau|)
334 DO iii = 1, num_x_nodes
335 y_values(iii) = exp(-x_values(iii)*tau)
336 END DO
337
338 DO jjj = 1, num_integ_points
339 omega = omega_tj(jjj)
340 DO iii = 1, num_x_nodes
341 x_value = x_values(iii)
342 mat_a(iii, jjj) = cos(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2)
343 END DO
344 END DO
345 END IF
346
347 ! Singular value decomposition of mat_A
348 CALL dgesdd('A', num_x_nodes, num_integ_points, mat_a, num_x_nodes, sing_values, mat_u, num_x_nodes, &
349 mat_sinvvsinvt, num_x_nodes, work, lwork, iwork, info)
350 cpassert(info == 0)
351
352 ! integration weights = V Sigma U^T y
353 DO jjj = 1, num_integ_points
354 DO iii = 1, num_integ_points
355 mat_sinvvsinvsigma(iii, jjj) = mat_sinvvsinvt(jjj, iii)*sing_values(jjj) &
356 /(regularization**2 + sing_values(jjj)**2)
357 END DO
358 END DO
359
360 CALL dgemm('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_u, num_x_nodes, y_values, num_x_nodes, &
361 0.0_dp, vec_uty, num_x_nodes)
362 CALL dgemm('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_sinvvsinvsigma, num_integ_points, &
363 vec_uty, num_x_nodes, 0.0_dp, weight_work, num_integ_points)
364
365 weights(jquad, :) = weight_work(:)
366
367 max_error_tmp = 0.0_dp
368 DO iii = 1, num_x_nodes
369 func_val = 0.0_dp
370 IF (time_to_frequency) THEN
371 DO jjj = 1, num_integ_points
372 IF (sine_transform) THEN
373 func_val = func_val + weight_work(jjj)*sin(omega*tau_tj(jjj))*exp(-x_values(iii)*tau_tj(jjj))
374 ELSE
375 func_val = func_val + weight_work(jjj)*cos(omega*tau_tj(jjj))*exp(-x_values(iii)*tau_tj(jjj))
376 END IF
377 END DO
378 ELSE
379 DO jjj = 1, num_integ_points
380 omega = omega_tj(jjj)
381 func_val = func_val + weight_work(jjj)*cos(tau*omega)*2.0_dp*x_values(iii) &
382 /(x_values(iii)**2 + omega**2)
383 END DO
384 END IF
385 max_error_tmp = max(max_error_tmp, abs(y_values(iii) - func_val))
386 END DO
387 max_error = max(max_error, max_error_tmp)
388
389 END DO
390
391 DEALLOCATE (x_values, y_values, mat_a, weight_work, sing_values, mat_u, mat_sinvvsinvt, work, iwork, &
392 mat_sinvvsinvsigma, vec_uty)
393
394 CALL timestop(handle)
395
396 END SUBROUTINE fit_l_sq_weights
397
398! **************************************************************************************************
399!> \brief test the singular value decomposition for the computation of integration weights for the
400!> Fourier transform between time and frequency grid in cubic-scaling RPA
401!> \param nR ...
402!> \param iw ...
403! **************************************************************************************************
404 SUBROUTINE test_least_square_ft(nR, iw)
405 INTEGER, INTENT(IN) :: nr, iw
406
407 INTEGER :: ierr, ir, num_integ_points
408 REAL(kind=dp) :: max_error, multiplicator, rc, rc_max
409 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tau_tj, tj, x_tw
410 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: weights_cos_tf_t_to_w
411
412 rc_max = 1.0e+7
413
414 multiplicator = rc_max**(1.0_dp/(real(nr, kind=dp) - 1.0_dp))
415
416 DO num_integ_points = 1, 20
417
418 ALLOCATE (x_tw(2*num_integ_points))
419 x_tw = 0.0_dp
420 ALLOCATE (tau_tj(num_integ_points))
421 tau_tj = 0.0_dp
422 ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
423 weights_cos_tf_t_to_w = 0.0_dp
424 ALLOCATE (tj(num_integ_points))
425 tj = 0.0_dp
426
427 DO ir = 0, nr - 1
428
429 rc = 2.0_dp*multiplicator**ir
430
431 ierr = 0
432 CALL get_rpa_minimax_coeff(num_integ_points, rc, x_tw, ierr, print_warning=.false.)
433
434 tj(:) = x_tw(1:num_integ_points)
435
436 x_tw = 0.0_dp
437
438 CALL get_exp_minimax_coeff(num_integ_points, rc, x_tw)
439
440 tau_tj(:) = x_tw(1:num_integ_points)/2.0_dp
441
442 CALL fit_l_sq_weights(cosine_time_to_frequency, num_integ_points, tau_tj, &
443 weights_cos_tf_t_to_w, tj, 1.0_dp, rc, max_error, 200, 0.0_dp)
444
445 IF (iw > 0) THEN
446 WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, rc, max_error
447 END IF
448
449 END DO
450
451 DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tj)
452
453 END DO
454
455 END SUBROUTINE test_least_square_ft
456
457END MODULE time_frequency_grids
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Interface to the Greenx library.
subroutine, public greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, regularization_minimax, imaginary_time, time_weights_at_zero_frequency, frequency, frequency_weights, cosine_time_to_frequency_weights, cosine_frequency_to_time_weights, sine_time_to_frequency_weights, ierr)
Get a minimax grid from GreenX when it is available.
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
subroutine, public get_exp_minimax_coeff_gw(k, e_range, aw)
...
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition minimax_exp.F:29
subroutine, public get_exp_minimax_coeff(k, rc, aw, mm_error, which_coeffs)
Get best minimax approximation for given input parameters. Automatically chooses the most exact set o...
Routines to calculate the minimax coefficients for approximating 1/x as 1/x ~ 1/pi SUM_{i}...
Definition minimax_rpa.F:14
subroutine, public get_rpa_minimax_coeff_larger_grid(k, e_range, aw)
...
subroutine, public get_rpa_minimax_coeff(k, e_range, aw, ierr, print_warning)
The a_i and w_i coefficient are stored in aw such that the first 1:K elements correspond to a_i and t...
Definition minimax_rpa.F:41
Definition and construction of time/frequency grids for correlation methods.
integer, parameter cosine_time_to_frequency
subroutine, public build_minimax_time_frequency_grid(num_points, energy_min, energy_max, regularization, num_points_per_magnitude, grid, build_frequency, build_time, build_transforms, build_sine, time_scaling, time_weight_scaling, max_fit_error, print_warning, unit_nr, prefer_external_backend, used_external_backend)
Build a minimax time/frequency grid through the common backend boundary.
subroutine, public build_clenshaw_grid(num_points, grid)
Build a Clenshaw-Curtis frequency grid.
subroutine, public test_least_square_ft(nr, iw)
test the singular value decomposition for the computation of integration weights for the Fourier tran...
subroutine, public time_frequency_grid_release(grid)
Release all data owned by a time_frequency_grid_type object.