(git:ed6f26b)
Loading...
Searching...
No Matches
rpa_main.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines to calculate RI-RPA energy
10!> \par History
11!> 06.2012 created [Mauro Del Ben]
12!> 04.2015 GW routines added [Jan Wilhelm]
13!> 10.2015 Cubic-scaling RPA routines added [Jan Wilhelm]
14!> 10.2018 Cubic-scaling SOS-MP2 added [Frederick Stein]
15!> 03.2019 Refactoring [Frederick Stein]
16! **************************************************************************************************
17MODULE rpa_main
18 USE bibliography, ONLY: &
19 bates2013, delben2013, delben2015, freeman1977, gruneis2009, ren2011, ren2013, &
20 wilhelm2016a, wilhelm2016b, wilhelm2017, wilhelm2018, cite_reference
21 USE bse_main, ONLY: start_bse_calculation
22 USE cp_blacs_env, ONLY: cp_blacs_env_create, &
23 cp_blacs_env_release, &
24 cp_blacs_env_type
25 USE cp_cfm_types, ONLY: cp_cfm_type
26 USE cp_dbcsr_api, ONLY: dbcsr_add, &
27 dbcsr_clear, &
28 dbcsr_get_info, &
29 dbcsr_p_type, &
30 dbcsr_type
31 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
32 USE cp_fm_struct, ONLY: cp_fm_struct_create, &
33 cp_fm_struct_release, &
34 cp_fm_struct_type
35 USE cp_fm_types, ONLY: cp_fm_create, &
36 cp_fm_release, &
37 cp_fm_set_all, &
38 cp_fm_to_fm, &
39 cp_fm_type
40 USE dbt_api, ONLY: dbt_type
41 USE dgemm_counter_types, ONLY: dgemm_counter_init, &
42 dgemm_counter_type, &
43 dgemm_counter_write
44 USE group_dist_types, ONLY: create_group_dist, &
45 get_group_dist, &
46 group_dist_d1_type, &
47 maxsize, &
48 release_group_dist
49 USE hfx_types, ONLY: block_ind_type, &
50 hfx_compression_type
51 USE input_constants, ONLY: rpa_exchange_axk, &
52 rpa_exchange_none, &
53 rpa_exchange_sosex, &
54 sigma_none, &
55 wfc_mm_style_gemm
56 USE kinds, ONLY: dp, &
57 int_8
58 USE kpoint_types, ONLY: kpoint_type
59 USE machine, ONLY: m_flush, &
60 m_memory
61 USE mathconstants, ONLY: pi, &
62 z_zero
63 USE message_passing, ONLY: mp_comm_type, &
64 mp_para_env_release, &
65 mp_para_env_type
66 USE minimax_exp, ONLY: check_exp_minimax_range
67 USE mp2_grids, ONLY: get_clenshaw_grid, &
68 get_minimax_grid, &
69 init_greenx_grids
70 USE mp2_laplace, ONLY: sos_mp2_postprocessing
71 USE mp2_ri_grad_util, ONLY: array2fm
72 USE mp2_types, ONLY: mp2_type, &
73 three_dim_real_array, &
74 two_dim_int_array, &
75 two_dim_real_array
76 USE qs_environment_types, ONLY: get_qs_env, &
77 qs_environment_type
78 USE rpa_exchange, ONLY: rpa_exchange_needed_mem, &
79 rpa_exchange_work_type
80 USE rpa_grad, ONLY: rpa_grad_copy_q, &
81 rpa_grad_create, &
82 rpa_grad_finalize, &
83 rpa_grad_matrix_operations, &
84 rpa_grad_needed_mem, &
85 rpa_grad_type
86 USE rpa_gw, ONLY: allocate_matrices_gw, &
87 allocate_matrices_gw_im_time, &
88 compute_gw_self_energy, &
89 compute_qp_energies, &
90 compute_w_cubic_gw, &
91 deallocate_matrices_gw, &
92 deallocate_matrices_gw_im_time, &
93 get_fermi_level_offset
94 USE rpa_gw_ic, ONLY: calculate_ic_correction
95 USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_mos, &
96 invert_eps_compute_w_and_erpa_kp
97 USE rpa_im_time, ONLY: compute_mat_p_omega, &
98 zero_mat_p_omega
99 USE rpa_im_time_force_methods, ONLY: calc_laplace_loop_forces, &
100 calc_post_loop_forces, &
101 calc_rpa_loop_forces, &
102 init_im_time_forces, &
103 keep_initial_quad
104 USE rpa_im_time_force_types, ONLY: im_time_force_release, &
105 im_time_force_type
106 USE rpa_sigma_functional, ONLY: finalize_rpa_sigma, &
107 rpa_sigma_create, &
108 rpa_sigma_matrix_spectral, &
109 rpa_sigma_type
110 USE rpa_util, ONLY: q_trace_and_add_unit_matrix, &
111 alloc_im_time, &
112 calc_mat_q, &
113 compute_erpa_by_freq_int, &
114 contract_p_omega_with_mat_l, &
115 dealloc_im_time, &
116 remove_scaling_factor_rpa
117 USE util, ONLY: get_limit
118#if defined (__GREENX)
119 USE gx_minimax, ONLY: gx_minimax_grid
120#endif
121
122#include "./base/base_uses.f90"
123
124 IMPLICIT NONE
125
126 PRIVATE
127
128 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_main'
129
130 PUBLIC :: rpa_ri_compute_en
131
132CONTAINS
133
134! **************************************************************************************************
135!> \brief ...
136!> \param qs_env ...
137!> \param Erpa ...
138!> \param mp2_env ...
139!> \param BIb_C ...
140!> \param BIb_C_gw ...
141!> \param BIb_C_bse_ij ...
142!> \param BIb_C_bse_ab ...
143!> \param para_env ...
144!> \param para_env_sub ...
145!> \param color_sub ...
146!> \param gd_array ...
147!> \param gd_B_virtual ...
148!> \param gd_B_all ...
149!> \param gd_B_occ_bse ...
150!> \param gd_B_virt_bse ...
151!> \param mo_coeff ...
152!> \param fm_matrix_PQ ...
153!> \param fm_matrix_L_kpoints ...
154!> \param fm_matrix_Minv_L_kpoints ...
155!> \param fm_matrix_Minv ...
156!> \param fm_matrix_Minv_Vtrunc_Minv ...
157!> \param kpoints ...
158!> \param Eigenval ...
159!> \param nmo ...
160!> \param homo ...
161!> \param dimen_RI ...
162!> \param dimen_RI_red ...
163!> \param gw_corr_lev_occ ...
164!> \param gw_corr_lev_virt ...
165!> \param bse_lev_virt ...
166!> \param unit_nr ...
167!> \param do_ri_sos_laplace_mp2 ...
168!> \param my_do_gw ...
169!> \param do_im_time ...
170!> \param do_bse ...
171!> \param matrix_s ...
172!> \param mat_munu ...
173!> \param mat_P_global ...
174!> \param t_3c_M ...
175!> \param t_3c_O ...
176!> \param t_3c_O_compressed ...
177!> \param t_3c_O_ind ...
178!> \param starts_array_mc ...
179!> \param ends_array_mc ...
180!> \param starts_array_mc_block ...
181!> \param ends_array_mc_block ...
182!> \param calc_forces ...
183! **************************************************************************************************
184 SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, &
185 para_env, para_env_sub, color_sub, &
186 gd_array, gd_B_virtual, gd_B_all, gd_B_occ_bse, gd_B_virt_bse, &
187 mo_coeff, fm_matrix_PQ, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
188 fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, kpoints, &
189 Eigenval, nmo, homo, dimen_RI, dimen_RI_red, gw_corr_lev_occ, gw_corr_lev_virt, &
190 bse_lev_virt, &
191 unit_nr, do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_bse, matrix_s, &
192 mat_munu, mat_P_global, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
193 starts_array_mc, ends_array_mc, &
194 starts_array_mc_block, ends_array_mc_block, calc_forces)
195
196 TYPE(qs_environment_type), POINTER :: qs_env
197 REAL(kind=dp), INTENT(OUT) :: erpa
198 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
199 TYPE(three_dim_real_array), DIMENSION(:), &
200 INTENT(INOUT) :: bib_c, bib_c_gw
201 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
202 INTENT(INOUT) :: bib_c_bse_ij, bib_c_bse_ab
203 TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
204 INTEGER, INTENT(INOUT) :: color_sub
205 TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_array
206 TYPE(group_dist_d1_type), DIMENSION(:), &
207 INTENT(INOUT) :: gd_b_virtual
208 TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_b_all, gd_b_occ_bse, gd_b_virt_bse
209 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
210 TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_pq
211 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_l_kpoints, &
212 fm_matrix_minv_l_kpoints, &
213 fm_matrix_minv, &
214 fm_matrix_minv_vtrunc_minv
215 TYPE(kpoint_type), POINTER :: kpoints
216 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
217 INTENT(INOUT) :: eigenval
218 INTEGER, INTENT(IN) :: nmo
219 INTEGER, DIMENSION(:), INTENT(IN) :: homo
220 INTEGER, INTENT(IN) :: dimen_ri, dimen_ri_red
221 INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt
222 INTEGER, INTENT(IN) :: bse_lev_virt, unit_nr
223 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, my_do_gw, &
224 do_im_time, do_bse
225 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
226 TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
227 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_p_global
228 TYPE(dbt_type) :: t_3c_m
229 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_o
230 TYPE(hfx_compression_type), ALLOCATABLE, &
231 DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_o_compressed
232 TYPE(block_ind_type), ALLOCATABLE, &
233 DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_o_ind
234 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
235 starts_array_mc_block, &
236 ends_array_mc_block
237 LOGICAL, INTENT(IN) :: calc_forces
238
239 CHARACTER(LEN=*), PARAMETER :: routinen = 'rpa_ri_compute_en'
240
241 INTEGER :: best_integ_group_size, best_num_integ_point, color_rpa_group, dimen_homo_square, &
242 dimen_nm_gw, dimen_virt_square, handle, handle2, handle3, ierr, iib, &
243 input_num_integ_groups, integ_group_size, ispin, jjb, min_integ_group_size, &
244 my_ab_comb_bse_end, my_ab_comb_bse_size, my_ab_comb_bse_start, my_group_l_end, &
245 my_group_l_size, my_group_l_start, my_ij_comb_bse_end, my_ij_comb_bse_size, &
246 my_ij_comb_bse_start, my_nm_gw_end, my_nm_gw_size, my_nm_gw_start, ncol_block_mat, &
247 ngroup, nrow_block_mat, nspins, num_integ_group, num_integ_points, pos_integ_group
248 INTEGER(KIND=int_8) :: mem
249 INTEGER, ALLOCATABLE, DIMENSION(:) :: dimen_ia, my_ia_end, my_ia_size, &
250 my_ia_start, virtual
251 LOGICAL :: do_kpoints_from_gamma, do_minimax_quad, &
252 my_open_shell, skip_integ_group_opt
253 REAL(kind=dp) :: allowed_memory, avail_mem, e_range, emax, emin, mem_for_iak, mem_for_qk, &
254 mem_min, mem_per_group, mem_per_rank, mem_per_repl, mem_real
255 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: eigenval_kp
256 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_q, fm_mat_q_gemm, fm_mat_s, &
257 fm_mat_s_gw
258 TYPE(cp_fm_type), DIMENSION(1) :: fm_mat_r_gw, fm_mat_s_ab_bse, &
259 fm_mat_s_ij_bse
260 TYPE(mp_para_env_type), POINTER :: para_env_rpa
261 TYPE(two_dim_real_array), ALLOCATABLE, &
262 DIMENSION(:) :: bib_c_2d, bib_c_2d_gw
263 TYPE(two_dim_real_array), DIMENSION(1) :: bib_c_2d_bse_ab, bib_c_2d_bse_ij
264
265 CALL timeset(routinen, handle)
266
267 CALL cite_reference(delben2013)
268 CALL cite_reference(delben2015)
269
270 IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_axk) THEN
271 CALL cite_reference(bates2013)
272 ELSE IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_sosex) THEN
273 CALL cite_reference(freeman1977)
274 CALL cite_reference(gruneis2009)
275 END IF
276 IF (mp2_env%ri_rpa%do_rse) THEN
277 CALL cite_reference(ren2011)
278 CALL cite_reference(ren2013)
279 END IF
280
281 IF (my_do_gw) THEN
282 CALL cite_reference(wilhelm2016a)
283 CALL cite_reference(wilhelm2017)
284 CALL cite_reference(wilhelm2018)
285 END IF
286
287 IF (do_im_time) THEN
288 CALL cite_reference(wilhelm2016b)
289 END IF
290
291 nspins = SIZE(homo)
292 my_open_shell = (nspins == 2)
293 ALLOCATE (virtual(nspins), dimen_ia(nspins), my_ia_end(nspins), my_ia_start(nspins), my_ia_size(nspins))
294 virtual(:) = nmo - homo(:)
295 dimen_ia(:) = virtual(:)*homo(:)
296
297 ALLOCATE (eigenval_kp(nmo, 1, nspins))
298 eigenval_kp(:, 1, :) = eigenval(:, :)
299
300 IF (do_im_time) mp2_env%ri_rpa%minimax_quad = .true.
301 do_minimax_quad = mp2_env%ri_rpa%minimax_quad
302
303 IF (do_ri_sos_laplace_mp2) THEN
304 num_integ_points = mp2_env%ri_laplace%n_quadrature
305 input_num_integ_groups = mp2_env%ri_laplace%num_integ_groups
306
307 ! check the range for the minimax approximation
308 e_range = mp2_env%e_range
309 IF (mp2_env%e_range <= 1.0_dp .OR. mp2_env%e_gap <= 0.0_dp) THEN
310 emin = huge(dp)
311 emax = 0.0_dp
312 DO ispin = 1, nspins
313 IF (homo(ispin) > 0) THEN
314 emin = min(emin, 2.0_dp*(eigenval(homo(ispin) + 1, ispin) - eigenval(homo(ispin), ispin)))
315 emax = max(emax, 2.0_dp*(maxval(eigenval(:, ispin)) - minval(eigenval(:, ispin))))
316 END IF
317 END DO
318 e_range = emax/emin
319 END IF
320 IF (e_range < 2.0_dp) e_range = 2.0_dp
321 ierr = 0
322 CALL check_exp_minimax_range(num_integ_points, e_range, ierr)
323 IF (ierr /= 0) THEN
324 jjb = num_integ_points - 1
325 DO iib = 1, jjb
326 num_integ_points = num_integ_points - 1
327 ierr = 0
328 CALL check_exp_minimax_range(num_integ_points, e_range, ierr)
329 IF (ierr == 0) EXIT
330 END DO
331 END IF
332 cpassert(num_integ_points >= 1)
333 ELSE
334 num_integ_points = mp2_env%ri_rpa%rpa_num_quad_points
335 input_num_integ_groups = mp2_env%ri_rpa%rpa_num_integ_groups
336 IF (my_do_gw .AND. do_minimax_quad) THEN
337 IF (num_integ_points > 34) THEN
338 IF (unit_nr > 0) &
339 CALL cp_warn(__location__, &
340 "The required number of quadrature point exceeds the maximum possible in the "// &
341 "Minimax quadrature scheme. The number of quadrature point has been reset to 30.")
342 num_integ_points = 30
343 END IF
344 ELSE
345 IF (do_minimax_quad .AND. num_integ_points > 20) THEN
346 IF (unit_nr > 0) &
347 CALL cp_warn(__location__, &
348 "The required number of quadrature point exceeds the maximum possible in the "// &
349 "Minimax quadrature scheme. The number of quadrature point has been reset to 20.")
350 num_integ_points = 20
351 END IF
352 END IF
353 END IF
354 allowed_memory = mp2_env%mp2_memory
355
356 CALL get_group_dist(gd_array, color_sub, my_group_l_start, my_group_l_end, my_group_l_size)
357
358 ngroup = para_env%num_pe/para_env_sub%num_pe
359
360 ! for imaginary time or periodic GW or BSE, we use all processors for a single frequency/time point
361 IF (do_im_time .OR. mp2_env%ri_g0w0%do_periodic .OR. do_bse) THEN
362
363 integ_group_size = ngroup
364 best_num_integ_point = num_integ_points
365
366 ELSE
367
368 ! Calculate available memory and create integral group according to that
369 ! mem_for_iaK is the memory needed for storing the 3 centre integrals
370 mem_for_iak = real(sum(dimen_ia), kind=dp)*dimen_ri_red*8.0_dp/(1024_dp**2)
371 mem_for_qk = real(dimen_ri_red, kind=dp)*nspins*dimen_ri_red*8.0_dp/(1024_dp**2)
372
373 CALL m_memory(mem)
374 mem_real = (mem + 1024*1024 - 1)/(1024*1024)
375 CALL para_env%min(mem_real)
376
377 mem_per_rank = 0.0_dp
378
379 ! B_ia_P
380 mem_per_repl = mem_for_iak
381 ! Q (regular and for dgemm)
382 mem_per_repl = mem_per_repl + 2.0_dp*mem_for_qk
383
384 IF (calc_forces) CALL rpa_grad_needed_mem(homo, virtual, dimen_ri_red, mem_per_rank, mem_per_repl, do_ri_sos_laplace_mp2)
385 CALL rpa_exchange_needed_mem(mp2_env, homo, virtual, dimen_ri_red, para_env, mem_per_rank, mem_per_repl)
386
387 mem_min = mem_per_repl/para_env%num_pe + mem_per_rank
388
389 IF (unit_nr > 0) THEN
390 WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum required memory per MPI process:', mem_min, ' MiB'
391 WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Available memory per MPI process:', mem_real, ' MiB'
392 END IF
393
394 ! Use only the allowed amount of memory
395 mem_real = min(mem_real, allowed_memory)
396 ! For the memory estimate, we require the amount of required memory per replication group and the available memory
397 mem_real = mem_real - mem_per_rank
398
399 mem_per_group = mem_real*para_env_sub%num_pe
400
401 ! here we try to find the best rpa/laplace group size
402 skip_integ_group_opt = .false.
403
404 ! Check the input number of integration groups
405 IF (input_num_integ_groups > 0) THEN
406 IF (mod(num_integ_points, input_num_integ_groups) == 0) THEN
407 IF (mod(ngroup, input_num_integ_groups) == 0) THEN
408 best_integ_group_size = ngroup/input_num_integ_groups
409 best_num_integ_point = num_integ_points/input_num_integ_groups
410 skip_integ_group_opt = .true.
411 ELSE
412 IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Total number of groups not multiple of NUM_INTEG_GROUPS'
413 END IF
414 ELSE
415 IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'NUM_QUAD_POINTS not multiple of NUM_INTEG_GROUPS'
416 END IF
417 END IF
418
419 IF (.NOT. skip_integ_group_opt) THEN
420 best_integ_group_size = ngroup
421 best_num_integ_point = num_integ_points
422
423 min_integ_group_size = max(1, ngroup/num_integ_points)
424
425 integ_group_size = min_integ_group_size - 1
426 DO iib = min_integ_group_size + 1, ngroup
427 integ_group_size = integ_group_size + 1
428
429 ! check that the ngroup is a multiple of integ_group_size
430 IF (mod(ngroup, integ_group_size) /= 0) cycle
431
432 ! check for memory
433 avail_mem = integ_group_size*mem_per_group
434 IF (avail_mem < mem_per_repl) cycle
435
436 ! check that the integration groups have the same size
437 num_integ_group = ngroup/integ_group_size
438 IF (mod(num_integ_points, num_integ_group) /= 0) cycle
439
440 best_num_integ_point = num_integ_points/num_integ_group
441 best_integ_group_size = integ_group_size
442
443 EXIT
444
445 END DO
446 END IF
447
448 integ_group_size = best_integ_group_size
449
450 END IF
451
452 IF (unit_nr > 0 .AND. .NOT. do_im_time) THEN
453 IF (do_ri_sos_laplace_mp2) THEN
454 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
455 "RI_INFO| Group size for laplace numerical integration:", integ_group_size*para_env_sub%num_pe
456 WRITE (unit=unit_nr, fmt="(T3,A)") &
457 "INTEG_INFO| MINIMAX approximation"
458 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
459 "INTEG_INFO| Number of integration points:", num_integ_points
460 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
461 "INTEG_INFO| Number of integration points per Laplace group:", best_num_integ_point
462 ELSE
463 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
464 "RI_INFO| Group size for frequency integration:", integ_group_size*para_env_sub%num_pe
465 IF (do_minimax_quad) THEN
466 WRITE (unit=unit_nr, fmt="(T3,A)") &
467 "INTEG_INFO| MINIMAX quadrature"
468 ELSE
469 WRITE (unit=unit_nr, fmt="(T3,A)") &
470 "INTEG_INFO| Clenshaw-Curtius quadrature"
471 END IF
472 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
473 "INTEG_INFO| Number of integration points:", num_integ_points
474 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
475 "INTEG_INFO| Number of integration points per RPA group:", best_num_integ_point
476 END IF
477 CALL m_flush(unit_nr)
478 END IF
479
480 num_integ_group = ngroup/integ_group_size
481
482 pos_integ_group = mod(color_sub, integ_group_size)
483 color_rpa_group = color_sub/integ_group_size
484
485 CALL timeset(routinen//"_reorder", handle2)
486
487 ! not necessary for imaginary time
488
489 ALLOCATE (bib_c_2d(nspins))
490
491 IF (.NOT. do_im_time) THEN
492
493 ! reorder the local data in such a way to help the next stage of matrix creation
494 ! now the data inside the group are divided into a ia x K matrix
495 DO ispin = 1, nspins
496 CALL calculate_bib_c_2d(bib_c_2d(ispin)%array, bib_c(ispin)%array, para_env_sub, dimen_ia(ispin), &
497 homo(ispin), virtual(ispin), gd_b_virtual(ispin), &
498 my_ia_size(ispin), my_ia_start(ispin), my_ia_end(ispin), my_group_l_size)
499
500 DEALLOCATE (bib_c(ispin)%array)
501 CALL release_group_dist(gd_b_virtual(ispin))
502
503 END DO
504
505 ! in the GW case, BIb_C_2D_gw is an nm x K matrix, with n: number of corr GW levels, m=nmo
506 IF (my_do_gw) THEN
507 ALLOCATE (bib_c_2d_gw(nspins))
508
509 CALL timeset(routinen//"_reorder_gw", handle3)
510
511 dimen_nm_gw = nmo*(gw_corr_lev_occ(1) + gw_corr_lev_virt(1))
512
513 ! The same for open shell
514 DO ispin = 1, nspins
515 CALL calculate_bib_c_2d(bib_c_2d_gw(ispin)%array, bib_c_gw(ispin)%array, para_env_sub, dimen_nm_gw, &
516 gw_corr_lev_occ(ispin) + gw_corr_lev_virt(ispin), nmo, gd_b_all, &
517 my_nm_gw_size, my_nm_gw_start, my_nm_gw_end, my_group_l_size)
518 DEALLOCATE (bib_c_gw(ispin)%array)
519 END DO
520
521 CALL release_group_dist(gd_b_all)
522
523 CALL timestop(handle3)
524
525 END IF
526 END IF
527
528 IF (do_bse) THEN
529
530 CALL timeset(routinen//"_reorder_bse1", handle3)
531
532 dimen_homo_square = homo(1)**2
533 ! We do not implement an explicit bse_lev_occ different to homo here, because the small number of occupied levels
534 ! does not critically influence the memory
535 CALL calculate_bib_c_2d(bib_c_2d_bse_ij(1)%array, bib_c_bse_ij, para_env_sub, dimen_homo_square, &
536 homo(1), homo(1), gd_b_occ_bse, &
537 my_ij_comb_bse_size, my_ij_comb_bse_start, my_ij_comb_bse_end, my_group_l_size)
538
539 DEALLOCATE (bib_c_bse_ij)
540 CALL release_group_dist(gd_b_occ_bse)
541
542 CALL timestop(handle3)
543
544 CALL timeset(routinen//"_reorder_bse2", handle3)
545
546 dimen_virt_square = bse_lev_virt**2
547
548 CALL calculate_bib_c_2d(bib_c_2d_bse_ab(1)%array, bib_c_bse_ab, para_env_sub, dimen_virt_square, &
549 bse_lev_virt, bse_lev_virt, gd_b_virt_bse, &
550 my_ab_comb_bse_size, my_ab_comb_bse_start, my_ab_comb_bse_end, my_group_l_size)
551
552 DEALLOCATE (bib_c_bse_ab)
553 CALL release_group_dist(gd_b_virt_bse)
554
555 CALL timestop(handle3)
556
557 END IF
558
559 CALL timestop(handle2)
560
561 IF (num_integ_group > 1) THEN
562 ALLOCATE (para_env_rpa)
563 CALL para_env_rpa%from_split(para_env, color_rpa_group)
564 ELSE
565 para_env_rpa => para_env
566 END IF
567
568 ! now create the matrices needed for the calculation, Q, S and G
569 ! Q and G will have omega dependence
570
571 IF (do_im_time) THEN
572 ALLOCATE (fm_mat_q(nspins), fm_mat_q_gemm(1), fm_mat_s(1))
573 ELSE
574 ALLOCATE (fm_mat_q(nspins), fm_mat_q_gemm(nspins), fm_mat_s(nspins))
575 END IF
576
577 CALL create_integ_mat(bib_c_2d, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
578 dimen_ri_red, dimen_ia, color_rpa_group, &
579 mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
580 my_ia_size, my_ia_start, my_ia_end, &
581 my_group_l_size, my_group_l_start, my_group_l_end, &
582 para_env_rpa, fm_mat_s, nrow_block_mat, ncol_block_mat, &
583 dimen_ia_for_block_size=dimen_ia(1), &
584 do_im_time=do_im_time, fm_mat_q_gemm=fm_mat_q_gemm, fm_mat_q=fm_mat_q, qs_env=qs_env)
585
586 DEALLOCATE (bib_c_2d, my_ia_end, my_ia_size, my_ia_start)
587
588 ! for GW, we need other matrix fm_mat_S, always allocate the container to prevent crying compilers
589 ALLOCATE (fm_mat_s_gw(nspins))
590 IF (my_do_gw .AND. .NOT. do_im_time) THEN
591
592 CALL create_integ_mat(bib_c_2d_gw, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
593 dimen_ri_red, [dimen_nm_gw, dimen_nm_gw], color_rpa_group, &
594 mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
595 [my_nm_gw_size, my_nm_gw_size], [my_nm_gw_start, my_nm_gw_start], [my_nm_gw_end, my_nm_gw_end], &
596 my_group_l_size, my_group_l_start, my_group_l_end, &
597 para_env_rpa, fm_mat_s_gw, nrow_block_mat, ncol_block_mat, &
598 fm_mat_q(1)%matrix_struct%context, fm_mat_q(1)%matrix_struct%context, &
599 fm_mat_q=fm_mat_r_gw)
600 DEALLOCATE (bib_c_2d_gw)
601
602 END IF
603
604 ! for Bethe-Salpeter, we need other matrix fm_mat_S
605 IF (do_bse) THEN
606 CALL create_integ_mat(bib_c_2d_bse_ij, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
607 dimen_ri_red, [dimen_homo_square], color_rpa_group, &
608 mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
609 [my_ij_comb_bse_size], [my_ij_comb_bse_start], [my_ij_comb_bse_end], &
610 my_group_l_size, my_group_l_start, my_group_l_end, &
611 para_env_rpa, fm_mat_s_ij_bse, nrow_block_mat, ncol_block_mat, &
612 fm_mat_q(1)%matrix_struct%context, fm_mat_q(1)%matrix_struct%context)
613
614 CALL create_integ_mat(bib_c_2d_bse_ab, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
615 dimen_ri_red, [dimen_virt_square], color_rpa_group, &
616 mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
617 [my_ab_comb_bse_size], [my_ab_comb_bse_start], [my_ab_comb_bse_end], &
618 my_group_l_size, my_group_l_start, my_group_l_end, &
619 para_env_rpa, fm_mat_s_ab_bse, nrow_block_mat, ncol_block_mat, &
620 fm_mat_q(1)%matrix_struct%context, fm_mat_q(1)%matrix_struct%context)
621
622 END IF
623
624 do_kpoints_from_gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
625 IF (do_kpoints_from_gamma) THEN
626 CALL get_bandstruc_and_k_dependent_mos(qs_env, eigenval_kp)
627 END IF
628
629 ! Now start the RPA calculation
630 ! fm_mo_coeff_occ, fm_mo_coeff_virt will be deallocated here
631 CALL rpa_num_int(qs_env, erpa, mp2_env, para_env, para_env_rpa, para_env_sub, unit_nr, &
632 homo, virtual, dimen_ri, dimen_ri_red, dimen_ia, dimen_nm_gw, &
633 eigenval_kp, num_integ_points, num_integ_group, color_rpa_group, &
634 fm_matrix_pq, fm_mat_s, fm_mat_q_gemm, fm_mat_q, fm_mat_s_gw, fm_mat_r_gw(1), &
635 fm_mat_s_ij_bse(1), fm_mat_s_ab_bse(1), &
636 my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
637 bse_lev_virt, &
638 do_minimax_quad, &
639 do_im_time, mo_coeff, &
640 fm_matrix_l_kpoints, fm_matrix_minv_l_kpoints, &
641 fm_matrix_minv, fm_matrix_minv_vtrunc_minv, mat_munu, mat_p_global, &
642 t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, &
643 starts_array_mc, ends_array_mc, &
644 starts_array_mc_block, ends_array_mc_block, &
645 matrix_s, do_kpoints_from_gamma, kpoints, gd_array, color_sub, &
646 do_ri_sos_laplace_mp2=do_ri_sos_laplace_mp2, calc_forces=calc_forces)
647
648 CALL release_group_dist(gd_array)
649
650 IF (num_integ_group > 1) CALL mp_para_env_release(para_env_rpa)
651
652 IF (.NOT. do_im_time) THEN
653 CALL cp_fm_release(fm_mat_q_gemm)
654 CALL cp_fm_release(fm_mat_s)
655 END IF
656 CALL cp_fm_release(fm_mat_q)
657
658 IF (my_do_gw .AND. .NOT. do_im_time) THEN
659 CALL cp_fm_release(fm_mat_s_gw)
660 CALL cp_fm_release(fm_mat_r_gw(1))
661 END IF
662
663 IF (do_bse) THEN
664 CALL cp_fm_release(fm_mat_s_ij_bse(1))
665 CALL cp_fm_release(fm_mat_s_ab_bse(1))
666 END IF
667
668 CALL timestop(handle)
669
670 END SUBROUTINE rpa_ri_compute_en
671
672! **************************************************************************************************
673!> \brief reorder the local data in such a way to help the next stage of matrix creation;
674!> now the data inside the group are divided into a ia x K matrix (BIb_C_2D);
675!> Subroutine created to avoid massive double coding
676!> \param BIb_C_2D ...
677!> \param BIb_C ...
678!> \param para_env_sub ...
679!> \param dimen_ia ...
680!> \param homo ...
681!> \param virtual ...
682!> \param gd_B_virtual ...
683!> \param my_ia_size ...
684!> \param my_ia_start ...
685!> \param my_ia_end ...
686!> \param my_group_L_size ...
687!> \author Jan Wilhelm, 03/2015
688! **************************************************************************************************
689 SUBROUTINE calculate_bib_c_2d(BIb_C_2D, BIb_C, para_env_sub, dimen_ia, homo, virtual, &
690 gd_B_virtual, &
691 my_ia_size, my_ia_start, my_ia_end, my_group_L_size)
692
693 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
694 INTENT(OUT) :: bib_c_2d
695 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
696 INTENT(IN) :: bib_c
697 TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
698 INTEGER, INTENT(IN) :: dimen_ia, homo, virtual
699 TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_b_virtual
700 INTEGER :: my_ia_size, my_ia_start, my_ia_end, &
701 my_group_l_size
702
703 INTEGER, PARAMETER :: occ_chunk = 128
704
705 INTEGER :: ia_global, iib, itmp(2), jjb, my_b_size, my_b_virtual_start, occ_high, occ_low, &
706 proc_receive, proc_send, proc_shift, rec_b_size, rec_b_virtual_end, rec_b_virtual_start
707 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), TARGET :: bib_c_rec_1d
708 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: bib_c_rec
709
710 itmp = get_limit(dimen_ia, para_env_sub%num_pe, para_env_sub%mepos)
711 my_ia_start = itmp(1)
712 my_ia_end = itmp(2)
713 my_ia_size = my_ia_end - my_ia_start + 1
714
715 CALL get_group_dist(gd_b_virtual, para_env_sub%mepos, sizes=my_b_size, starts=my_b_virtual_start)
716
717 ! reorder data
718 ALLOCATE (bib_c_2d(my_group_l_size, my_ia_size))
719
720!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,ia_global) &
721!$OMP SHARED(homo,my_B_size,virtual,my_B_virtual_start,my_ia_start,my_ia_end,BIb_C,BIb_C_2D,&
722!$OMP my_group_L_size)
723 DO iib = 1, homo
724 DO jjb = 1, my_b_size
725 ia_global = (iib - 1)*virtual + my_b_virtual_start + jjb - 1
726 IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
727 bib_c_2d(1:my_group_l_size, ia_global - my_ia_start + 1) = bib_c(1:my_group_l_size, jjb, iib)
728 END IF
729 END DO
730 END DO
731
732 IF (para_env_sub%num_pe > 1) THEN
733 ALLOCATE (bib_c_rec_1d(int(my_group_l_size, int_8)*maxsize(gd_b_virtual)*min(homo, occ_chunk)))
734 DO proc_shift = 1, para_env_sub%num_pe - 1
735 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
736 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
737
738 CALL get_group_dist(gd_b_virtual, proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
739
740 ! do this in chunks to avoid high memory overhead
741 DO occ_low = 1, homo, occ_chunk
742 occ_high = min(homo, occ_low + occ_chunk - 1)
743 bib_c_rec(1:my_group_l_size, 1:rec_b_size, 1:occ_high - occ_low + 1) => &
744 bib_c_rec_1d(1:int(my_group_l_size, int_8)*rec_b_size*(occ_high - occ_low + 1))
745 CALL para_env_sub%sendrecv(bib_c(:, :, occ_low:occ_high), proc_send, &
746 bib_c_rec(:, :, 1:occ_high - occ_low + 1), proc_receive)
747!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,ia_global) &
748!$OMP SHARED(occ_low,occ_high,rec_B_size,virtual,rec_B_virtual_start,my_ia_start,my_ia_end,BIb_C_rec,BIb_C_2D,&
749!$OMP my_group_L_size)
750 DO iib = occ_low, occ_high
751 DO jjb = 1, rec_b_size
752 ia_global = (iib - 1)*virtual + rec_b_virtual_start + jjb - 1
753 IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
754 bib_c_2d(1:my_group_l_size, ia_global - my_ia_start + 1) = bib_c_rec(1:my_group_l_size, jjb, iib - occ_low + 1)
755 END IF
756 END DO
757 END DO
758 END DO
759
760 END DO
761 DEALLOCATE (bib_c_rec_1d)
762 END IF
763
764 END SUBROUTINE calculate_bib_c_2d
765
766! **************************************************************************************************
767!> \brief ...
768!> \param BIb_C_2D ...
769!> \param para_env ...
770!> \param para_env_sub ...
771!> \param color_sub ...
772!> \param ngroup ...
773!> \param integ_group_size ...
774!> \param dimen_RI ...
775!> \param dimen_ia ...
776!> \param color_rpa_group ...
777!> \param ext_row_block_size ...
778!> \param ext_col_block_size ...
779!> \param unit_nr ...
780!> \param my_ia_size ...
781!> \param my_ia_start ...
782!> \param my_ia_end ...
783!> \param my_group_L_size ...
784!> \param my_group_L_start ...
785!> \param my_group_L_end ...
786!> \param para_env_RPA ...
787!> \param fm_mat_S ...
788!> \param nrow_block_mat ...
789!> \param ncol_block_mat ...
790!> \param blacs_env_ext ...
791!> \param blacs_env_ext_S ...
792!> \param dimen_ia_for_block_size ...
793!> \param do_im_time ...
794!> \param fm_mat_Q_gemm ...
795!> \param fm_mat_Q ...
796!> \param qs_env ...
797! **************************************************************************************************
798 SUBROUTINE create_integ_mat(BIb_C_2D, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
799 dimen_RI, dimen_ia, color_rpa_group, &
800 ext_row_block_size, ext_col_block_size, unit_nr, &
801 my_ia_size, my_ia_start, my_ia_end, &
802 my_group_L_size, my_group_L_start, my_group_L_end, &
803 para_env_RPA, fm_mat_S, nrow_block_mat, ncol_block_mat, &
804 blacs_env_ext, blacs_env_ext_S, dimen_ia_for_block_size, &
805 do_im_time, fm_mat_Q_gemm, fm_mat_Q, qs_env)
806
807 TYPE(two_dim_real_array), DIMENSION(:), &
808 INTENT(INOUT) :: bib_c_2d
809 TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_sub
810 INTEGER, INTENT(IN) :: color_sub, ngroup, integ_group_size, &
811 dimen_ri
812 INTEGER, DIMENSION(:), INTENT(IN) :: dimen_ia
813 INTEGER, INTENT(IN) :: color_rpa_group, ext_row_block_size, &
814 ext_col_block_size, unit_nr
815 INTEGER, DIMENSION(:), INTENT(IN) :: my_ia_size, my_ia_start, my_ia_end
816 INTEGER, INTENT(IN) :: my_group_l_size, my_group_l_start, &
817 my_group_l_end
818 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_rpa
819 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_s
820 INTEGER, INTENT(INOUT) :: nrow_block_mat, ncol_block_mat
821 TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: blacs_env_ext, blacs_env_ext_s
822 INTEGER, INTENT(IN), OPTIONAL :: dimen_ia_for_block_size
823 LOGICAL, INTENT(IN), OPTIONAL :: do_im_time
824 TYPE(cp_fm_type), DIMENSION(:), OPTIONAL :: fm_mat_q_gemm, fm_mat_q
825 TYPE(qs_environment_type), INTENT(IN), OPTIONAL, &
826 POINTER :: qs_env
827
828 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_integ_mat'
829
830 INTEGER :: col_row_proc_ratio, grid_2d(2), handle, &
831 iproc, iproc_col, iproc_row, ispin, &
832 mepos_in_rpa_group
833 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: group_grid_2_mepos
834 LOGICAL :: my_blacs_ext, my_blacs_s_ext, &
835 my_do_im_time
836 TYPE(cp_blacs_env_type), POINTER :: blacs_env, blacs_env_q
837 TYPE(cp_fm_struct_type), POINTER :: fm_struct
838 TYPE(group_dist_d1_type) :: gd_ia, gd_l
839
840 CALL timeset(routinen, handle)
841
842 cpassert(PRESENT(blacs_env_ext) .OR. PRESENT(dimen_ia_for_block_size))
843
844 my_blacs_ext = .false.
845 IF (PRESENT(blacs_env_ext)) my_blacs_ext = .true.
846
847 my_blacs_s_ext = .false.
848 IF (PRESENT(blacs_env_ext_s)) my_blacs_s_ext = .true.
849
850 my_do_im_time = .false.
851 IF (PRESENT(do_im_time)) my_do_im_time = do_im_time
852
853 NULLIFY (blacs_env)
854 ! create the RPA blacs env
855 IF (my_blacs_s_ext) THEN
856 blacs_env => blacs_env_ext_s
857 ELSE
858 IF (para_env_rpa%num_pe > 1) THEN
859 col_row_proc_ratio = max(1, dimen_ia_for_block_size/dimen_ri)
860
861 iproc_col = min(max(int(sqrt(real(para_env_rpa%num_pe*col_row_proc_ratio, kind=dp))), 1), para_env_rpa%num_pe) + 1
862 DO iproc = 1, para_env_rpa%num_pe
863 iproc_col = iproc_col - 1
864 IF (mod(para_env_rpa%num_pe, iproc_col) == 0) EXIT
865 END DO
866
867 iproc_row = para_env_rpa%num_pe/iproc_col
868 grid_2d(1) = iproc_row
869 grid_2d(2) = iproc_col
870 ELSE
871 grid_2d = 1
872 END IF
873 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_rpa, grid_2d=grid_2d)
874
875 IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
876 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
877 "MATRIX_INFO| Number row processes:", grid_2d(1)
878 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
879 "MATRIX_INFO| Number column processes:", grid_2d(2)
880 END IF
881
882 ! define the block_size for the row
883 IF (ext_row_block_size > 0) THEN
884 nrow_block_mat = ext_row_block_size
885 ELSE
886 nrow_block_mat = max(1, dimen_ri/grid_2d(1)/2)
887 END IF
888
889 ! define the block_size for the column
890 IF (ext_col_block_size > 0) THEN
891 ncol_block_mat = ext_col_block_size
892 ELSE
893 ncol_block_mat = max(1, dimen_ia_for_block_size/grid_2d(2)/2)
894 END IF
895
896 IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
897 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
898 "MATRIX_INFO| Row block size:", nrow_block_mat
899 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
900 "MATRIX_INFO| Column block size:", ncol_block_mat
901 END IF
902 END IF
903
904 IF (.NOT. my_do_im_time) THEN
905 DO ispin = 1, SIZE(bib_c_2d)
906 NULLIFY (fm_struct)
907 IF (my_blacs_ext) THEN
908 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_ri, &
909 ncol_global=dimen_ia(ispin), para_env=para_env_rpa)
910 ELSE
911 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_ri, &
912 ncol_global=dimen_ia(ispin), para_env=para_env_rpa, &
913 nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.true.)
914
915 END IF ! external blacs_env
916
917 CALL create_group_dist(gd_ia, my_ia_start(ispin), my_ia_end(ispin), my_ia_size(ispin), para_env_rpa)
918 CALL create_group_dist(gd_l, my_group_l_start, my_group_l_end, my_group_l_size, para_env_rpa)
919
920 ! create the info array
921
922 mepos_in_rpa_group = mod(color_sub, integ_group_size)
923 ALLOCATE (group_grid_2_mepos(0:integ_group_size - 1, 0:para_env_sub%num_pe - 1))
924 group_grid_2_mepos = 0
925 group_grid_2_mepos(mepos_in_rpa_group, para_env_sub%mepos) = para_env_rpa%mepos
926 CALL para_env_rpa%sum(group_grid_2_mepos)
927
928 CALL array2fm(bib_c_2d(ispin)%array, fm_struct, my_group_l_start, my_group_l_end, &
929 my_ia_start(ispin), my_ia_end(ispin), gd_l, gd_ia, &
930 group_grid_2_mepos, ngroup, para_env_sub%num_pe, fm_mat_s(ispin), &
931 integ_group_size, color_rpa_group)
932
933 DEALLOCATE (group_grid_2_mepos)
934 CALL cp_fm_struct_release(fm_struct)
935
936 ! deallocate the info array
937 CALL release_group_dist(gd_l)
938 CALL release_group_dist(gd_ia)
939
940 ! sum the local data across processes belonging to different RPA group.
941 IF (para_env_rpa%num_pe /= para_env%num_pe) THEN
942 block
943 TYPE(mp_comm_type) :: comm_exchange
944 comm_exchange = fm_mat_s(ispin)%matrix_struct%context%interconnect(para_env)
945 CALL comm_exchange%sum(fm_mat_s(ispin)%local_data)
946 CALL comm_exchange%free()
947 END block
948 END IF
949 END DO
950 END IF
951
952 IF (PRESENT(fm_mat_q_gemm) .AND. .NOT. my_do_im_time) THEN
953 ! create the Q matrix dimen_RIxdimen_RI where the result of the mat-mat-mult will be stored
954 NULLIFY (fm_struct)
955 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_ri, &
956 ncol_global=dimen_ri, para_env=para_env_rpa, &
957 nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.true.)
958 DO ispin = 1, SIZE(fm_mat_q_gemm)
959 CALL cp_fm_create(fm_mat_q_gemm(ispin), fm_struct, name="fm_mat_Q_gemm")
960 END DO
961 CALL cp_fm_struct_release(fm_struct)
962 END IF
963
964 IF (PRESENT(fm_mat_q)) THEN
965 NULLIFY (blacs_env_q)
966 IF (my_blacs_ext) THEN
967 blacs_env_q => blacs_env_ext
968 ELSE IF (para_env_rpa%num_pe == para_env%num_pe .AND. PRESENT(qs_env)) THEN
969 CALL get_qs_env(qs_env, blacs_env=blacs_env_q)
970 ELSE
971 CALL cp_blacs_env_create(blacs_env=blacs_env_q, para_env=para_env_rpa)
972 END IF
973 NULLIFY (fm_struct)
974 CALL cp_fm_struct_create(fm_struct, context=blacs_env_q, nrow_global=dimen_ri, &
975 ncol_global=dimen_ri, para_env=para_env_rpa)
976 DO ispin = 1, SIZE(fm_mat_q)
977 CALL cp_fm_create(fm_mat_q(ispin), fm_struct, name="fm_mat_Q")
978 END DO
979
980 CALL cp_fm_struct_release(fm_struct)
981
982 IF (.NOT. (my_blacs_ext .OR. (para_env_rpa%num_pe == para_env%num_pe .AND. PRESENT(qs_env)))) &
983 CALL cp_blacs_env_release(blacs_env_q)
984 END IF
985
986 ! release blacs_env
987 IF (.NOT. my_blacs_s_ext) THEN
988 CALL cp_blacs_env_release(blacs_env)
989 ELSE
990 NULLIFY (blacs_env)
991 END IF
992
993 CALL timestop(handle)
994
995 END SUBROUTINE create_integ_mat
996
997! **************************************************************************************************
998!> \brief ...
999!> \param qs_env ...
1000!> \param Erpa ...
1001!> \param mp2_env ...
1002!> \param para_env ...
1003!> \param para_env_RPA ...
1004!> \param para_env_sub ...
1005!> \param unit_nr ...
1006!> \param homo ...
1007!> \param virtual ...
1008!> \param dimen_RI ...
1009!> \param dimen_RI_red ...
1010!> \param dimen_ia ...
1011!> \param dimen_nm_gw ...
1012!> \param Eigenval ...
1013!> \param num_integ_points ...
1014!> \param num_integ_group ...
1015!> \param color_rpa_group ...
1016!> \param fm_matrix_PQ ...
1017!> \param fm_mat_S ...
1018!> \param fm_mat_Q_gemm ...
1019!> \param fm_mat_Q ...
1020!> \param fm_mat_S_gw ...
1021!> \param fm_mat_R_gw ...
1022!> \param fm_mat_S_ij_bse ...
1023!> \param fm_mat_S_ab_bse ...
1024!> \param my_do_gw ...
1025!> \param do_bse ...
1026!> \param gw_corr_lev_occ ...
1027!> \param gw_corr_lev_virt ...
1028!> \param bse_lev_virt ...
1029!> \param do_minimax_quad ...
1030!> \param do_im_time ...
1031!> \param mo_coeff ...
1032!> \param fm_matrix_L_kpoints ...
1033!> \param fm_matrix_Minv_L_kpoints ...
1034!> \param fm_matrix_Minv ...
1035!> \param fm_matrix_Minv_Vtrunc_Minv ...
1036!> \param mat_munu ...
1037!> \param mat_P_global ...
1038!> \param t_3c_M ...
1039!> \param t_3c_O ...
1040!> \param t_3c_O_compressed ...
1041!> \param t_3c_O_ind ...
1042!> \param starts_array_mc ...
1043!> \param ends_array_mc ...
1044!> \param starts_array_mc_block ...
1045!> \param ends_array_mc_block ...
1046!> \param matrix_s ...
1047!> \param do_kpoints_from_Gamma ...
1048!> \param kpoints ...
1049!> \param gd_array ...
1050!> \param color_sub ...
1051!> \param do_ri_sos_laplace_mp2 ...
1052!> \param calc_forces ...
1053! **************************************************************************************************
1054 SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_sub, unit_nr, &
1055 homo, virtual, dimen_RI, dimen_RI_red, dimen_ia, dimen_nm_gw, &
1056 Eigenval, num_integ_points, num_integ_group, color_rpa_group, &
1057 fm_matrix_PQ, fm_mat_S, fm_mat_Q_gemm, fm_mat_Q, fm_mat_S_gw, fm_mat_R_gw, &
1058 fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
1059 my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
1060 bse_lev_virt, &
1061 do_minimax_quad, do_im_time, mo_coeff, &
1062 fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
1063 fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, mat_munu, mat_P_global, &
1064 t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1065 starts_array_mc, ends_array_mc, &
1066 starts_array_mc_block, ends_array_mc_block, &
1067 matrix_s, do_kpoints_from_Gamma, kpoints, gd_array, color_sub, &
1068 do_ri_sos_laplace_mp2, calc_forces)
1069
1070 TYPE(qs_environment_type), POINTER :: qs_env
1071 REAL(kind=dp), INTENT(OUT) :: erpa
1072 TYPE(mp2_type) :: mp2_env
1073 TYPE(mp_para_env_type), POINTER :: para_env, para_env_rpa, para_env_sub
1074 INTEGER, INTENT(IN) :: unit_nr
1075 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
1076 INTEGER, INTENT(IN) :: dimen_ri, dimen_ri_red
1077 INTEGER, DIMENSION(:), INTENT(IN) :: dimen_ia
1078 INTEGER, INTENT(IN) :: dimen_nm_gw
1079 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1080 INTENT(INOUT) :: eigenval
1081 INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, &
1082 color_rpa_group
1083 TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_pq
1084 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_s
1085 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_q_gemm, fm_mat_q, fm_mat_s_gw
1086 TYPE(cp_fm_type), INTENT(IN) :: fm_mat_r_gw, fm_mat_s_ij_bse, &
1087 fm_mat_s_ab_bse
1088 LOGICAL, INTENT(IN) :: my_do_gw, do_bse
1089 INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt
1090 INTEGER, INTENT(IN) :: bse_lev_virt
1091 LOGICAL, INTENT(IN) :: do_minimax_quad, do_im_time
1092 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1093 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_l_kpoints, &
1094 fm_matrix_minv_l_kpoints, &
1095 fm_matrix_minv, &
1096 fm_matrix_minv_vtrunc_minv
1097 TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
1098 TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_p_global
1099 TYPE(dbt_type), INTENT(INOUT) :: t_3c_m
1100 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
1101 INTENT(INOUT) :: t_3c_o
1102 TYPE(hfx_compression_type), ALLOCATABLE, &
1103 DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_o_compressed
1104 TYPE(block_ind_type), ALLOCATABLE, &
1105 DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_o_ind
1106 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
1107 starts_array_mc_block, &
1108 ends_array_mc_block
1109 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1110 LOGICAL :: do_kpoints_from_gamma
1111 TYPE(kpoint_type), POINTER :: kpoints
1112 TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1113 INTEGER, INTENT(IN) :: color_sub
1114 LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, calc_forces
1115
1116 CHARACTER(LEN=*), PARAMETER :: routinen = 'rpa_num_int'
1117
1118 COMPLEX(KIND=dp), ALLOCATABLE, &
1119 DIMENSION(:, :, :, :) :: vec_sigma_c_gw
1120 INTEGER :: count_ev_sc_gw, cut_memory, group_size_p, gw_corr_lev_tot, handle, handle3, &
1121 i, ikp_local, ispin, iter_evgw, iter_sc_gw0, j, jquad, min_bsize, mm_style, nkp, &
1122 nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, pspin, qspin, &
1123 size_p
1124#if defined (__GREENX)
1125 INTEGER :: gi, ierr
1126#endif
1127 INTEGER(int_8) :: dbcsr_nflop
1128 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c
1129 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c
1130 INTEGER, DIMENSION(:), POINTER :: col_blk_size, prim_blk_sizes, &
1131 ri_blk_sizes
1132 LOGICAL :: do_apply_ic_corr_to_gw, do_gw_im_time, do_ic_model, do_kpoints_cubic_rpa, &
1133 do_periodic, do_print, do_ri_sigma_x, exit_ev_gw, first_cycle, &
1134 first_cycle_periodic_correction, my_open_shell, print_ic_values
1135 LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: has_mat_p_blocks
1136 REAL(kind=dp) :: a_scaling, alpha, dbcsr_time, e_exchange, e_exchange_corr, &
1137 eps_filter, eps_filter_im_time, ext_scaling, fermi_level_offset, fermi_level_offset_input, &
1138 my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old
1139#if defined (__GREENX)
1140 REAL(kind=dp) :: cosft_duality_error_greenx, emax, emin, &
1141 max_errors_greenx(3)
1142#endif
1143 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, e_fermi, tau_tj, tau_wj, tj, &
1144 trace_qomega, vec_omega_fit_gw, wj, &
1145 wkp_w
1146 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_w_gw, weights_cos_tf_t_to_w, &
1147 weights_cos_tf_w_to_t, &
1148 weights_sin_tf_t_to_w
1149 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: eigenval_last, eigenval_scf, &
1150 vec_sigma_x_gw
1151 TYPE(cp_cfm_type) :: cfm_mat_q
1152 TYPE(cp_fm_type) :: fm_mat_q_static_bse_gemm, fm_mat_ri_global_work, fm_mat_s_ia_bse, &
1153 fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, &
1154 fm_scaled_dm_virt_tau
1155 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_s_gw_work, fm_mat_w, &
1156 fm_mo_coeff_occ, fm_mo_coeff_virt
1157 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_l_kpoints, fm_mat_minv_l_kpoints
1158 TYPE(dbcsr_p_type) :: mat_dm, mat_l, mat_m_p_munu_occ, &
1159 mat_m_p_munu_virt, mat_minvvminv
1160 TYPE(dbcsr_p_type), ALLOCATABLE, &
1161 DIMENSION(:, :, :) :: mat_p_omega
1162 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_im_mo_mo, &
1163 matrix_berry_re_mo_mo
1164 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_p_omega_kp
1165 TYPE(dbcsr_type), POINTER :: mat_w, mat_work
1166 TYPE(dbt_type) :: t_3c_overl_int_ao_mo
1167 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_3c_overl_int_gw_ao, &
1168 t_3c_overl_int_gw_ri, &
1169 t_3c_overl_nnp_ic, &
1170 t_3c_overl_nnp_ic_reflected
1171 TYPE(dgemm_counter_type) :: dgemm_counter
1172 TYPE(hfx_compression_type), ALLOCATABLE, &
1173 DIMENSION(:) :: t_3c_o_mo_compressed
1174 TYPE(im_time_force_type) :: force_data
1175 TYPE(rpa_exchange_work_type) :: exchange_work
1176 TYPE(rpa_grad_type) :: rpa_grad
1177 TYPE(rpa_sigma_type) :: rpa_sigma
1178 TYPE(two_dim_int_array), ALLOCATABLE, DIMENSION(:) :: t_3c_o_mo_ind
1179
1180 CALL timeset(routinen, handle)
1181
1182 nspins = SIZE(homo)
1183 nmo = homo(1) + virtual(1)
1184
1185 my_open_shell = (nspins == 2)
1186
1187 do_gw_im_time = my_do_gw .AND. do_im_time
1188 do_ri_sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x
1189 do_ic_model = mp2_env%ri_g0w0%do_ic_model
1190 print_ic_values = mp2_env%ri_g0w0%print_ic_values
1191 do_periodic = mp2_env%ri_g0w0%do_periodic
1192 do_kpoints_cubic_rpa = mp2_env%ri_rpa_im_time%do_im_time_kpoints
1193
1194 ! For SOS-MP2 only gemm is implemented
1195 mm_style = wfc_mm_style_gemm
1196 IF (.NOT. do_ri_sos_laplace_mp2) mm_style = mp2_env%ri_rpa%mm_style
1197
1198 IF (my_do_gw) THEN
1199 ext_scaling = 0.2_dp
1200 omega_max_fit = mp2_env%ri_g0w0%omega_max_fit
1201 fermi_level_offset_input = mp2_env%ri_g0w0%fermi_level_offset
1202 iter_evgw = mp2_env%ri_g0w0%iter_evGW
1203 iter_sc_gw0 = mp2_env%ri_g0w0%iter_sc_GW0
1204 IF ((.NOT. do_im_time)) THEN
1205 IF (iter_sc_gw0 .NE. 1 .AND. iter_evgw .NE. 1) cpabort("Mixed scGW0/evGW not implemented.")
1206 ! in case of scGW0 with the N^4 algorithm, we use the evGW code but use the DFT eigenvalues for W
1207 IF (iter_sc_gw0 .NE. 1) iter_evgw = iter_sc_gw0
1208 END IF
1209 ELSE
1210 ext_scaling = 0.0_dp
1211 iter_evgw = 1
1212 iter_sc_gw0 = 1
1213 END IF
1214
1215 IF (do_kpoints_cubic_rpa .AND. do_ri_sos_laplace_mp2) THEN
1216 cpabort("RI-SOS-Laplace-MP2 with k-point-sampling is not implemented.")
1217 END IF
1218
1219 do_apply_ic_corr_to_gw = .false.
1220 IF (mp2_env%ri_g0w0%ic_corr_list(1)%array(1) > 0.0_dp) do_apply_ic_corr_to_gw = .true.
1221
1222 IF (do_im_time) THEN
1223 cpassert(do_minimax_quad .OR. do_ri_sos_laplace_mp2)
1224
1225 group_size_p = mp2_env%ri_rpa_im_time%group_size_P
1226 cut_memory = mp2_env%ri_rpa_im_time%cut_memory
1227 eps_filter = mp2_env%ri_rpa_im_time%eps_filter
1228 eps_filter_im_time = mp2_env%ri_rpa_im_time%eps_filter* &
1229 mp2_env%ri_rpa_im_time%eps_filter_factor
1230
1231 min_bsize = mp2_env%ri_rpa_im_time%min_bsize
1232
1233 CALL alloc_im_time(qs_env, para_env, dimen_ri, dimen_ri_red, &
1234 num_integ_points, nspins, fm_mat_q(1), fm_mo_coeff_occ, fm_mo_coeff_virt, &
1235 fm_matrix_minv_l_kpoints, fm_matrix_l_kpoints, mat_p_global, &
1236 t_3c_o, matrix_s, kpoints, eps_filter_im_time, &
1237 cut_memory, nkp, num_cells_dm, num_3c_repl, &
1238 size_p, ikp_local, &
1239 index_to_cell_3c, &
1240 cell_to_index_3c, &
1241 col_blk_size, &
1242 do_ic_model, do_kpoints_cubic_rpa, &
1243 do_kpoints_from_gamma, do_ri_sigma_x, my_open_shell, &
1244 has_mat_p_blocks, wkp_w, &
1245 cfm_mat_q, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, &
1246 fm_mat_ri_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, &
1247 fm_mo_coeff_virt_scaled, mat_dm, mat_l, mat_m_p_munu_occ, mat_m_p_munu_virt, &
1248 mat_minvvminv, mat_p_omega, mat_p_omega_kp, mat_work, mo_coeff, &
1249 fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, homo, nmo)
1250
1251 IF (calc_forces) CALL init_im_time_forces(force_data, fm_matrix_pq, t_3c_m, unit_nr, mp2_env, qs_env)
1252
1253 IF (my_do_gw) THEN
1254
1255 CALL dbcsr_get_info(mat_p_global%matrix, &
1256 row_blk_size=ri_blk_sizes)
1257
1258 CALL dbcsr_get_info(matrix_s(1)%matrix, &
1259 row_blk_size=prim_blk_sizes)
1260
1261 gw_corr_lev_tot = gw_corr_lev_occ(1) + gw_corr_lev_virt(1)
1262
1263 IF (.NOT. do_kpoints_cubic_rpa) THEN
1264 CALL allocate_matrices_gw_im_time(gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, &
1265 num_integ_points, unit_nr, &
1266 ri_blk_sizes, do_ic_model, &
1267 para_env, fm_mat_w, fm_mat_q(1), &
1268 mo_coeff, &
1269 t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, &
1270 t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, &
1271 starts_array_mc, ends_array_mc, &
1272 t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, &
1273 matrix_s, mat_w, t_3c_o, &
1274 t_3c_o_compressed, t_3c_o_ind, &
1275 qs_env)
1276
1277 END IF
1278 END IF
1279
1280 END IF
1281 IF (do_ic_model) THEN
1282 ! image charge model only implemented for cubic scaling GW
1283 cpassert(do_gw_im_time)
1284 cpassert(.NOT. do_periodic)
1285 IF (cut_memory .NE. 1) cpabort("For IC, use MEMORY_CUT 1 in the LOW_SCALING section.")
1286 END IF
1287
1288 ALLOCATE (e_fermi(nspins))
1289 IF (do_minimax_quad .OR. do_ri_sos_laplace_mp2) THEN
1290 do_print = .NOT. do_ic_model
1291#if defined (__GREENX)
1292 CALL init_greenx_grids(qs_env, para_env, homo, eigenval, do_ri_sos_laplace_mp2, &
1293 do_kpoints_cubic_rpa, emin, emax, e_fermi(1))
1294 CALL gx_minimax_grid(num_integ_points, emin, emax, tau_tj, tau_wj, tj, wj, &
1295 weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
1296 max_errors_greenx, cosft_duality_error_greenx, ierr, &
1297 bare_cos_sin_weights=.true., &
1298 regularization=qs_env%mp2_env%ri_g0w0%regularization_minimax)
1299 ! Factor 4 is hard-coded in the RPA weights in the internal CP2K minimax routines
1300 wj(:) = wj(:)*4.0_dp
1301 IF (ierr == 0) THEN
1302 IF (unit_nr > 0) THEN
1303 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
1304 "GREENX MINIMAX_INFO| Number of integration points:", num_integ_points
1305 WRITE (unit=unit_nr, fmt="(T3,A,T61,3F20.7)") &
1306 "GREENX MINIMAX_INFO| Gap (Emin):", emin
1307 WRITE (unit=unit_nr, fmt="(T3,A,T61,3F20.7)") &
1308 "GREENX MINIMAX_INFO| Maximum eigenvalue difference (Emax):", emax
1309 WRITE (unit=unit_nr, fmt="(T3,A,T61,3F20.4)") &
1310 "GREENX MINIMAX_INFO| Energy range (Emax/Emin):", emax/emin
1311 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
1312 "GREENX MINIMAX_INFO| Frequency grid (scaled):", "Weights", "Abscissas"
1313 DO gi = 1, num_integ_points
1314 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") wj(gi), tj(gi)
1315 END DO
1316 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
1317 "GREENX MINIMAX_INFO| Time grid (scaled):", "Weights", "Abscissas"
1318 DO gi = 1, num_integ_points
1319 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") tau_wj(gi), tau_tj(gi)
1320 END DO
1321 CALL m_flush(unit_nr)
1322 END IF
1323 ELSE
1324 IF (unit_nr > 0) THEN
1325 WRITE (unit=unit_nr, fmt="(T3,A,T75)") &
1326 "GREENX MINIMAX_INFO| Grid not available, use internal CP2K grids"
1327 CALL m_flush(unit_nr)
1328 END IF
1329 IF (ALLOCATED(tau_tj)) THEN
1330 DEALLOCATE (tau_tj)
1331 END IF
1332 IF (ALLOCATED(tau_wj)) THEN
1333 DEALLOCATE (tau_wj)
1334 END IF
1335 IF (ALLOCATED(tj)) THEN
1336 DEALLOCATE (tj)
1337 END IF
1338 IF (ALLOCATED(wj)) THEN
1339 DEALLOCATE (wj)
1340 END IF
1341 IF (ALLOCATED(weights_cos_tf_t_to_w)) THEN
1342 DEALLOCATE (weights_cos_tf_t_to_w)
1343 END IF
1344 IF (ALLOCATED(weights_cos_tf_w_to_t)) THEN
1345 DEALLOCATE (weights_cos_tf_w_to_t)
1346 END IF
1347 IF (ALLOCATED(weights_sin_tf_t_to_w)) THEN
1348 DEALLOCATE (weights_sin_tf_t_to_w)
1349 END IF
1350 CALL get_minimax_grid(para_env, unit_nr, homo, eigenval, num_integ_points, do_im_time, &
1351 do_ri_sos_laplace_mp2, do_print, &
1352 tau_tj, tau_wj, qs_env, do_gw_im_time, &
1353 do_kpoints_cubic_rpa, e_fermi(1), tj, wj, &
1354 weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
1355 qs_env%mp2_env%ri_g0w0%regularization_minimax)
1356 END IF
1357#else
1358 CALL get_minimax_grid(para_env, unit_nr, homo, eigenval, num_integ_points, do_im_time, &
1359 do_ri_sos_laplace_mp2, do_print, &
1360 tau_tj, tau_wj, qs_env, do_gw_im_time, &
1361 do_kpoints_cubic_rpa, e_fermi(1), tj, wj, &
1362 weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
1363 qs_env%mp2_env%ri_g0w0%regularization_minimax)
1364#endif
1365
1366 !For sos_laplace_mp2 and low-scaling RPA, potentially need to store/retrieve the initial weights
1367 IF (qs_env%mp2_env%ri_rpa_im_time%keep_quad) THEN
1368 CALL keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, &
1369 weights_cos_tf_w_to_t, do_ri_sos_laplace_mp2, do_im_time, &
1370 num_integ_points, unit_nr, qs_env)
1371 END IF
1372 ELSE
1373 IF (calc_forces) cpabort("Forces with Clenshaw-Curtis grid not implemented.")
1374 CALL get_clenshaw_grid(para_env, para_env_rpa, unit_nr, homo, virtual, eigenval, num_integ_points, &
1375 num_integ_group, color_rpa_group, fm_mat_s, my_do_gw, &
1376 ext_scaling, a_scaling, tj, wj)
1377 END IF
1378
1379 ! This array is needed for RPA
1380 IF (.NOT. do_ri_sos_laplace_mp2) THEN
1381 ALLOCATE (trace_qomega(dimen_ri_red))
1382 END IF
1383
1384 IF (do_ri_sos_laplace_mp2 .AND. .NOT. do_im_time) THEN
1385 alpha = 1.0_dp
1386 ELSE IF (my_open_shell .OR. do_ri_sos_laplace_mp2) THEN
1387 alpha = 2.0_dp
1388 ELSE
1389 alpha = 4.0_dp
1390 END IF
1391 IF (my_do_gw) THEN
1392 CALL allocate_matrices_gw(vec_sigma_c_gw, color_rpa_group, dimen_nm_gw, &
1393 gw_corr_lev_occ, gw_corr_lev_virt, homo, &
1394 nmo, num_integ_group, num_integ_points, unit_nr, &
1395 gw_corr_lev_tot, num_fit_points, omega_max_fit, &
1396 do_minimax_quad, do_periodic, do_ri_sigma_x,.NOT. do_im_time, &
1397 first_cycle_periodic_correction, &
1398 a_scaling, eigenval, tj, vec_omega_fit_gw, vec_sigma_x_gw, &
1399 delta_corr, eigenval_last, eigenval_scf, vec_w_gw, &
1400 fm_mat_s_gw, fm_mat_s_gw_work, &
1401 para_env, mp2_env, kpoints, nkp, nkp_self_energy, &
1402 do_kpoints_cubic_rpa, do_kpoints_from_gamma)
1403
1404 IF (do_bse) THEN
1405
1406 CALL cp_fm_create(fm_mat_q_static_bse_gemm, fm_mat_q_gemm(1)%matrix_struct)
1407 CALL cp_fm_to_fm(fm_mat_q_gemm(1), fm_mat_q_static_bse_gemm)
1408 CALL cp_fm_set_all(fm_mat_q_static_bse_gemm, 0.0_dp)
1409
1410 END IF
1411
1412 END IF
1413
1414 IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_create(rpa_grad, fm_mat_q(1), &
1415 fm_mat_s, homo, virtual, mp2_env, eigenval(:, 1, :), &
1416 unit_nr, do_ri_sos_laplace_mp2)
1417!doubtful
1418 IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) &
1419 CALL exchange_work%create(qs_env, para_env_sub, mat_munu, dimen_ri_red, &
1420 fm_mat_s, fm_mat_q(1), fm_mat_q_gemm(1), homo, virtual)
1421 erpa = 0.0_dp
1422 IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) e_exchange = 0.0_dp
1423 first_cycle = .true.
1424 omega_old = 0.0_dp
1425 CALL dgemm_counter_init(dgemm_counter, unit_nr, mp2_env%ri_rpa%print_dgemm_info)
1426
1427 DO count_ev_sc_gw = 1, iter_evgw
1428 dbcsr_time = 0.0_dp
1429 dbcsr_nflop = 0
1430
1431 IF (do_ic_model) cycle
1432
1433 ! reset some values, important when doing eigenvalue self-consistent GW
1434 IF (my_do_gw) THEN
1435 erpa = 0.0_dp
1436 vec_sigma_c_gw = z_zero
1437 first_cycle = .true.
1438 END IF
1439
1440 ! calculate Q_PQ(it)
1441 IF (do_im_time) THEN ! not using Imaginary time
1442
1443 IF (.NOT. do_kpoints_cubic_rpa) THEN
1444 DO ispin = 1, nspins
1445 e_fermi(ispin) = (eigenval(homo(ispin), 1, ispin) + eigenval(homo(ispin) + 1, 1, ispin))*0.5_dp
1446 END DO
1447 END IF
1448
1449 tau = 0.0_dp
1450 tau_old = 0.0_dp
1451
1452 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(/T3,A,T66,i15)") &
1453 "MEMORY_INFO| Memory cut:", cut_memory
1454 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,ES15.2)") &
1455 "SPARSITY_INFO| Eps filter for M virt/occ tensors:", eps_filter
1456 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,ES15.2)") &
1457 "SPARSITY_INFO| Eps filter for P matrix:", eps_filter_im_time
1458 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,i15)") &
1459 "SPARSITY_INFO| Minimum tensor block size:", min_bsize
1460
1461 ! for evGW, we have to ensure that mat_P_omega is zero
1462 CALL zero_mat_p_omega(mat_p_omega(:, :, 1))
1463
1464 ! compute the matrix Q(it) and Fourier transform it directly to mat_P_omega(iw)
1465 CALL compute_mat_p_omega(mat_p_omega(:, :, 1), fm_scaled_dm_occ_tau, &
1466 fm_scaled_dm_virt_tau, fm_mo_coeff_occ(1), fm_mo_coeff_virt(1), &
1467 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1468 mat_p_global, matrix_s, 1, &
1469 t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, &
1470 starts_array_mc, ends_array_mc, &
1471 starts_array_mc_block, ends_array_mc_block, &
1472 weights_cos_tf_t_to_w, tj, tau_tj, e_fermi(1), eps_filter, alpha, &
1473 eps_filter_im_time, eigenval(:, 1, 1), nmo, &
1474 num_integ_points, cut_memory, &
1475 unit_nr, mp2_env, para_env, &
1476 qs_env, do_kpoints_from_gamma, &
1477 index_to_cell_3c, cell_to_index_3c, &
1478 has_mat_p_blocks, do_ri_sos_laplace_mp2, &
1479 dbcsr_time, dbcsr_nflop)
1480
1481 ! the same for open shell, use fm_mo_coeff_occ_beta and fm_mo_coeff_virt_beta
1482 IF (my_open_shell) THEN
1483 CALL zero_mat_p_omega(mat_p_omega(:, :, 2))
1484 CALL compute_mat_p_omega(mat_p_omega(:, :, 2), fm_scaled_dm_occ_tau, &
1485 fm_scaled_dm_virt_tau, fm_mo_coeff_occ(2), &
1486 fm_mo_coeff_virt(2), &
1487 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1488 mat_p_global, matrix_s, 2, &
1489 t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, &
1490 starts_array_mc, ends_array_mc, &
1491 starts_array_mc_block, ends_array_mc_block, &
1492 weights_cos_tf_t_to_w, tj, tau_tj, e_fermi(2), eps_filter, alpha, &
1493 eps_filter_im_time, eigenval(:, 1, 2), nmo, &
1494 num_integ_points, cut_memory, &
1495 unit_nr, mp2_env, para_env, &
1496 qs_env, do_kpoints_from_gamma, &
1497 index_to_cell_3c, cell_to_index_3c, &
1498 has_mat_p_blocks, do_ri_sos_laplace_mp2, &
1499 dbcsr_time, dbcsr_nflop)
1500
1501 !For RPA, we sum up the P matrices. If no force needed, can clean-up the beta spin one
1502 IF (.NOT. do_ri_sos_laplace_mp2) THEN
1503 DO j = 1, SIZE(mat_p_omega, 2)
1504 DO i = 1, SIZE(mat_p_omega, 1)
1505 CALL dbcsr_add(mat_p_omega(i, j, 1)%matrix, mat_p_omega(i, j, 2)%matrix, 1.0_dp, 1.0_dp)
1506 IF (.NOT. calc_forces) CALL dbcsr_clear(mat_p_omega(i, j, 2)%matrix)
1507 END DO
1508 END DO
1509 END IF
1510 END IF ! my_open_shell
1511
1512 END IF ! do im time
1513
1514 IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1515 CALL rpa_sigma_create(rpa_sigma, mp2_env%ri_rpa%sigma_param, fm_mat_q(1), unit_nr, para_env)
1516 END IF
1517
1518 DO jquad = 1, num_integ_points
1519 IF (modulo(jquad, num_integ_group) /= color_rpa_group) cycle
1520
1521 CALL timeset(routinen//"_RPA_matrix_operations", handle3)
1522
1523 IF (do_ri_sos_laplace_mp2) THEN
1524 omega = tau_tj(jquad)
1525 ELSE
1526 IF (do_minimax_quad) THEN
1527 omega = tj(jquad)
1528 ELSE
1529 omega = a_scaling/tan(tj(jquad))
1530 END IF
1531 END IF ! do_ri_sos_laplace_mp2
1532
1533 IF (do_im_time) THEN
1534 ! in case we do imag time, we already calculated fm_mat_Q by a Fourier transform from im. time
1535
1536 IF (.NOT. (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma)) THEN
1537
1538 DO ispin = 1, SIZE(mat_p_omega, 3)
1539 CALL contract_p_omega_with_mat_l(mat_p_omega(jquad, 1, ispin)%matrix, mat_l%matrix, mat_work, &
1540 eps_filter_im_time, fm_mat_work, dimen_ri, dimen_ri_red, &
1541 fm_mat_minv_l_kpoints(1, 1), fm_mat_q(ispin))
1542 END DO
1543 END IF
1544
1545 ELSE
1546 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3, A, 1X, I3, 1X, A, 1X, I3)") &
1547 "INTEG_INFO| Started with Integration point", jquad, "of", num_integ_points
1548
1549 IF (first_cycle .AND. count_ev_sc_gw > 1) THEN
1550 IF (iter_sc_gw0 == 1) THEN
1551 DO ispin = 1, nspins
1552 CALL remove_scaling_factor_rpa(fm_mat_s(ispin), virtual(ispin), &
1553 eigenval_last(:, 1, ispin), homo(ispin), omega_old)
1554 END DO
1555 ELSE
1556 DO ispin = 1, nspins
1557 CALL remove_scaling_factor_rpa(fm_mat_s(ispin), virtual(ispin), &
1558 eigenval_scf(:, 1, ispin), homo(ispin), omega_old)
1559 END DO
1560 END IF
1561 END IF
1562
1563 IF (iter_sc_gw0 > 1) THEN
1564 DO ispin = 1, nspins
1565 CALL calc_mat_q(fm_mat_s(ispin), do_ri_sos_laplace_mp2, first_cycle, virtual(ispin), &
1566 eigenval_scf(:, 1, ispin), homo(ispin), omega, omega_old, jquad, mm_style, &
1567 dimen_ri_red, dimen_ia(ispin), alpha, fm_mat_q(ispin), &
1568 fm_mat_q_gemm(ispin), do_bse, fm_mat_q_static_bse_gemm, dgemm_counter, &
1569 num_integ_points, count_ev_sc_gw)
1570 END DO
1571
1572 ! For SOS-MP2 we need both matrices separately
1573 IF (.NOT. do_ri_sos_laplace_mp2) THEN
1574 DO ispin = 2, nspins
1575 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_q(1), beta=1.0_dp, matrix_b=fm_mat_q(ispin))
1576 END DO
1577 END IF
1578 ELSE
1579 DO ispin = 1, nspins
1580 CALL calc_mat_q(fm_mat_s(ispin), do_ri_sos_laplace_mp2, first_cycle, virtual(ispin), &
1581 eigenval(:, 1, ispin), homo(ispin), omega, omega_old, jquad, mm_style, &
1582 dimen_ri_red, dimen_ia(ispin), alpha, fm_mat_q(ispin), &
1583 fm_mat_q_gemm(ispin), do_bse, fm_mat_q_static_bse_gemm, dgemm_counter, &
1584 num_integ_points, count_ev_sc_gw)
1585 END DO
1586
1587 ! For SOS-MP2 we need both matrices separately
1588 IF (.NOT. do_ri_sos_laplace_mp2) THEN
1589 DO ispin = 2, nspins
1590 CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_mat_q(1), beta=1.0_dp, matrix_b=fm_mat_q(ispin))
1591 END DO
1592 END IF
1593
1594 END IF
1595
1596 END IF ! im time
1597
1598 ! Calculate RPA exchange energy correction
1599 IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
1600 e_exchange_corr = 0.0_dp
1601 CALL exchange_work%compute(fm_mat_q(1), eigenval(:, 1, :), fm_mat_s, omega, e_exchange_corr, mp2_env)
1602
1603 ! Evaluate the final exchange energy correction
1604 e_exchange = e_exchange + e_exchange_corr*wj(jquad)
1605 END IF
1606
1607 ! for developing Sigma functional closed and open shell are taken cared for
1608 IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1609 CALL rpa_sigma_matrix_spectral(rpa_sigma, fm_mat_q(1), wj(jquad), para_env_rpa)
1610 END IF
1611
1612 IF (do_ri_sos_laplace_mp2) THEN
1613
1614 CALL sos_mp2_postprocessing(fm_mat_q, erpa, tau_wj(jquad))
1615
1616 IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
1617 fm_mat_q, fm_mat_q_gemm, dgemm_counter, fm_mat_s, omega, homo, virtual, &
1618 eigenval(:, 1, :), tau_wj(jquad), unit_nr)
1619 ELSE
1620 IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_copy_q(fm_mat_q(1), rpa_grad)
1621
1622 CALL q_trace_and_add_unit_matrix(dimen_ri_red, trace_qomega, fm_mat_q(1))
1623
1624 IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
1625 CALL invert_eps_compute_w_and_erpa_kp(dimen_ri, num_integ_points, jquad, nkp, count_ev_sc_gw, para_env, &
1626 erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, &
1627 wkp_w, do_gw_im_time, do_ri_sigma_x, do_kpoints_from_gamma, &
1628 cfm_mat_q, ikp_local, &
1629 mat_p_omega(:, :, 1), mat_p_omega_kp, qs_env, eps_filter_im_time, unit_nr, &
1630 kpoints, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, &
1631 fm_mat_w, fm_mat_ri_global_work, mat_minvvminv, &
1632 fm_matrix_minv, fm_matrix_minv_vtrunc_minv)
1633 ELSE
1634 CALL compute_erpa_by_freq_int(dimen_ri_red, trace_qomega, fm_mat_q(1), para_env_rpa, erpa, wj(jquad))
1635 END IF
1636
1637 IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
1638 fm_mat_q, fm_mat_q_gemm, dgemm_counter, fm_mat_s, omega, homo, virtual, &
1639 eigenval(:, 1, :), wj(jquad), unit_nr)
1640 END IF ! do_ri_sos_laplace_mp2
1641
1642 ! save omega and reset the first_cycle flag
1643 first_cycle = .false.
1644 omega_old = omega
1645
1646 CALL timestop(handle3)
1647
1648 IF (my_do_gw) THEN
1649
1650 CALL get_fermi_level_offset(fermi_level_offset, fermi_level_offset_input, eigenval(:, 1, :), homo)
1651
1652 ! do_im_time = TRUE means low-scaling calculation
1653 IF (do_im_time) THEN
1654 ! only for molecules
1655 IF (.NOT. (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma)) THEN
1656 CALL compute_w_cubic_gw(fm_mat_w, fm_mat_q(1), fm_mat_work, dimen_ri, fm_mat_minv_l_kpoints, num_integ_points, &
1657 tj, tau_tj, weights_cos_tf_w_to_t, jquad, omega)
1658 END IF
1659 ELSE
1660 CALL compute_gw_self_energy(vec_sigma_c_gw, dimen_nm_gw, dimen_ri_red, gw_corr_lev_occ, &
1661 gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, &
1662 do_im_time, do_periodic, first_cycle_periodic_correction, &
1663 fermi_level_offset, &
1664 omega, eigenval(:, 1, :), delta_corr, vec_omega_fit_gw, vec_w_gw, wj, &
1665 fm_mat_q(1), fm_mat_r_gw, fm_mat_s_gw, &
1666 fm_mat_s_gw_work, mo_coeff(1), para_env, &
1667 para_env_rpa, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, &
1668 kpoints, qs_env, mp2_env)
1669 END IF
1670 END IF
1671
1672 IF (unit_nr > 0) CALL m_flush(unit_nr)
1673 CALL para_env%sync() ! sync to see output
1674
1675 END DO ! jquad
1676
1677 IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1678 CALL finalize_rpa_sigma(rpa_sigma, unit_nr, mp2_env%ri_rpa%e_sigma_corr, para_env, do_minimax_quad)
1679 IF (do_minimax_quad) mp2_env%ri_rpa%e_sigma_corr = mp2_env%ri_rpa%e_sigma_corr/2.0_dp
1680 CALL para_env%sum(mp2_env%ri_rpa%e_sigma_corr)
1681 END IF
1682
1683 CALL para_env%sum(erpa)
1684
1685 IF (.NOT. (do_ri_sos_laplace_mp2)) THEN
1686 erpa = erpa/(pi*2.0_dp)
1687 IF (do_minimax_quad) erpa = erpa/2.0_dp
1688 END IF
1689
1690 IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
1691 CALL para_env%sum(e_exchange)
1692 e_exchange = e_exchange/(pi*2.0_dp)
1693 IF (do_minimax_quad) e_exchange = e_exchange/2.0_dp
1694 mp2_env%ri_rpa%ener_exchange = e_exchange
1695 END IF
1696
1697 IF (calc_forces .AND. do_ri_sos_laplace_mp2 .AND. do_im_time) THEN
1698 IF (my_open_shell) THEN
1699 pspin = 1
1700 qspin = 2
1701 CALL calc_laplace_loop_forces(force_data, mat_p_omega(:, 1, :), t_3c_m, t_3c_o(1, 1), &
1702 t_3c_o_compressed(1, 1, :), t_3c_o_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1703 fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1704 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1705 starts_array_mc, ends_array_mc, starts_array_mc_block, &
1706 ends_array_mc_block, num_integ_points, nmo, eigenval(:, 1, :), &
1707 tau_tj, tau_wj, cut_memory, pspin, qspin, my_open_shell, &
1708 unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1709 pspin = 2
1710 qspin = 1
1711 CALL calc_laplace_loop_forces(force_data, mat_p_omega(:, 1, :), t_3c_m, t_3c_o(1, 1), &
1712 t_3c_o_compressed(1, 1, :), t_3c_o_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1713 fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1714 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1715 starts_array_mc, ends_array_mc, starts_array_mc_block, &
1716 ends_array_mc_block, num_integ_points, nmo, eigenval(:, 1, :), &
1717 tau_tj, tau_wj, cut_memory, pspin, qspin, my_open_shell, &
1718 unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1719
1720 ELSE
1721 pspin = 1
1722 qspin = 1
1723 CALL calc_laplace_loop_forces(force_data, mat_p_omega(:, 1, :), t_3c_m, t_3c_o(1, 1), &
1724 t_3c_o_compressed(1, 1, :), t_3c_o_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1725 fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1726 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1727 starts_array_mc, ends_array_mc, starts_array_mc_block, &
1728 ends_array_mc_block, num_integ_points, nmo, eigenval(:, 1, :), &
1729 tau_tj, tau_wj, cut_memory, pspin, qspin, my_open_shell, &
1730 unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1731 END IF
1732 CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1733 END IF !laplace SOS-MP2
1734
1735 IF (calc_forces .AND. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) THEN
1736 DO ispin = 1, nspins
1737 CALL calc_rpa_loop_forces(force_data, mat_p_omega(:, 1, :), t_3c_m, t_3c_o(1, 1), &
1738 t_3c_o_compressed(1, 1, :), t_3c_o_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1739 fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1740 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1741 starts_array_mc, ends_array_mc, starts_array_mc_block, &
1742 ends_array_mc_block, num_integ_points, nmo, eigenval(:, 1, :), &
1743 e_fermi(ispin), weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, tj, &
1744 wj, tau_tj, cut_memory, ispin, my_open_shell, unit_nr, dbcsr_time, &
1745 dbcsr_nflop, mp2_env, qs_env)
1746 END DO
1747 CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1748 END IF
1749
1750 IF (do_im_time) THEN
1751
1752 my_flop_rate = real(dbcsr_nflop, dp)/(1.0e09_dp*dbcsr_time)
1753 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(/T3,A,T73,ES8.2)") &
1754 "PERFORMANCE| DBCSR total number of flops:", real(dbcsr_nflop*para_env%num_pe, dp)
1755 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.2)") &
1756 "PERFORMANCE| DBCSR total execution time:", dbcsr_time
1757 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.2)") &
1758 "PERFORMANCE| DBCSR flop rate (Gflops / MPI rank):", my_flop_rate
1759
1760 ELSE
1761
1762 CALL dgemm_counter_write(dgemm_counter, para_env)
1763
1764 END IF
1765
1766 ! GW: for low-scaling calculation: Compute self-energy Sigma(i*tau), Sigma(i*omega)
1767 ! for low-scaling and ordinary-scaling: analytic continuation from Sigma(iw) -> Sigma(w)
1768 ! and correction of quasiparticle energies e_n^GW
1769 IF (my_do_gw) THEN
1770
1771 CALL compute_qp_energies(vec_sigma_c_gw, count_ev_sc_gw, gw_corr_lev_occ, &
1772 gw_corr_lev_tot, gw_corr_lev_virt, homo, &
1773 nmo, num_fit_points, num_integ_points, &
1774 unit_nr, do_apply_ic_corr_to_gw, do_im_time, &
1775 do_periodic, do_ri_sigma_x, first_cycle_periodic_correction, &
1776 e_fermi, eps_filter, fermi_level_offset, &
1777 delta_corr, eigenval, &
1778 eigenval_last, eigenval_scf, iter_sc_gw0, exit_ev_gw, tau_tj, tj, &
1779 vec_omega_fit_gw, vec_sigma_x_gw, mp2_env%ri_g0w0%ic_corr_list, &
1780 weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, &
1781 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_mo_coeff_occ, &
1782 fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
1783 mo_coeff(1), fm_mat_w, para_env, para_env_rpa, mat_dm, mat_minvvminv, &
1784 t_3c_o, t_3c_m, t_3c_overl_int_ao_mo, t_3c_o_compressed, t_3c_o_mo_compressed, &
1785 t_3c_o_ind, t_3c_o_mo_ind, &
1786 t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, &
1787 matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, mat_w, matrix_s, &
1788 kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_rpa, &
1789 starts_array_mc, ends_array_mc)
1790
1791 ! if HOMO-LUMO gap differs by less than mp2_env%ri_g0w0%eps_ev_sc_iter, exit ev sc GW loop
1792 IF (exit_ev_gw) EXIT
1793
1794 END IF ! my_do_gw if
1795
1796 END DO ! evGW loop
1797
1798 IF (do_ic_model) THEN
1799
1800 IF (my_open_shell) THEN
1801
1802 CALL calculate_ic_correction(eigenval(:, 1, 1), mat_minvvminv%matrix, &
1803 t_3c_overl_nnp_ic(1), t_3c_overl_nnp_ic_reflected(1), &
1804 gw_corr_lev_tot, &
1805 gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), unit_nr, &
1806 print_ic_values, para_env, do_alpha=.true.)
1807
1808 CALL calculate_ic_correction(eigenval(:, 1, 2), mat_minvvminv%matrix, &
1809 t_3c_overl_nnp_ic(2), t_3c_overl_nnp_ic_reflected(2), &
1810 gw_corr_lev_tot, &
1811 gw_corr_lev_occ(2), gw_corr_lev_virt(2), homo(2), unit_nr, &
1812 print_ic_values, para_env, do_beta=.true.)
1813
1814 ELSE
1815
1816 CALL calculate_ic_correction(eigenval(:, 1, 1), mat_minvvminv%matrix, &
1817 t_3c_overl_nnp_ic(1), t_3c_overl_nnp_ic_reflected(1), &
1818 gw_corr_lev_tot, &
1819 gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), unit_nr, &
1820 print_ic_values, para_env)
1821
1822 END IF
1823
1824 END IF
1825
1826 ! postprocessing after GW for Bethe-Salpeter
1827 IF (do_bse) THEN
1828 ! Check used GW flavor; in Case of evGW we use W0 for BSE
1829 ! Use environment variable, since local iter_evGW is overwritten if evGW0 is invoked
1830 IF (mp2_env%ri_g0w0%iter_evGW > 1) THEN
1831 IF (unit_nr > 0) THEN
1832 CALL cp_warn(__location__, &
1833 "BSE@evGW applies W0, i.e. screening with DFT energies to the BSE!")
1834 END IF
1835 END IF
1836 ! Create a copy of fm_mat_S for usage in BSE
1837 CALL cp_fm_create(fm_mat_s_ia_bse, fm_mat_s(1)%matrix_struct)
1838 CALL cp_fm_to_fm(fm_mat_s(1), fm_mat_s_ia_bse)
1839 ! Remove energy/frequency factor from 3c-Integral for BSE
1840 IF (iter_sc_gw0 == 1) THEN
1841 CALL remove_scaling_factor_rpa(fm_mat_s_ia_bse, virtual(1), &
1842 eigenval_last(:, 1, 1), homo(1), omega)
1843 ELSE
1844 CALL remove_scaling_factor_rpa(fm_mat_s_ia_bse, virtual(1), &
1845 eigenval_scf(:, 1, 1), homo(1), omega)
1846 END IF
1847 ! Main routine for all BSE postprocessing
1848 CALL start_bse_calculation(fm_mat_s_ia_bse, fm_mat_s_ij_bse, fm_mat_s_ab_bse, &
1849 fm_mat_q_static_bse_gemm, &
1850 eigenval, eigenval_scf, &
1851 homo, virtual, dimen_ri, dimen_ri_red, bse_lev_virt, &
1852 gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
1853 ! Release BSE-copy of fm_mat_S
1854 CALL cp_fm_release(fm_mat_s_ia_bse)
1855 END IF
1856
1857 IF (my_do_gw) THEN
1858 CALL deallocate_matrices_gw(fm_mat_s_gw_work, vec_w_gw, vec_sigma_c_gw, vec_omega_fit_gw, &
1859 mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw, &
1860 eigenval_last, eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
1861 kpoints, vec_sigma_x_gw,.NOT. do_im_time)
1862 END IF
1863
1864 IF (do_im_time) THEN
1865
1866 CALL dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, &
1867 fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, index_to_cell_3c, &
1868 cell_to_index_3c, do_ic_model, &
1869 do_kpoints_cubic_rpa, do_kpoints_from_gamma, do_ri_sigma_x, &
1870 has_mat_p_blocks, &
1871 wkp_w, cfm_mat_q, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, &
1872 fm_matrix_minv, fm_matrix_minv_vtrunc_minv, fm_mat_ri_global_work, fm_mat_work, &
1873 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_l, &
1874 mat_minvvminv, mat_p_omega, mat_p_omega_kp, &
1875 t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, mat_work, qs_env)
1876
1877 IF (my_do_gw) THEN
1878 CALL deallocate_matrices_gw_im_time(weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, do_ic_model, &
1879 do_kpoints_cubic_rpa, fm_mat_w, &
1880 t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, &
1881 t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, &
1882 t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, &
1883 mat_w, qs_env)
1884 END IF
1885
1886 END IF
1887
1888 IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) CALL exchange_work%release()
1889
1890 IF (.NOT. do_ri_sos_laplace_mp2) THEN
1891 DEALLOCATE (tj)
1892 DEALLOCATE (wj)
1893 DEALLOCATE (trace_qomega)
1894 END IF
1895
1896 IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
1897 DEALLOCATE (tau_tj)
1898 DEALLOCATE (tau_wj)
1899 END IF
1900
1901 IF (do_im_time .AND. calc_forces) THEN
1902 CALL im_time_force_release(force_data)
1903 END IF
1904
1905 IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_finalize(rpa_grad, mp2_env, para_env_sub, para_env, &
1906 qs_env, gd_array, color_sub, do_ri_sos_laplace_mp2, &
1907 homo, virtual)
1908
1909 CALL timestop(handle)
1910
1911 END SUBROUTINE rpa_num_int
1912
1913END MODULE rpa_main
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
group_dist_types::create_group_dist
Definition group_dist_types.F:39
group_dist_types::get_group_dist
Definition group_dist_types.F:47
group_dist_types::release_group_dist
Definition group_dist_types.F:53
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::wilhelm2016b
integer, save, public wilhelm2016b
Definition bibliography.F:36
bibliography::delben2015
integer, save, public delben2015
Definition bibliography.F:36
bibliography::wilhelm2016a
integer, save, public wilhelm2016a
Definition bibliography.F:36
bibliography::wilhelm2018
integer, save, public wilhelm2018
Definition bibliography.F:36
bibliography::delben2013
integer, save, public delben2013
Definition bibliography.F:36
bibliography::ren2013
integer, save, public ren2013
Definition bibliography.F:36
bibliography::wilhelm2017
integer, save, public wilhelm2017
Definition bibliography.F:36
bibliography::ren2011
integer, save, public ren2011
Definition bibliography.F:36
bibliography::gruneis2009
integer, save, public gruneis2009
Definition bibliography.F:36
bibliography::freeman1977
integer, save, public freeman1977
Definition bibliography.F:36
bibliography::bates2013
integer, save, public bates2013
Definition bibliography.F:36
bse_main
Main routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_main.F:14
bse_main::start_bse_calculation
subroutine, public start_bse_calculation(fm_mat_s_ia_bse, fm_mat_s_ij_bse, fm_mat_s_ab_bse, fm_mat_q_static_bse_gemm, eigenval, eigenval_scf, homo, virtual, dimen_ri, dimen_ri_red, bse_lev_virt, gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
Main subroutine managing BSE calculations.
Definition bse_main.F:83
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:794
cp_dbcsr_api::dbcsr_clear
subroutine, public dbcsr_clear(matrix)
...
Definition cp_dbcsr_api.F:302
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:167
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:528
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
dgemm_counter_types
Counters to determine the performance of parallel DGEMMs.
Definition dgemm_counter_types.F:14
dgemm_counter_types::dgemm_counter_init
elemental subroutine, public dgemm_counter_init(dgemm_counter, unit_nr, print_info)
Initialize a dgemm_counter.
Definition dgemm_counter_types.F:49
dgemm_counter_types::dgemm_counter_write
subroutine, public dgemm_counter_write(dgemm_counter, para_env)
calculate and print flop rates
Definition dgemm_counter_types.F:103
group_dist_types
Types to describe group distributions.
Definition group_dist_types.F:14
group_dist_types::maxsize
elemental integer function, public maxsize(this)
...
Definition group_dist_types.F:270
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::rpa_exchange_none
integer, parameter, public rpa_exchange_none
Definition input_constants.F:1070
input_constants::sigma_none
integer, parameter, public sigma_none
Definition input_constants.F:1334
input_constants::rpa_exchange_sosex
integer, parameter, public rpa_exchange_sosex
Definition input_constants.F:1070
input_constants::wfc_mm_style_gemm
integer, parameter, public wfc_mm_style_gemm
Definition input_constants.F:1059
input_constants::rpa_exchange_axk
integer, parameter, public rpa_exchange_axk
Definition input_constants.F:1070
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
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
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:440
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:130
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
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
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
Definition minimax_exp.F:29
minimax_exp::check_exp_minimax_range
subroutine, public check_exp_minimax_range(k, rc, ierr)
Check that a minimax approximation is available for given input k, Rc. ierr == 0: everything ok ierr ...
Definition minimax_exp.F:101
mp2_grids
Routines to calculate frequency and time grids (integration points and weights) for correlation metho...
Definition mp2_grids.F:14
mp2_grids::init_greenx_grids
subroutine, public init_greenx_grids(qs_env, para_env, homo, eigenval, do_ri_sos_laplace_mp2, do_kpoints_cubic_rpa, emin, emax, e_fermi)
returns minimal and maximal energy values for the E_range for the minimax grid selection
Definition mp2_grids.F:1545
mp2_grids::get_minimax_grid
subroutine, public get_minimax_grid(para_env, unit_nr, homo, eigenval, num_integ_points, do_im_time, do_ri_sos_laplace_mp2, do_print, tau_tj, tau_wj, qs_env, do_gw_im_time, do_kpoints_cubic_rpa, e_fermi, tj, wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, regularization)
...
Definition mp2_grids.F:75
mp2_grids::get_clenshaw_grid
subroutine, public get_clenshaw_grid(para_env, para_env_rpa, unit_nr, homo, virtual, eigenval, num_integ_points, num_integ_group, color_rpa_group, fm_mat_s, my_do_gw, ext_scaling, a_scaling, tj, wj)
...
Definition mp2_grids.F:319
mp2_laplace
Routines to calculate MP2 energy with laplace approach.
Definition mp2_laplace.F:13
mp2_laplace::sos_mp2_postprocessing
subroutine, public sos_mp2_postprocessing(fm_mat_q, erpa, tau_wjquad)
...
Definition mp2_laplace.F:81
mp2_ri_grad_util
Routines for calculating RI-MP2 gradients.
Definition mp2_ri_grad_util.F:13
mp2_ri_grad_util::array2fm
subroutine, public array2fm(mat2d, fm_struct, my_start_row, my_end_row, my_start_col, my_end_col, gd_row, gd_col, group_grid_2_mepos, ngroup_row, ngroup_col, fm_mat, integ_group_size, color_group, do_release_mat)
redistribute local part of array to fm
Definition mp2_ri_grad_util.F:498
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
rpa_exchange
Auxiliary routines needed for RPA-exchange given blacs_env to another.
Definition rpa_exchange.F:18
rpa_exchange::rpa_exchange_needed_mem
subroutine, public rpa_exchange_needed_mem(mp2_env, homo, virtual, dimen_ri, para_env, mem_per_rank, mem_per_repl)
...
Definition rpa_exchange.F:146
rpa_grad
Routines to calculate RI-RPA and SOS-MP2 gradients.
Definition rpa_grad.F:13
rpa_grad::rpa_grad_finalize
subroutine, public rpa_grad_finalize(rpa_grad, mp2_env, para_env_sub, para_env, qs_env, gd_array, color_sub, do_ri_sos_laplace_mp2, homo, virtual)
...
Definition rpa_grad.F:2029
rpa_grad::rpa_grad_copy_q
subroutine, public rpa_grad_copy_q(fm_mat_q, rpa_grad)
...
Definition rpa_grad.F:713
rpa_grad::rpa_grad_needed_mem
pure subroutine, public rpa_grad_needed_mem(homo, virtual, dimen_ri, mem_per_rank, mem_per_repl, do_ri_sos_laplace_mp2)
Calculates the necessary minimum memory for the Gradient code ion MiB.
Definition rpa_grad.F:122
rpa_grad::rpa_grad_create
subroutine, public rpa_grad_create(rpa_grad, fm_mat_q, fm_mat_s, homo, virtual, mp2_env, eigenval, unit_nr, do_ri_sos_laplace_mp2)
Creates the arrays of a rpa_grad_type.
Definition rpa_grad.F:168
rpa_grad::rpa_grad_matrix_operations
subroutine, public rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, fm_mat_q, fm_mat_q_gemm, dgemm_counter, fm_mat_s, omega, homo, virtual, eigenval, weight, unit_nr)
...
Definition rpa_grad.F:738
rpa_gw_ic
Routines to calculate image charge corrections.
Definition rpa_gw_ic.F:13
rpa_gw_ic::calculate_ic_correction
subroutine, public calculate_ic_correction(eigenval, mat_sinvvsinv, t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, unit_nr, print_ic_values, para_env, do_alpha, do_beta)
...
Definition rpa_gw_ic.F:58
rpa_gw_kpoints_util
Routines treating GW and RPA calculations with kpoints.
Definition rpa_gw_kpoints_util.F:13
rpa_gw_kpoints_util::invert_eps_compute_w_and_erpa_kp
subroutine, public invert_eps_compute_w_and_erpa_kp(dimen_ri, num_integ_points, jquad, nkp, count_ev_sc_gw, para_env, erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, wkp_w, do_gw_im_time, do_ri_sigma_x, do_kpoints_from_gamma, cfm_mat_q, ikp_local, mat_p_omega, mat_p_omega_kp, qs_env, eps_filter_im_time, unit_nr, kpoints, fm_mat_minv_l_kpoints, fm_matrix_l_kpoints, fm_mat_w, fm_mat_ri_global_work, mat_minvvminv, fm_matrix_minv, fm_matrix_minv_vtrunc_minv)
...
Definition rpa_gw_kpoints_util.F:135
rpa_gw_kpoints_util::get_bandstruc_and_k_dependent_mos
subroutine, public get_bandstruc_and_k_dependent_mos(qs_env, eigenval_kp)
...
Definition rpa_gw_kpoints_util.F:1428
rpa_gw
Routines for GW, continuous development [Jan Wilhelm].
Definition rpa_gw.F:14
rpa_gw::deallocate_matrices_gw
subroutine, public deallocate_matrices_gw(fm_mat_s_gw_work, vec_w_gw, vec_sigma_c_gw, vec_omega_fit_gw, vec_sigma_x_minus_vxc_gw, eigenval_last, eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, kpoints, vec_sigma_x_gw, my_do_gw)
...
Definition rpa_gw.F:615
rpa_gw::allocate_matrices_gw
subroutine, public allocate_matrices_gw(vec_sigma_c_gw, color_rpa_group, dimen_nm_gw, gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, num_integ_group, num_integ_points, unit_nr, gw_corr_lev_tot, num_fit_points, omega_max_fit, do_minimax_quad, do_periodic, do_ri_sigma_x, my_do_gw, first_cycle_periodic_correction, a_scaling, eigenval, tj, vec_omega_fit_gw, vec_sigma_x_gw, delta_corr, eigenval_last, eigenval_scf, vec_w_gw, fm_mat_s_gw, fm_mat_s_gw_work, para_env, mp2_env, kpoints, nkp, nkp_self_energy, do_kpoints_cubic_rpa, do_kpoints_from_gamma)
...
Definition rpa_gw.F:350
rpa_gw::compute_gw_self_energy
subroutine, public compute_gw_self_energy(vec_sigma_c_gw, dimen_nm_gw, dimen_ri, gw_corr_lev_occ, gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, do_im_time, do_periodic, first_cycle_periodic_correction, fermi_level_offset, omega, eigenval, delta_corr, vec_omega_fit_gw, vec_w_gw, wj, fm_mat_q, fm_mat_r_gw, fm_mat_s_gw, fm_mat_s_gw_work, mo_coeff, para_env, para_env_rpa, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, kpoints, qs_env, mp2_env)
...
Definition rpa_gw.F:809
rpa_gw::get_fermi_level_offset
subroutine, public get_fermi_level_offset(fermi_level_offset, fermi_level_offset_input, eigenval, homo)
...
Definition rpa_gw.F:913
rpa_gw::allocate_matrices_gw_im_time
subroutine, public allocate_matrices_gw_im_time(gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, num_integ_points, unit_nr, ri_blk_sizes, do_ic_model, para_env, fm_mat_w, fm_mat_q, mo_coeff, t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, starts_array_mc, ends_array_mc, t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, matrix_s, mat_w, t_3c_overl_int, t_3c_o_compressed, t_3c_o_ind, qs_env)
...
Definition rpa_gw.F:204
rpa_gw::deallocate_matrices_gw_im_time
subroutine, public deallocate_matrices_gw_im_time(weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, do_ic_model, do_kpoints_cubic_rpa, fm_mat_w, t_3c_overl_int_ao_mo, t_3c_o_mo_compressed, t_3c_o_mo_ind, t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, mat_w, qs_env)
...
Definition rpa_gw.F:690
rpa_gw::compute_qp_energies
subroutine, public compute_qp_energies(vec_sigma_c_gw, count_ev_sc_gw, gw_corr_lev_occ, gw_corr_lev_tot, gw_corr_lev_virt, homo, nmo, num_fit_points, num_integ_points, unit_nr, do_apply_ic_corr_to_gw, do_im_time, do_periodic, do_ri_sigma_x, first_cycle_periodic_correction, e_fermi, eps_filter, fermi_level_offset, delta_corr, eigenval, eigenval_last, eigenval_scf, iter_sc_gw0, exit_ev_gw, tau_tj, tj, vec_omega_fit_gw, vec_sigma_x_gw, ic_corr_list, weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, mo_coeff, fm_mat_w, para_env, para_env_rpa, mat_dm, mat_minvvminv, t_3c_o, t_3c_m, t_3c_overl_int_ao_mo, t_3c_o_compressed, t_3c_o_mo_compressed, t_3c_o_ind, t_3c_o_mo_ind, t_3c_overl_int_gw_ri, t_3c_overl_int_gw_ao, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, mat_w, matrix_s, kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_rpa, starts_array_mc, ends_array_mc)
...
Definition rpa_gw.F:1234
rpa_gw::compute_w_cubic_gw
subroutine, public compute_w_cubic_gw(fm_mat_w, fm_mat_q, fm_mat_work, dimen_ri, fm_mat_l, num_integ_points, tj, tau_tj, weights_cos_tf_w_to_t, jquad, omega)
...
Definition rpa_gw.F:955
rpa_im_time_force_methods
Routines needed for cubic-scaling RPA and SOS-Laplace-MP2 forces.
Definition rpa_im_time_force_methods.F:12
rpa_im_time_force_methods::calc_laplace_loop_forces
subroutine, public calc_laplace_loop_forces(force_data, mat_p_omega, t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, num_integ_points, nmo, eigenval, tau_tj, tau_wj, cut_memory, pspin, qspin, open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
Updates the cubic-scaling SOS-Laplace-MP2 contribution to the forces at each quadrature point.
Definition rpa_im_time_force_methods.F:669
rpa_im_time_force_methods::calc_post_loop_forces
subroutine, public calc_post_loop_forces(force_data, unit_nr, qs_env)
All the forces that can be calculated after the loop on the Laplace quaradture, using terms collected...
Definition rpa_im_time_force_methods.F:2213
rpa_im_time_force_methods::calc_rpa_loop_forces
subroutine, public calc_rpa_loop_forces(force_data, mat_p_omega, t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, num_integ_points, nmo, eigenval, e_fermi, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, tj, wj, tau_tj, cut_memory, ispin, open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
Updates the cubic-scaling RPA contribution to the forces at each quadrature point....
Definition rpa_im_time_force_methods.F:1220
rpa_im_time_force_methods::keep_initial_quad
subroutine, public keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, do_laplace, do_im_time, num_integ_points, unit_nr, qs_env)
Overwrites the "optimal" Laplace quadrature with that of the first step.
Definition rpa_im_time_force_methods.F:3931
rpa_im_time_force_methods::init_im_time_forces
subroutine, public init_im_time_forces(force_data, fm_matrix_pq, t_3c_m, unit_nr, mp2_env, qs_env)
Initializes and pre-calculates all needed tensors for the forces.
Definition rpa_im_time_force_methods.F:190
rpa_im_time_force_types
Types needed for cubic-scaling RPA and SOS-Laplace-MP2 forces.
Definition rpa_im_time_force_types.F:12
rpa_im_time_force_types::im_time_force_release
subroutine, public im_time_force_release(force_data)
Cleans everything up.
Definition rpa_im_time_force_types.F:91
rpa_im_time
Routines for low-scaling RPA/GW with imaginary time.
Definition rpa_im_time.F:13
rpa_im_time::zero_mat_p_omega
subroutine, public zero_mat_p_omega(mat_p_omega)
...
Definition rpa_im_time.F:615
rpa_im_time::compute_mat_p_omega
subroutine, public compute_mat_p_omega(mat_p_omega, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_p_global, matrix_s, ispin, t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, weights_cos_tf_t_to_w, tj, tau_tj, e_fermi, eps_filter, alpha, eps_filter_im_time, eigenval, nmo, num_integ_points, cut_memory, unit_nr, mp2_env, para_env, qs_env, do_kpoints_from_gamma, index_to_cell_3c, cell_to_index_3c, has_mat_p_blocks, do_ri_sos_laplace_mp2, dbcsr_time, dbcsr_nflop)
...
Definition rpa_im_time.F:139
rpa_main
Routines to calculate RI-RPA energy.
Definition rpa_main.F:17
rpa_main::rpa_ri_compute_en
subroutine, public rpa_ri_compute_en(qs_env, erpa, mp2_env, bib_c, bib_c_gw, bib_c_bse_ij, bib_c_bse_ab, para_env, para_env_sub, color_sub, gd_array, gd_b_virtual, gd_b_all, gd_b_occ_bse, gd_b_virt_bse, mo_coeff, fm_matrix_pq, fm_matrix_l_kpoints, fm_matrix_minv_l_kpoints, fm_matrix_minv, fm_matrix_minv_vtrunc_minv, kpoints, eigenval, nmo, homo, dimen_ri, dimen_ri_red, gw_corr_lev_occ, gw_corr_lev_virt, bse_lev_virt, unit_nr, do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_bse, matrix_s, mat_munu, mat_p_global, t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, calc_forces)
...
Definition rpa_main.F:195
rpa_sigma_functional
Routines to calculate RI-RPA energy and Sigma correction to the RPA energies using the cubic spline b...
Definition rpa_sigma_functional.F:13
rpa_sigma_functional::rpa_sigma_create
subroutine, public rpa_sigma_create(rpa_sigma, sigma_param, fm_mat_q, unit_nr, para_env)
... Collect the Q(w) (fm_mat_Q) matrix to create rpa_sigma a derived type variable....
Definition rpa_sigma_functional.F:64
rpa_sigma_functional::finalize_rpa_sigma
subroutine, public finalize_rpa_sigma(rpa_sigma, unit_nr, e_sigma_corr, para_env, do_minimax_quad)
... Save the calculated value of E_c correction to the global variable and memory clean.
Definition rpa_sigma_functional.F:195
rpa_sigma_functional::rpa_sigma_matrix_spectral
subroutine, public rpa_sigma_matrix_spectral(rpa_sigma, fm_mat_q, wj, para_env_rpa)
... Diagonalize and store the eigenvalues of fm_mat_Q in rpa_sigmasigma_eigenvalue.
Definition rpa_sigma_functional.F:134
rpa_util
Utility functions for RPA calculations.
Definition rpa_util.F:13
rpa_util::compute_erpa_by_freq_int
subroutine, public compute_erpa_by_freq_int(dimen_ri, trace_qomega, fm_mat_q, para_env_rpa, erpa, wjquad)
...
Definition rpa_util.F:929
rpa_util::alloc_im_time
subroutine, public alloc_im_time(qs_env, para_env, dimen_ri, dimen_ri_red, num_integ_points, nspins, fm_mat_q, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_matrix_minv_l_kpoints, fm_matrix_l_kpoints, mat_p_global, t_3c_o, matrix_s, kpoints, eps_filter_im_time, cut_memory, nkp, num_cells_dm, num_3c_repl, size_p, ikp_local, index_to_cell_3c, cell_to_index_3c, col_blk_size, do_ic_model, do_kpoints_cubic_rpa, do_kpoints_from_gamma, do_ri_sigma_x, my_open_shell, has_mat_p_blocks, wkp_w, cfm_mat_q, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, fm_mat_ri_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_l, mat_m_p_munu_occ, mat_m_p_munu_virt, mat_minvvminv, mat_p_omega, mat_p_omega_kp, mat_work, mo_coeff, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, homo, nmo)
...
Definition rpa_util.F:147
rpa_util::dealloc_im_time
subroutine, public dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, index_to_cell_3c, cell_to_index_3c, do_ic_model, do_kpoints_cubic_rpa, do_kpoints_from_gamma, do_ri_sigma_x, has_mat_p_blocks, wkp_w, cfm_mat_q, fm_mat_minv_l_kpoints, fm_mat_l_kpoints, fm_matrix_minv, fm_matrix_minv_vtrunc_minv, fm_mat_ri_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_l, mat_minvvminv, mat_p_omega, mat_p_omega_kp, t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, mat_work, qs_env)
...
Definition rpa_util.F:1116
rpa_util::q_trace_and_add_unit_matrix
subroutine, public q_trace_and_add_unit_matrix(dimen_ri, trace_qomega, fm_mat_q)
...
Definition rpa_util.F:878
rpa_util::calc_mat_q
subroutine, public calc_mat_q(fm_mat_s, do_ri_sos_laplace_mp2, first_cycle, virtual, eigenval, homo, omega, omega_old, jquad, mm_style, dimen_ri, dimen_ia, alpha, fm_mat_q, fm_mat_q_gemm, do_bse, fm_mat_q_static_bse_gemm, dgemm_counter, num_integ_points, count_ev_sc_gw)
...
Definition rpa_util.F:661
rpa_util::contract_p_omega_with_mat_l
subroutine, public contract_p_omega_with_mat_l(mat_p_omega, mat_l, mat_work, eps_filter_im_time, fm_mat_work, dimen_ri, dimen_ri_red, fm_mat_l, fm_mat_q)
...
Definition rpa_util.F:1266
rpa_util::remove_scaling_factor_rpa
subroutine, public remove_scaling_factor_rpa(fm_mat_s, virtual, eigenval_last, homo, omega_old)
...
Definition rpa_util.F:711
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
dgemm_counter_types::dgemm_counter_type
Definition dgemm_counter_types.F:32
group_dist_types::group_dist_d1_type
Definition group_dist_types.F:35
hfx_types::block_ind_type
Definition hfx_types.F:367
hfx_types::hfx_compression_type
Definition hfx_types.F:351
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:150
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
mp2_types::mp2_type
Definition mp2_types.F:345
mp2_types::three_dim_real_array
Definition mp2_types.F:87
mp2_types::two_dim_int_array
Definition mp2_types.F:75
mp2_types::two_dim_real_array
Definition mp2_types.F:83
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
rpa_exchange::rpa_exchange_work_type
Definition rpa_exchange.F:109
rpa_grad::rpa_grad_type
Definition rpa_grad.F:98
rpa_im_time_force_types::im_time_force_type
Definition rpa_im_time_force_types.F:37
rpa_sigma_functional::rpa_sigma_type
Definition rpa_sigma_functional.F:41