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