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