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 Definition and construction of time/frequency grids for correlation methods.
10 : !> \par History
11 : !> 05.2019 Refactored from rpa_ri_gpw [Frederick Stein]
12 : ! **************************************************************************************************
13 : MODULE time_frequency_grids
14 : USE greenx_interface, ONLY: greenx_get_minimax_grid
15 : USE kinds, ONLY: dp
16 : USE mathconstants, ONLY: pi,&
17 : z_zero
18 : USE minimax_exp, ONLY: get_exp_minimax_coeff
19 : USE minimax_exp_gw, ONLY: get_exp_minimax_coeff_gw
20 : USE minimax_rpa, ONLY: get_rpa_minimax_coeff,&
21 : get_rpa_minimax_coeff_larger_grid
22 : #include "./base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 :
28 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'time_frequency_grids'
29 :
30 : TYPE time_frequency_grid_type
31 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: frequency
32 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: frequency_weights
33 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: imaginary_time
34 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: time_weights_at_zero_frequency
35 : ! Complete discrete transforms, including cosine/sine factors.
36 : ! Matrix rows index output grid points, columns index input grid points.
37 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: cosine_time_to_frequency
38 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: cosine_frequency_to_time
39 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: sine_time_to_frequency
40 : END TYPE time_frequency_grid_type
41 :
42 : INTEGER, PARAMETER :: cosine_time_to_frequency = 1, &
43 : sine_time_to_frequency = 2, &
44 : cosine_frequency_to_time = 3
45 : INTEGER, PARAMETER, PUBLIC :: positive_time = 1, negative_time = 2
46 :
47 : PUBLIC :: build_clenshaw_grid, build_minimax_time_frequency_grid, test_least_square_ft, &
48 : time_frequency_grid_release, time_frequency_grid_type, transform_time_to_frequency
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Release all data owned by a time_frequency_grid_type object.
54 : !> \param grid The grid object to release.
55 : ! **************************************************************************************************
56 1130 : SUBROUTINE time_frequency_grid_release(grid)
57 :
58 : TYPE(time_frequency_grid_type), INTENT(INOUT) :: grid
59 :
60 1130 : IF (ALLOCATED(grid%frequency)) DEALLOCATE (grid%frequency)
61 1130 : IF (ALLOCATED(grid%frequency_weights)) DEALLOCATE (grid%frequency_weights)
62 1130 : IF (ALLOCATED(grid%imaginary_time)) DEALLOCATE (grid%imaginary_time)
63 1130 : IF (ALLOCATED(grid%time_weights_at_zero_frequency)) DEALLOCATE (grid%time_weights_at_zero_frequency)
64 1130 : IF (ALLOCATED(grid%cosine_time_to_frequency)) DEALLOCATE (grid%cosine_time_to_frequency)
65 1130 : IF (ALLOCATED(grid%cosine_frequency_to_time)) DEALLOCATE (grid%cosine_frequency_to_time)
66 1130 : IF (ALLOCATED(grid%sine_time_to_frequency)) DEALLOCATE (grid%sine_time_to_frequency)
67 :
68 1130 : END SUBROUTINE time_frequency_grid_release
69 :
70 : ! **************************************************************************************************
71 : !> \brief Transform real positive/negative-time samples to complex positive-frequency values.
72 : !> The even part uses the cosine transform and the odd part the sine transform.
73 : !> \param grid Valid grid with complete cosine/sine time-to-frequency matrices of matching dimensions.
74 : !> \param samples_tau Real samples indexed (state, time, branch), with branches
75 : !> positive_time and negative_time. The samples include their physical signs.
76 : !> \param values_omega Complex output indexed (state, frequency), evaluated at the leading
77 : !> SIZE(values_omega, 2) positive-frequency nodes of grid.
78 : ! **************************************************************************************************
79 680 : SUBROUTINE transform_time_to_frequency(grid, samples_tau, values_omega)
80 : TYPE(time_frequency_grid_type), INTENT(IN) :: grid
81 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: samples_tau
82 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: values_omega
83 :
84 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_time_to_frequency'
85 :
86 : INTEGER :: handle, ifreq, itime, nstate
87 : REAL(KIND=dp) :: weight_cos, weight_sin
88 680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: samples_even, samples_odd
89 :
90 680 : CALL timeset(routineN, handle)
91 :
92 680 : nstate = SIZE(samples_tau, 1)
93 680 : CPASSERT(SIZE(samples_tau, 3) == 2)
94 680 : CPASSERT(SIZE(values_omega, 1) == nstate)
95 680 : CPASSERT(SIZE(samples_tau, 2) == SIZE(grid%cosine_time_to_frequency, 2))
96 680 : CPASSERT(SIZE(values_omega, 2) <= SIZE(grid%cosine_time_to_frequency, 1))
97 :
98 2720 : ALLOCATE (samples_even(nstate), samples_odd(nstate))
99 32548 : values_omega = z_zero
100 7226 : DO itime = 1, SIZE(samples_tau, 2)
101 73796 : samples_even(:) = 0.5_dp*(samples_tau(:, itime, positive_time) + samples_tau(:, itime, negative_time))
102 73796 : samples_odd(:) = 0.5_dp*(samples_tau(:, itime, positive_time) - samples_tau(:, itime, negative_time))
103 39662 : DO ifreq = 1, SIZE(values_omega, 2)
104 32436 : weight_cos = grid%cosine_time_to_frequency(ifreq, itime)
105 32436 : weight_sin = grid%sine_time_to_frequency(ifreq, itime)
106 : values_omega(:, ifreq) = values_omega(:, ifreq) + &
107 423510 : CMPLX(weight_cos*samples_even, weight_sin*samples_odd, KIND=dp)
108 : END DO
109 : END DO
110 :
111 680 : CALL timestop(handle)
112 :
113 1360 : END SUBROUTINE transform_time_to_frequency
114 :
115 : ! **************************************************************************************************
116 : !> \brief Build a Clenshaw-Curtis frequency grid.
117 : !> \param num_points Number of integration points.
118 : !> \param grid Grid object to fill.
119 : ! **************************************************************************************************
120 128 : SUBROUTINE build_clenshaw_grid(num_points, grid)
121 :
122 : INTEGER, INTENT(IN) :: num_points
123 : TYPE(time_frequency_grid_type), INTENT(OUT) :: grid
124 :
125 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_clenshaw_grid'
126 :
127 : INTEGER :: handle, jquad
128 :
129 128 : CALL timeset(routineN, handle)
130 :
131 128 : CPASSERT(num_points > 0)
132 :
133 512 : ALLOCATE (grid%frequency(num_points), grid%frequency_weights(num_points))
134 :
135 5550 : DO jquad = 1, num_points - 1
136 5422 : grid%frequency(jquad) = jquad*pi/(2.0_dp*num_points)
137 5550 : grid%frequency_weights(jquad) = pi/(num_points*SIN(grid%frequency(jquad))**2)
138 : END DO
139 128 : grid%frequency(num_points) = pi/2.0_dp
140 : grid%frequency_weights(num_points) = &
141 128 : pi/(2.0_dp*num_points*SIN(grid%frequency(num_points))**2)
142 :
143 128 : CALL timestop(handle)
144 :
145 128 : END SUBROUTINE build_clenshaw_grid
146 :
147 : ! **************************************************************************************************
148 : !> \brief Build a minimax time/frequency grid through the common backend boundary.
149 : !> \param num_points Number of minimax points.
150 : !> \param energy_min Lower end of the physical energy interval.
151 : !> \param energy_max Upper end of the physical energy interval.
152 : !> \param regularization Regularization used for the fitted transform weights.
153 : !> \param num_points_per_magnitude Number of fitting points per decade.
154 : !> \param grid Grid object to fill.
155 : !> \param build_frequency Whether to construct the frequency components.
156 : !> \param build_time Whether to construct the imaginary-time components.
157 : !> \param build_transforms Whether to construct the complete time/frequency transform matrices.
158 : !> \param build_sine Whether to construct the sine time-to-frequency matrix.
159 : !> \param time_scaling Scaling applied to the imaginary-time abscissas.
160 : !> \param time_weight_scaling Scaling applied to the imaginary-time weights.
161 : !> \param max_fit_error Maximum fitting error across all requested transforms.
162 : !> \param print_warning Whether the minimax coefficient routine prints warnings.
163 : !> \param unit_nr Output unit used by an external backend.
164 : !> \param prefer_external_backend Whether to try an available external backend first.
165 : !> \param used_external_backend Whether the external backend supplied the grid.
166 : ! **************************************************************************************************
167 352 : SUBROUTINE build_minimax_time_frequency_grid(num_points, energy_min, energy_max, regularization, &
168 : num_points_per_magnitude, grid, build_frequency, build_time, &
169 : build_transforms, build_sine, time_scaling, time_weight_scaling, &
170 : max_fit_error, print_warning, unit_nr, prefer_external_backend, &
171 : used_external_backend)
172 :
173 : INTEGER, INTENT(IN) :: num_points
174 : REAL(KIND=dp), INTENT(IN) :: energy_min, energy_max, regularization
175 : INTEGER, INTENT(IN) :: num_points_per_magnitude
176 : TYPE(time_frequency_grid_type), INTENT(OUT) :: grid
177 : LOGICAL, INTENT(IN) :: build_frequency, build_time, &
178 : build_transforms, build_sine
179 : REAL(KIND=dp), INTENT(IN) :: time_scaling, time_weight_scaling
180 : REAL(KIND=dp), INTENT(OUT) :: max_fit_error
181 : LOGICAL, INTENT(IN) :: print_warning
182 : INTEGER, INTENT(IN) :: unit_nr
183 : LOGICAL, INTENT(IN) :: prefer_external_backend
184 : LOGICAL, INTENT(OUT), OPTIONAL :: used_external_backend
185 :
186 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_minimax_time_frequency_grid'
187 :
188 : INTEGER :: external_ierr, handle, ierr
189 : REAL(KIND=dp) :: e_range, max_error
190 352 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: points_and_weights
191 :
192 352 : CALL timeset(routineN, handle)
193 :
194 352 : IF (PRESENT(used_external_backend)) used_external_backend = .FALSE.
195 352 : IF (prefer_external_backend) THEN
196 : CALL greenx_get_minimax_grid(unit_nr, num_points, energy_min, energy_max, regularization, &
197 : grid%imaginary_time, grid%time_weights_at_zero_frequency, &
198 : grid%frequency, grid%frequency_weights, &
199 : grid%cosine_time_to_frequency, &
200 : grid%cosine_frequency_to_time, &
201 214 : grid%sine_time_to_frequency, external_ierr)
202 214 : IF (external_ierr == 0) THEN
203 78 : max_fit_error = 0.0_dp
204 78 : IF (PRESENT(used_external_backend)) used_external_backend = .TRUE.
205 78 : CALL timestop(handle)
206 : RETURN
207 : END IF
208 136 : CALL time_frequency_grid_release(grid)
209 : END IF
210 :
211 274 : CPASSERT(num_points > 0)
212 274 : CPASSERT(energy_min > 0.0_dp)
213 274 : CPASSERT(energy_max >= energy_min)
214 274 : CPASSERT(time_scaling > 0.0_dp)
215 274 : CPASSERT(time_weight_scaling > 0.0_dp)
216 274 : CPASSERT(.NOT. build_transforms .OR. (build_frequency .AND. build_time))
217 274 : CPASSERT(.NOT. build_sine .OR. build_transforms)
218 274 : max_error = 0.0_dp
219 :
220 274 : e_range = energy_max/energy_min
221 822 : ALLOCATE (points_and_weights(2*num_points))
222 :
223 274 : IF (build_frequency) THEN
224 210 : IF (num_points <= 20) THEN
225 210 : CALL get_rpa_minimax_coeff(num_points, e_range, points_and_weights, ierr, print_warning)
226 : ELSE
227 0 : CALL get_rpa_minimax_coeff_larger_grid(num_points, e_range, points_and_weights)
228 : END IF
229 :
230 630 : ALLOCATE (grid%frequency(num_points))
231 420 : ALLOCATE (grid%frequency_weights(num_points))
232 2554 : grid%frequency(:) = points_and_weights(1:num_points)*energy_min
233 2554 : grid%frequency_weights(:) = points_and_weights(num_points + 1:)*energy_min
234 210 : IF (num_points >= 26) grid%frequency_weights(:) = grid%frequency_weights(:)*4.0_dp
235 : END IF
236 :
237 274 : IF (build_time) THEN
238 244 : IF (num_points <= 20) THEN
239 244 : CALL get_exp_minimax_coeff(num_points, e_range, points_and_weights)
240 : ELSE
241 0 : CALL get_exp_minimax_coeff_gw(num_points, e_range, points_and_weights)
242 : END IF
243 :
244 732 : ALLOCATE (grid%imaginary_time(num_points))
245 488 : ALLOCATE (grid%time_weights_at_zero_frequency(num_points))
246 2712 : grid%imaginary_time(:) = points_and_weights(1:num_points)/time_scaling/energy_min
247 2712 : grid%time_weights_at_zero_frequency(:) = points_and_weights(num_points + 1:)/time_weight_scaling/energy_min
248 :
249 244 : IF (build_transforms) THEN
250 720 : ALLOCATE (grid%cosine_time_to_frequency(num_points, num_points))
251 : CALL fit_l_sq_weights(cosine_time_to_frequency, num_points, grid%imaginary_time, &
252 : grid%cosine_time_to_frequency, grid%frequency, energy_min, energy_max, &
253 180 : max_error, num_points_per_magnitude, regularization)
254 :
255 540 : ALLOCATE (grid%cosine_frequency_to_time(num_points, num_points))
256 : CALL fit_l_sq_weights(cosine_frequency_to_time, num_points, grid%imaginary_time, &
257 : grid%cosine_frequency_to_time, grid%frequency, energy_min, energy_max, &
258 180 : max_error, num_points_per_magnitude, regularization)
259 :
260 180 : IF (build_sine) THEN
261 444 : ALLOCATE (grid%sine_time_to_frequency(num_points, num_points))
262 : CALL fit_l_sq_weights(sine_time_to_frequency, num_points, grid%imaginary_time, &
263 : grid%sine_time_to_frequency, grid%frequency, energy_min, energy_max, &
264 148 : max_error, num_points_per_magnitude, regularization)
265 : END IF
266 : END IF
267 : END IF
268 :
269 274 : max_fit_error = max_error
270 :
271 274 : CALL timestop(handle)
272 :
273 626 : END SUBROUTINE build_minimax_time_frequency_grid
274 :
275 : ! **************************************************************************************************
276 : !> \brief Calculate a complete least-squares time/frequency transform matrix.
277 : !> \param transform_kind Type of transform to fit.
278 : !> \param num_integ_points Number of integration points.
279 : !> \param tau_tj Imaginary-time integration points.
280 : !> \param weights Complete transform coefficients, including cosine/sine factors.
281 : !> Rows index output grid points, columns index input grid points.
282 : !> \param omega_tj Frequency integration points.
283 : !> \param E_min Lower end of the fitting interval.
284 : !> \param E_max Upper end of the fitting interval.
285 : !> \param max_error Maximum fitting error.
286 : !> \param num_points_per_magnitude Number of fitting points per decade.
287 : !> \param regularization Regularization parameter for the pseudoinverse.
288 : ! **************************************************************************************************
289 508 : SUBROUTINE fit_l_sq_weights(transform_kind, num_integ_points, tau_tj, weights, omega_tj, E_min, E_max, &
290 : max_error, num_points_per_magnitude, regularization)
291 :
292 : INTEGER, INTENT(IN) :: transform_kind, num_integ_points
293 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tau_tj
294 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: weights
295 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: omega_tj
296 : REAL(KIND=dp), INTENT(IN) :: E_min, E_max
297 : REAL(KIND=dp), INTENT(OUT) :: max_error
298 : INTEGER, INTENT(IN) :: num_points_per_magnitude
299 : REAL(KIND=dp), INTENT(IN) :: regularization
300 :
301 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_l_sq_weights'
302 :
303 : INTEGER :: handle, iii, info, jjj, jquad, lwork, &
304 : num_x_nodes
305 508 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
306 : LOGICAL :: sine_transform, time_to_frequency
307 : REAL(KIND=dp) :: func_val, max_error_tmp, multiplicator, &
308 : omega, tau, x_value
309 508 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, vec_UTy, weight_work, work, &
310 508 : x_values, y_values
311 508 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, &
312 508 : mat_SinvVSinvT, mat_U
313 :
314 508 : CALL timeset(routineN, handle)
315 :
316 508 : CPASSERT(transform_kind >= cosine_time_to_frequency)
317 508 : CPASSERT(transform_kind <= cosine_frequency_to_time)
318 508 : CPASSERT(SIZE(tau_tj) == num_integ_points)
319 508 : CPASSERT(SIZE(omega_tj) == num_integ_points)
320 508 : CPASSERT(SIZE(weights, 1) == num_integ_points)
321 508 : CPASSERT(SIZE(weights, 2) == num_integ_points)
322 :
323 508 : time_to_frequency = transform_kind /= cosine_frequency_to_time
324 508 : sine_transform = transform_kind == sine_time_to_frequency
325 :
326 : ! take num_points_per_magnitude points per 10-interval
327 508 : num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude
328 :
329 : ! take at least as many x points as integration points to have clear
330 : ! input for the singular value decomposition
331 508 : num_x_nodes = MAX(num_x_nodes, num_integ_points)
332 :
333 : ALLOCATE (x_values(num_x_nodes), y_values(num_x_nodes), mat_A(num_x_nodes, num_integ_points), &
334 : weight_work(num_integ_points), sing_values(num_integ_points), mat_U(num_x_nodes, num_x_nodes), &
335 : mat_SinvVSinvT(num_x_nodes, num_integ_points), work(8*num_integ_points*num_integ_points + &
336 : 12*num_integ_points + 2*num_x_nodes), iwork(8*num_integ_points), &
337 11176 : mat_SinvVSinvSigma(num_integ_points, num_x_nodes), vec_UTy(num_x_nodes))
338 508 : mat_SinvVSinvSigma = 0.0_dp
339 :
340 : ! double the value nessary for 'A' to achieve good performance
341 508 : lwork = SIZE(work)
342 :
343 508 : multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
344 157308 : DO iii = 1, num_x_nodes
345 157308 : x_values(iii) = E_min*multiplicator**(iii - 1)
346 : END DO
347 :
348 508 : max_error = 0.0_dp
349 :
350 7044 : DO jquad = 1, num_integ_points
351 :
352 6536 : IF (time_to_frequency) THEN
353 4328 : omega = omega_tj(jquad)
354 :
355 4328 : IF (sine_transform) THEN
356 : ! y=2*omega/(x^2+omega^2)
357 562120 : DO iii = 1, num_x_nodes
358 562120 : y_values(iii) = 2.0_dp*omega/(x_values(iii)**2 + omega**2)
359 : END DO
360 : ELSE
361 : ! y=2*x/(x^2+omega^2)
362 593408 : DO iii = 1, num_x_nodes
363 593408 : y_values(iii) = 2.0_dp*x_values(iii)/(x_values(iii)**2 + omega**2)
364 : END DO
365 : END IF
366 :
367 76824 : DO jjj = 1, num_integ_points
368 17544824 : DO iii = 1, num_x_nodes
369 17540496 : IF (sine_transform) THEN
370 8689600 : mat_A(iii, jjj) = SIN(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
371 : ELSE
372 8778400 : mat_A(iii, jjj) = COS(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
373 : END IF
374 : END DO
375 : END DO
376 : ELSE
377 2208 : tau = tau_tj(jquad)
378 :
379 : ! y=exp(-x*|tau|)
380 593408 : DO iii = 1, num_x_nodes
381 593408 : y_values(iii) = EXP(-x_values(iii)*tau)
382 : END DO
383 :
384 38580 : DO jjj = 1, num_integ_points
385 36372 : omega = omega_tj(jjj)
386 8816980 : DO iii = 1, num_x_nodes
387 8778400 : x_value = x_values(iii)
388 8814772 : mat_A(iii, jjj) = COS(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2)
389 : END DO
390 : END DO
391 : END IF
392 :
393 : ! Singular value decomposition of mat_A
394 : CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
395 6536 : mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
396 6536 : CPASSERT(info == 0)
397 :
398 : ! integration weights = V Sigma U^T y
399 115404 : DO jjj = 1, num_integ_points
400 2129468 : DO iii = 1, num_integ_points
401 : mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)*sing_values(jjj) &
402 2122932 : /(regularization**2 + sing_values(jjj)**2)
403 : END DO
404 : END DO
405 :
406 : CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, &
407 6536 : 0.0_dp, vec_UTy, num_x_nodes)
408 : CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, &
409 6536 : vec_UTy, num_x_nodes, 0.0_dp, weight_work, num_integ_points)
410 :
411 6536 : IF (sine_transform) THEN
412 38244 : weights(jquad, :) = weight_work(:)*SIN(omega*tau_tj(:))
413 4416 : ELSE IF (time_to_frequency) THEN
414 38580 : weights(jquad, :) = weight_work(:)*COS(omega*tau_tj(:))
415 : ELSE
416 38580 : weights(jquad, :) = weight_work(:)*COS(tau*omega_tj(:))
417 : END IF
418 :
419 6536 : max_error_tmp = 0.0_dp
420 1748936 : DO iii = 1, num_x_nodes
421 1742400 : func_val = 0.0_dp
422 1742400 : IF (time_to_frequency) THEN
423 18619200 : DO jjj = 1, num_integ_points
424 18619200 : func_val = func_val + weights(jquad, jjj)*EXP(-x_values(iii)*tau_tj(jjj))
425 : END DO
426 : ELSE
427 9369600 : DO jjj = 1, num_integ_points
428 8778400 : omega = omega_tj(jjj)
429 : func_val = func_val + weights(jquad, jjj)*2.0_dp*x_values(iii) &
430 9369600 : /(x_values(iii)**2 + omega**2)
431 : END DO
432 : END IF
433 1748936 : max_error_tmp = MAX(max_error_tmp, ABS(y_values(iii) - func_val))
434 : END DO
435 7044 : max_error = MAX(max_error, max_error_tmp)
436 :
437 : END DO
438 :
439 508 : CALL timestop(handle)
440 :
441 1016 : END SUBROUTINE fit_l_sq_weights
442 :
443 : ! **************************************************************************************************
444 : !> \brief test the singular value decomposition for the computation of integration weights for the
445 : !> Fourier transform between time and frequency grid in cubic-scaling RPA
446 : !> \param nR ...
447 : !> \param iw ...
448 : ! **************************************************************************************************
449 0 : SUBROUTINE test_least_square_ft(nR, iw)
450 : INTEGER, INTENT(IN) :: nR, iw
451 :
452 : REAL(KIND=dp), PARAMETER :: Rc_max = 1.0E+07_dp
453 :
454 : INTEGER :: ierr, iR, num_integ_points
455 : REAL(KIND=dp) :: max_error, multiplicator, Rc
456 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tau_tj, tj, x_tw
457 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: weights_cos_tf_t_to_w
458 :
459 0 : multiplicator = Rc_max**(1.0_dp/(REAL(nR, KIND=dp) - 1.0_dp))
460 :
461 0 : DO num_integ_points = 1, 20
462 :
463 0 : ALLOCATE (x_tw(2*num_integ_points))
464 0 : x_tw = 0.0_dp
465 0 : ALLOCATE (tau_tj(num_integ_points))
466 0 : tau_tj = 0.0_dp
467 0 : ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
468 0 : weights_cos_tf_t_to_w = 0.0_dp
469 0 : ALLOCATE (tj(num_integ_points))
470 0 : tj = 0.0_dp
471 :
472 0 : DO iR = 0, nR - 1
473 :
474 0 : Rc = 2.0_dp*multiplicator**iR
475 :
476 0 : ierr = 0
477 0 : CALL get_rpa_minimax_coeff(num_integ_points, Rc, x_tw, ierr, print_warning=.FALSE.)
478 :
479 0 : tj(:) = x_tw(1:num_integ_points)
480 :
481 0 : x_tw = 0.0_dp
482 :
483 0 : CALL get_exp_minimax_coeff(num_integ_points, Rc, x_tw)
484 :
485 0 : tau_tj(:) = x_tw(1:num_integ_points)/2.0_dp
486 :
487 : CALL fit_l_sq_weights(cosine_time_to_frequency, num_integ_points, tau_tj, &
488 0 : weights_cos_tf_t_to_w, tj, 1.0_dp, Rc, max_error, 200, 0.0_dp)
489 :
490 0 : IF (iw > 0) THEN
491 0 : WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, Rc, max_error
492 : END IF
493 :
494 : END DO
495 :
496 0 : DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tj)
497 :
498 : END DO
499 :
500 0 : END SUBROUTINE test_least_square_ft
501 :
502 0 : END MODULE time_frequency_grids
|
