(git:1a29073)
mp2_grids.F
Go to the documentation of this file.
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  SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, &
72  do_im_time, do_ri_sos_laplace_mp2, do_print, tau_tj, tau_wj, qs_env, do_gw_im_time, &
73  do_kpoints_cubic_RPA, e_fermi, tj, wj, weights_cos_tf_t_to_w, &
74  weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, regularization)
75 
76  TYPE(mp_para_env_type), INTENT(IN) :: para_env
77  INTEGER, INTENT(IN) :: unit_nr
78  INTEGER, DIMENSION(:), INTENT(IN) :: homo
79  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: eigenval
80  INTEGER, INTENT(IN) :: num_integ_points
81  LOGICAL, INTENT(IN) :: do_im_time, do_ri_sos_laplace_mp2, &
82  do_print
83  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
84  INTENT(OUT) :: tau_tj, tau_wj
85  TYPE(qs_environment_type), POINTER :: qs_env
86  LOGICAL, INTENT(IN) :: do_gw_im_time, do_kpoints_cubic_rpa
87  REAL(kind=dp), INTENT(OUT) :: e_fermi
88  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
89  INTENT(INOUT) :: tj, wj
90  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
91  INTENT(OUT) :: weights_cos_tf_t_to_w, &
92  weights_cos_tf_w_to_t, &
93  weights_sin_tf_t_to_w
94  REAL(kind=dp), INTENT(IN), OPTIONAL :: regularization
95 
96  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_minimax_grid'
97  INTEGER, PARAMETER :: num_points_per_magnitude = 200
98 
99  INTEGER :: handle, ierr, ispin, jquad, nspins
100  LOGICAL :: my_do_kpoints, my_open_shell
101  REAL(kind=dp) :: emax, emin, max_error_min, my_e_range, &
102  my_regularization, scaling
103  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: x_tw
104  TYPE(section_vals_type), POINTER :: input
105 
106  CALL timeset(routinen, handle)
107 
108  ! Test for spin unrestricted
109  nspins = SIZE(homo)
110  my_open_shell = (nspins == 2)
111 
112  ! Test whether all necessary variables are available
113  my_do_kpoints = .false.
114  IF (.NOT. do_ri_sos_laplace_mp2) THEN
115  my_do_kpoints = do_kpoints_cubic_rpa
116  END IF
117 
118  my_regularization = 0.0_dp
119  IF (PRESENT(regularization)) THEN
120  my_regularization = regularization
121  END IF
122 
123  IF (my_do_kpoints) THEN
124  CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, emin, emax, e_fermi)
125  my_e_range = emax/emin
126  ELSE
127  IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN
128  emin = huge(dp)
129  emax = 0.0_dp
130  DO ispin = 1, nspins
131  IF (homo(ispin) > 0) THEN
132  emin = min(emin, eigenval(homo(ispin) + 1, 1, ispin) - eigenval(homo(ispin), 1, ispin))
133  emax = max(emax, maxval(eigenval(:, :, ispin)) - minval(eigenval(:, :, ispin)))
134  END IF
135  END DO
136  my_e_range = emax/emin
137  qs_env%mp2_env%e_range = my_e_range
138  qs_env%mp2_env%e_gap = emin
139 
140  CALL get_qs_env(qs_env, input=input)
141  CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=my_e_range)
142  CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_GAP", r_val=emin)
143  ELSE
144  my_e_range = qs_env%mp2_env%E_range
145  emin = qs_env%mp2_env%E_gap
146  emax = emin*my_e_range
147  END IF
148  END IF
149 
150  IF (num_integ_points > 20 .AND. my_e_range < 100.0_dp) THEN
151  IF (unit_nr > 0) &
152  CALL cp_warn(__location__, &
153  "You requested a large minimax grid (> 20 points) for a small minimax range R (R < 100). "// &
154  "That may lead to numerical "// &
155  "instabilities when computing minimax grid weights. You can prevent small ranges by choosing "// &
156  "a larger basis set with higher angular momenta or alternatively using all-electron calculations.")
157  END IF
158 
159  IF (.NOT. do_ri_sos_laplace_mp2) THEN
160  ALLOCATE (x_tw(2*num_integ_points))
161  x_tw = 0.0_dp
162  ierr = 0
163  IF (num_integ_points .LE. 20) THEN
164  CALL get_rpa_minimax_coeff(num_integ_points, my_e_range, x_tw, ierr)
165  ELSE
166  CALL get_rpa_minimax_coeff_larger_grid(num_integ_points, my_e_range, x_tw)
167  END IF
168 
169  ALLOCATE (tj(num_integ_points))
170  tj = 0.0_dp
171 
172  ALLOCATE (wj(num_integ_points))
173  wj = 0.0_dp
174 
175  DO jquad = 1, num_integ_points
176  tj(jquad) = x_tw(jquad)
177  wj(jquad) = x_tw(jquad + num_integ_points)
178  END DO
179 
180  DEALLOCATE (x_tw)
181 
182  IF (unit_nr > 0 .AND. do_print) THEN
183  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
184  "MINIMAX_INFO| Number of integration points:", num_integ_points
185  WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.4)") &
186  "MINIMAX_INFO| Gap for the minimax approximation:", emin
187  WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.4)") &
188  "MINIMAX_INFO| Range for the minimax approximation:", my_e_range
189  WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") "MINIMAX_INFO| Minimax parameters:", "Weights", "Abscissas"
190  DO jquad = 1, num_integ_points
191  WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") wj(jquad), tj(jquad)
192  END DO
193  CALL m_flush(unit_nr)
194  END IF
195 
196  ! scale the minimax parameters
197  tj(:) = tj(:)*emin
198  wj(:) = wj(:)*emin
199  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  emin = emin*2.0_dp
204  END IF
205 
206  ! set up the minimax time grid
207  IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
208 
209  ALLOCATE (x_tw(2*num_integ_points))
210  x_tw = 0.0_dp
211 
212  IF (num_integ_points .LE. 20) THEN
213  CALL get_exp_minimax_coeff(num_integ_points, my_e_range, x_tw)
214  ELSE
215  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  scaling = 2.0_dp
220  IF (do_ri_sos_laplace_mp2) scaling = 1.0_dp
221 
222  ALLOCATE (tau_tj(0:num_integ_points))
223  tau_tj = 0.0_dp
224 
225  ALLOCATE (tau_wj(num_integ_points))
226  tau_wj = 0.0_dp
227 
228  DO jquad = 1, num_integ_points
229  tau_tj(jquad) = x_tw(jquad)/scaling
230  tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling
231  END DO
232 
233  DEALLOCATE (x_tw)
234 
235  IF (unit_nr > 0 .AND. do_print) THEN
236  WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.4)") &
237  "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  "MINIMAX_INFO| Gap:", emin
241  WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
242  "MINIMAX_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas"
243  DO jquad = 1, num_integ_points
244  WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad)
245  END DO
246  CALL m_flush(unit_nr)
247  END IF
248 
249  ! scale grid from [1,R] to [Emin,Emax]
250  tau_tj(:) = tau_tj(:)/emin
251  tau_wj(:) = tau_wj(:)/emin
252 
253  IF (.NOT. do_ri_sos_laplace_mp2) THEN
254  ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
255  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  my_regularization)
260 
261  ! get the weights for the cosine transform W^c(iw) -> W^c(it)
262  ALLOCATE (weights_cos_tf_w_to_t(num_integ_points, num_integ_points))
263  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  my_regularization)
268 
269  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  ALLOCATE (weights_sin_tf_t_to_w(num_integ_points, num_integ_points))
273  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  my_regularization)
278 
279  IF (unit_nr > 0) THEN
280  WRITE (unit=unit_nr, fmt="(T3,A,T66,ES15.2)") &
281  "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  CALL timestop(handle)
290 
291  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  SUBROUTINE get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, &
312  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  CALL timeset(routinen, handle)
334 
335  nspins = SIZE(homo)
336  my_open_shell = (nspins == 2)
337 
338  ! Now, start to prepare the different grid
339  ALLOCATE (tj(num_integ_points))
340  tj = 0.0_dp
341 
342  ALLOCATE (wj(num_integ_points))
343  wj = 0.0_dp
344 
345  DO jquad = 1, num_integ_points - 1
346  tj(jquad) = jquad*pi/(2.0_dp*num_integ_points)
347  wj(jquad) = pi/(num_integ_points*sin(tj(jquad))**2)
348  END DO
349  tj(num_integ_points) = pi/2.0_dp
350  wj(num_integ_points) = pi/(2.0_dp*num_integ_points*sin(tj(num_integ_points))**2)
351 
352  IF (my_do_gw .AND. ext_scaling > 0.0_dp) THEN
353  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  tj, wj, fm_mat_s)
358  END IF
359 
360  IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.5)') 'INTEG_INFO| Scaling parameter:', a_scaling
361 
362  wj(:) = wj(:)*a_scaling
363 
364  CALL timestop(handle)
365 
366  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  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cottj, d_ia, d_ia_beta, iaia_ri, &
404  iaia_ri_beta, m_ia, m_ia_beta
405  TYPE(mp_para_env_type), POINTER :: para_env_col, para_env_col_beta
406 
407  CALL timeset(routinen, handle)
408 
409  nspins = SIZE(homo)
410  my_open_shell = (nspins == 2)
411 
412  eps = 1.0e-10_dp
413 
414  ALLOCATE (cottj(num_integ_points))
415 
416  ! calculate the cotangent of the abscissa tj
417  DO jquad = 1, num_integ_points
418  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  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  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  d_ia_beta, iaia_ri_beta, m_ia_beta, fm_mat_s(2), para_env_col_beta)
428 
429  right_term_ref = right_term_ref + right_term_ref_beta
430  END IF
431 
432  ! bcast the result
433  IF (para_env%mepos == 0) THEN
434  CALL para_env%bcast(right_term_ref, 0)
435  ELSE
436  right_term_ref = 0.0_dp
437  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  conv_param = 100.0_dp*epsilon(right_term_ref)
443  step = 0.5_dp
444  a_low = 0.0_dp
445  a_high = step
446  right_term = -right_term_ref
447  DO icycle = 1, num_integ_points*2
448  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  para_env, para_env_col, para_env_col_beta)
454  left_term = left_term/4.0_dp/pi*a_scaling
455 
456  IF (abs(left_term) > abs(right_term) .OR. abs(left_term + right_term) <= conv_param) EXIT
457  a_low = a_high
458  a_high = a_high + step
459 
460  END DO
461 
462  IF (abs(left_term + right_term) >= conv_param) THEN
463  IF (a_scaling >= 2*num_integ_points*step) THEN
464  a_scaling = 1.0_dp
465  ELSE
466 
467  DO icycle = 1, num_integ_points*2
468  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  para_env, para_env_col, para_env_col_beta)
474  left_term = left_term/4.0_dp/pi*a_scaling
475 
476  IF (abs(left_term) > abs(right_term)) THEN
477  a_high = a_scaling
478  ELSE
479  a_low = a_scaling
480  END IF
481 
482  IF (abs(a_high - a_low) < 1.0e-5_dp) EXIT
483 
484  END DO
485 
486  END IF
487  END IF
488 
489  a_scaling_ext = a_scaling
490  CALL para_env%bcast(a_scaling_ext, 0)
491 
492  DEALLOCATE (cottj)
493  DEALLOCATE (iaia_ri)
494  DEALLOCATE (d_ia)
495  DEALLOCATE (m_ia)
496  CALL mp_para_env_release(para_env_col)
497 
498  IF (my_open_shell) THEN
499  DEALLOCATE (iaia_ri_beta)
500  DEALLOCATE (d_ia_beta)
501  DEALLOCATE (m_ia_beta)
502  CALL mp_para_env_release(para_env_col_beta)
503  END IF
504 
505  CALL timestop(handle)
506 
507  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  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  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  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  col_indices=col_indices)
555 
556  ! allocate the local buffer of iaia_RI integrals (dp kind)
557  ALLOCATE (iaia_ri_dp(ncol_local))
558  iaia_ri_dp = 0.0_dp
559 
560  ! 2) perform the local multiplication SUM_K (ia|K)*(ia|K)
561  DO iib = 1, ncol_local
562  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  color_col = fm_mat_s%matrix_struct%context%mepos(2)
582  ALLOCATE (para_env_col)
583  CALL para_env_col%from_split(para_env_rpa, color_col)
584 
585  CALL para_env_col%sum(iaia_ri_dp)
586 
587  ! convert the iaia_RI_dp into double-double precision
588  ALLOCATE (iaia_ri(ncol_local))
589  DO iib = 1, ncol_local
590  iaia_ri(iib) = iaia_ri_dp(iib)
591  END DO
592  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  ALLOCATE (d_ia(ncol_local))
598 
599  ALLOCATE (m_ia(ncol_local))
600 
601  DO iib = 1, ncol_local
602  i_global = col_indices(iib)
603 
604  iocc = max(1, i_global - 1)/virtual + 1
605  avirt = i_global - (iocc - 1)*virtual
606  eigen_diff = eigenval(avirt + homo) - eigenval(iocc)
607 
608  d_ia(iib) = eigen_diff
609  END DO
610 
611  DO iib = 1, ncol_local
612  m_ia(iib) = d_ia(iib)*d_ia(iib) + 2.0_dp*d_ia(iib)*iaia_ri(iib)
613  END DO
614 
615  right_term_ref = 0.0_dp
616  DO iib = 1, ncol_local
617  right_term_ref = right_term_ref + (sqrt(m_ia(iib)) - d_ia(iib) - iaia_ri(iib))
618  END DO
619  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  color_row = fm_mat_s%matrix_struct%context%mepos(1)
623  ALLOCATE (para_env_row)
624  CALL para_env_row%from_split(para_env_rpa, color_row)
625 
626  ! allocate communication array for rows
627  CALL para_env_row%sum(right_term_ref)
628 
629  CALL mp_para_env_release(para_env_row)
630 
631  CALL timestop(handle)
632 
633  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  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  left_term = 0.0_dp
676  first_deriv = 0.0_dp
677  left_term_beta = 0.0_dp
678  first_deriv_beta = 0.0_dp
679  DO jquad = 1, num_integ_points
680  ! parallelize over integration points
681  IF (modulo(jquad, num_integ_group) /= color_rpa_group) cycle
682  omega = a_scaling*cottj(jquad)
683 
684  DO iib = 1, ncol_local
685  ! parallelize over ia elements in the para_env_row group
686  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  (m_ia(iib) - d_ia(iib)**2)/(omega**2 + d_ia(iib)**2))
691  first_deriv = first_deriv + wj(jquad)*cottj(jquad)**2* &
692  ((-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  IF (my_open_shell) THEN
696  DO iib = 1, ncol_local_beta
697  ! parallelize over ia elements in the para_env_row group
698  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  (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  ((-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  CALL para_env%sum(left_term)
713  CALL para_env%sum(first_deriv)
714 
715  IF (my_open_shell) THEN
716  CALL para_env%sum(left_term_beta)
717  CALL para_env%sum(first_deriv_beta)
718 
719  left_term = left_term + left_term_beta
720  first_deriv = first_deriv + first_deriv_beta
721  END IF
722 
723  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
758  REAL(kind=dp) :: multiplicator, omega
759  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_uty, work, &
760  x_values, y_values
761  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_a, mat_sinvvsinvsigma, &
762  mat_sinvvsinvt, mat_u
763 
764  CALL timeset(routinen, handle)
765 
766  ! take num_points_per_magnitude points per 10-interval
767  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  num_x_nodes = max(num_x_nodes, num_integ_points)
772 
773  ALLOCATE (x_values(num_x_nodes))
774  x_values = 0.0_dp
775  ALLOCATE (y_values(num_x_nodes))
776  y_values = 0.0_dp
777  ALLOCATE (mat_a(num_x_nodes, num_integ_points))
778  mat_a = 0.0_dp
779  ALLOCATE (tau_wj_work(num_integ_points))
780  tau_wj_work = 0.0_dp
781  ALLOCATE (sing_values(num_integ_points))
782  sing_values = 0.0_dp
783  ALLOCATE (mat_u(num_x_nodes, num_x_nodes))
784  mat_u = 0.0_dp
785  ALLOCATE (mat_sinvvsinvt(num_x_nodes, num_integ_points))
786 
787  mat_sinvvsinvt = 0.0_dp
788  ! double the value nessary for 'A' to achieve good performance
789  lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
790  ALLOCATE (work(lwork))
791  work = 0.0_dp
792  ALLOCATE (iwork(8*num_integ_points))
793  iwork = 0
794  ALLOCATE (mat_sinvvsinvsigma(num_integ_points, num_x_nodes))
795  mat_sinvvsinvsigma = 0.0_dp
796  ALLOCATE (vec_uty(num_x_nodes))
797  vec_uty = 0.0_dp
798 
799  max_error = 0.0_dp
800 
801  ! loop over all omega frequency points
802  DO jquad = 1, num_integ_points
803 
804  ! set the x-values logarithmically in the interval [Emin,Emax]
805  multiplicator = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
806  DO iii = 1, num_x_nodes
807  x_values(iii) = e_min*multiplicator**(iii - 1)
808  END DO
809 
810  omega = omega_tj(jquad)
811 
812  ! y=2x/(x^2+omega_k^2)
813  DO iii = 1, num_x_nodes
814  y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2 + omega**2)
815  END DO
816 
817  ! calculate mat_A
818  DO jjj = 1, num_integ_points
819  DO iii = 1, num_x_nodes
820  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  mat_sinvvsinvt, num_x_nodes, work, lwork, iwork, info)
827 
828  cpassert(info == 0)
829 
830  ! integration weights = V Sigma U^T y
831  ! 1) V*Sigma
832  DO jjj = 1, num_integ_points
833  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  /(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  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  num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
847 
848  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  y_values, num_integ_points, num_x_nodes)
852 
853  END DO ! jquad
854 
855  DEALLOCATE (x_values, y_values, mat_a, tau_wj_work, sing_values, mat_u, mat_sinvvsinvt, &
856  work, iwork, mat_sinvvsinvsigma, vec_uty)
857 
858  CALL timestop(handle)
859 
860  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
894  REAL(kind=dp) :: chi2_min_jquad, multiplicator, omega
895  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, tau_wj_work, vec_uty, work, &
896  work_array, x_values, y_values
897  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_a, mat_sinvvsinvsigma, &
898  mat_sinvvsinvt, mat_u
899 
900  CALL timeset(routinen, handle)
901 
902  ! take num_points_per_magnitude points per 10-interval
903  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  num_x_nodes = max(num_x_nodes, num_integ_points)
908 
909  ALLOCATE (x_values(num_x_nodes))
910  x_values = 0.0_dp
911  ALLOCATE (y_values(num_x_nodes))
912  y_values = 0.0_dp
913  ALLOCATE (mat_a(num_x_nodes, num_integ_points))
914  mat_a = 0.0_dp
915  ALLOCATE (tau_wj_work(num_integ_points))
916  tau_wj_work = 0.0_dp
917  ALLOCATE (work_array(2*num_integ_points))
918  work_array = 0.0_dp
919  ALLOCATE (sing_values(num_integ_points))
920  sing_values = 0.0_dp
921  ALLOCATE (mat_u(num_x_nodes, num_x_nodes))
922  mat_u = 0.0_dp
923  ALLOCATE (mat_sinvvsinvt(num_x_nodes, num_integ_points))
924 
925  mat_sinvvsinvt = 0.0_dp
926  ! double the value nessary for 'A' to achieve good performance
927  lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
928  ALLOCATE (work(lwork))
929  work = 0.0_dp
930  ALLOCATE (iwork(8*num_integ_points))
931  iwork = 0
932  ALLOCATE (mat_sinvvsinvsigma(num_integ_points, num_x_nodes))
933  mat_sinvvsinvsigma = 0.0_dp
934  ALLOCATE (vec_uty(num_x_nodes))
935  vec_uty = 0.0_dp
936 
937  max_error = 0.0_dp
938 
939  ! loop over all omega frequency points
940  DO jquad = 1, num_integ_points
941 
942  chi2_min_jquad = 100.0_dp
943 
944  ! set the x-values logarithmically in the interval [Emin,Emax]
945  multiplicator = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
946  DO iii = 1, num_x_nodes
947  x_values(iii) = e_min*multiplicator**(iii - 1)
948  END DO
949 
950  omega = omega_tj(jquad)
951 
952  ! y=2x/(x^2+omega_k^2)
953  DO iii = 1, num_x_nodes
954 ! y_values(iii) = 2.0_dp*x_values(iii)/((x_values(iii))**2+omega**2)
955  y_values(iii) = 2.0_dp*omega/((x_values(iii))**2 + omega**2)
956  END DO
957 
958  ! calculate mat_A
959  DO jjj = 1, num_integ_points
960  DO iii = 1, num_x_nodes
961  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  mat_sinvvsinvt, num_x_nodes, work, lwork, iwork, info)
968 
969  cpassert(info == 0)
970 
971  ! integration weights = V Sigma U^T y
972  ! 1) V*Sigma
973  DO jjj = 1, num_integ_points
974  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  /(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  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  num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
988 
989  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  y_values, num_integ_points, num_x_nodes)
993 
994  END DO ! jquad
995 
996  DEALLOCATE (x_values, y_values, mat_a, tau_wj_work, work_array, sing_values, mat_u, mat_sinvvsinvt, &
997  work, iwork, mat_sinvvsinvsigma, vec_uty)
998 
999  CALL timestop(handle)
1000 
1001  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  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  max_error_tmp = 0.0_dp
1026 
1027  DO kkk = 1, num_x_nodes
1028 
1029  func_val = 0.0_dp
1030 
1031  CALL eval_fit_func_tau_grid_cosine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
1032 
1033  IF (abs(y_values(kkk) - func_val) > max_error_tmp) THEN
1034  max_error_tmp = abs(y_values(kkk) - func_val)
1035  func_val_temp = func_val
1036  END IF
1037 
1038  END DO
1039 
1040  IF (max_error_tmp > max_error) THEN
1041 
1042  max_error = max_error_tmp
1043 
1044  END IF
1045 
1046  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  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  func_val = 0.0_dp
1069 
1070  DO iii = 1, num_integ_points
1071 
1072  ! calculate value of the fit function
1073  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  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  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  func_val = 0.0_dp
1100 
1101  DO iii = 1, num_integ_points
1102 
1103  ! calculate value of the fit function
1104  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  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  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  max_error_tmp = 0.0_dp
1133 
1134  DO kkk = 1, num_x_nodes
1135 
1136  func_val = 0.0_dp
1137 
1138  CALL eval_fit_func_tau_grid_sine(func_val, x_values(kkk), num_integ_points, tau_tj, tau_wj_work, omega)
1139 
1140  IF (abs(y_values(kkk) - func_val) > max_error_tmp) THEN
1141  max_error_tmp = abs(y_values(kkk) - func_val)
1142  func_val_temp = func_val
1143  END IF
1144 
1145  END DO
1146 
1147  IF (max_error_tmp > max_error) THEN
1148 
1149  max_error = max_error_tmp
1150 
1151  END IF
1152 
1153  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tau_tj, tau_wj, tj, wj, x_tw
1167  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: weights_cos_tf_t_to_w
1168 
1169  rc_max = 1.0e+7
1170 
1171  multiplicator = rc_max**(1.0_dp/(real(nr, kind=dp) - 1.0_dp))
1172 
1173  DO num_integ_points = 1, 20
1174 
1175  ALLOCATE (x_tw(2*num_integ_points))
1176  x_tw = 0.0_dp
1177  ALLOCATE (tau_tj(num_integ_points))
1178  tau_tj = 0.0_dp
1179  ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
1180  weights_cos_tf_t_to_w = 0.0_dp
1181  ALLOCATE (tau_wj(num_integ_points))
1182  tau_wj = 0.0_dp
1183  ALLOCATE (tj(num_integ_points))
1184  tj = 0.0_dp
1185  ALLOCATE (wj(num_integ_points))
1186  wj = 0.0_dp
1187 
1188  DO ir = 0, nr - 1
1189 
1190  rc = 2.0_dp*multiplicator**ir
1191 
1192  ierr = 0
1193  CALL get_rpa_minimax_coeff(num_integ_points, rc, x_tw, ierr, print_warning=.false.)
1194 
1195  DO jquad = 1, num_integ_points
1196  tj(jquad) = x_tw(jquad)
1197  wj(jquad) = x_tw(jquad + num_integ_points)
1198  END DO
1199 
1200  x_tw = 0.0_dp
1201 
1202  CALL get_exp_minimax_coeff(num_integ_points, rc, x_tw)
1203 
1204  DO jquad = 1, num_integ_points
1205  tau_tj(jquad) = x_tw(jquad)/2.0_dp
1206  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  1.0_dp, rc, max_error, 200, 0.0_dp)
1212 
1213  IF (iw > 0) THEN
1214  WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, rc, max_error
1215  END IF
1216 
1217  END DO
1218 
1219  DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tau_wj, wj, tj)
1220 
1221  END DO
1222 
1223  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
1257  REAL(kind=dp) :: chi2_min_jquad, multiplicator, omega, &
1258  tau, x_value
1259  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: omega_wj_work, sing_values, vec_uty, &
1260  work, work_array, x_values, y_values
1261  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_a, mat_sinvvsinvsigma, &
1262  mat_sinvvsinvt, mat_u
1263 
1264  CALL timeset(routinen, handle)
1265 
1266  ! take num_points_per_magnitude points per 10-interval
1267  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  num_x_nodes = max(num_x_nodes, num_integ_points)
1272 
1273  ALLOCATE (x_values(num_x_nodes))
1274  x_values = 0.0_dp
1275  ALLOCATE (y_values(num_x_nodes))
1276  y_values = 0.0_dp
1277  ALLOCATE (mat_a(num_x_nodes, num_integ_points))
1278  mat_a = 0.0_dp
1279  ALLOCATE (omega_wj_work(num_integ_points))
1280  omega_wj_work = 0.0_dp
1281  ALLOCATE (work_array(2*num_integ_points))
1282  work_array = 0.0_dp
1283  ALLOCATE (sing_values(num_integ_points))
1284  sing_values = 0.0_dp
1285  ALLOCATE (mat_u(num_x_nodes, num_x_nodes))
1286  mat_u = 0.0_dp
1287  ALLOCATE (mat_sinvvsinvt(num_x_nodes, num_integ_points))
1288 
1289  mat_sinvvsinvt = 0.0_dp
1290  ! double the value nessary for 'A' to achieve good performance
1291  lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
1292  ALLOCATE (work(lwork))
1293  work = 0.0_dp
1294  ALLOCATE (iwork(8*num_integ_points))
1295  iwork = 0
1296  ALLOCATE (mat_sinvvsinvsigma(num_integ_points, num_x_nodes))
1297  mat_sinvvsinvsigma = 0.0_dp
1298  ALLOCATE (vec_uty(num_x_nodes))
1299  vec_uty = 0.0_dp
1300 
1301  ! set the x-values logarithmically in the interval [Emin,Emax]
1302  multiplicator = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
1303  DO iii = 1, num_x_nodes
1304  x_values(iii) = e_min*multiplicator**(iii - 1)
1305  END DO
1306 
1307  max_error = 0.0_dp
1308 
1309  ! loop over all tau time points
1310  DO jquad = 1, num_integ_points
1311 
1312  chi2_min_jquad = 100.0_dp
1313 
1314  tau = tau_tj(jquad)
1315 
1316  ! y=exp(-x*|tau_k|)
1317  DO iii = 1, num_x_nodes
1318  y_values(iii) = exp(-x_values(iii)*tau)
1319  END DO
1320 
1321  ! calculate mat_A
1322  DO jjj = 1, num_integ_points
1323  DO iii = 1, num_x_nodes
1324  omega = omega_tj(jjj)
1325  x_value = x_values(iii)
1326  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  mat_sinvvsinvt, num_x_nodes, work, lwork, iwork, info)
1333 
1334  cpassert(info == 0)
1335 
1336  ! integration weights = V Sigma U^T y
1337  ! 1) V*Sigma
1338  DO jjj = 1, num_integ_points
1339  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  /(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  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  num_x_nodes, 0.0_dp, omega_wj_work, num_integ_points)
1353 
1354  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  y_values, num_integ_points, num_x_nodes)
1358 
1359  END DO ! jquad
1360 
1361  DEALLOCATE (x_values, y_values, mat_a, omega_wj_work, work_array, sing_values, mat_u, mat_sinvvsinvt, &
1362  work, iwork, mat_sinvvsinvsigma, vec_uty)
1363 
1364  CALL timestop(handle)
1365 
1366  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  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  CALL timeset(routinen, handle)
1395 
1396  max_error_tmp = 0.0_dp
1397 
1398  DO kkk = 1, num_x_nodes
1399 
1400  func_val = 0.0_dp
1401 
1402  CALL eval_fit_func_omega_grid_cosine(func_val, x_values(kkk), num_integ_points, omega_tj, omega_wj_work, tau)
1403 
1404  IF (abs(y_values(kkk) - func_val) > max_error_tmp) THEN
1405  max_error_tmp = abs(y_values(kkk) - func_val)
1406  func_val_temp = func_val
1407  END IF
1408 
1409  END DO
1410 
1411  IF (max_error_tmp > max_error) THEN
1412 
1413  max_error = max_error_tmp
1414 
1415  END IF
1416 
1417  CALL timestop(handle)
1418 
1419  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  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  func_val = 0.0_dp
1442 
1443  DO iii = 1, num_integ_points
1444 
1445  ! calculate value of the fit function
1446  omega = omega_tj(iii)
1447  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  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  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  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  CALL timeset(routinen, handle)
1480 
1481  CALL get_qs_env(qs_env, &
1482  kpoints=kpoint)
1483 
1484  mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
1485  CALL get_mo_set(mo_set, nmo=nmo)
1486 
1487  CALL get_kpoint_info(kpoint, kp_range=kp_range)
1488  kplocal = kp_range(2) - kp_range(1) + 1
1489 
1490  gap = 1000.0_dp
1491  max_eig_diff = 0.0_dp
1492  e_homo = -1000.0_dp
1493  e_lumo = 1000.0_dp
1494 
1495  DO ikpgr = 1, kplocal
1496  kp => kpoint%kp_env(ikpgr)%kpoint_env
1497  nspin = SIZE(kp%mos, 2)
1498  DO ispin = 1, nspin
1499  mo_set => kp%mos(1, ispin)
1500  CALL get_mo_set(mo_set, eigenvalues=eigenvalues, homo=homo)
1501  e_homo_temp = eigenvalues(homo)
1502  e_lumo_temp = eigenvalues(homo + 1)
1503 
1504  IF (e_homo_temp > e_homo) e_homo = e_homo_temp
1505  IF (e_lumo_temp < e_lumo) e_lumo = e_lumo_temp
1506  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  tmp(1) = e_homo
1514  tmp(2) = -e_lumo
1515  tmp(3) = max_eig_diff
1516  CALL para_env%max(tmp)
1517 
1518  gap = -tmp(2) - tmp(1)
1519  e_fermi = (tmp(1) - tmp(2))*0.5_dp
1520  max_eig_diff = tmp(3)
1521 
1522  CALL timestop(handle)
1523 
1524  END SUBROUTINE
1525 
1526 END MODULE mp2_grids
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition: grid_cpu_task_list.c:195
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition: input_section_types.F:1216
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
message_passing::mp_para_env_release
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Definition: message_passing.F:2027
minimax_exp_gw
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition: minimax_exp_gw.F:17
minimax_exp_gw::get_exp_minimax_coeff_gw
subroutine, public get_exp_minimax_coeff_gw(k, E_range, aw)
...
Definition: minimax_exp_gw.F:38
minimax_exp
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition: minimax_exp.F:29
minimax_exp::get_exp_minimax_coeff
subroutine, public get_exp_minimax_coeff(k, Rc, aw, mm_error, which_coeffs)
Get best minimax approximation for given input parameters. Automatically chooses the most exact set o...
Definition: minimax_exp.F:127
minimax_rpa
Routines to calculate the minimax coefficients for approximating 1/x as 1/x ~ 1/pi SUM_{i}^{K} w_i x^...
Definition: minimax_rpa.F:14
minimax_rpa::get_rpa_minimax_coeff_larger_grid
subroutine, public get_rpa_minimax_coeff_larger_grid(k, E_range, aw)
...
Definition: minimax_rpa.F:10122
minimax_rpa::get_rpa_minimax_coeff
subroutine, public get_rpa_minimax_coeff(k, E_range, aw, ierr, print_warning)
The a_i and w_i coefficient are stored in aw such that the first 1:K elements correspond to a_i and t...
Definition: minimax_rpa.F:41
mp2_grids
Routines to calculate frequency and time grids (integration points and weights) for correlation metho...
Definition: mp2_grids.F:14
mp2_grids::get_l_sq_wghts_sin_tf_t_to_w
subroutine, public get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, omega_tj, E_min, E_max, max_error, num_points_per_magnitude, regularization)
Calculate integration weights for the tau grid (in dependency of the omega node)
Definition: mp2_grids.F:876
mp2_grids::get_clenshaw_grid
subroutine, public get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, num_integ_group, color_rpa_group, fm_mat_S, my_do_gw, ext_scaling, a_scaling, tj, wj)
...
Definition: mp2_grids.F:314
mp2_grids::test_least_square_ft
subroutine, public test_least_square_ft(nR, iw)
test the singular value decomposition for the computation of integration weights for the Fourier tran...
Definition: mp2_grids.F:1162
mp2_grids::get_l_sq_wghts_cos_tf_w_to_t
subroutine, public get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, omega_tj, E_min, E_max, max_error, num_points_per_magnitude, regularization)
...
Definition: mp2_grids.F:1239
mp2_grids::get_minimax_grid
subroutine, public get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, do_ri_sos_laplace_mp2, do_print, tau_tj, tau_wj, qs_env, do_gw_im_time, do_kpoints_cubic_RPA, e_fermi, tj, wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, regularization)
...
Definition: mp2_grids.F:75
mp2_grids::get_l_sq_wghts_cos_tf_t_to_w
subroutine, public get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, omega_tj, E_min, E_max, max_error, num_points_per_magnitude, regularization)
Calculate integration weights for the tau grid (in dependency of the omega node)
Definition: mp2_grids.F:740
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397