Line data Source code
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 : ! **************************************************************************************************
14 : MODULE 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 :
46 : CONTAINS
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 126 : 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 126 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x_tw
104 : TYPE(section_vals_type), POINTER :: input
105 :
106 126 : CALL timeset(routineN, handle)
107 :
108 : ! Test for spin unrestricted
109 126 : nspins = SIZE(homo)
110 126 : my_open_shell = (nspins == 2)
111 :
112 : ! Test whether all necessary variables are available
113 126 : my_do_kpoints = .FALSE.
114 126 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
115 68 : my_do_kpoints = do_kpoints_cubic_RPA
116 : END IF
117 :
118 126 : my_regularization = 0.0_dp
119 126 : IF (PRESENT(regularization)) THEN
120 126 : my_regularization = regularization
121 : END IF
122 :
123 126 : IF (my_do_kpoints) THEN
124 0 : CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi)
125 0 : my_E_Range = Emax/Emin
126 : ELSE
127 126 : IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN
128 0 : Emin = HUGE(dp)
129 0 : Emax = 0.0_dp
130 0 : DO ispin = 1, nspins
131 0 : IF (homo(ispin) > 0) THEN
132 0 : Emin = MIN(Emin, Eigenval(homo(ispin) + 1, 1, ispin) - Eigenval(homo(ispin), 1, ispin))
133 0 : Emax = MAX(Emax, MAXVAL(Eigenval(:, :, ispin)) - MINVAL(Eigenval(:, :, ispin)))
134 : END IF
135 : END DO
136 0 : my_E_Range = Emax/Emin
137 0 : qs_env%mp2_env%e_range = my_e_range
138 0 : qs_env%mp2_env%e_gap = Emin
139 :
140 0 : CALL get_qs_env(qs_env, input=input)
141 0 : CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=my_e_range)
142 0 : CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_GAP", r_val=emin)
143 : ELSE
144 126 : my_E_range = qs_env%mp2_env%E_range
145 126 : Emin = qs_env%mp2_env%E_gap
146 126 : Emax = Emin*my_E_range
147 : END IF
148 : END IF
149 :
150 126 : IF (num_integ_points > 20 .AND. my_E_Range < 100.0_dp) THEN
151 0 : 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 0 : "a larger basis set with higher angular momenta or alternatively using all-electron calculations.")
157 : END IF
158 :
159 126 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
160 204 : ALLOCATE (x_tw(2*num_integ_points))
161 584 : x_tw = 0.0_dp
162 68 : ierr = 0
163 68 : IF (num_integ_points .LE. 20) THEN
164 68 : CALL get_rpa_minimax_coeff(num_integ_points, my_E_Range, x_tw, ierr)
165 : ELSE
166 0 : CALL get_rpa_minimax_coeff_larger_grid(num_integ_points, my_E_Range, x_tw)
167 : END IF
168 :
169 204 : ALLOCATE (tj(num_integ_points))
170 326 : tj = 0.0_dp
171 :
172 204 : ALLOCATE (wj(num_integ_points))
173 326 : wj = 0.0_dp
174 :
175 326 : DO jquad = 1, num_integ_points
176 258 : tj(jquad) = x_tw(jquad)
177 326 : 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 68 : IF (num_integ_points >= 26) THEN
182 0 : wj(:) = wj(:)*4.0_dp
183 : END IF
184 :
185 68 : DEALLOCATE (x_tw)
186 :
187 68 : IF (unit_nr > 0 .AND. do_print) THEN
188 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
189 33 : "MINIMAX_INFO| Number of integration points:", num_integ_points
190 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
191 33 : "MINIMAX_INFO| Gap for the minimax approximation:", Emin
192 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
193 33 : "MINIMAX_INFO| Range for the minimax approximation:", my_E_Range
194 33 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") "MINIMAX_INFO| Minimax parameters:", "Weights", "Abscissas"
195 157 : DO jquad = 1, num_integ_points
196 157 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(jquad), tj(jquad)
197 : END DO
198 33 : CALL m_flush(unit_nr)
199 : END IF
200 :
201 : ! scale the minimax parameters
202 326 : tj(:) = tj(:)*Emin
203 326 : 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 58 : Emin = Emin*2.0_dp
209 : END IF
210 :
211 : ! set up the minimax time grid
212 126 : IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
213 :
214 300 : ALLOCATE (x_tw(2*num_integ_points))
215 836 : x_tw = 0.0_dp
216 :
217 100 : IF (num_integ_points .LE. 20) THEN
218 100 : CALL get_exp_minimax_coeff(num_integ_points, my_E_Range, x_tw)
219 : ELSE
220 0 : 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 100 : scaling = 2.0_dp
225 100 : IF (do_ri_sos_laplace_mp2) scaling = 1.0_dp
226 :
227 300 : ALLOCATE (tau_tj(num_integ_points))
228 468 : tau_tj = 0.0_dp
229 :
230 300 : ALLOCATE (tau_wj(num_integ_points))
231 468 : tau_wj = 0.0_dp
232 :
233 468 : DO jquad = 1, num_integ_points
234 368 : tau_tj(jquad) = x_tw(jquad)/scaling
235 468 : tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling
236 : END DO
237 :
238 100 : DEALLOCATE (x_tw)
239 :
240 100 : IF (unit_nr > 0 .AND. do_print) THEN
241 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
242 49 : "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 49 : "MINIMAX_INFO| Gap:", Emin
246 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
247 49 : "MINIMAX_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas"
248 228 : DO jquad = 1, num_integ_points
249 228 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad)
250 : END DO
251 49 : CALL m_flush(unit_nr)
252 : END IF
253 :
254 : ! scale grid from [1,R] to [Emin,Emax]
255 468 : tau_tj(:) = tau_tj(:)/Emin
256 468 : tau_wj(:) = tau_wj(:)/Emin
257 :
258 100 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
259 168 : ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
260 678 : 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 42 : my_regularization)
265 :
266 : ! get the weights for the cosine transform W^c(iw) -> W^c(it)
267 168 : ALLOCATE (weights_cos_tf_w_to_t(num_integ_points, num_integ_points))
268 678 : 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 42 : my_regularization)
273 :
274 42 : 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 40 : ALLOCATE (weights_sin_tf_t_to_w(num_integ_points, num_integ_points))
278 310 : 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 10 : my_regularization)
283 :
284 10 : IF (unit_nr > 0) THEN
285 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
286 5 : "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 126 : CALL timestop(handle)
295 :
296 126 : 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 94 : SUBROUTINE get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, &
317 94 : 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 94 : CALL timeset(routineN, handle)
339 :
340 94 : nspins = SIZE(homo)
341 94 : my_open_shell = (nspins == 2)
342 :
343 : ! Now, start to prepare the different grid
344 282 : ALLOCATE (tj(num_integ_points))
345 3684 : tj = 0.0_dp
346 :
347 188 : ALLOCATE (wj(num_integ_points))
348 3684 : wj = 0.0_dp
349 :
350 3590 : DO jquad = 1, num_integ_points - 1
351 3496 : tj(jquad) = jquad*pi/(2.0_dp*num_integ_points)
352 3590 : wj(jquad) = pi/(num_integ_points*SIN(tj(jquad))**2)
353 : END DO
354 94 : tj(num_integ_points) = pi/2.0_dp
355 94 : wj(num_integ_points) = pi/(2.0_dp*num_integ_points*SIN(tj(num_integ_points))**2)
356 :
357 94 : IF (my_do_gw .AND. ext_scaling > 0.0_dp) THEN
358 58 : 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 36 : tj, wj, fm_mat_S)
363 : END IF
364 :
365 94 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.5)') 'INTEG_INFO| Scaling parameter:', a_scaling
366 :
367 3684 : wj(:) = wj(:)*a_scaling
368 :
369 94 : CALL timestop(handle)
370 :
371 94 : 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 36 : 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 36 : 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 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cottj, D_ia, D_ia_beta, iaia_RI, &
409 36 : iaia_RI_beta, M_ia, M_ia_beta
410 : TYPE(mp_para_env_type), POINTER :: para_env_col, para_env_col_beta
411 :
412 36 : CALL timeset(routineN, handle)
413 :
414 36 : nspins = SIZE(homo)
415 36 : my_open_shell = (nspins == 2)
416 :
417 36 : eps = 1.0E-10_dp
418 :
419 108 : ALLOCATE (cottj(num_integ_points))
420 :
421 : ! calculate the cotangent of the abscissa tj
422 406 : DO jquad = 1, num_integ_points
423 406 : 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 36 : 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 36 : 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 6 : D_ia_beta, iaia_RI_beta, M_ia_beta, fm_mat_S(2), para_env_col_beta)
433 :
434 6 : right_term_ref = right_term_ref + right_term_ref_beta
435 : END IF
436 :
437 : ! bcast the result
438 36 : IF (para_env%mepos == 0) THEN
439 18 : CALL para_env%bcast(right_term_ref, 0)
440 : ELSE
441 18 : right_term_ref = 0.0_dp
442 18 : 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 36 : conv_param = 100.0_dp*EPSILON(right_term_ref)
448 36 : step = 0.5_dp
449 36 : a_low = 0.0_dp
450 36 : a_high = step
451 36 : right_term = -right_term_ref
452 98 : DO icycle = 1, num_integ_points*2
453 90 : 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 90 : para_env, para_env_col, para_env_col_beta)
459 90 : left_term = left_term/4.0_dp/pi*a_scaling
460 :
461 90 : IF (ABS(left_term) > ABS(right_term) .OR. ABS(left_term + right_term) <= conv_param) EXIT
462 62 : a_low = a_high
463 98 : a_high = a_high + step
464 :
465 : END DO
466 :
467 36 : IF (ABS(left_term + right_term) >= conv_param) THEN
468 32 : IF (a_scaling >= 2*num_integ_points*step) THEN
469 10 : a_scaling = 1.0_dp
470 : ELSE
471 :
472 340 : DO icycle = 1, num_integ_points*2
473 336 : 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 336 : para_env, para_env_col, para_env_col_beta)
479 336 : left_term = left_term/4.0_dp/pi*a_scaling
480 :
481 336 : IF (ABS(left_term) > ABS(right_term)) THEN
482 : a_high = a_scaling
483 : ELSE
484 186 : a_low = a_scaling
485 : END IF
486 :
487 340 : IF (ABS(a_high - a_low) < 1.0e-5_dp) EXIT
488 :
489 : END DO
490 :
491 : END IF
492 : END IF
493 :
494 36 : a_scaling_ext = a_scaling
495 36 : CALL para_env%bcast(a_scaling_ext, 0)
496 :
497 36 : DEALLOCATE (cottj)
498 36 : DEALLOCATE (iaia_RI)
499 36 : DEALLOCATE (D_ia)
500 36 : DEALLOCATE (M_ia)
501 36 : CALL mp_para_env_release(para_env_col)
502 :
503 36 : IF (my_open_shell) THEN
504 6 : DEALLOCATE (iaia_RI_beta)
505 6 : DEALLOCATE (D_ia_beta)
506 6 : DEALLOCATE (M_ia_beta)
507 6 : CALL mp_para_env_release(para_env_col_beta)
508 : END IF
509 :
510 36 : CALL timestop(handle)
511 :
512 72 : 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 42 : 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 42 : 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 42 : 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 42 : col_indices=col_indices)
560 :
561 : ! allocate the local buffer of iaia_RI integrals (dp kind)
562 124 : ALLOCATE (iaia_RI_dp(ncol_local))
563 2722 : iaia_RI_dp = 0.0_dp
564 :
565 : ! 2) perform the local multiplication SUM_K (ia|K)*(ia|K)
566 2722 : DO iiB = 1, ncol_local
567 186900 : 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 42 : color_col = fm_mat_S%matrix_struct%context%mepos(2)
587 42 : ALLOCATE (para_env_col)
588 42 : CALL para_env_col%from_split(para_env_RPA, color_col)
589 :
590 42 : CALL para_env_col%sum(iaia_RI_dp)
591 :
592 : ! convert the iaia_RI_dp into double-double precision
593 124 : ALLOCATE (iaia_RI(ncol_local))
594 2722 : DO iiB = 1, ncol_local
595 2722 : iaia_RI(iiB) = iaia_RI_dp(iiB)
596 : END DO
597 42 : 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 124 : ALLOCATE (D_ia(ncol_local))
603 :
604 82 : ALLOCATE (M_ia(ncol_local))
605 :
606 2722 : DO iiB = 1, ncol_local
607 2680 : i_global = col_indices(iiB)
608 :
609 2680 : iocc = MAX(1, i_global - 1)/virtual + 1
610 2680 : avirt = i_global - (iocc - 1)*virtual
611 2680 : eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
612 :
613 2722 : D_ia(iiB) = eigen_diff
614 : END DO
615 :
616 2722 : DO iiB = 1, ncol_local
617 2722 : M_ia(iiB) = D_ia(iiB)*D_ia(iiB) + 2.0_dp*D_ia(iiB)*iaia_RI(iiB)
618 : END DO
619 :
620 42 : right_term_ref = 0.0_dp
621 2722 : DO iiB = 1, ncol_local
622 2722 : right_term_ref = right_term_ref + (SQRT(M_ia(iiB)) - D_ia(iiB) - iaia_RI(iiB))
623 : END DO
624 42 : 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 42 : color_row = fm_mat_S%matrix_struct%context%mepos(1)
628 42 : ALLOCATE (para_env_row)
629 42 : CALL para_env_row%from_split(para_env_RPA, color_row)
630 :
631 : ! allocate communication array for rows
632 42 : CALL para_env_row%sum(right_term_ref)
633 :
634 42 : CALL mp_para_env_release(para_env_row)
635 :
636 42 : CALL timestop(handle)
637 :
638 42 : 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 426 : 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 426 : left_term = 0.0_dp
681 426 : first_deriv = 0.0_dp
682 426 : left_term_beta = 0.0_dp
683 426 : first_deriv_beta = 0.0_dp
684 4370 : DO jquad = 1, num_integ_points
685 : ! parallelize over integration points
686 3944 : IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE
687 2172 : omega = a_scaling*cottj(jquad)
688 :
689 161964 : DO iiB = 1, ncol_local
690 : ! parallelize over ia elements in the para_env_row group
691 159792 : 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 144592 : (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2))
696 : first_deriv = first_deriv + wj(jquad)*cottj(jquad)**2* &
697 161964 : ((-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 2598 : IF (my_open_shell) THEN
701 14170 : DO iiB = 1, ncol_local_beta
702 : ! parallelize over ia elements in the para_env_row group
703 13860 : 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 13860 : (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 14170 : ((-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 426 : CALL para_env%sum(left_term)
718 426 : CALL para_env%sum(first_deriv)
719 :
720 426 : IF (my_open_shell) THEN
721 68 : CALL para_env%sum(left_term_beta)
722 68 : CALL para_env%sum(first_deriv_beta)
723 :
724 68 : left_term = left_term + left_term_beta
725 68 : first_deriv = first_deriv + first_deriv_beta
726 : END IF
727 :
728 426 : 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 70 : 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 70 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
763 : REAL(KIND=dp) :: multiplicator, omega
764 70 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_UTy, work, &
765 : x_values, y_values
766 70 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, &
767 70 : mat_SinvVSinvT, mat_U
768 :
769 70 : CALL timeset(routineN, handle)
770 :
771 : ! take num_points_per_magnitude points per 10-interval
772 70 : 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 70 : num_x_nodes = MAX(num_x_nodes, num_integ_points)
777 :
778 210 : ALLOCATE (x_values(num_x_nodes))
779 24070 : x_values = 0.0_dp
780 140 : ALLOCATE (y_values(num_x_nodes))
781 24070 : y_values = 0.0_dp
782 280 : ALLOCATE (mat_A(num_x_nodes, num_integ_points))
783 150976 : mat_A = 0.0_dp
784 210 : ALLOCATE (tau_wj_work(num_integ_points))
785 576 : tau_wj_work = 0.0_dp
786 140 : ALLOCATE (sing_values(num_integ_points))
787 576 : sing_values = 0.0_dp
788 280 : ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
789 9144070 : mat_U = 0.0_dp
790 210 : ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
791 :
792 150976 : mat_SinvVSinvT = 0.0_dp
793 : ! double the value nessary for 'A' to achieve good performance
794 70 : lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
795 210 : ALLOCATE (work(lwork))
796 105230 : work = 0.0_dp
797 210 : ALLOCATE (iwork(8*num_integ_points))
798 4118 : iwork = 0
799 210 : ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
800 174470 : mat_SinvVSinvSigma = 0.0_dp
801 140 : ALLOCATE (vec_UTy(num_x_nodes))
802 24070 : vec_UTy = 0.0_dp
803 :
804 70 : max_error = 0.0_dp
805 :
806 : ! loop over all omega frequency points
807 576 : DO jquad = 1, num_integ_points
808 :
809 : ! set the x-values logarithmically in the interval [Emin,Emax]
810 506 : multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
811 150906 : DO iii = 1, num_x_nodes
812 150906 : x_values(iii) = E_min*multiplicator**(iii - 1)
813 : END DO
814 :
815 506 : omega = omega_tj(jquad)
816 :
817 : ! y=2x/(x^2+omega_k^2)
818 150906 : DO iii = 1, num_x_nodes
819 150906 : y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2 + omega**2)
820 : END DO
821 :
822 : ! calculate mat_A
823 6892 : DO jjj = 1, num_integ_points
824 1590892 : DO iii = 1, num_x_nodes
825 1590386 : 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 506 : mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
832 :
833 506 : CPASSERT(info == 0)
834 :
835 : ! integration weights = V Sigma U^T y
836 : ! 1) V*Sigma
837 6892 : DO jjj = 1, num_integ_points
838 114582 : 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 114076 : /(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 506 : 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 506 : num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
852 :
853 6892 : 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 576 : y_values, num_integ_points, num_x_nodes)
857 :
858 : END DO ! jquad
859 :
860 0 : DEALLOCATE (x_values, y_values, mat_A, tau_wj_work, sing_values, mat_U, mat_SinvVSinvT, &
861 70 : work, iwork, mat_SinvVSinvSigma, vec_UTy)
862 :
863 70 : CALL timestop(handle)
864 :
865 70 : 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 38 : 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 38 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
899 : REAL(KIND=dp) :: chi2_min_jquad, multiplicator, omega
900 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_UTy, work, &
901 38 : work_array, x_values, y_values
902 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, &
903 38 : mat_SinvVSinvT, mat_U
904 :
905 38 : CALL timeset(routineN, handle)
906 :
907 : ! take num_points_per_magnitude points per 10-interval
908 38 : 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 38 : num_x_nodes = MAX(num_x_nodes, num_integ_points)
913 :
914 114 : ALLOCATE (x_values(num_x_nodes))
915 12838 : x_values = 0.0_dp
916 76 : ALLOCATE (y_values(num_x_nodes))
917 12838 : y_values = 0.0_dp
918 152 : ALLOCATE (mat_A(num_x_nodes, num_integ_points))
919 119656 : mat_A = 0.0_dp
920 114 : ALLOCATE (tau_wj_work(num_integ_points))
921 456 : tau_wj_work = 0.0_dp
922 114 : ALLOCATE (work_array(2*num_integ_points))
923 874 : work_array = 0.0_dp
924 76 : ALLOCATE (sing_values(num_integ_points))
925 456 : sing_values = 0.0_dp
926 152 : ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
927 4972838 : mat_U = 0.0_dp
928 114 : ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
929 :
930 119656 : mat_SinvVSinvT = 0.0_dp
931 : ! double the value nessary for 'A' to achieve good performance
932 38 : lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
933 114 : ALLOCATE (work(lwork))
934 79758 : work = 0.0_dp
935 114 : ALLOCATE (iwork(8*num_integ_points))
936 3382 : iwork = 0
937 114 : ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
938 132038 : mat_SinvVSinvSigma = 0.0_dp
939 76 : ALLOCATE (vec_UTy(num_x_nodes))
940 12838 : vec_UTy = 0.0_dp
941 :
942 38 : max_error = 0.0_dp
943 :
944 : ! loop over all omega frequency points
945 456 : DO jquad = 1, num_integ_points
946 :
947 418 : chi2_min_jquad = 100.0_dp
948 :
949 : ! set the x-values logarithmically in the interval [Emin,Emax]
950 418 : multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
951 119618 : DO iii = 1, num_x_nodes
952 119618 : x_values(iii) = E_min*multiplicator**(iii - 1)
953 : END DO
954 :
955 418 : omega = omega_tj(jquad)
956 :
957 : ! y=2x/(x^2+omega_k^2)
958 119618 : DO iii = 1, num_x_nodes
959 : ! y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2+omega**2)
960 119618 : y_values(iii) = 2.0_dp*omega/((x_values(iii))**2 + omega**2)
961 : END DO
962 :
963 : ! calculate mat_A
964 6556 : DO jjj = 1, num_integ_points
965 1501756 : DO iii = 1, num_x_nodes
966 1501338 : 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 418 : mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
973 :
974 418 : CPASSERT(info == 0)
975 :
976 : ! integration weights = V Sigma U^T y
977 : ! 1) V*Sigma
978 6556 : DO jjj = 1, num_integ_points
979 113534 : 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 113116 : /(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 418 : 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 418 : num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
993 :
994 6556 : 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 456 : y_values, num_integ_points, num_x_nodes)
998 :
999 : END DO ! jquad
1000 :
1001 0 : DEALLOCATE (x_values, y_values, mat_A, tau_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, &
1002 38 : work, iwork, mat_SinvVSinvSigma, vec_UTy)
1003 :
1004 38 : CALL timestop(handle)
1005 :
1006 38 : 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 506 : 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 506 : max_error_tmp = 0.0_dp
1031 :
1032 150906 : DO kkk = 1, num_x_nodes
1033 :
1034 : func_val = 0.0_dp
1035 :
1036 150400 : CALL eval_fit_func_tau_grid_cosine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
1037 :
1038 150906 : IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
1039 2486 : max_error_tmp = ABS(y_values(kkk) - func_val)
1040 2486 : func_val_temp = func_val
1041 : END IF
1042 :
1043 : END DO
1044 :
1045 506 : IF (max_error_tmp > max_error) THEN
1046 :
1047 124 : max_error = max_error_tmp
1048 :
1049 : END IF
1050 :
1051 506 : 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 150400 : 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 150400 : func_val = 0.0_dp
1074 :
1075 1734400 : DO iii = 1, num_integ_points
1076 :
1077 : ! calculate value of the fit function
1078 1734400 : 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 150400 : 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 119200 : 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 119200 : func_val = 0.0_dp
1105 :
1106 1614400 : DO iii = 1, num_integ_points
1107 :
1108 : ! calculate value of the fit function
1109 1614400 : 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 119200 : 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 418 : 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 418 : max_error_tmp = 0.0_dp
1138 :
1139 119618 : DO kkk = 1, num_x_nodes
1140 :
1141 119200 : func_val = 0.0_dp
1142 :
1143 119200 : CALL eval_fit_func_tau_grid_sine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
1144 :
1145 119618 : IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
1146 2410 : max_error_tmp = ABS(y_values(kkk) - func_val)
1147 2410 : func_val_temp = func_val
1148 : END IF
1149 :
1150 : END DO
1151 :
1152 418 : IF (max_error_tmp > max_error) THEN
1153 :
1154 42 : max_error = max_error_tmp
1155 :
1156 : END IF
1157 :
1158 418 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tau_tj, tau_wj, tj, wj, x_tw
1172 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: weights_cos_tf_t_to_w
1173 :
1174 0 : Rc_max = 1.0E+7
1175 :
1176 0 : multiplicator = Rc_max**(1.0_dp/(REAL(nR, KIND=dp) - 1.0_dp))
1177 :
1178 0 : DO num_integ_points = 1, 20
1179 :
1180 0 : ALLOCATE (x_tw(2*num_integ_points))
1181 0 : x_tw = 0.0_dp
1182 0 : ALLOCATE (tau_tj(num_integ_points))
1183 0 : tau_tj = 0.0_dp
1184 0 : ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
1185 0 : weights_cos_tf_t_to_w = 0.0_dp
1186 0 : ALLOCATE (tau_wj(num_integ_points))
1187 0 : tau_wj = 0.0_dp
1188 0 : ALLOCATE (tj(num_integ_points))
1189 0 : tj = 0.0_dp
1190 0 : ALLOCATE (wj(num_integ_points))
1191 0 : wj = 0.0_dp
1192 :
1193 0 : DO iR = 0, nR - 1
1194 :
1195 0 : Rc = 2.0_dp*multiplicator**iR
1196 :
1197 0 : ierr = 0
1198 0 : CALL get_rpa_minimax_coeff(num_integ_points, Rc, x_tw, ierr, print_warning=.FALSE.)
1199 :
1200 0 : DO jquad = 1, num_integ_points
1201 0 : tj(jquad) = x_tw(jquad)
1202 0 : wj(jquad) = x_tw(jquad + num_integ_points)
1203 : END DO
1204 :
1205 0 : x_tw = 0.0_dp
1206 :
1207 0 : CALL get_exp_minimax_coeff(num_integ_points, Rc, x_tw)
1208 :
1209 0 : DO jquad = 1, num_integ_points
1210 0 : tau_tj(jquad) = x_tw(jquad)/2.0_dp
1211 0 : 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 0 : 1.0_dp, Rc, max_error, 200, 0.0_dp)
1217 :
1218 0 : IF (iw > 0) THEN
1219 0 : WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, Rc, max_error
1220 : END IF
1221 :
1222 : END DO
1223 :
1224 0 : DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tau_wj, wj, tj)
1225 :
1226 : END DO
1227 :
1228 0 : 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 70 : 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 70 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1262 : REAL(KIND=dp) :: chi2_min_jquad, multiplicator, omega, &
1263 : tau, x_value
1264 70 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: omega_wj_work, sing_values, vec_UTy, &
1265 70 : work, work_array, x_values, y_values
1266 70 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, &
1267 70 : mat_SinvVSinvT, mat_U
1268 :
1269 70 : CALL timeset(routineN, handle)
1270 :
1271 : ! take num_points_per_magnitude points per 10-interval
1272 70 : 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 70 : num_x_nodes = MAX(num_x_nodes, num_integ_points)
1277 :
1278 210 : ALLOCATE (x_values(num_x_nodes))
1279 24070 : x_values = 0.0_dp
1280 140 : ALLOCATE (y_values(num_x_nodes))
1281 24070 : y_values = 0.0_dp
1282 280 : ALLOCATE (mat_A(num_x_nodes, num_integ_points))
1283 150976 : mat_A = 0.0_dp
1284 210 : ALLOCATE (omega_wj_work(num_integ_points))
1285 576 : omega_wj_work = 0.0_dp
1286 210 : ALLOCATE (work_array(2*num_integ_points))
1287 1082 : work_array = 0.0_dp
1288 140 : ALLOCATE (sing_values(num_integ_points))
1289 576 : sing_values = 0.0_dp
1290 280 : ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
1291 9144070 : mat_U = 0.0_dp
1292 210 : ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
1293 :
1294 150976 : mat_SinvVSinvT = 0.0_dp
1295 : ! double the value nessary for 'A' to achieve good performance
1296 70 : lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
1297 210 : ALLOCATE (work(lwork))
1298 105230 : work = 0.0_dp
1299 210 : ALLOCATE (iwork(8*num_integ_points))
1300 4118 : iwork = 0
1301 210 : ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
1302 174470 : mat_SinvVSinvSigma = 0.0_dp
1303 140 : ALLOCATE (vec_UTy(num_x_nodes))
1304 24070 : vec_UTy = 0.0_dp
1305 :
1306 : ! set the x-values logarithmically in the interval [Emin,Emax]
1307 70 : multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
1308 24070 : DO iii = 1, num_x_nodes
1309 24070 : x_values(iii) = E_min*multiplicator**(iii - 1)
1310 : END DO
1311 :
1312 70 : max_error = 0.0_dp
1313 :
1314 : ! loop over all tau time points
1315 576 : DO jquad = 1, num_integ_points
1316 :
1317 506 : chi2_min_jquad = 100.0_dp
1318 :
1319 506 : tau = tau_tj(jquad)
1320 :
1321 : ! y=exp(-x*|tau_k|)
1322 150906 : DO iii = 1, num_x_nodes
1323 150906 : y_values(iii) = EXP(-x_values(iii)*tau)
1324 : END DO
1325 :
1326 : ! calculate mat_A
1327 6892 : DO jjj = 1, num_integ_points
1328 1590892 : DO iii = 1, num_x_nodes
1329 1584000 : omega = omega_tj(jjj)
1330 1584000 : x_value = x_values(iii)
1331 1590386 : 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 506 : mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
1338 :
1339 506 : CPASSERT(info == 0)
1340 :
1341 : ! integration weights = V Sigma U^T y
1342 : ! 1) V*Sigma
1343 6892 : DO jjj = 1, num_integ_points
1344 114582 : 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 114076 : /(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 506 : 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 506 : num_x_nodes, 0.0_dp, omega_wj_work, num_integ_points)
1358 :
1359 6892 : 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 576 : y_values, num_integ_points, num_x_nodes)
1363 :
1364 : END DO ! jquad
1365 :
1366 0 : DEALLOCATE (x_values, y_values, mat_A, omega_wj_work, work_array, sing_values, mat_U, mat_SinvVSinvT, &
1367 70 : work, iwork, mat_SinvVSinvSigma, vec_UTy)
1368 :
1369 70 : CALL timestop(handle)
1370 :
1371 70 : 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 506 : 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 506 : CALL timeset(routineN, handle)
1400 :
1401 506 : max_error_tmp = 0.0_dp
1402 :
1403 150906 : DO kkk = 1, num_x_nodes
1404 :
1405 : func_val = 0.0_dp
1406 :
1407 150400 : CALL eval_fit_func_omega_grid_cosine(func_val, x_values(kkk), num_integ_points, omega_tj, omega_wj_work, tau)
1408 :
1409 150906 : IF (ABS(y_values(kkk) - func_val) > max_error_tmp) THEN
1410 4228 : max_error_tmp = ABS(y_values(kkk) - func_val)
1411 4228 : func_val_temp = func_val
1412 : END IF
1413 :
1414 : END DO
1415 :
1416 506 : IF (max_error_tmp > max_error) THEN
1417 :
1418 140 : max_error = max_error_tmp
1419 :
1420 : END IF
1421 :
1422 506 : CALL timestop(handle)
1423 :
1424 506 : 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 150400 : 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 150400 : func_val = 0.0_dp
1447 :
1448 1734400 : DO iii = 1, num_integ_points
1449 :
1450 : ! calculate value of the fit function
1451 1584000 : omega = omega_tj(iii)
1452 1734400 : 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 150400 : 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 8 : 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 4 : 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 4 : CALL timeset(routineN, handle)
1485 :
1486 : CALL get_qs_env(qs_env, &
1487 4 : kpoints=kpoint)
1488 :
1489 4 : mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1490 4 : CALL get_mo_set(mo_set, nmo=nmo)
1491 :
1492 4 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
1493 4 : kplocal = kp_range(2) - kp_range(1) + 1
1494 :
1495 4 : gap = 1000.0_dp
1496 4 : max_eig_diff = 0.0_dp
1497 4 : e_homo = -1000.0_dp
1498 4 : e_lumo = 1000.0_dp
1499 :
1500 12 : DO ikpgr = 1, kplocal
1501 8 : kp => kpoint%kp_env(ikpgr)%kpoint_env
1502 8 : nspin = SIZE(kp%mos, 2)
1503 20 : DO ispin = 1, nspin
1504 8 : mo_set => kp%mos(1, ispin)
1505 8 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, homo=homo)
1506 8 : e_homo_temp = eigenvalues(homo)
1507 8 : e_lumo_temp = eigenvalues(homo + 1)
1508 :
1509 8 : IF (e_homo_temp > e_homo) e_homo = e_homo_temp
1510 8 : IF (e_lumo_temp < e_lumo) e_lumo = e_lumo_temp
1511 24 : 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 4 : tmp(1) = e_homo
1519 4 : tmp(2) = -e_lumo
1520 4 : tmp(3) = max_eig_diff
1521 4 : CALL para_env%max(tmp)
1522 :
1523 4 : gap = -tmp(2) - tmp(1)
1524 4 : e_fermi = (tmp(1) - tmp(2))*0.5_dp
1525 4 : max_eig_diff = tmp(3)
1526 :
1527 4 : CALL timestop(handle)
1528 :
1529 4 : 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 200 : 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 200 : CALL timeset(routineN, handle)
1562 : ! Test for spin unrestricted
1563 200 : nspins = SIZE(homo)
1564 200 : my_open_shell = (nspins == 2)
1565 :
1566 : ! Test whether all necessary variables are available
1567 200 : my_do_kpoints = .FALSE.
1568 200 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1569 142 : my_do_kpoints = do_kpoints_cubic_RPA
1570 : END IF
1571 :
1572 142 : IF (my_do_kpoints) THEN
1573 4 : CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi)
1574 : my_E_Range = Emax/Emin
1575 : ELSE
1576 196 : IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN
1577 150 : Emin = HUGE(dp)
1578 150 : Emax = 0.0_dp
1579 336 : DO ispin = 1, nspins
1580 336 : IF (homo(ispin) > 0) THEN
1581 182 : Emin = MIN(Emin, Eigenval(homo(ispin) + 1, 1, ispin) - Eigenval(homo(ispin), 1, ispin))
1582 14640 : Emax = MAX(Emax, MAXVAL(Eigenval(:, :, ispin)) - MINVAL(Eigenval(:, :, ispin)))
1583 : END IF
1584 : END DO
1585 150 : my_E_Range = Emax/Emin
1586 150 : qs_env%mp2_env%e_range = my_e_range
1587 150 : qs_env%mp2_env%e_gap = Emin
1588 :
1589 150 : CALL get_qs_env(qs_env, input=input)
1590 150 : CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=my_e_range)
1591 150 : 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 46 : Emin = qs_env%mp2_env%E_gap
1595 46 : 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 200 : IF (do_ri_sos_laplace_mp2) THEN
1603 58 : Emin = Emin*2.0_dp
1604 : END IF
1605 :
1606 200 : CALL timestop(handle)
1607 200 : END SUBROUTINE
1608 :
1609 : END MODULE mp2_grids
|
