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 Routines to calculate frequency and time grids (integration points and weights)
10 : !> for correlation methods
11 : !> \par History
12 : !> 05.2019 Refactored from rpa_ri_gpw [Frederick Stein]
13 : ! **************************************************************************************************
14 : MODULE mp2_grids
15 : USE cp_fm_types, ONLY: cp_fm_get_info,&
16 : cp_fm_type
17 : USE greenx_interface, ONLY: greenx_get_minimax_grid
18 : USE input_section_types, ONLY: section_vals_type,&
19 : section_vals_val_set
20 : USE kinds, ONLY: dp
21 : USE kpoint_types, ONLY: get_kpoint_info,&
22 : kpoint_env_type,&
23 : kpoint_type
24 : USE machine, ONLY: m_flush
25 : USE mathconstants, ONLY: pi
26 : USE message_passing, ONLY: mp_para_env_release,&
27 : mp_para_env_type
28 : USE minimax_exp, ONLY: get_exp_minimax_coeff
29 : USE minimax_exp_gw, ONLY: get_exp_minimax_coeff_gw
30 : USE minimax_rpa, ONLY: get_rpa_minimax_coeff,&
31 : get_rpa_minimax_coeff_larger_grid
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_environment_type
34 : USE qs_mo_types, ONLY: get_mo_set,&
35 : mo_set_type
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 :
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_grids'
43 :
44 : PUBLIC :: get_minimax_grid, get_clenshaw_grid, test_least_square_ft, get_l_sq_wghts_cos_tf_t_to_w, &
45 : get_l_sq_wghts_cos_tf_w_to_t, get_l_sq_wghts_sin_tf_t_to_w
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief ...
51 : !> \param para_env ...
52 : !> \param unit_nr ...
53 : !> \param homo ...
54 : !> \param Eigenval ...
55 : !> \param num_integ_points ...
56 : !> \param do_im_time ...
57 : !> \param do_ri_sos_laplace_mp2 ...
58 : !> \param do_print ...
59 : !> \param tau_tj ...
60 : !> \param tau_wj ...
61 : !> \param qs_env ...
62 : !> \param do_gw_im_time ...
63 : !> \param do_kpoints_cubic_RPA ...
64 : !> \param e_fermi ...
65 : !> \param tj ...
66 : !> \param wj ...
67 : !> \param weights_cos_tf_t_to_w ...
68 : !> \param weights_cos_tf_w_to_t ...
69 : !> \param weights_sin_tf_t_to_w ...
70 : !> \param regularization ...
71 : ! **************************************************************************************************
72 204 : SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, &
73 : do_im_time, do_ri_sos_laplace_mp2, do_print, tau_tj, tau_wj, qs_env, do_gw_im_time, &
74 : do_kpoints_cubic_RPA, e_fermi, tj, wj, weights_cos_tf_t_to_w, &
75 : weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, regularization)
76 :
77 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
78 : INTEGER, INTENT(IN) :: unit_nr
79 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
80 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval
81 : INTEGER, INTENT(IN) :: num_integ_points
82 : LOGICAL, INTENT(IN) :: do_im_time, do_ri_sos_laplace_mp2, &
83 : do_print
84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
85 : INTENT(OUT) :: tau_tj, tau_wj
86 : TYPE(qs_environment_type), POINTER :: qs_env
87 : LOGICAL, INTENT(IN) :: do_gw_im_time, do_kpoints_cubic_RPA
88 : REAL(KIND=dp), INTENT(OUT) :: e_fermi
89 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
90 : INTENT(OUT) :: tj, wj
91 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
92 : INTENT(OUT) :: weights_cos_tf_t_to_w, &
93 : weights_cos_tf_w_to_t, &
94 : weights_sin_tf_t_to_w
95 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: regularization
96 :
97 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_minimax_grid'
98 : INTEGER, PARAMETER :: num_points_per_magnitude = 200
99 :
100 : INTEGER :: handle, ierr, jquad, nspins
101 : LOGICAL :: my_do_kpoints
102 : REAL(KIND=dp) :: E_Range, Emax, Emin, max_error_min, &
103 : my_regularization, scaling
104 204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x_tw
105 :
106 204 : CALL timeset(routineN, handle)
107 :
108 : CALL determine_energy_range(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, &
109 204 : do_kpoints_cubic_RPA, Emin, Emax, e_range, e_fermi)
110 :
111 : CALL greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, &
112 : tau_tj, tau_wj, qs_env%mp2_env%ri_g0w0%regularization_minimax, &
113 : tj, wj, weights_cos_tf_t_to_w, &
114 204 : weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, ierr)
115 :
116 : ! Shortcut if Greenx was available and successful
117 204 : IF (ierr == 0) THEN
118 74 : CALL timestop(handle)
119 : RETURN
120 : END IF
121 :
122 : ! Test for spin unrestricted
123 130 : nspins = SIZE(homo)
124 :
125 : ! Test whether all necessary variables are available
126 130 : my_do_kpoints = .FALSE.
127 130 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
128 130 : my_do_kpoints = do_kpoints_cubic_RPA
129 : END IF
130 :
131 : my_regularization = 0.0_dp
132 130 : IF (PRESENT(regularization)) THEN
133 130 : my_regularization = regularization
134 :
135 130 : IF (num_integ_points > 20 .AND. e_range < 100.0_dp) THEN
136 0 : IF (unit_nr > 0) &
137 : CALL cp_warn(__LOCATION__, &
138 : "You requested a large minimax grid (> 20 points) for a small minimax range R (R < 100). "// &
139 : "That may lead to numerical "// &
140 : "instabilities when computing minimax grid weights. You can prevent small ranges by choosing "// &
141 0 : "a larger basis set with higher angular momenta or alternatively using all-electron calculations.")
142 : END IF
143 :
144 130 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
145 216 : ALLOCATE (x_tw(2*num_integ_points))
146 620 : x_tw = 0.0_dp
147 72 : ierr = 0
148 72 : IF (num_integ_points .LE. 20) THEN
149 72 : CALL get_rpa_minimax_coeff(num_integ_points, e_range, x_tw, ierr)
150 : ELSE
151 0 : CALL get_rpa_minimax_coeff_larger_grid(num_integ_points, e_range, x_tw)
152 : END IF
153 :
154 216 : ALLOCATE (tj(num_integ_points))
155 346 : tj = 0.0_dp
156 :
157 144 : ALLOCATE (wj(num_integ_points))
158 346 : wj = 0.0_dp
159 :
160 346 : DO jquad = 1, num_integ_points
161 274 : tj(jquad) = x_tw(jquad)
162 346 : wj(jquad) = x_tw(jquad + num_integ_points)
163 : END DO
164 :
165 : ! for the smaller grids, the factor of 4 is included in get_rpa_minimax_coeff for wj
166 72 : IF (num_integ_points >= 26) THEN
167 0 : wj(:) = wj(:)*4.0_dp
168 : END IF
169 :
170 72 : DEALLOCATE (x_tw)
171 :
172 72 : IF (unit_nr > 0 .AND. do_print) THEN
173 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
174 35 : "MINIMAX_INFO| Number of integration points:", num_integ_points
175 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
176 35 : "MINIMAX_INFO| Gap for the minimax approximation:", Emin
177 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
178 35 : "MINIMAX_INFO| Range for the minimax approximation:", e_range
179 35 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") "MINIMAX_INFO| Minimax parameters:", "Weights", "Abscissas"
180 167 : DO jquad = 1, num_integ_points
181 167 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(jquad), tj(jquad)
182 : END DO
183 35 : CALL m_flush(unit_nr)
184 : END IF
185 :
186 : ! scale the minimax parameters
187 346 : tj(:) = tj(:)*Emin
188 346 : wj(:) = wj(:)*Emin
189 : END IF
190 :
191 : ! set up the minimax time grid
192 130 : IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
193 :
194 300 : ALLOCATE (x_tw(2*num_integ_points))
195 836 : x_tw = 0.0_dp
196 :
197 100 : IF (num_integ_points .LE. 20) THEN
198 100 : CALL get_exp_minimax_coeff(num_integ_points, e_range, x_tw)
199 : ELSE
200 0 : CALL get_exp_minimax_coeff_gw(num_integ_points, e_range, x_tw)
201 : END IF
202 :
203 : ! For RPA we include already a factor of two (see later steps)
204 100 : scaling = 2.0_dp
205 100 : IF (do_ri_sos_laplace_mp2) scaling = 1.0_dp
206 :
207 300 : ALLOCATE (tau_tj(num_integ_points))
208 468 : tau_tj = 0.0_dp
209 :
210 200 : ALLOCATE (tau_wj(num_integ_points))
211 468 : tau_wj = 0.0_dp
212 :
213 468 : DO jquad = 1, num_integ_points
214 368 : tau_tj(jquad) = x_tw(jquad)/scaling
215 468 : tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling
216 : END DO
217 :
218 100 : DEALLOCATE (x_tw)
219 :
220 100 : IF (unit_nr > 0 .AND. do_print) THEN
221 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
222 49 : "MINIMAX_INFO| Range for the minimax approximation:", e_range
223 : ! For testing the gap
224 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
225 49 : "MINIMAX_INFO| Gap:", Emin
226 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
227 49 : "MINIMAX_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas"
228 228 : DO jquad = 1, num_integ_points
229 228 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad)
230 : END DO
231 49 : CALL m_flush(unit_nr)
232 : END IF
233 :
234 : ! scale grid from [1,R] to [Emin,Emax]
235 468 : tau_tj(:) = tau_tj(:)/Emin
236 468 : tau_wj(:) = tau_wj(:)/Emin
237 :
238 100 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
239 168 : ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
240 678 : weights_cos_tf_t_to_w = 0.0_dp
241 :
242 : CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, tj, &
243 : Emin, Emax, max_error_min, num_points_per_magnitude, &
244 42 : my_regularization)
245 :
246 : ! get the weights for the cosine transform W^c(iw) -> W^c(it)
247 126 : ALLOCATE (weights_cos_tf_w_to_t(num_integ_points, num_integ_points))
248 678 : weights_cos_tf_w_to_t = 0.0_dp
249 :
250 : CALL get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, tj, &
251 : Emin, Emax, max_error_min, num_points_per_magnitude, &
252 42 : my_regularization)
253 :
254 42 : IF (do_gw_im_time) THEN
255 :
256 : ! get the weights for the sine transform Sigma^sin(it) -> Sigma^sin(iw) (PRB 94, 165109 (2016), Eq. 71)
257 30 : ALLOCATE (weights_sin_tf_t_to_w(num_integ_points, num_integ_points))
258 310 : weights_sin_tf_t_to_w = 0.0_dp
259 :
260 : CALL get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, tj, &
261 : Emin, Emax, max_error_min, num_points_per_magnitude, &
262 10 : my_regularization)
263 :
264 10 : IF (unit_nr > 0) THEN
265 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
266 5 : "MINIMAX_INFO| Maximum deviation of the imag. time fit:", max_error_min
267 : END IF
268 : END IF
269 :
270 : END IF
271 :
272 : END IF
273 : END IF
274 :
275 130 : CALL timestop(handle)
276 :
277 130 : END SUBROUTINE get_minimax_grid
278 :
279 : ! **************************************************************************************************
280 : !> \brief ...
281 : !> \param para_env ...
282 : !> \param para_env_RPA ...
283 : !> \param unit_nr ...
284 : !> \param homo ...
285 : !> \param virtual ...
286 : !> \param Eigenval ...
287 : !> \param num_integ_points ...
288 : !> \param num_integ_group ...
289 : !> \param color_rpa_group ...
290 : !> \param fm_mat_S ...
291 : !> \param my_do_gw ...
292 : !> \param ext_scaling ...
293 : !> \param a_scaling ...
294 : !> \param tj ...
295 : !> \param wj ...
296 : ! **************************************************************************************************
297 96 : SUBROUTINE get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, &
298 96 : num_integ_group, color_rpa_group, fm_mat_S, my_do_gw, &
299 : ext_scaling, a_scaling, tj, wj)
300 :
301 : TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_RPA
302 : INTEGER, INTENT(IN) :: unit_nr
303 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
304 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval
305 : INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, &
306 : color_rpa_group
307 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
308 : LOGICAL, INTENT(IN) :: my_do_gw
309 : REAL(KIND=dp), INTENT(IN) :: ext_scaling
310 : REAL(KIND=dp), INTENT(OUT) :: a_scaling
311 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
312 : INTENT(OUT) :: tj, wj
313 :
314 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_clenshaw_grid'
315 :
316 : INTEGER :: handle, jquad, nspins
317 : LOGICAL :: my_open_shell
318 :
319 96 : CALL timeset(routineN, handle)
320 :
321 96 : nspins = SIZE(homo)
322 96 : my_open_shell = (nspins == 2)
323 :
324 : ! Now, start to prepare the different grid
325 288 : ALLOCATE (tj(num_integ_points))
326 3766 : tj = 0.0_dp
327 :
328 192 : ALLOCATE (wj(num_integ_points))
329 3766 : wj = 0.0_dp
330 :
331 3670 : DO jquad = 1, num_integ_points - 1
332 3574 : tj(jquad) = jquad*pi/(2.0_dp*num_integ_points)
333 3670 : wj(jquad) = pi/(num_integ_points*SIN(tj(jquad))**2)
334 : END DO
335 96 : tj(num_integ_points) = pi/2.0_dp
336 96 : wj(num_integ_points) = pi/(2.0_dp*num_integ_points*SIN(tj(num_integ_points))**2)
337 :
338 96 : IF (my_do_gw .AND. ext_scaling > 0.0_dp) THEN
339 58 : a_scaling = ext_scaling
340 : ELSE
341 : CALL calc_scaling_factor(a_scaling, para_env, para_env_RPA, homo, virtual, Eigenval, &
342 : num_integ_points, num_integ_group, color_rpa_group, &
343 38 : tj, wj, fm_mat_S)
344 : END IF
345 :
346 96 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.5)') 'INTEG_INFO| Scaling parameter:', a_scaling
347 :
348 3766 : wj(:) = wj(:)*a_scaling
349 :
350 96 : CALL timestop(handle)
351 :
352 96 : END SUBROUTINE get_clenshaw_grid
353 :
354 : ! **************************************************************************************************
355 : !> \brief ...
356 : !> \param a_scaling_ext ...
357 : !> \param para_env ...
358 : !> \param para_env_RPA ...
359 : !> \param homo ...
360 : !> \param virtual ...
361 : !> \param Eigenval ...
362 : !> \param num_integ_points ...
363 : !> \param num_integ_group ...
364 : !> \param color_rpa_group ...
365 : !> \param tj_ext ...
366 : !> \param wj_ext ...
367 : !> \param fm_mat_S ...
368 : ! **************************************************************************************************
369 38 : SUBROUTINE calc_scaling_factor(a_scaling_ext, para_env, para_env_RPA, homo, virtual, Eigenval, &
370 : num_integ_points, num_integ_group, color_rpa_group, &
371 38 : tj_ext, wj_ext, fm_mat_S)
372 : REAL(KIND=dp), INTENT(OUT) :: a_scaling_ext
373 : TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_RPA
374 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
375 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval
376 : INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, &
377 : color_rpa_group
378 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
379 : INTENT(IN) :: tj_ext, wj_ext
380 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
381 :
382 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_scaling_factor'
383 :
384 : INTEGER :: handle, icycle, jquad, ncol_local, &
385 : ncol_local_beta, nspins
386 : LOGICAL :: my_open_shell
387 : REAL(KIND=dp) :: a_high, a_low, a_scaling, conv_param, eps, first_deriv, left_term, &
388 : right_term, right_term_ref, right_term_ref_beta, step
389 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cottj, D_ia, D_ia_beta, iaia_RI, &
390 38 : iaia_RI_beta, M_ia, M_ia_beta
391 : TYPE(mp_para_env_type), POINTER :: para_env_col, para_env_col_beta
392 :
393 38 : CALL timeset(routineN, handle)
394 :
395 38 : nspins = SIZE(homo)
396 38 : my_open_shell = (nspins == 2)
397 :
398 38 : eps = 1.0E-10_dp
399 :
400 114 : ALLOCATE (cottj(num_integ_points))
401 :
402 : ! calculate the cotangent of the abscissa tj
403 488 : DO jquad = 1, num_integ_points
404 488 : cottj(jquad) = 1.0_dp/TAN(tj_ext(jquad))
405 : END DO
406 :
407 : CALL calc_ia_ia_integrals(para_env_RPA, homo(1), virtual(1), ncol_local, right_term_ref, Eigenval(:, 1, 1), &
408 38 : D_ia, iaia_RI, M_ia, fm_mat_S(1), para_env_col)
409 :
410 : ! In the open shell case do point 1-2-3 for the beta spin
411 38 : IF (my_open_shell) THEN
412 : CALL calc_ia_ia_integrals(para_env_RPA, homo(2), virtual(2), ncol_local_beta, right_term_ref_beta, Eigenval(:, 1, 2), &
413 8 : D_ia_beta, iaia_RI_beta, M_ia_beta, fm_mat_S(2), para_env_col_beta)
414 :
415 8 : right_term_ref = right_term_ref + right_term_ref_beta
416 : END IF
417 :
418 : ! bcast the result
419 38 : IF (para_env%mepos == 0) THEN
420 19 : CALL para_env%bcast(right_term_ref, 0)
421 : ELSE
422 19 : right_term_ref = 0.0_dp
423 19 : CALL para_env%bcast(right_term_ref, 0)
424 : END IF
425 :
426 : ! 5) start iteration for solving the non-linear equation by bisection
427 : ! find limit, here step=0.5 seems a good compromise
428 38 : conv_param = 100.0_dp*EPSILON(right_term_ref)
429 38 : step = 0.5_dp
430 38 : a_low = 0.0_dp
431 38 : a_high = step
432 38 : right_term = -right_term_ref
433 100 : DO icycle = 1, num_integ_points*2
434 92 : a_scaling = a_high
435 :
436 : CALL calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
437 : M_ia, cottj, wj_ext, D_ia, D_ia_beta, M_ia_beta, &
438 : ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
439 92 : para_env, para_env_col, para_env_col_beta)
440 92 : left_term = left_term/4.0_dp/pi*a_scaling
441 :
442 92 : IF (ABS(left_term) > ABS(right_term) .OR. ABS(left_term + right_term) <= conv_param) EXIT
443 62 : a_low = a_high
444 100 : a_high = a_high + step
445 :
446 : END DO
447 :
448 38 : IF (ABS(left_term + right_term) >= conv_param) THEN
449 32 : IF (a_scaling >= 2*num_integ_points*step) THEN
450 10 : a_scaling = 1.0_dp
451 : ELSE
452 :
453 340 : DO icycle = 1, num_integ_points*2
454 336 : a_scaling = (a_low + a_high)/2.0_dp
455 :
456 : CALL calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
457 : M_ia, cottj, wj_ext, D_ia, D_ia_beta, M_ia_beta, &
458 : ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
459 336 : para_env, para_env_col, para_env_col_beta)
460 336 : left_term = left_term/4.0_dp/pi*a_scaling
461 :
462 336 : IF (ABS(left_term) > ABS(right_term)) THEN
463 : a_high = a_scaling
464 : ELSE
465 186 : a_low = a_scaling
466 : END IF
467 :
468 340 : IF (ABS(a_high - a_low) < 1.0e-5_dp) EXIT
469 :
470 : END DO
471 :
472 : END IF
473 : END IF
474 :
475 38 : a_scaling_ext = a_scaling
476 38 : CALL para_env%bcast(a_scaling_ext, 0)
477 :
478 38 : DEALLOCATE (cottj)
479 38 : DEALLOCATE (iaia_RI)
480 38 : DEALLOCATE (D_ia)
481 38 : DEALLOCATE (M_ia)
482 38 : CALL mp_para_env_release(para_env_col)
483 :
484 38 : IF (my_open_shell) THEN
485 8 : DEALLOCATE (iaia_RI_beta)
486 8 : DEALLOCATE (D_ia_beta)
487 8 : DEALLOCATE (M_ia_beta)
488 8 : CALL mp_para_env_release(para_env_col_beta)
489 : END IF
490 :
491 38 : CALL timestop(handle)
492 :
493 76 : END SUBROUTINE calc_scaling_factor
494 :
495 : ! **************************************************************************************************
496 : !> \brief ...
497 : !> \param para_env_RPA ...
498 : !> \param homo ...
499 : !> \param virtual ...
500 : !> \param ncol_local ...
501 : !> \param right_term_ref ...
502 : !> \param Eigenval ...
503 : !> \param D_ia ...
504 : !> \param iaia_RI ...
505 : !> \param M_ia ...
506 : !> \param fm_mat_S ...
507 : !> \param para_env_col ...
508 : ! **************************************************************************************************
509 46 : SUBROUTINE calc_ia_ia_integrals(para_env_RPA, homo, virtual, ncol_local, right_term_ref, Eigenval, &
510 : D_ia, iaia_RI, M_ia, fm_mat_S, para_env_col)
511 :
512 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA
513 : INTEGER, INTENT(IN) :: homo, virtual
514 : INTEGER, INTENT(OUT) :: ncol_local
515 : REAL(KIND=dp), INTENT(OUT) :: right_term_ref
516 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
517 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
518 : INTENT(OUT) :: D_ia, iaia_RI, M_ia
519 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
520 : TYPE(mp_para_env_type), POINTER :: para_env_col
521 :
522 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_ia_ia_integrals'
523 :
524 : INTEGER :: avirt, color_col, color_row, handle, &
525 : i_global, iiB, iocc, nrow_local
526 46 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
527 : REAL(KIND=dp) :: eigen_diff
528 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: iaia_RI_dp
529 : TYPE(mp_para_env_type), POINTER :: para_env_row
530 :
531 46 : CALL timeset(routineN, handle)
532 :
533 : ! calculate the (ia|ia) RI integrals
534 : ! ----------------------------------
535 : ! 1) get info fm_mat_S
536 : CALL cp_fm_get_info(matrix=fm_mat_S, &
537 : nrow_local=nrow_local, &
538 : ncol_local=ncol_local, &
539 : row_indices=row_indices, &
540 46 : col_indices=col_indices)
541 :
542 : ! allocate the local buffer of iaia_RI integrals (dp kind)
543 136 : ALLOCATE (iaia_RI_dp(ncol_local))
544 2764 : iaia_RI_dp = 0.0_dp
545 :
546 : ! 2) perform the local multiplication SUM_K (ia|K)*(ia|K)
547 2764 : DO iiB = 1, ncol_local
548 188006 : iaia_RI_dp(iiB) = iaia_RI_dp(iiB) + DOT_PRODUCT(fm_mat_S%local_data(:, iiB), fm_mat_S%local_data(:, iiB))
549 : END DO
550 :
551 : ! 3) sum the result with the processes of the RPA_group having the same columns
552 : ! _______ia______ _
553 : ! | | | | | | |
554 : ! --> | 1 | 5 | 9 | 13| SUM --> | |
555 : ! |___|__ |___|___| |_|
556 : ! | | | | | | |
557 : ! --> | 2 | 6 | 10| 14| SUM --> | |
558 : ! K |___|___|___|___| |_| (ia|ia)_RI
559 : ! | | | | | | |
560 : ! --> | 3 | 7 | 11| 15| SUM --> | |
561 : ! |___|___|___|___| |_|
562 : ! | | | | | | |
563 : ! --> | 4 | 8 | 12| 16| SUM --> | |
564 : ! |___|___|___|___| |_|
565 : !
566 :
567 46 : color_col = fm_mat_S%matrix_struct%context%mepos(2)
568 46 : ALLOCATE (para_env_col)
569 46 : CALL para_env_col%from_split(para_env_RPA, color_col)
570 :
571 46 : CALL para_env_col%sum(iaia_RI_dp)
572 :
573 : ! convert the iaia_RI_dp into double-double precision
574 136 : ALLOCATE (iaia_RI(ncol_local))
575 2764 : DO iiB = 1, ncol_local
576 2764 : iaia_RI(iiB) = iaia_RI_dp(iiB)
577 : END DO
578 46 : DEALLOCATE (iaia_RI_dp)
579 :
580 : ! 4) calculate the right hand term, D_ia is the matrix containing the
581 : ! orbital energy differences, M_ia is the diagonal of the full RPA 'excitation'
582 : ! matrix
583 136 : ALLOCATE (D_ia(ncol_local))
584 :
585 90 : ALLOCATE (M_ia(ncol_local))
586 :
587 2764 : DO iiB = 1, ncol_local
588 2718 : i_global = col_indices(iiB)
589 :
590 2718 : iocc = MAX(1, i_global - 1)/virtual + 1
591 2718 : avirt = i_global - (iocc - 1)*virtual
592 2718 : eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
593 :
594 2764 : D_ia(iiB) = eigen_diff
595 : END DO
596 :
597 2764 : DO iiB = 1, ncol_local
598 2764 : M_ia(iiB) = D_ia(iiB)*D_ia(iiB) + 2.0_dp*D_ia(iiB)*iaia_RI(iiB)
599 : END DO
600 :
601 46 : right_term_ref = 0.0_dp
602 2764 : DO iiB = 1, ncol_local
603 2764 : right_term_ref = right_term_ref + (SQRT(M_ia(iiB)) - D_ia(iiB) - iaia_RI(iiB))
604 : END DO
605 46 : right_term_ref = right_term_ref/2.0_dp
606 :
607 : ! sum the result with the processes of the RPA_group having the same row
608 46 : color_row = fm_mat_S%matrix_struct%context%mepos(1)
609 46 : ALLOCATE (para_env_row)
610 46 : CALL para_env_row%from_split(para_env_RPA, color_row)
611 :
612 : ! allocate communication array for rows
613 46 : CALL para_env_row%sum(right_term_ref)
614 :
615 46 : CALL mp_para_env_release(para_env_row)
616 :
617 46 : CALL timestop(handle)
618 :
619 46 : END SUBROUTINE calc_ia_ia_integrals
620 :
621 : ! **************************************************************************************************
622 : !> \brief ...
623 : !> \param a_scaling ...
624 : !> \param left_term ...
625 : !> \param first_deriv ...
626 : !> \param num_integ_points ...
627 : !> \param my_open_shell ...
628 : !> \param M_ia ...
629 : !> \param cottj ...
630 : !> \param wj ...
631 : !> \param D_ia ...
632 : !> \param D_ia_beta ...
633 : !> \param M_ia_beta ...
634 : !> \param ncol_local ...
635 : !> \param ncol_local_beta ...
636 : !> \param num_integ_group ...
637 : !> \param color_rpa_group ...
638 : !> \param para_env ...
639 : !> \param para_env_col ...
640 : !> \param para_env_col_beta ...
641 : ! **************************************************************************************************
642 428 : SUBROUTINE calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
643 : M_ia, cottj, wj, D_ia, D_ia_beta, M_ia_beta, &
644 : ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
645 : para_env, para_env_col, para_env_col_beta)
646 : REAL(KIND=dp), INTENT(IN) :: a_scaling
647 : REAL(KIND=dp), INTENT(INOUT) :: left_term, first_deriv
648 : INTEGER, INTENT(IN) :: num_integ_points
649 : LOGICAL, INTENT(IN) :: my_open_shell
650 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
651 : INTENT(IN) :: M_ia, cottj, wj, D_ia, D_ia_beta, &
652 : M_ia_beta
653 : INTEGER, INTENT(IN) :: ncol_local, ncol_local_beta, &
654 : num_integ_group, color_rpa_group
655 : TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_col
656 : TYPE(mp_para_env_type), POINTER :: para_env_col_beta
657 :
658 : INTEGER :: iiB, jquad
659 : REAL(KIND=dp) :: first_deriv_beta, left_term_beta, omega
660 :
661 428 : left_term = 0.0_dp
662 428 : first_deriv = 0.0_dp
663 428 : left_term_beta = 0.0_dp
664 428 : first_deriv_beta = 0.0_dp
665 4452 : DO jquad = 1, num_integ_points
666 : ! parallelize over integration points
667 4024 : IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE
668 2212 : omega = a_scaling*cottj(jquad)
669 :
670 162484 : DO iiB = 1, ncol_local
671 : ! parallelize over ia elements in the para_env_row group
672 160272 : IF (MODULO(iiB, para_env_col%num_pe) /= para_env_col%mepos) CYCLE
673 : ! calculate left_term
674 : left_term = left_term + wj(jquad)* &
675 : (LOG(1.0_dp + (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2)) - &
676 145072 : (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2))
677 : first_deriv = first_deriv + wj(jquad)*cottj(jquad)**2* &
678 162484 : ((-M_ia(iiB) + D_ia(iiB)**2)**2/((omega**2 + D_ia(iiB)**2)**2*(omega**2 + M_ia(iiB))))
679 : END DO
680 :
681 2640 : IF (my_open_shell) THEN
682 14490 : DO iiB = 1, ncol_local_beta
683 : ! parallelize over ia elements in the para_env_row group
684 14140 : IF (MODULO(iiB, para_env_col_beta%num_pe) /= para_env_col_beta%mepos) CYCLE
685 : ! calculate left_term
686 : left_term_beta = left_term_beta + wj(jquad)* &
687 : (LOG(1.0_dp + (M_ia_beta(iiB) - D_ia_beta(iiB)**2)/(omega**2 + D_ia_beta(iiB)**2)) - &
688 14140 : (M_ia_beta(iiB) - D_ia_beta(iiB)**2)/(omega**2 + D_ia_beta(iiB)**2))
689 : first_deriv_beta = &
690 : first_deriv_beta + wj(jquad)*cottj(jquad)**2* &
691 14490 : ((-M_ia_beta(iiB) + D_ia_beta(iiB)**2)**2/((omega**2 + D_ia_beta(iiB)**2)**2*(omega**2 + M_ia_beta(iiB))))
692 : END DO
693 : END IF
694 :
695 : END DO
696 :
697 : ! sum the contribution from all proc, starting form the row group
698 428 : CALL para_env%sum(left_term)
699 428 : CALL para_env%sum(first_deriv)
700 :
701 428 : IF (my_open_shell) THEN
702 70 : CALL para_env%sum(left_term_beta)
703 70 : CALL para_env%sum(first_deriv_beta)
704 :
705 70 : left_term = left_term + left_term_beta
706 70 : first_deriv = first_deriv + first_deriv_beta
707 : END IF
708 :
709 428 : END SUBROUTINE calculate_objfunc
710 :
711 : ! **************************************************************************************************
712 : !> \brief Calculate integration weights for the tau grid (in dependency of the omega node)
713 : !> \param num_integ_points ...
714 : !> \param tau_tj ...
715 : !> \param weights_cos_tf_t_to_w ...
716 : !> \param omega_tj ...
717 : !> \param E_min ...
718 : !> \param E_max ...
719 : !> \param max_error ...
720 : !> \param num_points_per_magnitude ...
721 : !> \param regularization ...
722 : ! **************************************************************************************************
723 70 : SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, omega_tj, &
724 : E_min, E_max, max_error, num_points_per_magnitude, &
725 : regularization)
726 :
727 : INTEGER, INTENT(IN) :: num_integ_points
728 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
729 : INTENT(IN) :: tau_tj
730 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
731 : INTENT(INOUT) :: weights_cos_tf_t_to_w
732 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
733 : INTENT(IN) :: omega_tj
734 : REAL(KIND=dp), INTENT(IN) :: E_min, E_max
735 : REAL(KIND=dp), INTENT(INOUT) :: max_error
736 : INTEGER, INTENT(IN) :: num_points_per_magnitude
737 : REAL(KIND=dp), INTENT(IN) :: regularization
738 :
739 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_cos_tf_t_to_w'
740 :
741 : INTEGER :: handle, iii, info, jjj, jquad, lwork, &
742 : num_x_nodes
743 70 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
744 : REAL(KIND=dp) :: multiplicator, omega
745 70 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_UTy, work, &
746 : x_values, y_values
747 70 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, &
748 70 : mat_SinvVSinvT, mat_U
749 :
750 70 : CALL timeset(routineN, handle)
751 :
752 : ! take num_points_per_magnitude points per 10-interval
753 70 : num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude
754 :
755 : ! take at least as many x points as integration points to have clear
756 : ! input for the singular value decomposition
757 70 : num_x_nodes = MAX(num_x_nodes, num_integ_points)
758 :
759 210 : ALLOCATE (x_values(num_x_nodes))
760 24070 : x_values = 0.0_dp
761 140 : ALLOCATE (y_values(num_x_nodes))
762 24070 : y_values = 0.0_dp
763 280 : ALLOCATE (mat_A(num_x_nodes, num_integ_points))
764 150976 : mat_A = 0.0_dp
765 210 : ALLOCATE (tau_wj_work(num_integ_points))
766 576 : tau_wj_work = 0.0_dp
767 140 : ALLOCATE (sing_values(num_integ_points))
768 576 : sing_values = 0.0_dp
769 280 : ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
770 9144070 : mat_U = 0.0_dp
771 210 : ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
772 :
773 150976 : mat_SinvVSinvT = 0.0_dp
774 : ! double the value nessary for 'A' to achieve good performance
775 70 : lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
776 210 : ALLOCATE (work(lwork))
777 105230 : work = 0.0_dp
778 210 : ALLOCATE (iwork(8*num_integ_points))
779 4118 : iwork = 0
780 210 : ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
781 174470 : mat_SinvVSinvSigma = 0.0_dp
782 140 : ALLOCATE (vec_UTy(num_x_nodes))
783 24070 : vec_UTy = 0.0_dp
784 :
785 70 : max_error = 0.0_dp
786 :
787 : ! loop over all omega frequency points
788 576 : DO jquad = 1, num_integ_points
789 :
790 : ! set the x-values logarithmically in the interval [Emin,Emax]
791 506 : multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
792 150906 : DO iii = 1, num_x_nodes
793 150906 : x_values(iii) = E_min*multiplicator**(iii - 1)
794 : END DO
795 :
796 506 : omega = omega_tj(jquad)
797 :
798 : ! y=2x/(x^2+omega_k^2)
799 150906 : DO iii = 1, num_x_nodes
800 150906 : y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2 + omega**2)
801 : END DO
802 :
803 : ! calculate mat_A
804 6892 : DO jjj = 1, num_integ_points
805 1590892 : DO iii = 1, num_x_nodes
806 1590386 : mat_A(iii, jjj) = COS(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
807 : END DO
808 : END DO
809 :
810 : ! Singular value decomposition of mat_A
811 : CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
812 506 : mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
813 :
814 506 : CPASSERT(info == 0)
815 :
816 : ! integration weights = V Sigma U^T y
817 : ! 1) V*Sigma
818 6892 : DO jjj = 1, num_integ_points
819 114582 : DO iii = 1, num_integ_points
820 : ! mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj)
821 : mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)*sing_values(jjj) &
822 114076 : /(regularization**2 + sing_values(jjj)**2)
823 : END DO
824 : END DO
825 :
826 : ! 2) U^T y
827 : CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, &
828 506 : 0.0_dp, vec_UTy, num_x_nodes)
829 :
830 : ! 3) (V*Sigma) * (U^T y)
831 : CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, &
832 506 : num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
833 :
834 6892 : weights_cos_tf_t_to_w(jquad, :) = tau_wj_work(:)
835 :
836 : CALL calc_max_error_fit_tau_grid_with_cosine(max_error, omega, tau_tj, tau_wj_work, x_values, &
837 576 : y_values, num_integ_points, num_x_nodes)
838 :
839 : END DO ! jquad
840 :
841 0 : DEALLOCATE (x_values, y_values, mat_A, tau_wj_work, sing_values, mat_U, mat_SinvVSinvT, &
842 70 : work, iwork, mat_SinvVSinvSigma, vec_UTy)
843 :
844 70 : CALL timestop(handle)
845 :
846 70 : END SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w
847 :
848 : ! **************************************************************************************************
849 : !> \brief Calculate integration weights for the tau grid (in dependency of the omega node)
850 : !> \param num_integ_points ...
851 : !> \param tau_tj ...
852 : !> \param weights_sin_tf_t_to_w ...
853 : !> \param omega_tj ...
854 : !> \param E_min ...
855 : !> \param E_max ...
856 : !> \param max_error ...
857 : !> \param num_points_per_magnitude ...
858 : !> \param regularization ...
859 : ! **************************************************************************************************
860 38 : SUBROUTINE get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, omega_tj, &
861 : E_min, E_max, max_error, num_points_per_magnitude, regularization)
862 :
863 : INTEGER, INTENT(IN) :: num_integ_points
864 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
865 : INTENT(IN) :: tau_tj
866 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
867 : INTENT(INOUT) :: weights_sin_tf_t_to_w
868 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
869 : INTENT(IN) :: omega_tj
870 : REAL(KIND=dp), INTENT(IN) :: E_min, E_max
871 : REAL(KIND=dp), INTENT(OUT) :: max_error
872 : INTEGER, INTENT(IN) :: num_points_per_magnitude
873 : REAL(KIND=dp), INTENT(IN) :: regularization
874 :
875 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_sin_tf_t_to_w'
876 :
877 : INTEGER :: handle, iii, info, jjj, jquad, lwork, &
878 : num_x_nodes
879 38 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
880 : REAL(KIND=dp) :: chi2_min_jquad, multiplicator, omega
881 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_UTy, work, &
882 38 : work_array, x_values, y_values
883 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, &
884 38 : mat_SinvVSinvT, mat_U
885 :
886 38 : CALL timeset(routineN, handle)
887 :
888 : ! take num_points_per_magnitude points per 10-interval
889 38 : num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude
890 :
891 : ! take at least as many x points as integration points to have clear
892 : ! input for the singular value decomposition
893 38 : num_x_nodes = MAX(num_x_nodes, num_integ_points)
894 :
895 114 : ALLOCATE (x_values(num_x_nodes))
896 12838 : x_values = 0.0_dp
897 76 : ALLOCATE (y_values(num_x_nodes))
898 12838 : y_values = 0.0_dp
899 152 : ALLOCATE (mat_A(num_x_nodes, num_integ_points))
900 119656 : mat_A = 0.0_dp
901 114 : ALLOCATE (tau_wj_work(num_integ_points))
902 456 : tau_wj_work = 0.0_dp
903 114 : ALLOCATE (work_array(2*num_integ_points))
904 874 : work_array = 0.0_dp
905 76 : ALLOCATE (sing_values(num_integ_points))
906 456 : sing_values = 0.0_dp
907 152 : ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
908 4972838 : mat_U = 0.0_dp
909 114 : ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
910 :
911 119656 : mat_SinvVSinvT = 0.0_dp
912 : ! double the value nessary for 'A' to achieve good performance
913 38 : lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
914 114 : ALLOCATE (work(lwork))
915 79758 : work = 0.0_dp
916 114 : ALLOCATE (iwork(8*num_integ_points))
917 3382 : iwork = 0
918 114 : ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
919 132038 : mat_SinvVSinvSigma = 0.0_dp
920 76 : ALLOCATE (vec_UTy(num_x_nodes))
921 12838 : vec_UTy = 0.0_dp
922 :
923 38 : max_error = 0.0_dp
924 :
925 : ! loop over all omega frequency points
926 456 : DO jquad = 1, num_integ_points
927 :
928 418 : chi2_min_jquad = 100.0_dp
929 :
930 : ! set the x-values logarithmically in the interval [Emin,Emax]
931 418 : multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
932 119618 : DO iii = 1, num_x_nodes
933 119618 : x_values(iii) = E_min*multiplicator**(iii - 1)
934 : END DO
935 :
936 418 : omega = omega_tj(jquad)
937 :
938 : ! y=2x/(x^2+omega_k^2)
939 119618 : DO iii = 1, num_x_nodes
940 : ! y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2+omega**2)
941 119618 : y_values(iii) = 2.0_dp*omega/((x_values(iii))**2 + omega**2)
942 : END DO
943 :
944 : ! calculate mat_A
945 6556 : DO jjj = 1, num_integ_points
946 1501756 : DO iii = 1, num_x_nodes
947 1501338 : mat_A(iii, jjj) = SIN(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
948 : END DO
949 : END DO
950 :
951 : ! Singular value decomposition of mat_A
952 : CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
953 418 : mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
954 :
955 418 : CPASSERT(info == 0)
956 :
957 : ! integration weights = V Sigma U^T y
958 : ! 1) V*Sigma
959 6556 : DO jjj = 1, num_integ_points
960 113534 : DO iii = 1, num_integ_points
961 : ! mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj)
962 : mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)*sing_values(jjj) &
963 113116 : /(regularization**2 + sing_values(jjj)**2)
964 : END DO
965 : END DO
966 :
967 : ! 2) U^T y
968 : CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, &
969 418 : 0.0_dp, vec_UTy, num_x_nodes)
970 :
971 : ! 3) (V*Sigma) * (U^T y)
972 : CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, &
973 418 : num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
974 :
975 6556 : weights_sin_tf_t_to_w(jquad, :) = tau_wj_work(:)
976 :
977 : CALL calc_max_error_fit_tau_grid_with_sine(max_error, omega, tau_tj, tau_wj_work, x_values, &
978 456 : y_values, num_integ_points, num_x_nodes)
979 :
980 : END DO ! jquad
981 :
982 0 : DEALLOCATE (x_values, y_values, mat_A, tau_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, &
983 38 : work, iwork, mat_SinvVSinvSigma, vec_UTy)
984 :
985 38 : CALL timestop(handle)
986 :
987 38 : END SUBROUTINE get_l_sq_wghts_sin_tf_t_to_w
988 :
989 : ! **************************************************************************************************
990 : !> \brief ...
991 : !> \param max_error ...
992 : !> \param omega ...
993 : !> \param tau_tj ...
994 : !> \param tau_wj_work ...
995 : !> \param x_values ...
996 : !> \param y_values ...
997 : !> \param num_integ_points ...
998 : !> \param num_x_nodes ...
999 : ! **************************************************************************************************
1000 506 : PURE SUBROUTINE calc_max_error_fit_tau_grid_with_cosine(max_error, omega, tau_tj, tau_wj_work, x_values, &
1001 : y_values, num_integ_points, num_x_nodes)
1002 :
1003 : REAL(KIND=dp), INTENT(INOUT) :: max_error, omega
1004 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1005 : INTENT(IN) :: tau_tj, tau_wj_work, x_values, y_values
1006 : INTEGER, INTENT(IN) :: num_integ_points, num_x_nodes
1007 :
1008 : INTEGER :: kkk
1009 : REAL(KIND=dp) :: func_val, func_val_temp, max_error_tmp
1010 :
1011 506 : max_error_tmp = 0.0_dp
1012 :
1013 150906 : DO kkk = 1, num_x_nodes
1014 :
1015 : func_val = 0.0_dp
1016 :
1017 150400 : CALL eval_fit_func_tau_grid_cosine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
1018 :
1019 150906 : IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
1020 2494 : max_error_tmp = ABS(y_values(kkk) - func_val)
1021 2494 : func_val_temp = func_val
1022 : END IF
1023 :
1024 : END DO
1025 :
1026 506 : IF (max_error_tmp > max_error) THEN
1027 :
1028 144 : max_error = max_error_tmp
1029 :
1030 : END IF
1031 :
1032 506 : END SUBROUTINE calc_max_error_fit_tau_grid_with_cosine
1033 :
1034 : ! **************************************************************************************************
1035 : !> \brief Evaluate fit function when calculating tau grid for cosine transform
1036 : !> \param func_val ...
1037 : !> \param x_value ...
1038 : !> \param num_integ_points ...
1039 : !> \param tau_tj ...
1040 : !> \param tau_wj_work ...
1041 : !> \param omega ...
1042 : ! **************************************************************************************************
1043 150400 : PURE SUBROUTINE eval_fit_func_tau_grid_cosine(func_val, x_value, num_integ_points, tau_tj, tau_wj_work, omega)
1044 :
1045 : REAL(KIND=dp), INTENT(OUT) :: func_val
1046 : REAL(KIND=dp), INTENT(IN) :: x_value
1047 : INTEGER, INTENT(IN) :: num_integ_points
1048 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1049 : INTENT(IN) :: tau_tj, tau_wj_work
1050 : REAL(KIND=dp), INTENT(IN) :: omega
1051 :
1052 : INTEGER :: iii
1053 :
1054 150400 : func_val = 0.0_dp
1055 :
1056 1734400 : DO iii = 1, num_integ_points
1057 :
1058 : ! calculate value of the fit function
1059 1734400 : func_val = func_val + tau_wj_work(iii)*COS(omega*tau_tj(iii))*EXP(-x_value*tau_tj(iii))
1060 :
1061 : END DO
1062 :
1063 150400 : END SUBROUTINE eval_fit_func_tau_grid_cosine
1064 :
1065 : ! **************************************************************************************************
1066 : !> \brief Evaluate fit function when calculating tau grid for sine transform
1067 : !> \param func_val ...
1068 : !> \param x_value ...
1069 : !> \param num_integ_points ...
1070 : !> \param tau_tj ...
1071 : !> \param tau_wj_work ...
1072 : !> \param omega ...
1073 : ! **************************************************************************************************
1074 119200 : PURE SUBROUTINE eval_fit_func_tau_grid_sine(func_val, x_value, num_integ_points, tau_tj, tau_wj_work, omega)
1075 :
1076 : REAL(KIND=dp), INTENT(INOUT) :: func_val
1077 : REAL(KIND=dp), INTENT(IN) :: x_value
1078 : INTEGER, INTENT(in) :: num_integ_points
1079 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1080 : INTENT(IN) :: tau_tj, tau_wj_work
1081 : REAL(KIND=dp), INTENT(IN) :: omega
1082 :
1083 : INTEGER :: iii
1084 :
1085 119200 : func_val = 0.0_dp
1086 :
1087 1614400 : DO iii = 1, num_integ_points
1088 :
1089 : ! calculate value of the fit function
1090 1614400 : func_val = func_val + tau_wj_work(iii)*SIN(omega*tau_tj(iii))*EXP(-x_value*tau_tj(iii))
1091 :
1092 : END DO
1093 :
1094 119200 : END SUBROUTINE eval_fit_func_tau_grid_sine
1095 :
1096 : ! **************************************************************************************************
1097 : !> \brief ...
1098 : !> \param max_error ...
1099 : !> \param omega ...
1100 : !> \param tau_tj ...
1101 : !> \param tau_wj_work ...
1102 : !> \param x_values ...
1103 : !> \param y_values ...
1104 : !> \param num_integ_points ...
1105 : !> \param num_x_nodes ...
1106 : ! **************************************************************************************************
1107 418 : PURE SUBROUTINE calc_max_error_fit_tau_grid_with_sine(max_error, omega, tau_tj, tau_wj_work, x_values, &
1108 : y_values, num_integ_points, num_x_nodes)
1109 :
1110 : REAL(KIND=dp), INTENT(INOUT) :: max_error, omega
1111 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1112 : INTENT(IN) :: tau_tj, tau_wj_work, x_values, y_values
1113 : INTEGER, INTENT(IN) :: num_integ_points, num_x_nodes
1114 :
1115 : INTEGER :: kkk
1116 : REAL(KIND=dp) :: func_val, func_val_temp, max_error_tmp
1117 :
1118 418 : max_error_tmp = 0.0_dp
1119 :
1120 119618 : DO kkk = 1, num_x_nodes
1121 :
1122 119200 : func_val = 0.0_dp
1123 :
1124 119200 : CALL eval_fit_func_tau_grid_sine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
1125 :
1126 119618 : IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
1127 2410 : max_error_tmp = ABS(y_values(kkk) - func_val)
1128 2410 : func_val_temp = func_val
1129 : END IF
1130 :
1131 : END DO
1132 :
1133 418 : IF (max_error_tmp > max_error) THEN
1134 :
1135 42 : max_error = max_error_tmp
1136 :
1137 : END IF
1138 :
1139 418 : END SUBROUTINE calc_max_error_fit_tau_grid_with_sine
1140 :
1141 : ! **************************************************************************************************
1142 : !> \brief test the singular value decomposition for the computation of integration weights for the
1143 : !> Fourier transform between time and frequency grid in cubic-scaling RPA
1144 : !> \param nR ...
1145 : !> \param iw ...
1146 : ! **************************************************************************************************
1147 0 : SUBROUTINE test_least_square_ft(nR, iw)
1148 : INTEGER, INTENT(IN) :: nR, iw
1149 :
1150 : INTEGER :: ierr, iR, jquad, num_integ_points
1151 : REAL(KIND=dp) :: max_error, multiplicator, Rc, Rc_max
1152 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tau_tj, tau_wj, tj, wj, x_tw
1153 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: weights_cos_tf_t_to_w
1154 :
1155 0 : Rc_max = 1.0E+7
1156 :
1157 0 : multiplicator = Rc_max**(1.0_dp/(REAL(nR, KIND=dp) - 1.0_dp))
1158 :
1159 0 : DO num_integ_points = 1, 20
1160 :
1161 0 : ALLOCATE (x_tw(2*num_integ_points))
1162 0 : x_tw = 0.0_dp
1163 0 : ALLOCATE (tau_tj(num_integ_points))
1164 0 : tau_tj = 0.0_dp
1165 0 : ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
1166 0 : weights_cos_tf_t_to_w = 0.0_dp
1167 0 : ALLOCATE (tau_wj(num_integ_points))
1168 0 : tau_wj = 0.0_dp
1169 0 : ALLOCATE (tj(num_integ_points))
1170 0 : tj = 0.0_dp
1171 0 : ALLOCATE (wj(num_integ_points))
1172 0 : wj = 0.0_dp
1173 :
1174 0 : DO iR = 0, nR - 1
1175 :
1176 0 : Rc = 2.0_dp*multiplicator**iR
1177 :
1178 0 : ierr = 0
1179 0 : CALL get_rpa_minimax_coeff(num_integ_points, Rc, x_tw, ierr, print_warning=.FALSE.)
1180 :
1181 0 : DO jquad = 1, num_integ_points
1182 0 : tj(jquad) = x_tw(jquad)
1183 0 : wj(jquad) = x_tw(jquad + num_integ_points)
1184 : END DO
1185 :
1186 0 : x_tw = 0.0_dp
1187 :
1188 0 : CALL get_exp_minimax_coeff(num_integ_points, Rc, x_tw)
1189 :
1190 0 : DO jquad = 1, num_integ_points
1191 0 : tau_tj(jquad) = x_tw(jquad)/2.0_dp
1192 0 : tau_wj(jquad) = x_tw(jquad + num_integ_points)/2.0_dp
1193 : END DO
1194 :
1195 : CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, &
1196 : weights_cos_tf_t_to_w, tj, &
1197 0 : 1.0_dp, Rc, max_error, 200, 0.0_dp)
1198 :
1199 0 : IF (iw > 0) THEN
1200 0 : WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, Rc, max_error
1201 : END IF
1202 :
1203 : END DO
1204 :
1205 0 : DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tau_wj, wj, tj)
1206 :
1207 : END DO
1208 :
1209 0 : END SUBROUTINE test_least_square_ft
1210 :
1211 : ! **************************************************************************************************
1212 : !> \brief ...
1213 : !> \param num_integ_points ...
1214 : !> \param tau_tj ...
1215 : !> \param weights_cos_tf_w_to_t ...
1216 : !> \param omega_tj ...
1217 : !> \param E_min ...
1218 : !> \param E_max ...
1219 : !> \param max_error ...
1220 : !> \param num_points_per_magnitude ...
1221 : !> \param regularization ...
1222 : ! **************************************************************************************************
1223 70 : SUBROUTINE get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, omega_tj, &
1224 : E_min, E_max, max_error, num_points_per_magnitude, regularization)
1225 :
1226 : INTEGER, INTENT(IN) :: num_integ_points
1227 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1228 : INTENT(IN) :: tau_tj
1229 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1230 : INTENT(INOUT) :: weights_cos_tf_w_to_t
1231 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1232 : INTENT(IN) :: omega_tj
1233 : REAL(KIND=dp), INTENT(IN) :: E_min, E_max
1234 : REAL(KIND=dp), INTENT(INOUT) :: max_error
1235 : INTEGER, INTENT(IN) :: num_points_per_magnitude
1236 : REAL(KIND=dp), INTENT(IN) :: regularization
1237 :
1238 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_l_sq_wghts_cos_tf_w_to_t'
1239 :
1240 : INTEGER :: handle, iii, info, jjj, jquad, lwork, &
1241 : num_x_nodes
1242 70 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1243 : REAL(KIND=dp) :: chi2_min_jquad, multiplicator, omega, &
1244 : tau, x_value
1245 70 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: omega_wj_work, sing_values, vec_UTy, &
1246 70 : work, work_array, x_values, y_values
1247 70 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, &
1248 70 : mat_SinvVSinvT, mat_U
1249 :
1250 70 : CALL timeset(routineN, handle)
1251 :
1252 : ! take num_points_per_magnitude points per 10-interval
1253 70 : num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude
1254 :
1255 : ! take at least as many x points as integration points to have clear
1256 : ! input for the singular value decomposition
1257 70 : num_x_nodes = MAX(num_x_nodes, num_integ_points)
1258 :
1259 210 : ALLOCATE (x_values(num_x_nodes))
1260 24070 : x_values = 0.0_dp
1261 140 : ALLOCATE (y_values(num_x_nodes))
1262 24070 : y_values = 0.0_dp
1263 280 : ALLOCATE (mat_A(num_x_nodes, num_integ_points))
1264 150976 : mat_A = 0.0_dp
1265 210 : ALLOCATE (omega_wj_work(num_integ_points))
1266 576 : omega_wj_work = 0.0_dp
1267 210 : ALLOCATE (work_array(2*num_integ_points))
1268 1082 : work_array = 0.0_dp
1269 140 : ALLOCATE (sing_values(num_integ_points))
1270 576 : sing_values = 0.0_dp
1271 280 : ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
1272 9144070 : mat_U = 0.0_dp
1273 210 : ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
1274 :
1275 150976 : mat_SinvVSinvT = 0.0_dp
1276 : ! double the value nessary for 'A' to achieve good performance
1277 70 : lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
1278 210 : ALLOCATE (work(lwork))
1279 105230 : work = 0.0_dp
1280 210 : ALLOCATE (iwork(8*num_integ_points))
1281 4118 : iwork = 0
1282 210 : ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
1283 174470 : mat_SinvVSinvSigma = 0.0_dp
1284 140 : ALLOCATE (vec_UTy(num_x_nodes))
1285 24070 : vec_UTy = 0.0_dp
1286 :
1287 : ! set the x-values logarithmically in the interval [Emin,Emax]
1288 70 : multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
1289 24070 : DO iii = 1, num_x_nodes
1290 24070 : x_values(iii) = E_min*multiplicator**(iii - 1)
1291 : END DO
1292 :
1293 70 : max_error = 0.0_dp
1294 :
1295 : ! loop over all tau time points
1296 576 : DO jquad = 1, num_integ_points
1297 :
1298 506 : chi2_min_jquad = 100.0_dp
1299 :
1300 506 : tau = tau_tj(jquad)
1301 :
1302 : ! y=exp(-x*|tau_k|)
1303 150906 : DO iii = 1, num_x_nodes
1304 150906 : y_values(iii) = EXP(-x_values(iii)*tau)
1305 : END DO
1306 :
1307 : ! calculate mat_A
1308 6892 : DO jjj = 1, num_integ_points
1309 1590892 : DO iii = 1, num_x_nodes
1310 1584000 : omega = omega_tj(jjj)
1311 1584000 : x_value = x_values(iii)
1312 1590386 : mat_A(iii, jjj) = COS(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2)
1313 : END DO
1314 : END DO
1315 :
1316 : ! Singular value decomposition of mat_A
1317 : CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
1318 506 : mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
1319 :
1320 506 : CPASSERT(info == 0)
1321 :
1322 : ! integration weights = V Sigma U^T y
1323 : ! 1) V*Sigma
1324 6892 : DO jjj = 1, num_integ_points
1325 114582 : DO iii = 1, num_integ_points
1326 : ! mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj)
1327 : mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)*sing_values(jjj) &
1328 114076 : /(regularization**2 + sing_values(jjj)**2)
1329 : END DO
1330 : END DO
1331 :
1332 : ! 2) U^T y
1333 : CALL DGEMM('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, y_values, num_x_nodes, &
1334 506 : 0.0_dp, vec_UTy, num_x_nodes)
1335 :
1336 : ! 3) (V*Sigma) * (U^T y)
1337 : CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, &
1338 506 : num_x_nodes, 0.0_dp, omega_wj_work, num_integ_points)
1339 :
1340 6892 : weights_cos_tf_w_to_t(jquad, :) = omega_wj_work(:)
1341 :
1342 : CALL calc_max_error_fit_omega_grid_with_cosine(max_error, tau, omega_tj, omega_wj_work, x_values, &
1343 576 : y_values, num_integ_points, num_x_nodes)
1344 :
1345 : END DO ! jquad
1346 :
1347 0 : DEALLOCATE (x_values, y_values, mat_A, omega_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, &
1348 70 : work, iwork, mat_SinvVSinvSigma, vec_UTy)
1349 :
1350 70 : CALL timestop(handle)
1351 :
1352 70 : END SUBROUTINE get_l_sq_wghts_cos_tf_w_to_t
1353 :
1354 : ! **************************************************************************************************
1355 : !> \brief ...
1356 : !> \param max_error ...
1357 : !> \param tau ...
1358 : !> \param omega_tj ...
1359 : !> \param omega_wj_work ...
1360 : !> \param x_values ...
1361 : !> \param y_values ...
1362 : !> \param num_integ_points ...
1363 : !> \param num_x_nodes ...
1364 : ! **************************************************************************************************
1365 506 : SUBROUTINE calc_max_error_fit_omega_grid_with_cosine(max_error, tau, omega_tj, omega_wj_work, x_values, &
1366 : y_values, num_integ_points, num_x_nodes)
1367 :
1368 : REAL(KIND=dp), INTENT(INOUT) :: max_error
1369 : REAL(KIND=dp), INTENT(IN) :: tau
1370 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1371 : INTENT(IN) :: omega_tj, omega_wj_work, x_values, &
1372 : y_values
1373 : INTEGER, INTENT(IN) :: num_integ_points, num_x_nodes
1374 :
1375 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_max_error_fit_omega_grid_with_cosine'
1376 :
1377 : INTEGER :: handle, kkk
1378 : REAL(KIND=dp) :: func_val, func_val_temp, max_error_tmp
1379 :
1380 506 : CALL timeset(routineN, handle)
1381 :
1382 506 : max_error_tmp = 0.0_dp
1383 :
1384 150906 : DO kkk = 1, num_x_nodes
1385 :
1386 : func_val = 0.0_dp
1387 :
1388 150400 : CALL eval_fit_func_omega_grid_cosine(func_val, x_values(kkk), num_integ_points, omega_tj, omega_wj_work, tau)
1389 :
1390 150906 : IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
1391 4244 : max_error_tmp = ABS(y_values(kkk) - func_val)
1392 4244 : func_val_temp = func_val
1393 : END IF
1394 :
1395 : END DO
1396 :
1397 506 : IF (max_error_tmp > max_error) THEN
1398 :
1399 132 : max_error = max_error_tmp
1400 :
1401 : END IF
1402 :
1403 506 : CALL timestop(handle)
1404 :
1405 506 : END SUBROUTINE calc_max_error_fit_omega_grid_with_cosine
1406 :
1407 : ! **************************************************************************************************
1408 : !> \brief ...
1409 : !> \param func_val ...
1410 : !> \param x_value ...
1411 : !> \param num_integ_points ...
1412 : !> \param omega_tj ...
1413 : !> \param omega_wj_work ...
1414 : !> \param tau ...
1415 : ! **************************************************************************************************
1416 150400 : PURE SUBROUTINE eval_fit_func_omega_grid_cosine(func_val, x_value, num_integ_points, omega_tj, omega_wj_work, tau)
1417 : REAL(KIND=dp), INTENT(OUT) :: func_val
1418 : REAL(KIND=dp), INTENT(IN) :: x_value
1419 : INTEGER, INTENT(IN) :: num_integ_points
1420 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1421 : INTENT(IN) :: omega_tj, omega_wj_work
1422 : REAL(KIND=dp), INTENT(IN) :: tau
1423 :
1424 : INTEGER :: iii
1425 : REAL(KIND=dp) :: omega
1426 :
1427 150400 : func_val = 0.0_dp
1428 :
1429 1734400 : DO iii = 1, num_integ_points
1430 :
1431 : ! calculate value of the fit function
1432 1584000 : omega = omega_tj(iii)
1433 1734400 : func_val = func_val + omega_wj_work(iii)*COS(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2)
1434 :
1435 : END DO
1436 :
1437 150400 : END SUBROUTINE eval_fit_func_omega_grid_cosine
1438 :
1439 : ! **************************************************************************************************
1440 : !> \brief ...
1441 : !> \param qs_env ...
1442 : !> \param para_env ...
1443 : !> \param gap ...
1444 : !> \param max_eig_diff ...
1445 : !> \param e_fermi ...
1446 : ! **************************************************************************************************
1447 8 : SUBROUTINE gap_and_max_eig_diff_kpoints(qs_env, para_env, gap, max_eig_diff, e_fermi)
1448 :
1449 : TYPE(qs_environment_type), POINTER :: qs_env
1450 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1451 : REAL(KIND=dp), INTENT(OUT) :: gap, max_eig_diff, e_fermi
1452 :
1453 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gap_and_max_eig_diff_kpoints'
1454 :
1455 : INTEGER :: handle, homo, ikpgr, ispin, kplocal, &
1456 : nmo, nspin
1457 : INTEGER, DIMENSION(2) :: kp_range
1458 : REAL(KIND=dp) :: e_homo, e_homo_temp, e_lumo, e_lumo_temp
1459 : REAL(KIND=dp), DIMENSION(3) :: tmp
1460 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1461 : TYPE(kpoint_env_type), POINTER :: kp
1462 : TYPE(kpoint_type), POINTER :: kpoint
1463 : TYPE(mo_set_type), POINTER :: mo_set
1464 :
1465 4 : CALL timeset(routineN, handle)
1466 :
1467 : CALL get_qs_env(qs_env, &
1468 4 : kpoints=kpoint)
1469 :
1470 4 : mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1471 4 : CALL get_mo_set(mo_set, nmo=nmo)
1472 :
1473 4 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1474 4 : kplocal = kp_range(2) - kp_range(1) + 1
1475 :
1476 4 : gap = 1000.0_dp
1477 4 : max_eig_diff = 0.0_dp
1478 4 : e_homo = -1000.0_dp
1479 4 : e_lumo = 1000.0_dp
1480 :
1481 12 : DO ikpgr = 1, kplocal
1482 8 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1483 8 : nspin = SIZE(kp%mos, 2)
1484 20 : DO ispin = 1, nspin
1485 8 : mo_set => kp%mos(1, ispin)
1486 8 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, homo=homo)
1487 8 : e_homo_temp = eigenvalues(homo)
1488 8 : e_lumo_temp = eigenvalues(homo + 1)
1489 :
1490 8 : IF (e_homo_temp > e_homo) e_homo = e_homo_temp
1491 8 : IF (e_lumo_temp < e_lumo) e_lumo = e_lumo_temp
1492 24 : IF (eigenvalues(nmo) - eigenvalues(1) > max_eig_diff) max_eig_diff = eigenvalues(nmo) - eigenvalues(1)
1493 :
1494 : END DO
1495 : END DO
1496 :
1497 : ! Collect all three numbers in an array
1498 : ! Reverse sign of lumo to reduce number of MPI calls
1499 4 : tmp(1) = e_homo
1500 4 : tmp(2) = -e_lumo
1501 4 : tmp(3) = max_eig_diff
1502 4 : CALL para_env%max(tmp)
1503 :
1504 4 : gap = -tmp(2) - tmp(1)
1505 4 : e_fermi = (tmp(1) - tmp(2))*0.5_dp
1506 4 : max_eig_diff = tmp(3)
1507 :
1508 4 : CALL timestop(handle)
1509 :
1510 4 : END SUBROUTINE
1511 :
1512 : ! **************************************************************************************************
1513 : !> \brief returns minimal and maximal energy values for the E_range for the minimax grid selection
1514 : !> \param qs_env ...
1515 : !> \param para_env ...
1516 : !> \param homo index of the homo level for the respective spin channel
1517 : !> \param Eigenval eigenvalues
1518 : !> \param do_ri_sos_laplace_mp2 flag for SOS-MP2
1519 : !> \param do_kpoints_cubic_RPA flag for cubic-scaling RPA with k-points
1520 : !> \param Emin minimal eigenvalue difference (gap of the system)
1521 : !> \param Emax maximal eigenvalue difference
1522 : !> \param e_range ...
1523 : !> \param e_fermi Fermi level
1524 : ! **************************************************************************************************
1525 204 : SUBROUTINE determine_energy_range(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, &
1526 : do_kpoints_cubic_RPA, Emin, Emax, e_range, e_fermi)
1527 :
1528 : TYPE(qs_environment_type), POINTER :: qs_env
1529 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1530 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
1531 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval
1532 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, &
1533 : do_kpoints_cubic_RPA
1534 : REAL(KIND=dp), INTENT(OUT) :: Emin, Emax, e_range, e_fermi
1535 :
1536 : CHARACTER(LEN=*), PARAMETER :: routineN = 'determine_energy_range'
1537 :
1538 : INTEGER :: handle, ispin, nspins
1539 : LOGICAL :: my_do_kpoints
1540 : TYPE(section_vals_type), POINTER :: input
1541 :
1542 204 : CALL timeset(routineN, handle)
1543 : ! Test for spin unrestricted
1544 204 : nspins = SIZE(homo)
1545 :
1546 : ! Test whether all necessary variables are available
1547 204 : my_do_kpoints = .FALSE.
1548 204 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1549 146 : my_do_kpoints = do_kpoints_cubic_RPA
1550 : END IF
1551 :
1552 146 : IF (my_do_kpoints) THEN
1553 4 : CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi)
1554 4 : E_Range = Emax/Emin
1555 : ELSE
1556 200 : IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN
1557 152 : Emin = HUGE(dp)
1558 152 : Emax = 0.0_dp
1559 340 : DO ispin = 1, nspins
1560 340 : IF (homo(ispin) > 0) THEN
1561 184 : Emin = MIN(Emin, Eigenval(homo(ispin) + 1, 1, ispin) - Eigenval(homo(ispin), 1, ispin))
1562 14732 : Emax = MAX(Emax, MAXVAL(Eigenval(:, :, ispin)) - MINVAL(Eigenval(:, :, ispin)))
1563 : END IF
1564 : END DO
1565 152 : E_Range = Emax/Emin
1566 152 : qs_env%mp2_env%e_range = e_range
1567 152 : qs_env%mp2_env%e_gap = Emin
1568 :
1569 152 : CALL get_qs_env(qs_env, input=input)
1570 152 : CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=e_range)
1571 152 : CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_GAP", r_val=emin)
1572 : ELSE
1573 48 : E_range = qs_env%mp2_env%E_range
1574 48 : Emin = qs_env%mp2_env%E_gap
1575 48 : Emax = Emin*E_range
1576 : END IF
1577 : END IF
1578 :
1579 : ! When we perform SOS-MP2, we need an additional factor of 2 for the energies (compare with mp2_laplace.F)
1580 : ! We do not need weights etc. for the cosine transform
1581 : ! We do not scale Emax because it is not needed for SOS-MP2
1582 204 : IF (do_ri_sos_laplace_mp2) THEN
1583 58 : Emin = Emin*2.0_dp
1584 58 : Emax = Emax*2.0_dp
1585 : END IF
1586 :
1587 204 : CALL timestop(handle)
1588 204 : END SUBROUTINE
1589 :
1590 : END MODULE mp2_grids
|