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