(git:fc9bb57)
Loading...
Searching...
No Matches
mp2_grids.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief 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! **************************************************************************************************
14MODULE 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
47CONTAINS
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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: x_tw
105
106 CALL timeset(routinen, handle)
107
108 CALL determine_energy_range(qs_env, para_env, homo, eigenval, do_ri_sos_laplace_mp2, &
109 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 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 "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 weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, ierr)
124
125! Shortcut if Greenx was available and successful
126 IF (ierr == 0) THEN
127 CALL timestop(handle)
128 RETURN
129 END IF
130
131 ! Test for spin unrestricted
132 nspins = SIZE(homo)
133
134 ! Test whether all necessary variables are available
135 my_do_kpoints = .false.
136 IF (.NOT. do_ri_sos_laplace_mp2) THEN
137 my_do_kpoints = do_kpoints_cubic_rpa
138 END IF
139
140 my_regularization = 0.0_dp
141 IF (PRESENT(regularization)) THEN
142 my_regularization = regularization
143
144 IF (num_integ_points > 20 .AND. e_range < 100.0_dp) THEN
145 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 "a larger basis set with higher angular momenta or alternatively using all-electron calculations.")
151 END IF
152 END IF
153
154 IF (.NOT. do_ri_sos_laplace_mp2) THEN
155 ALLOCATE (x_tw(2*num_integ_points))
156 x_tw = 0.0_dp
157 ierr = 0
158 IF (num_integ_points <= 20) THEN
159 CALL get_rpa_minimax_coeff(num_integ_points, e_range, x_tw, ierr)
160 ELSE
161 CALL get_rpa_minimax_coeff_larger_grid(num_integ_points, e_range, x_tw)
162 END IF
163
164 ALLOCATE (tj(num_integ_points))
165 tj = 0.0_dp
166
167 ALLOCATE (wj(num_integ_points))
168 wj = 0.0_dp
169
170 DO jquad = 1, num_integ_points
171 tj(jquad) = x_tw(jquad)
172 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 IF (num_integ_points >= 26) THEN
177 wj(:) = wj(:)*4.0_dp
178 END IF
179
180 DEALLOCATE (x_tw)
181
182 IF (unit_nr > 0 .AND. do_print) THEN
183 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
184 "MINIMAX_INFO| Number of integration points:", num_integ_points
185 WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.4)") &
186 "MINIMAX_INFO| Gap for the minimax approximation:", emin
187 WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.4)") &
188 "MINIMAX_INFO| Range for the minimax approximation:", e_range
189 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") "MINIMAX_INFO| Minimax parameters:", "Weights", "Abscissas"
190 DO jquad = 1, num_integ_points
191 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") wj(jquad), tj(jquad)
192 END DO
193 CALL m_flush(unit_nr)
194 END IF
195
196 ! scale the minimax parameters
197 tj(:) = tj(:)*emin
198 wj(:) = wj(:)*emin
199 END IF
200
201 ! set up the minimax time grid
202 IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
203
204 ALLOCATE (x_tw(2*num_integ_points))
205 x_tw = 0.0_dp
206
207 IF (num_integ_points <= 20) THEN
208 CALL get_exp_minimax_coeff(num_integ_points, e_range, x_tw)
209 ELSE
210 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 scaling = 2.0_dp
215 IF (do_ri_sos_laplace_mp2) scaling = 1.0_dp
216
217 ALLOCATE (tau_tj(num_integ_points))
218 tau_tj = 0.0_dp
219
220 ALLOCATE (tau_wj(num_integ_points))
221 tau_wj = 0.0_dp
222
223 DO jquad = 1, num_integ_points
224 tau_tj(jquad) = x_tw(jquad)/scaling
225 tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling
226 END DO
227
228 DEALLOCATE (x_tw)
229
230 IF (unit_nr > 0 .AND. do_print) THEN
231 WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.4)") &
232 "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 "MINIMAX_INFO| Gap:", emin
236 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
237 "MINIMAX_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas"
238 DO jquad = 1, num_integ_points
239 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad)
240 END DO
241 CALL m_flush(unit_nr)
242 END IF
243
244 ! scale grid from [1,R] to [Emin,Emax]
245 tau_tj(:) = tau_tj(:)/emin
246 tau_wj(:) = tau_wj(:)/emin
247
248 IF (.NOT. do_ri_sos_laplace_mp2) THEN
249 ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
250 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 my_regularization)
255
256 ! get the weights for the cosine transform W^c(iw) -> W^c(it)
257 ALLOCATE (weights_cos_tf_w_to_t(num_integ_points, num_integ_points))
258 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 my_regularization)
263
264 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 ALLOCATE (weights_sin_tf_t_to_w(num_integ_points, num_integ_points))
268 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 my_regularization)
273
274 IF (unit_nr > 0) THEN
275 WRITE (unit=unit_nr, fmt="(T3,A,T66,ES15.2)") &
276 "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 CALL timestop(handle)
286
287 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 SUBROUTINE get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, &
308 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 CALL timeset(routinen, handle)
330
331 nspins = SIZE(homo)
332 my_open_shell = (nspins == 2)
333
334 ! Now, start to prepare the different grid
335 ALLOCATE (tj(num_integ_points))
336 tj = 0.0_dp
337
338 ALLOCATE (wj(num_integ_points))
339 wj = 0.0_dp
340
341 DO jquad = 1, num_integ_points - 1
342 tj(jquad) = jquad*pi/(2.0_dp*num_integ_points)
343 wj(jquad) = pi/(num_integ_points*sin(tj(jquad))**2)
344 END DO
345 tj(num_integ_points) = pi/2.0_dp
346 wj(num_integ_points) = pi/(2.0_dp*num_integ_points*sin(tj(num_integ_points))**2)
347
348 IF (my_do_gw .AND. ext_scaling > 0.0_dp) THEN
349 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 tj, wj, fm_mat_s)
354 END IF
355
356 IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.5)') 'INTEG_INFO| Scaling parameter:', a_scaling
357
358 wj(:) = wj(:)*a_scaling
359
360 CALL timestop(handle)
361
362 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cottj, d_ia, d_ia_beta, iaia_ri, &
400 iaia_ri_beta, m_ia, m_ia_beta
401 TYPE(mp_para_env_type), POINTER :: para_env_col, para_env_col_beta
402
403 CALL timeset(routinen, handle)
404
405 nspins = SIZE(homo)
406 my_open_shell = (nspins == 2)
407
408 eps = 1.0e-10_dp
409
410 ALLOCATE (cottj(num_integ_points))
411
412 ! calculate the cotangent of the abscissa tj
413 DO jquad = 1, num_integ_points
414 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 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 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 d_ia_beta, iaia_ri_beta, m_ia_beta, fm_mat_s(2), para_env_col_beta)
424
425 right_term_ref = right_term_ref + right_term_ref_beta
426 END IF
427
428 ! bcast the result
429 IF (para_env%mepos == 0) THEN
430 CALL para_env%bcast(right_term_ref, 0)
431 ELSE
432 right_term_ref = 0.0_dp
433 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 conv_param = 100.0_dp*epsilon(right_term_ref)
439 step = 0.5_dp
440 a_low = 0.0_dp
441 a_high = step
442 right_term = -right_term_ref
443 DO icycle = 1, num_integ_points*2
444 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 para_env, para_env_col, para_env_col_beta)
450 left_term = left_term/4.0_dp/pi*a_scaling
451
452 IF (abs(left_term) > abs(right_term) .OR. abs(left_term + right_term) <= conv_param) EXIT
453 a_low = a_high
454 a_high = a_high + step
455
456 END DO
457
458 IF (abs(left_term + right_term) >= conv_param) THEN
459 IF (a_scaling >= 2*num_integ_points*step) THEN
460 a_scaling = 1.0_dp
461 ELSE
462
463 DO icycle = 1, num_integ_points*2
464 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 para_env, para_env_col, para_env_col_beta)
470 left_term = left_term/4.0_dp/pi*a_scaling
471
472 IF (abs(left_term) > abs(right_term)) THEN
473 a_high = a_scaling
474 ELSE
475 a_low = a_scaling
476 END IF
477
478 IF (abs(a_high - a_low) < 1.0e-5_dp) EXIT
479
480 END DO
481
482 END IF
483 END IF
484
485 a_scaling_ext = a_scaling
486 CALL para_env%bcast(a_scaling_ext, 0)
487
488 DEALLOCATE (cottj)
489 DEALLOCATE (iaia_ri)
490 DEALLOCATE (d_ia)
491 DEALLOCATE (m_ia)
492 CALL mp_para_env_release(para_env_col)
493
494 IF (my_open_shell) THEN
495 DEALLOCATE (iaia_ri_beta)
496 DEALLOCATE (d_ia_beta)
497 DEALLOCATE (m_ia_beta)
498 CALL mp_para_env_release(para_env_col_beta)
499 END IF
500
501 CALL timestop(handle)
502
503 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 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 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 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 col_indices=col_indices)
551
552 ! allocate the local buffer of iaia_RI integrals (dp kind)
553 ALLOCATE (iaia_ri_dp(ncol_local))
554 iaia_ri_dp = 0.0_dp
555
556 ! 2) perform the local multiplication SUM_K (ia|K)*(ia|K)
557 DO iib = 1, ncol_local
558 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 color_col = fm_mat_s%matrix_struct%context%mepos(2)
578 ALLOCATE (para_env_col)
579 CALL para_env_col%from_split(para_env_rpa, color_col)
580
581 CALL para_env_col%sum(iaia_ri_dp)
582
583 ! convert the iaia_RI_dp into double-double precision
584 ALLOCATE (iaia_ri(ncol_local))
585 DO iib = 1, ncol_local
586 iaia_ri(iib) = iaia_ri_dp(iib)
587 END DO
588 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 ALLOCATE (d_ia(ncol_local))
594
595 ALLOCATE (m_ia(ncol_local))
596
597 DO iib = 1, ncol_local
598 i_global = col_indices(iib)
599
600 iocc = max(1, i_global - 1)/virtual + 1
601 avirt = i_global - (iocc - 1)*virtual
602 eigen_diff = eigenval(avirt + homo) - eigenval(iocc)
603
604 d_ia(iib) = eigen_diff
605 END DO
606
607 DO iib = 1, ncol_local
608 m_ia(iib) = d_ia(iib)*d_ia(iib) + 2.0_dp*d_ia(iib)*iaia_ri(iib)
609 END DO
610
611 right_term_ref = 0.0_dp
612 DO iib = 1, ncol_local
613 right_term_ref = right_term_ref + (sqrt(m_ia(iib)) - d_ia(iib) - iaia_ri(iib))
614 END DO
615 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 color_row = fm_mat_s%matrix_struct%context%mepos(1)
619 ALLOCATE (para_env_row)
620 CALL para_env_row%from_split(para_env_rpa, color_row)
621
622 ! allocate communication array for rows
623 CALL para_env_row%sum(right_term_ref)
624
625 CALL mp_para_env_release(para_env_row)
626
627 CALL timestop(handle)
628
629 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 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 left_term = 0.0_dp
672 first_deriv = 0.0_dp
673 left_term_beta = 0.0_dp
674 first_deriv_beta = 0.0_dp
675 DO jquad = 1, num_integ_points
676 ! parallelize over integration points
677 IF (modulo(jquad, num_integ_group) /= color_rpa_group) cycle
678 omega = a_scaling*cottj(jquad)
679
680 DO iib = 1, ncol_local
681 ! parallelize over ia elements in the para_env_row group
682 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 (m_ia(iib) - d_ia(iib)**2)/(omega**2 + d_ia(iib)**2))
687 first_deriv = first_deriv + wj(jquad)*cottj(jquad)**2* &
688 ((-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 IF (my_open_shell) THEN
692 DO iib = 1, ncol_local_beta
693 ! parallelize over ia elements in the para_env_row group
694 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 (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 ((-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 CALL para_env%sum(left_term)
709 CALL para_env%sum(first_deriv)
710
711 IF (my_open_shell) THEN
712 CALL para_env%sum(left_term_beta)
713 CALL para_env%sum(first_deriv_beta)
714
715 left_term = left_term + left_term_beta
716 first_deriv = first_deriv + first_deriv_beta
717 END IF
718
719 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
754 REAL(kind=dp) :: multiplicator, omega
755 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_uty, work, &
756 x_values, y_values
757 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_a, mat_sinvvsinvsigma, &
758 mat_sinvvsinvt, mat_u
759
760 CALL timeset(routinen, handle)
761
762 ! take num_points_per_magnitude points per 10-interval
763 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 num_x_nodes = max(num_x_nodes, num_integ_points)
768
769 ALLOCATE (x_values(num_x_nodes))
770 x_values = 0.0_dp
771 ALLOCATE (y_values(num_x_nodes))
772 y_values = 0.0_dp
773 ALLOCATE (mat_a(num_x_nodes, num_integ_points))
774 mat_a = 0.0_dp
775 ALLOCATE (tau_wj_work(num_integ_points))
776 tau_wj_work = 0.0_dp
777 ALLOCATE (sing_values(num_integ_points))
778 sing_values = 0.0_dp
779 ALLOCATE (mat_u(num_x_nodes, num_x_nodes))
780 mat_u = 0.0_dp
781 ALLOCATE (mat_sinvvsinvt(num_x_nodes, num_integ_points))
782
783 mat_sinvvsinvt = 0.0_dp
784 ! double the value nessary for 'A' to achieve good performance
785 lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
786 ALLOCATE (work(lwork))
787 work = 0.0_dp
788 ALLOCATE (iwork(8*num_integ_points))
789 iwork = 0
790 ALLOCATE (mat_sinvvsinvsigma(num_integ_points, num_x_nodes))
791 mat_sinvvsinvsigma = 0.0_dp
792 ALLOCATE (vec_uty(num_x_nodes))
793 vec_uty = 0.0_dp
794
795 max_error = 0.0_dp
796
797 ! loop over all omega frequency points
798 DO jquad = 1, num_integ_points
799
800 ! set the x-values logarithmically in the interval [Emin,Emax]
801 multiplicator = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
802 DO iii = 1, num_x_nodes
803 x_values(iii) = e_min*multiplicator**(iii - 1)
804 END DO
805
806 omega = omega_tj(jquad)
807
808 ! y=2x/(x^2+omega_k^2)
809 DO iii = 1, num_x_nodes
810 y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2 + omega**2)
811 END DO
812
813 ! calculate mat_A
814 DO jjj = 1, num_integ_points
815 DO iii = 1, num_x_nodes
816 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 mat_sinvvsinvt, num_x_nodes, work, lwork, iwork, info)
823
824 cpassert(info == 0)
825
826 ! integration weights = V Sigma U^T y
827 ! 1) V*Sigma
828 DO jjj = 1, num_integ_points
829 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 /(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 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 num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
843
844 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 y_values, num_integ_points, num_x_nodes)
848
849 END DO ! jquad
850
851 DEALLOCATE (x_values, y_values, mat_a, tau_wj_work, sing_values, mat_u, mat_sinvvsinvt, &
852 work, iwork, mat_sinvvsinvsigma, vec_uty)
853
854 CALL timestop(handle)
855
856 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
890 REAL(kind=dp) :: chi2_min_jquad, multiplicator, omega
891 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_uty, work, &
892 work_array, x_values, y_values
893 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_a, mat_sinvvsinvsigma, &
894 mat_sinvvsinvt, mat_u
895
896 CALL timeset(routinen, handle)
897
898 ! take num_points_per_magnitude points per 10-interval
899 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 num_x_nodes = max(num_x_nodes, num_integ_points)
904
905 ALLOCATE (x_values(num_x_nodes))
906 x_values = 0.0_dp
907 ALLOCATE (y_values(num_x_nodes))
908 y_values = 0.0_dp
909 ALLOCATE (mat_a(num_x_nodes, num_integ_points))
910 mat_a = 0.0_dp
911 ALLOCATE (tau_wj_work(num_integ_points))
912 tau_wj_work = 0.0_dp
913 ALLOCATE (work_array(2*num_integ_points))
914 work_array = 0.0_dp
915 ALLOCATE (sing_values(num_integ_points))
916 sing_values = 0.0_dp
917 ALLOCATE (mat_u(num_x_nodes, num_x_nodes))
918 mat_u = 0.0_dp
919 ALLOCATE (mat_sinvvsinvt(num_x_nodes, num_integ_points))
920
921 mat_sinvvsinvt = 0.0_dp
922 ! double the value nessary for 'A' to achieve good performance
923 lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
924 ALLOCATE (work(lwork))
925 work = 0.0_dp
926 ALLOCATE (iwork(8*num_integ_points))
927 iwork = 0
928 ALLOCATE (mat_sinvvsinvsigma(num_integ_points, num_x_nodes))
929 mat_sinvvsinvsigma = 0.0_dp
930 ALLOCATE (vec_uty(num_x_nodes))
931 vec_uty = 0.0_dp
932
933 max_error = 0.0_dp
934
935 ! loop over all omega frequency points
936 DO jquad = 1, num_integ_points
937
938 chi2_min_jquad = 100.0_dp
939
940 ! set the x-values logarithmically in the interval [Emin,Emax]
941 multiplicator = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
942 DO iii = 1, num_x_nodes
943 x_values(iii) = e_min*multiplicator**(iii - 1)
944 END DO
945
946 omega = omega_tj(jquad)
947
948 ! y=2x/(x^2+omega_k^2)
949 DO iii = 1, num_x_nodes
950! y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2+omega**2)
951 y_values(iii) = 2.0_dp*omega/((x_values(iii))**2 + omega**2)
952 END DO
953
954 ! calculate mat_A
955 DO jjj = 1, num_integ_points
956 DO iii = 1, num_x_nodes
957 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 mat_sinvvsinvt, num_x_nodes, work, lwork, iwork, info)
964
965 cpassert(info == 0)
966
967 ! integration weights = V Sigma U^T y
968 ! 1) V*Sigma
969 DO jjj = 1, num_integ_points
970 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 /(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 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 num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
984
985 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 y_values, num_integ_points, num_x_nodes)
989
990 END DO ! jquad
991
992 DEALLOCATE (x_values, y_values, mat_a, tau_wj_work, work_array, sing_values, mat_u, mat_sinvvsinvt, &
993 work, iwork, mat_sinvvsinvsigma, vec_uty)
994
995 CALL timestop(handle)
996
997 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 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 max_error_tmp = 0.0_dp
1022
1023 DO kkk = 1, num_x_nodes
1024
1025 func_val = 0.0_dp
1026
1027 CALL eval_fit_func_tau_grid_cosine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
1028
1029 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 IF (max_error_tmp > max_error) THEN
1037
1038 max_error = max_error_tmp
1039
1040 END IF
1041
1042 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 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 func_val = 0.0_dp
1065
1066 DO iii = 1, num_integ_points
1067
1068 ! calculate value of the fit function
1069 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 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 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 func_val = 0.0_dp
1096
1097 DO iii = 1, num_integ_points
1098
1099 ! calculate value of the fit function
1100 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 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 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 max_error_tmp = 0.0_dp
1129
1130 DO kkk = 1, num_x_nodes
1131
1132 func_val = 0.0_dp
1133
1134 CALL eval_fit_func_tau_grid_sine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
1135
1136 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 IF (max_error_tmp > max_error) THEN
1144
1145 max_error = max_error_tmp
1146
1147 END IF
1148
1149 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tau_tj, tau_wj, tj, wj, x_tw
1163 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: weights_cos_tf_t_to_w
1164
1165 rc_max = 1.0e+7
1166
1167 multiplicator = rc_max**(1.0_dp/(real(nr, kind=dp) - 1.0_dp))
1168
1169 DO num_integ_points = 1, 20
1170
1171 ALLOCATE (x_tw(2*num_integ_points))
1172 x_tw = 0.0_dp
1173 ALLOCATE (tau_tj(num_integ_points))
1174 tau_tj = 0.0_dp
1175 ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
1176 weights_cos_tf_t_to_w = 0.0_dp
1177 ALLOCATE (tau_wj(num_integ_points))
1178 tau_wj = 0.0_dp
1179 ALLOCATE (tj(num_integ_points))
1180 tj = 0.0_dp
1181 ALLOCATE (wj(num_integ_points))
1182 wj = 0.0_dp
1183
1184 DO ir = 0, nr - 1
1185
1186 rc = 2.0_dp*multiplicator**ir
1187
1188 ierr = 0
1189 CALL get_rpa_minimax_coeff(num_integ_points, rc, x_tw, ierr, print_warning=.false.)
1190
1191 DO jquad = 1, num_integ_points
1192 tj(jquad) = x_tw(jquad)
1193 wj(jquad) = x_tw(jquad + num_integ_points)
1194 END DO
1195
1196 x_tw = 0.0_dp
1197
1198 CALL get_exp_minimax_coeff(num_integ_points, rc, x_tw)
1199
1200 DO jquad = 1, num_integ_points
1201 tau_tj(jquad) = x_tw(jquad)/2.0_dp
1202 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 1.0_dp, rc, max_error, 200, 0.0_dp)
1208
1209 IF (iw > 0) THEN
1210 WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, rc, max_error
1211 END IF
1212
1213 END DO
1214
1215 DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tau_wj, wj, tj)
1216
1217 END DO
1218
1219 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1253 REAL(kind=dp) :: chi2_min_jquad, multiplicator, omega, &
1254 tau, x_value
1255 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: omega_wj_work, sing_values, vec_uty, &
1256 work, work_array, x_values, y_values
1257 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_a, mat_sinvvsinvsigma, &
1258 mat_sinvvsinvt, mat_u
1259
1260 CALL timeset(routinen, handle)
1261
1262 ! take num_points_per_magnitude points per 10-interval
1263 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 num_x_nodes = max(num_x_nodes, num_integ_points)
1268
1269 ALLOCATE (x_values(num_x_nodes))
1270 x_values = 0.0_dp
1271 ALLOCATE (y_values(num_x_nodes))
1272 y_values = 0.0_dp
1273 ALLOCATE (mat_a(num_x_nodes, num_integ_points))
1274 mat_a = 0.0_dp
1275 ALLOCATE (omega_wj_work(num_integ_points))
1276 omega_wj_work = 0.0_dp
1277 ALLOCATE (work_array(2*num_integ_points))
1278 work_array = 0.0_dp
1279 ALLOCATE (sing_values(num_integ_points))
1280 sing_values = 0.0_dp
1281 ALLOCATE (mat_u(num_x_nodes, num_x_nodes))
1282 mat_u = 0.0_dp
1283 ALLOCATE (mat_sinvvsinvt(num_x_nodes, num_integ_points))
1284
1285 mat_sinvvsinvt = 0.0_dp
1286 ! double the value nessary for 'A' to achieve good performance
1287 lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
1288 ALLOCATE (work(lwork))
1289 work = 0.0_dp
1290 ALLOCATE (iwork(8*num_integ_points))
1291 iwork = 0
1292 ALLOCATE (mat_sinvvsinvsigma(num_integ_points, num_x_nodes))
1293 mat_sinvvsinvsigma = 0.0_dp
1294 ALLOCATE (vec_uty(num_x_nodes))
1295 vec_uty = 0.0_dp
1296
1297 ! set the x-values logarithmically in the interval [Emin,Emax]
1298 multiplicator = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
1299 DO iii = 1, num_x_nodes
1300 x_values(iii) = e_min*multiplicator**(iii - 1)
1301 END DO
1302
1303 max_error = 0.0_dp
1304
1305 ! loop over all tau time points
1306 DO jquad = 1, num_integ_points
1307
1308 chi2_min_jquad = 100.0_dp
1309
1310 tau = tau_tj(jquad)
1311
1312 ! y=exp(-x*|tau_k|)
1313 DO iii = 1, num_x_nodes
1314 y_values(iii) = exp(-x_values(iii)*tau)
1315 END DO
1316
1317 ! calculate mat_A
1318 DO jjj = 1, num_integ_points
1319 DO iii = 1, num_x_nodes
1320 omega = omega_tj(jjj)
1321 x_value = x_values(iii)
1322 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 mat_sinvvsinvt, num_x_nodes, work, lwork, iwork, info)
1329
1330 cpassert(info == 0)
1331
1332 ! integration weights = V Sigma U^T y
1333 ! 1) V*Sigma
1334 DO jjj = 1, num_integ_points
1335 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 /(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 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 num_x_nodes, 0.0_dp, omega_wj_work, num_integ_points)
1349
1350 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 y_values, num_integ_points, num_x_nodes)
1354
1355 END DO ! jquad
1356
1357 DEALLOCATE (x_values, y_values, mat_a, omega_wj_work, work_array, sing_values, mat_u, mat_sinvvsinvt, &
1358 work, iwork, mat_sinvvsinvsigma, vec_uty)
1359
1360 CALL timestop(handle)
1361
1362 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 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 CALL timeset(routinen, handle)
1391
1392 max_error_tmp = 0.0_dp
1393
1394 DO kkk = 1, num_x_nodes
1395
1396 func_val = 0.0_dp
1397
1398 CALL eval_fit_func_omega_grid_cosine(func_val, x_values(kkk), num_integ_points, omega_tj, omega_wj_work, tau)
1399
1400 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 IF (max_error_tmp > max_error) THEN
1408
1409 max_error = max_error_tmp
1410
1411 END IF
1412
1413 CALL timestop(handle)
1414
1415 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 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 func_val = 0.0_dp
1438
1439 DO iii = 1, num_integ_points
1440
1441 ! calculate value of the fit function
1442 omega = omega_tj(iii)
1443 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 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 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 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 CALL timeset(routinen, handle)
1476
1477 CALL get_qs_env(qs_env, &
1478 kpoints=kpoint)
1479
1480 mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1481 CALL get_mo_set(mo_set, nmo=nmo)
1482
1483 CALL get_kpoint_info(kpoint, kp_range=kp_range)
1484 kplocal = kp_range(2) - kp_range(1) + 1
1485
1486 gap = 1000.0_dp
1487 max_eig_diff = 0.0_dp
1488 e_homo = -1000.0_dp
1489 e_lumo = 1000.0_dp
1490
1491 DO ikpgr = 1, kplocal
1492 kp => kpoint%kp_env(ikpgr)%kpoint_env
1493 nspin = SIZE(kp%mos, 2)
1494 DO ispin = 1, nspin
1495 mo_set => kp%mos(1, ispin)
1496 CALL get_mo_set(mo_set, eigenvalues=eigenvalues, homo=homo)
1497 e_homo_temp = eigenvalues(homo)
1498 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 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 tmp(1) = e_homo
1510 tmp(2) = -e_lumo
1511 tmp(3) = max_eig_diff
1512 CALL para_env%max(tmp)
1513
1514 gap = -tmp(2) - tmp(1)
1515 e_fermi = (tmp(1) - tmp(2))*0.5_dp
1516 max_eig_diff = tmp(3)
1517
1518 CALL timestop(handle)
1519
1520 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 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 CALL timeset(routinen, handle)
1553 ! Test for spin unrestricted
1554 nspins = SIZE(homo)
1555
1556 ! Test whether all necessary variables are available
1557 my_do_kpoints = .false.
1558 IF (.NOT. do_ri_sos_laplace_mp2) THEN
1559 my_do_kpoints = do_kpoints_cubic_rpa
1560 END IF
1561
1562 IF (my_do_kpoints) THEN
1563 CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, emin, emax, e_fermi)
1564 e_range = emax/emin
1565 ELSE
1566 IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN
1567 emin = huge(dp)
1568 emax = 0.0_dp
1569 DO ispin = 1, nspins
1570 IF (homo(ispin) > 0) THEN
1571 emin = min(emin, eigenval(homo(ispin) + 1, 1, ispin) - eigenval(homo(ispin), 1, ispin))
1572 emax = max(emax, maxval(eigenval(:, :, ispin)) - minval(eigenval(:, :, ispin)))
1573 END IF
1574 END DO
1575 e_range = emax/emin
1576 qs_env%mp2_env%e_range = e_range
1577 qs_env%mp2_env%e_gap = emin
1578
1579 CALL get_qs_env(qs_env, input=input)
1580 CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=e_range)
1581 CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_GAP", r_val=emin)
1582 ELSE
1583 e_range = qs_env%mp2_env%E_range
1584 emin = qs_env%mp2_env%E_gap
1585 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 IF (do_ri_sos_laplace_mp2) THEN
1593 emin = emin*2.0_dp
1594 emax = emax*2.0_dp
1595 END IF
1596
1597 CALL timestop(handle)
1598 END SUBROUTINE determine_energy_range
1599
1600END MODULE mp2_grids
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
greenx_interface
Interface to the Greenx library.
Definition greenx_interface.F:13
greenx_interface::greenx_get_minimax_grid
subroutine, public greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, tau_tj, tau_wj, regularization_minimax, tj, wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, ierr)
...
Definition greenx_interface.F:319
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition input_section_types.F:1233
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:544
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
message_passing::mp_para_env_release
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Definition message_passing.F:2162
minimax_exp_gw
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition minimax_exp_gw.F:17
minimax_exp_gw::get_exp_minimax_coeff_gw
subroutine, public get_exp_minimax_coeff_gw(k, e_range, aw)
...
Definition minimax_exp_gw.F:38
minimax_exp
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition minimax_exp.F:29
minimax_exp::get_exp_minimax_coeff
subroutine, public get_exp_minimax_coeff(k, rc, aw, mm_error, which_coeffs)
Get best minimax approximation for given input parameters. Automatically chooses the most exact set o...
Definition minimax_exp.F:127
minimax_rpa
Routines to calculate the minimax coefficients for approximating 1/x as 1/x ~ 1/pi SUM_{i}^{K} w_i x^...
Definition minimax_rpa.F:14
minimax_rpa::get_rpa_minimax_coeff_larger_grid
subroutine, public get_rpa_minimax_coeff_larger_grid(k, e_range, aw)
...
Definition minimax_rpa.F:10119
minimax_rpa::get_rpa_minimax_coeff
subroutine, public get_rpa_minimax_coeff(k, e_range, aw, ierr, print_warning)
The a_i and w_i coefficient are stored in aw such that the first 1:K elements correspond to a_i and t...
Definition minimax_rpa.F:41
mp2_grids
Routines to calculate frequency and time grids (integration points and weights) for correlation metho...
Definition mp2_grids.F:14
mp2_grids::get_minimax_grid
subroutine, public get_minimax_grid(para_env, unit_nr, homo, eigenval, num_integ_points, do_im_time, do_ri_sos_laplace_mp2, do_print, tau_tj, tau_wj, qs_env, do_gw_im_time, do_kpoints_cubic_rpa, e_fermi, tj, wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, regularization)
...
Definition mp2_grids.F:76
mp2_grids::get_l_sq_wghts_cos_tf_w_to_t
subroutine, public get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, omega_tj, e_min, e_max, max_error, num_points_per_magnitude, regularization)
...
Definition mp2_grids.F:1235
mp2_grids::get_l_sq_wghts_cos_tf_t_to_w
subroutine, public get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, omega_tj, e_min, e_max, max_error, num_points_per_magnitude, regularization)
Calculate integration weights for the tau grid (in dependency of the omega node)
Definition mp2_grids.F:736
mp2_grids::get_l_sq_wghts_sin_tf_t_to_w
subroutine, public get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, omega_tj, e_min, e_max, max_error, num_points_per_magnitude, regularization)
Calculate integration weights for the tau grid (in dependency of the omega node)
Definition mp2_grids.F:872
mp2_grids::test_least_square_ft
subroutine, public test_least_square_ft(nr, iw)
test the singular value decomposition for the computation of integration weights for the Fourier tran...
Definition mp2_grids.F:1158
mp2_grids::get_clenshaw_grid
subroutine, public get_clenshaw_grid(para_env, para_env_rpa, unit_nr, homo, virtual, eigenval, num_integ_points, num_integ_group, color_rpa_group, fm_mat_s, my_do_gw, ext_scaling, a_scaling, tj, wj)
...
Definition mp2_grids.F:310
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:539
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_env_type
Keeps information about a specific k-point.
Definition kpoint_types.F:90
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:187
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47