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