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