Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to calculate 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 : ! **************************************************************************************************
17 : MODULE 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 :
132 : CONTAINS
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 1470 : 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 294 : gd_array, gd_B_virtual, gd_B_all, gd_B_occ_bse, gd_B_virt_bse, &
187 294 : 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 588 : 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 588 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dimen_ia, my_ia_end, my_ia_size, &
250 294 : 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 294 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_kp
256 588 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_Q, fm_mat_Q_gemm, fm_mat_S, &
257 294 : fm_mat_S_gw
258 882 : TYPE(cp_fm_type), DIMENSION(1) :: fm_mat_R_gw, fm_mat_S_ab_bse, &
259 588 : fm_mat_S_ij_bse
260 : TYPE(mp_para_env_type), POINTER :: para_env_RPA
261 : TYPE(two_dim_real_array), ALLOCATABLE, &
262 588 : DIMENSION(:) :: BIb_C_2D, BIb_C_2D_gw
263 2940 : TYPE(two_dim_real_array), DIMENSION(1) :: BIb_C_2D_bse_ab, BIb_C_2D_bse_ij
264 :
265 294 : CALL timeset(routineN, handle)
266 :
267 294 : CALL cite_reference(DelBen2013)
268 294 : CALL cite_reference(DelBen2015)
269 :
270 294 : IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_axk) THEN
271 10 : CALL cite_reference(Bates2013)
272 284 : ELSE IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_sosex) THEN
273 2 : CALL cite_reference(Freeman1977)
274 2 : CALL cite_reference(Gruneis2009)
275 : END IF
276 294 : IF (mp2_env%ri_rpa%do_rse) THEN
277 4 : CALL cite_reference(Ren2011)
278 4 : CALL cite_reference(Ren2013)
279 : END IF
280 :
281 294 : IF (my_do_gw) THEN
282 104 : CALL cite_reference(Wilhelm2016a)
283 104 : CALL cite_reference(Wilhelm2017)
284 104 : CALL cite_reference(Wilhelm2018)
285 : END IF
286 :
287 294 : IF (do_im_time) THEN
288 134 : CALL cite_reference(Wilhelm2016b)
289 : END IF
290 :
291 294 : nspins = SIZE(homo)
292 294 : my_open_shell = (nspins == 2)
293 2058 : ALLOCATE (virtual(nspins), dimen_ia(nspins), my_ia_end(nspins), my_ia_start(nspins), my_ia_size(nspins))
294 644 : virtual(:) = nmo - homo(:)
295 644 : dimen_ia(:) = virtual(:)*homo(:)
296 :
297 1176 : ALLOCATE (Eigenval_kp(nmo, 1, nspins))
298 8842 : Eigenval_kp(:, 1, :) = Eigenval(:, :)
299 :
300 294 : IF (do_im_time) mp2_env%ri_rpa%minimax_quad = .TRUE.
301 294 : do_minimax_quad = mp2_env%ri_rpa%minimax_quad
302 :
303 294 : IF (do_ri_sos_laplace_mp2) THEN
304 58 : num_integ_points = mp2_env%ri_laplace%n_quadrature
305 58 : input_num_integ_groups = mp2_env%ri_laplace%num_integ_groups
306 :
307 : ! check the range for the minimax approximation
308 58 : E_Range = mp2_env%e_range
309 58 : 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 88 : DO ispin = 1, nspins
313 88 : IF (homo(ispin) > 0) THEN
314 50 : Emin = MIN(Emin, 2.0_dp*(Eigenval(homo(ispin) + 1, ispin) - Eigenval(homo(ispin), ispin)))
315 2740 : Emax = MAX(Emax, 2.0_dp*(MAXVAL(Eigenval(:, ispin)) - MINVAL(Eigenval(:, ispin))))
316 : END IF
317 : END DO
318 38 : E_Range = Emax/Emin
319 : END IF
320 58 : IF (E_Range < 2.0_dp) E_Range = 2.0_dp
321 : ierr = 0
322 58 : CALL check_exp_minimax_range(num_integ_points, E_Range, ierr)
323 58 : IF (ierr /= 0) THEN
324 : jjB = num_integ_points - 1
325 0 : DO iiB = 1, jjB
326 0 : num_integ_points = num_integ_points - 1
327 : ierr = 0
328 0 : CALL check_exp_minimax_range(num_integ_points, E_Range, ierr)
329 0 : IF (ierr == 0) EXIT
330 : END DO
331 : END IF
332 58 : CPASSERT(num_integ_points >= 1)
333 : ELSE
334 236 : num_integ_points = mp2_env%ri_rpa%rpa_num_quad_points
335 236 : input_num_integ_groups = mp2_env%ri_rpa%rpa_num_integ_groups
336 236 : IF (my_do_gw .AND. do_minimax_quad) THEN
337 46 : IF (num_integ_points > 34) THEN
338 0 : IF (unit_nr > 0) &
339 : CALL cp_warn(__LOCATION__, &
340 : "The required number of quadrature point exceeds the maximum possible in the "// &
341 0 : "Minimax quadrature scheme. The number of quadrature point has been reset to 30.")
342 0 : num_integ_points = 30
343 : END IF
344 : ELSE
345 190 : IF (do_minimax_quad .AND. num_integ_points > 20) THEN
346 0 : IF (unit_nr > 0) &
347 : CALL cp_warn(__LOCATION__, &
348 : "The required number of quadrature point exceeds the maximum possible in the "// &
349 0 : "Minimax quadrature scheme. The number of quadrature point has been reset to 20.")
350 0 : num_integ_points = 20
351 : END IF
352 : END IF
353 : END IF
354 294 : allowed_memory = mp2_env%mp2_memory
355 :
356 294 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
357 :
358 294 : 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 294 : IF (do_im_time .OR. mp2_env%ri_g0w0%do_periodic .OR. do_bse) THEN
362 :
363 166 : integ_group_size = ngroup
364 166 : 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 284 : mem_for_iaK = REAL(SUM(dimen_ia), KIND=dp)*dimen_RI_red*8.0_dp/(1024_dp**2)
371 128 : mem_for_QK = REAL(dimen_RI_red, KIND=dp)*nspins*dimen_RI_red*8.0_dp/(1024_dp**2)
372 :
373 128 : CALL m_memory(mem)
374 128 : mem_real = (mem + 1024*1024 - 1)/(1024*1024)
375 128 : CALL para_env%min(mem_real)
376 :
377 128 : 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 128 : mem_per_repl = mem_per_repl + 2.0_dp*mem_for_QK
383 :
384 128 : 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 128 : CALL rpa_exchange_needed_mem(mp2_env, homo, virtual, dimen_RI_red, para_env, mem_per_rank, mem_per_repl)
386 :
387 128 : mem_min = mem_per_repl/para_env%num_pe + mem_per_rank
388 :
389 128 : IF (unit_nr > 0) THEN
390 64 : WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum required memory per MPI process:', mem_min, ' MiB'
391 64 : 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 128 : 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 128 : mem_real = mem_real - mem_per_rank
398 :
399 128 : mem_per_group = mem_real*para_env_sub%num_pe
400 :
401 : ! here we try to find the best rpa/laplace group size
402 128 : skip_integ_group_opt = .FALSE.
403 :
404 : ! Check the input number of integration groups
405 128 : IF (input_num_integ_groups > 0) THEN
406 2 : IF (MOD(num_integ_points, input_num_integ_groups) == 0) THEN
407 0 : IF (MOD(ngroup, input_num_integ_groups) == 0) THEN
408 0 : best_integ_group_size = ngroup/input_num_integ_groups
409 0 : best_num_integ_point = num_integ_points/input_num_integ_groups
410 0 : skip_integ_group_opt = .TRUE.
411 : ELSE
412 0 : 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 2 : 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 128 : IF (.NOT. skip_integ_group_opt) THEN
420 128 : best_integ_group_size = ngroup
421 128 : best_num_integ_point = num_integ_points
422 :
423 128 : min_integ_group_size = MAX(1, ngroup/num_integ_points)
424 :
425 128 : integ_group_size = min_integ_group_size - 1
426 142 : DO iiB = min_integ_group_size + 1, ngroup
427 106 : integ_group_size = integ_group_size + 1
428 :
429 : ! check that the ngroup is a multiple of integ_group_size
430 106 : IF (MOD(ngroup, integ_group_size) /= 0) CYCLE
431 :
432 : ! check for memory
433 106 : avail_mem = integ_group_size*mem_per_group
434 106 : IF (avail_mem < mem_per_repl) CYCLE
435 :
436 : ! check that the integration groups have the same size
437 106 : num_integ_group = ngroup/integ_group_size
438 106 : IF (MOD(num_integ_points, num_integ_group) /= 0) CYCLE
439 :
440 92 : best_num_integ_point = num_integ_points/num_integ_group
441 92 : best_integ_group_size = integ_group_size
442 :
443 142 : EXIT
444 :
445 : END DO
446 : END IF
447 :
448 128 : integ_group_size = best_integ_group_size
449 :
450 : END IF
451 :
452 294 : IF (unit_nr > 0 .AND. .NOT. do_im_time) THEN
453 80 : IF (do_ri_sos_laplace_mp2) THEN
454 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
455 14 : "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 14 : "INTEG_INFO| MINIMAX approximation"
458 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
459 14 : "INTEG_INFO| Number of integration points:", num_integ_points
460 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
461 14 : "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 66 : "RI_INFO| Group size for frequency integration:", integ_group_size*para_env_sub%num_pe
465 66 : IF (do_minimax_quad) THEN
466 : WRITE (UNIT=unit_nr, FMT="(T3,A)") &
467 19 : "INTEG_INFO| MINIMAX quadrature"
468 : ELSE
469 : WRITE (UNIT=unit_nr, FMT="(T3,A)") &
470 47 : "INTEG_INFO| Clenshaw-Curtius quadrature"
471 : END IF
472 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
473 66 : "INTEG_INFO| Number of integration points:", num_integ_points
474 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
475 66 : "INTEG_INFO| Number of integration points per RPA group:", best_num_integ_point
476 : END IF
477 80 : CALL m_flush(unit_nr)
478 : END IF
479 :
480 294 : num_integ_group = ngroup/integ_group_size
481 :
482 294 : pos_integ_group = MOD(color_sub, integ_group_size)
483 294 : color_rpa_group = color_sub/integ_group_size
484 :
485 294 : CALL timeset(routineN//"_reorder", handle2)
486 :
487 : ! not necessary for imaginary time
488 :
489 1232 : ALLOCATE (BIb_C_2D(nspins))
490 :
491 294 : 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 348 : 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 188 : my_ia_size(ispin), my_ia_start(ispin), my_ia_end(ispin), my_group_L_size)
499 :
500 188 : DEALLOCATE (BIb_C(ispin)%array)
501 536 : 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 160 : IF (my_do_gw) THEN
507 178 : ALLOCATE (BIb_C_2D_gw(nspins))
508 :
509 58 : CALL timeset(routineN//"_reorder_gw", handle3)
510 :
511 58 : dimen_nm_gw = nmo*(gw_corr_lev_occ(1) + gw_corr_lev_virt(1))
512 :
513 : ! The same for open shell
514 120 : 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 62 : my_nm_gw_size, my_nm_gw_start, my_nm_gw_end, my_group_L_size)
518 120 : DEALLOCATE (BIb_C_gw(ispin)%array)
519 : END DO
520 :
521 58 : CALL release_group_dist(gd_B_all)
522 :
523 116 : CALL timestop(handle3)
524 :
525 : END IF
526 : END IF
527 :
528 354 : IF (do_bse) THEN
529 :
530 30 : CALL timeset(routineN//"_reorder_bse1", handle3)
531 :
532 30 : 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 30 : my_ij_comb_bse_size, my_ij_comb_bse_start, my_ij_comb_bse_end, my_group_L_size)
538 :
539 30 : DEALLOCATE (BIb_C_bse_ij)
540 30 : CALL release_group_dist(gd_B_occ_bse)
541 :
542 30 : CALL timestop(handle3)
543 :
544 30 : CALL timeset(routineN//"_reorder_bse2", handle3)
545 :
546 30 : 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 30 : my_ab_comb_bse_size, my_ab_comb_bse_start, my_ab_comb_bse_end, my_group_L_size)
551 :
552 30 : DEALLOCATE (BIb_C_bse_ab)
553 30 : CALL release_group_dist(gd_B_virt_bse)
554 :
555 30 : CALL timestop(handle3)
556 :
557 : END IF
558 :
559 294 : CALL timestop(handle2)
560 :
561 294 : IF (num_integ_group > 1) THEN
562 92 : ALLOCATE (para_env_RPA)
563 92 : CALL para_env_RPA%from_split(para_env, color_rpa_group)
564 : ELSE
565 202 : 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 294 : IF (do_im_time) THEN
572 832 : ALLOCATE (fm_mat_Q(nspins), fm_mat_Q_gemm(1), fm_mat_S(1))
573 : ELSE
574 1364 : 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 294 : 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 644 : 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 1232 : ALLOCATE (fm_mat_S_gw(nspins))
590 294 : 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 522 : fm_mat_Q=fm_mat_R_gw)
600 120 : DEALLOCATE (BIb_C_2D_gw)
601 :
602 : END IF
603 :
604 : ! for Bethe-Salpeter, we need other matrix fm_mat_S
605 294 : 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 150 : 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 150 : fm_mat_Q(1)%matrix_struct%context, fm_mat_Q(1)%matrix_struct%context)
621 :
622 : END IF
623 :
624 294 : do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
625 294 : IF (do_kpoints_from_Gamma) THEN
626 16 : 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 294 : do_ri_sos_laplace_mp2=do_ri_sos_laplace_mp2, calc_forces=calc_forces)
647 :
648 294 : CALL release_group_dist(gd_array)
649 :
650 294 : IF (num_integ_group > 1) CALL mp_para_env_release(para_env_RPA)
651 :
652 294 : IF (.NOT. do_im_time) THEN
653 160 : CALL cp_fm_release(fm_mat_Q_gemm)
654 160 : CALL cp_fm_release(fm_mat_S)
655 : END IF
656 294 : CALL cp_fm_release(fm_mat_Q)
657 :
658 294 : IF (my_do_gw .AND. .NOT. do_im_time) THEN
659 58 : CALL cp_fm_release(fm_mat_S_gw)
660 58 : CALL cp_fm_release(fm_mat_R_gw(1))
661 : END IF
662 :
663 294 : IF (do_bse) THEN
664 30 : CALL cp_fm_release(fm_mat_S_ij_bse(1))
665 30 : CALL cp_fm_release(fm_mat_S_ab_bse(1))
666 : END IF
667 :
668 294 : CALL timestop(handle)
669 :
670 1176 : 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 310 : 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 310 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: BIb_C_rec_1D
708 310 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: BIb_C_rec
709 :
710 310 : itmp = get_limit(dimen_ia, para_env_sub%num_pe, para_env_sub%mepos)
711 310 : my_ia_start = itmp(1)
712 310 : my_ia_end = itmp(2)
713 310 : my_ia_size = my_ia_end - my_ia_start + 1
714 :
715 310 : 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 1234 : 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 310 : !$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 310 : IF (para_env_sub%num_pe > 1) THEN
733 30 : ALLOCATE (BIb_C_rec_1D(INT(my_group_L_size, int_8)*maxsize(gd_B_virtual)*MIN(homo, occ_chunk)))
734 20 : DO proc_shift = 1, para_env_sub%num_pe - 1
735 10 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
736 10 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
737 :
738 10 : 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 20 : DO occ_low = 1, homo, occ_chunk
742 10 : 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 10 : 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 31970 : 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 10 : !$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 10 : DEALLOCATE (BIb_C_rec_1D)
762 : END IF
763 :
764 310 : 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 412 : SUBROUTINE create_integ_mat(BIb_C_2D, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
799 412 : dimen_RI, dimen_ia, color_rpa_group, &
800 : ext_row_block_size, ext_col_block_size, unit_nr, &
801 412 : my_ia_size, my_ia_start, my_ia_end, &
802 : my_group_L_size, my_group_L_start, my_group_L_end, &
803 412 : 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 412 : 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 412 : 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 412 : TYPE(group_dist_d1_type) :: gd_ia, gd_L
839 :
840 412 : CALL timeset(routineN, handle)
841 :
842 412 : CPASSERT(PRESENT(blacs_env_ext) .OR. PRESENT(dimen_ia_for_block_size))
843 :
844 412 : my_blacs_ext = .FALSE.
845 412 : IF (PRESENT(blacs_env_ext)) my_blacs_ext = .TRUE.
846 :
847 412 : my_blacs_S_ext = .FALSE.
848 412 : IF (PRESENT(blacs_env_ext_S)) my_blacs_S_ext = .TRUE.
849 :
850 412 : my_do_im_time = .FALSE.
851 412 : IF (PRESENT(do_im_time)) my_do_im_time = do_im_time
852 :
853 412 : NULLIFY (blacs_env)
854 : ! create the RPA blacs env
855 412 : IF (my_blacs_S_ext) THEN
856 118 : blacs_env => blacs_env_ext_S
857 : ELSE
858 294 : IF (para_env_RPA%num_pe > 1) THEN
859 202 : col_row_proc_ratio = MAX(1, dimen_ia_for_block_size/dimen_RI)
860 :
861 202 : 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 202 : DO iproc = 1, para_env_RPA%num_pe
863 202 : iproc_col = iproc_col - 1
864 202 : IF (MOD(para_env_RPA%num_pe, iproc_col) == 0) EXIT
865 : END DO
866 :
867 202 : iproc_row = para_env_RPA%num_pe/iproc_col
868 202 : grid_2D(1) = iproc_row
869 202 : grid_2D(2) = iproc_col
870 : ELSE
871 276 : grid_2D = 1
872 : END IF
873 294 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_RPA, grid_2d=grid_2D)
874 :
875 294 : IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
876 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
877 80 : "MATRIX_INFO| Number row processes:", grid_2D(1)
878 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
879 80 : "MATRIX_INFO| Number column processes:", grid_2D(2)
880 : END IF
881 :
882 : ! define the block_size for the row
883 294 : IF (ext_row_block_size > 0) THEN
884 0 : nrow_block_mat = ext_row_block_size
885 : ELSE
886 294 : nrow_block_mat = MAX(1, dimen_RI/grid_2D(1)/2)
887 : END IF
888 :
889 : ! define the block_size for the column
890 294 : IF (ext_col_block_size > 0) THEN
891 0 : ncol_block_mat = ext_col_block_size
892 : ELSE
893 294 : ncol_block_mat = MAX(1, dimen_ia_for_block_size/grid_2D(2)/2)
894 : END IF
895 :
896 294 : IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
897 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
898 80 : "MATRIX_INFO| Row block size:", nrow_block_mat
899 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
900 80 : "MATRIX_INFO| Column block size:", ncol_block_mat
901 : END IF
902 : END IF
903 :
904 345 : IF (.NOT. my_do_im_time) THEN
905 588 : DO ispin = 1, SIZE(BIb_C_2D)
906 310 : NULLIFY (fm_struct)
907 310 : IF (my_blacs_ext) THEN
908 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
909 122 : 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 188 : nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.TRUE.)
914 :
915 : END IF ! external blacs_env
916 :
917 310 : CALL create_group_dist(gd_ia, my_ia_start(ispin), my_ia_end(ispin), my_ia_size(ispin), para_env_RPA)
918 310 : 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 310 : mepos_in_RPA_group = MOD(color_sub, integ_group_size)
923 1240 : ALLOCATE (group_grid_2_mepos(0:integ_group_size - 1, 0:para_env_sub%num_pe - 1))
924 1106 : group_grid_2_mepos = 0
925 310 : group_grid_2_mepos(mepos_in_RPA_group, para_env_sub%mepos) = para_env_RPA%mepos
926 310 : 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 310 : integ_group_size, color_rpa_group)
932 :
933 310 : DEALLOCATE (group_grid_2_mepos)
934 310 : CALL cp_fm_struct_release(fm_struct)
935 :
936 : ! deallocate the info array
937 310 : CALL release_group_dist(gd_L)
938 310 : CALL release_group_dist(gd_ia)
939 :
940 : ! sum the local data across processes belonging to different RPA group.
941 588 : IF (para_env_RPA%num_pe /= para_env%num_pe) THEN
942 : BLOCK
943 : TYPE(mp_comm_type) :: comm_exchange
944 144 : comm_exchange = fm_mat_S(ispin)%matrix_struct%context%interconnect(para_env)
945 144 : CALL comm_exchange%sum(fm_mat_S(ispin)%local_data)
946 288 : CALL comm_exchange%free()
947 : END BLOCK
948 : END IF
949 : END DO
950 : END IF
951 :
952 412 : 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 160 : 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 160 : nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.TRUE.)
958 348 : DO ispin = 1, SIZE(fm_mat_Q_gemm)
959 348 : CALL cp_fm_create(fm_mat_Q_gemm(ispin), fm_struct, name="fm_mat_Q_gemm")
960 : END DO
961 160 : CALL cp_fm_struct_release(fm_struct)
962 : END IF
963 :
964 412 : IF (PRESENT(fm_mat_Q)) THEN
965 352 : NULLIFY (blacs_env_Q)
966 352 : IF (my_blacs_ext) THEN
967 58 : blacs_env_Q => blacs_env_ext
968 294 : ELSE IF (para_env_RPA%num_pe == para_env%num_pe .AND. PRESENT(qs_env)) THEN
969 202 : CALL get_qs_env(qs_env, blacs_env=blacs_env_Q)
970 : ELSE
971 92 : CALL cp_blacs_env_create(blacs_env=blacs_env_Q, para_env=para_env_RPA)
972 : END IF
973 352 : NULLIFY (fm_struct)
974 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_Q, nrow_global=dimen_RI, &
975 352 : ncol_global=dimen_RI, para_env=para_env_RPA)
976 760 : DO ispin = 1, SIZE(fm_mat_Q)
977 760 : CALL cp_fm_create(fm_mat_Q(ispin), fm_struct, name="fm_mat_Q")
978 : END DO
979 :
980 352 : CALL cp_fm_struct_release(fm_struct)
981 :
982 352 : IF (.NOT. (my_blacs_ext .OR. (para_env_RPA%num_pe == para_env%num_pe .AND. PRESENT(qs_env)))) &
983 92 : CALL cp_blacs_env_release(blacs_env_Q)
984 : END IF
985 :
986 : ! release blacs_env
987 412 : IF (.NOT. my_blacs_S_ext) THEN
988 294 : CALL cp_blacs_env_release(blacs_env)
989 : ELSE
990 118 : NULLIFY (blacs_env)
991 : END IF
992 :
993 412 : CALL timestop(handle)
994 :
995 412 : 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 294 : SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_sub, unit_nr, &
1055 294 : homo, virtual, dimen_RI, dimen_RI_red, dimen_ia, dimen_nm_gw, &
1056 : Eigenval, num_integ_points, num_integ_group, color_rpa_group, &
1057 588 : 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 294 : my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
1060 : bse_lev_virt, &
1061 294 : 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 294 : 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 294 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c
1129 294 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c
1130 588 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, prim_blk_sizes, &
1131 294 : 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 294 : 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 588 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, e_fermi, tau_tj, tau_wj, tj, &
1144 294 : trace_Qomega, vec_omega_fit_gw, wj, &
1145 294 : wkp_W
1146 294 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_W_gw, weights_cos_tf_t_to_w, &
1147 294 : weights_cos_tf_w_to_t, &
1148 294 : weights_sin_tf_t_to_w
1149 294 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_last, Eigenval_scf, &
1150 294 : 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 294 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_S_gw_work, fm_mat_W, &
1156 294 : fm_mo_coeff_occ, fm_mo_coeff_virt
1157 294 : 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 294 : DIMENSION(:, :, :) :: mat_P_omega
1162 294 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_im_mo_mo, &
1163 294 : matrix_berry_re_mo_mo
1164 294 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_kp
1165 : TYPE(dbcsr_type), POINTER :: mat_W, mat_work
1166 2058 : TYPE(dbt_type) :: t_3c_overl_int_ao_mo
1167 294 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_3c_overl_int_gw_AO, &
1168 294 : t_3c_overl_int_gw_RI, &
1169 294 : t_3c_overl_nnP_ic, &
1170 294 : t_3c_overl_nnP_ic_reflected
1171 : TYPE(dgemm_counter_type) :: dgemm_counter
1172 : TYPE(hfx_compression_type), ALLOCATABLE, &
1173 294 : DIMENSION(:) :: t_3c_O_mo_compressed
1174 21756 : TYPE(im_time_force_type) :: force_data
1175 294 : TYPE(rpa_exchange_work_type) :: exchange_work
1176 1470 : TYPE(rpa_grad_type) :: rpa_grad
1177 294 : TYPE(rpa_sigma_type) :: rpa_sigma
1178 294 : TYPE(two_dim_int_array), ALLOCATABLE, DIMENSION(:) :: t_3c_O_mo_ind
1179 :
1180 294 : CALL timeset(routineN, handle)
1181 :
1182 294 : nspins = SIZE(homo)
1183 294 : nmo = homo(1) + virtual(1)
1184 :
1185 294 : my_open_shell = (nspins == 2)
1186 :
1187 294 : do_gw_im_time = my_do_gw .AND. do_im_time
1188 294 : do_ri_Sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x
1189 294 : do_ic_model = mp2_env%ri_g0w0%do_ic_model
1190 294 : print_ic_values = mp2_env%ri_g0w0%print_ic_values
1191 294 : do_periodic = mp2_env%ri_g0w0%do_periodic
1192 294 : do_kpoints_cubic_RPA = mp2_env%ri_rpa_im_time%do_im_time_kpoints
1193 :
1194 : ! For SOS-MP2 only gemm is implemented
1195 294 : mm_style = wfc_mm_style_gemm
1196 294 : IF (.NOT. do_ri_sos_laplace_mp2) mm_style = mp2_env%ri_rpa%mm_style
1197 :
1198 294 : IF (my_do_gw) THEN
1199 104 : ext_scaling = 0.2_dp
1200 104 : omega_max_fit = mp2_env%ri_g0w0%omega_max_fit
1201 104 : fermi_level_offset_input = mp2_env%ri_g0w0%fermi_level_offset
1202 104 : iter_evGW = mp2_env%ri_g0w0%iter_evGW
1203 104 : iter_sc_GW0 = mp2_env%ri_g0w0%iter_sc_GW0
1204 104 : IF ((.NOT. do_im_time)) THEN
1205 58 : 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 58 : IF (iter_sc_GW0 .NE. 1) iter_evGW = iter_sc_GW0
1208 : END IF
1209 : ELSE
1210 190 : ext_scaling = 0.0_dp
1211 190 : iter_evGW = 1
1212 190 : iter_sc_GW0 = 1
1213 : END IF
1214 :
1215 294 : IF (do_kpoints_cubic_RPA .AND. do_ri_sos_laplace_mp2) THEN
1216 0 : CPABORT("RI-SOS-Laplace-MP2 with k-point-sampling is not implemented.")
1217 : END IF
1218 :
1219 294 : do_apply_ic_corr_to_gw = .FALSE.
1220 294 : IF (mp2_env%ri_g0w0%ic_corr_list(1)%array(1) > 0.0_dp) do_apply_ic_corr_to_gw = .TRUE.
1221 :
1222 294 : IF (do_im_time) THEN
1223 134 : CPASSERT(do_minimax_quad .OR. do_ri_sos_laplace_mp2)
1224 :
1225 134 : group_size_P = mp2_env%ri_rpa_im_time%group_size_P
1226 134 : cut_memory = mp2_env%ri_rpa_im_time%cut_memory
1227 134 : eps_filter = mp2_env%ri_rpa_im_time%eps_filter
1228 : eps_filter_im_time = mp2_env%ri_rpa_im_time%eps_filter* &
1229 134 : mp2_env%ri_rpa_im_time%eps_filter_factor
1230 :
1231 134 : 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 134 : fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, homo, nmo)
1250 :
1251 134 : IF (calc_forces) CALL init_im_time_forces(force_data, fm_matrix_PQ, t_3c_M, unit_nr, mp2_env, qs_env)
1252 :
1253 134 : IF (my_do_gw) THEN
1254 :
1255 : CALL dbcsr_get_info(mat_P_global%matrix, &
1256 46 : row_blk_size=RI_blk_sizes)
1257 :
1258 : CALL dbcsr_get_info(matrix_s(1)%matrix, &
1259 46 : row_blk_size=prim_blk_sizes)
1260 :
1261 46 : gw_corr_lev_tot = gw_corr_lev_occ(1) + gw_corr_lev_virt(1)
1262 :
1263 46 : 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 46 : qs_env)
1276 :
1277 : END IF
1278 : END IF
1279 :
1280 : END IF
1281 294 : IF (do_ic_model) THEN
1282 : ! image charge model only implemented for cubic scaling GW
1283 2 : CPASSERT(do_gw_im_time)
1284 2 : CPASSERT(.NOT. do_periodic)
1285 2 : IF (cut_memory .NE. 1) CPABORT("For IC, use MEMORY_CUT 1 in the LOW_SCALING section.")
1286 : END IF
1287 :
1288 882 : ALLOCATE (e_fermi(nspins))
1289 294 : IF (do_minimax_quad .OR. do_ri_sos_laplace_mp2) THEN
1290 200 : 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 200 : 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 200 : 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 1464 : wj(:) = wj(:)*4.0_dp
1301 200 : IF (ierr == 0) THEN
1302 74 : IF (unit_nr > 0) THEN
1303 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
1304 37 : "GREENX MINIMAX_INFO| Number of integration points:", num_integ_points
1305 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") &
1306 37 : "GREENX MINIMAX_INFO| Gap (Emin):", Emin
1307 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") &
1308 37 : "GREENX MINIMAX_INFO| Maximum eigenvalue difference (Emax):", Emax
1309 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.4)") &
1310 37 : "GREENX MINIMAX_INFO| Energy range (Emax/Emin):", Emax/Emin
1311 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
1312 37 : "GREENX MINIMAX_INFO| Frequency grid (scaled):", "Weights", "Abscissas"
1313 425 : DO gi = 1, num_integ_points
1314 425 : 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 37 : "GREENX MINIMAX_INFO| Time grid (scaled):", "Weights", "Abscissas"
1318 425 : DO gi = 1, num_integ_points
1319 425 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(gi), tau_tj(gi)
1320 : END DO
1321 37 : CALL m_flush(unit_nr)
1322 : END IF
1323 : ELSE
1324 126 : IF (unit_nr > 0) THEN
1325 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75)") &
1326 63 : "GREENX MINIMAX_INFO| Grid not available, use internal CP2K grids"
1327 63 : CALL m_flush(unit_nr)
1328 : END IF
1329 126 : IF (ALLOCATED(tau_tj)) THEN
1330 126 : DEALLOCATE (tau_tj)
1331 : END IF
1332 126 : IF (ALLOCATED(tau_wj)) THEN
1333 126 : DEALLOCATE (tau_wj)
1334 : END IF
1335 126 : IF (ALLOCATED(tj)) THEN
1336 126 : DEALLOCATE (tj)
1337 : END IF
1338 126 : IF (ALLOCATED(wj)) THEN
1339 126 : DEALLOCATE (wj)
1340 : END IF
1341 126 : IF (ALLOCATED(weights_cos_tf_t_to_w)) THEN
1342 0 : DEALLOCATE (weights_cos_tf_t_to_w)
1343 : END IF
1344 126 : IF (ALLOCATED(weights_cos_tf_w_to_t)) THEN
1345 0 : DEALLOCATE (weights_cos_tf_w_to_t)
1346 : END IF
1347 126 : IF (ALLOCATED(weights_sin_tf_t_to_w)) THEN
1348 0 : 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 126 : 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 200 : 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 200 : num_integ_points, unit_nr, qs_env)
1371 : END IF
1372 : ELSE
1373 94 : 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 94 : ext_scaling, a_scaling, tj, wj)
1377 : END IF
1378 :
1379 : ! This array is needed for RPA
1380 294 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1381 708 : ALLOCATE (trace_Qomega(dimen_RI_red))
1382 : END IF
1383 :
1384 294 : IF (do_ri_sos_laplace_mp2 .AND. .NOT. do_im_time) THEN
1385 28 : alpha = 1.0_dp
1386 266 : ELSE IF (my_open_shell .OR. do_ri_sos_laplace_mp2) THEN
1387 70 : alpha = 2.0_dp
1388 : ELSE
1389 196 : alpha = 4.0_dp
1390 : END IF
1391 294 : 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 104 : do_kpoints_cubic_RPA, do_kpoints_from_Gamma)
1403 :
1404 104 : IF (do_bse) THEN
1405 :
1406 30 : CALL cp_fm_create(fm_mat_Q_static_bse_gemm, fm_mat_Q_gemm(1)%matrix_struct)
1407 30 : CALL cp_fm_to_fm(fm_mat_Q_gemm(1), fm_mat_Q_static_bse_gemm)
1408 30 : CALL cp_fm_set_all(fm_mat_Q_static_bse_gemm, 0.0_dp)
1409 :
1410 : END IF
1411 :
1412 : END IF
1413 :
1414 294 : 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 40 : unit_nr, do_ri_sos_laplace_mp2)
1417 : !doubtful
1418 294 : 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 132 : fm_mat_S, fm_mat_Q(1), fm_mat_Q_gemm(1), homo, virtual)
1421 294 : Erpa = 0.0_dp
1422 294 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) e_exchange = 0.0_dp
1423 294 : first_cycle = .TRUE.
1424 294 : omega_old = 0.0_dp
1425 294 : CALL dgemm_counter_init(dgemm_counter, unit_nr, mp2_env%ri_rpa%print_dgemm_info)
1426 :
1427 698 : DO count_ev_sc_GW = 1, iter_evGW
1428 424 : dbcsr_time = 0.0_dp
1429 424 : dbcsr_nflop = 0
1430 :
1431 424 : IF (do_ic_model) CYCLE
1432 :
1433 : ! reset some values, important when doing eigenvalue self-consistent GW
1434 422 : IF (my_do_gw) THEN
1435 232 : Erpa = 0.0_dp
1436 156488 : vec_Sigma_c_gw = z_zero
1437 232 : first_cycle = .TRUE.
1438 : END IF
1439 :
1440 : ! calculate Q_PQ(it)
1441 422 : IF (do_im_time) THEN ! not using Imaginary time
1442 :
1443 146 : IF (.NOT. do_kpoints_cubic_RPA) THEN
1444 312 : DO ispin = 1, nspins
1445 312 : 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 146 : tau = 0.0_dp
1450 146 : tau_old = 0.0_dp
1451 :
1452 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(/T3,A,T66,i15)") &
1453 73 : "MEMORY_INFO| Memory cut:", cut_memory
1454 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
1455 73 : "SPARSITY_INFO| Eps filter for M virt/occ tensors:", eps_filter
1456 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
1457 73 : "SPARSITY_INFO| Eps filter for P matrix:", eps_filter_im_time
1458 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,i15)") &
1459 73 : "SPARSITY_INFO| Minimum tensor block size:", min_bsize
1460 :
1461 : ! for evGW, we have to ensure that mat_P_omega is zero
1462 146 : 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 146 : dbcsr_time, dbcsr_nflop)
1480 :
1481 : ! the same for open shell, use fm_mo_coeff_occ_beta and fm_mo_coeff_virt_beta
1482 146 : IF (my_open_shell) THEN
1483 28 : 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 28 : 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 28 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1503 90 : DO j = 1, SIZE(mat_P_omega, 2)
1504 598 : DO i = 1, SIZE(mat_P_omega, 1)
1505 508 : CALL dbcsr_add(mat_P_omega(i, j, 1)%matrix, mat_P_omega(i, j, 2)%matrix, 1.0_dp, 1.0_dp)
1506 578 : 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 422 : IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1515 10 : 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 12544 : DO jquad = 1, num_integ_points
1519 12122 : IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE
1520 :
1521 11482 : CALL timeset(routineN//"_RPA_matrix_operations", handle3)
1522 :
1523 11482 : IF (do_ri_sos_laplace_mp2) THEN
1524 206 : omega = tau_tj(jquad)
1525 : ELSE
1526 11276 : IF (do_minimax_quad) THEN
1527 1186 : omega = tj(jquad)
1528 : ELSE
1529 10090 : omega = a_scaling/TAN(tj(jquad))
1530 : END IF
1531 : END IF ! do_ri_sos_laplace_mp2
1532 :
1533 11482 : 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 1168 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
1537 :
1538 2324 : 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 2324 : fm_mat_Minv_L_kpoints(1, 1), fm_mat_Q(ispin))
1542 : END DO
1543 : END IF
1544 :
1545 : ELSE
1546 10314 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A, 1X, I3, 1X, A, 1X, I3)") &
1547 5157 : "INTEG_INFO| Started with Integration point", jquad, "of", num_integ_points
1548 :
1549 10314 : IF (first_cycle .AND. count_ev_sc_gw > 1) THEN
1550 116 : IF (iter_sc_gw0 == 1) THEN
1551 124 : DO ispin = 1, nspins
1552 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virtual(ispin), &
1553 124 : Eigenval_last(:, 1, ispin), homo(ispin), omega_old)
1554 : END DO
1555 : ELSE
1556 116 : DO ispin = 1, nspins
1557 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virtual(ispin), &
1558 116 : Eigenval_scf(:, 1, ispin), homo(ispin), omega_old)
1559 : END DO
1560 : END IF
1561 : END IF
1562 :
1563 10314 : IF (iter_sc_GW0 > 1) THEN
1564 12140 : 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 12140 : num_integ_points, count_ev_sc_GW)
1570 : END DO
1571 :
1572 : ! For SOS-MP2 we need both matrices separately
1573 6070 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1574 6070 : DO ispin = 2, nspins
1575 6070 : 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 8652 : 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 8652 : num_integ_points, count_ev_sc_GW)
1585 : END DO
1586 :
1587 : ! For SOS-MP2 we need both matrices separately
1588 4244 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1589 4280 : DO ispin = 2, nspins
1590 4280 : 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 11482 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
1600 12 : e_exchange_corr = 0.0_dp
1601 12 : 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 12 : 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 11482 : IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1609 30 : CALL rpa_sigma_matrix_spectral(rpa_sigma, fm_mat_Q(1), wj(jquad), para_env_RPA)
1610 : END IF
1611 :
1612 11482 : IF (do_ri_sos_laplace_mp2) THEN
1613 :
1614 206 : CALL SOS_MP2_postprocessing(fm_mat_Q, Erpa, tau_wj(jquad))
1615 :
1616 206 : 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 60 : Eigenval(:, 1, :), tau_wj(jquad), unit_nr)
1619 : ELSE
1620 11276 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_copy_Q(fm_mat_Q(1), rpa_grad)
1621 :
1622 11276 : CALL Q_trace_and_add_unit_matrix(dimen_RI_red, trace_Qomega, fm_mat_Q(1))
1623 :
1624 11276 : 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 120 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv)
1633 : ELSE
1634 11156 : 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 11276 : 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 48 : 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 11482 : first_cycle = .FALSE.
1644 11482 : omega_old = omega
1645 :
1646 11482 : CALL timestop(handle3)
1647 :
1648 11482 : IF (my_do_gw) THEN
1649 :
1650 10708 : 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 10708 : IF (do_im_time) THEN
1654 : ! only for molecules
1655 818 : 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 722 : 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 9890 : kpoints, qs_env, mp2_env)
1669 : END IF
1670 : END IF
1671 :
1672 11482 : IF (unit_nr > 0) CALL m_flush(unit_nr)
1673 24026 : CALL para_env%sync() ! sync to see output
1674 :
1675 : END DO ! jquad
1676 :
1677 422 : IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1678 10 : CALL finalize_rpa_sigma(rpa_sigma, unit_nr, mp2_env%ri_rpa%e_sigma_corr, para_env, do_minimax_quad)
1679 10 : IF (do_minimax_quad) mp2_env%ri_rpa%e_sigma_corr = mp2_env%ri_rpa%e_sigma_corr/2.0_dp
1680 10 : CALL para_env%sum(mp2_env%ri_rpa%e_sigma_corr)
1681 : END IF
1682 :
1683 422 : CALL para_env%sum(Erpa)
1684 :
1685 422 : IF (.NOT. (do_ri_sos_laplace_mp2)) THEN
1686 364 : Erpa = Erpa/(pi*2.0_dp)
1687 364 : IF (do_minimax_quad) Erpa = Erpa/2.0_dp
1688 : END IF
1689 :
1690 422 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
1691 12 : CALL para_env%sum(E_exchange)
1692 12 : E_exchange = E_exchange/(pi*2.0_dp)
1693 12 : IF (do_minimax_quad) E_exchange = E_exchange/2.0_dp
1694 12 : mp2_env%ri_rpa%ener_exchange = E_exchange
1695 : END IF
1696 :
1697 422 : IF (calc_forces .AND. do_ri_sos_laplace_mp2 .AND. do_im_time) THEN
1698 22 : IF (my_open_shell) THEN
1699 4 : Pspin = 1
1700 4 : 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 4 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1709 4 : Pspin = 2
1710 4 : 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 4 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1719 :
1720 : ELSE
1721 18 : Pspin = 1
1722 18 : 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 18 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1731 : END IF
1732 22 : CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1733 : END IF !laplace SOS-MP2
1734 :
1735 422 : IF (calc_forces .AND. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) THEN
1736 64 : 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 64 : dbcsr_nflop, mp2_env, qs_env)
1746 : END DO
1747 28 : CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1748 : END IF
1749 :
1750 422 : IF (do_im_time) THEN
1751 :
1752 146 : my_flop_rate = REAL(dbcsr_nflop, dp)/(1.0E09_dp*dbcsr_time)
1753 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(/T3,A,T73,ES8.2)") &
1754 73 : "PERFORMANCE| DBCSR total number of flops:", REAL(dbcsr_nflop*para_env%num_pe, dp)
1755 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.2)") &
1756 73 : "PERFORMANCE| DBCSR total execution time:", dbcsr_time
1757 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.2)") &
1758 73 : "PERFORMANCE| DBCSR flop rate (Gflops / MPI rank):", my_flop_rate
1759 :
1760 : ELSE
1761 :
1762 276 : 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 716 : 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 232 : 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 232 : IF (exit_ev_gw) EXIT
1793 :
1794 : END IF ! my_do_gw if
1795 :
1796 : END DO ! evGW loop
1797 :
1798 294 : IF (do_ic_model) THEN
1799 :
1800 2 : 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 0 : 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 0 : 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 2 : print_ic_values, para_env)
1821 :
1822 : END IF
1823 :
1824 : END IF
1825 :
1826 : ! postprocessing after GW for Bethe-Salpeter
1827 294 : 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 30 : IF (mp2_env%ri_g0w0%iter_evGW > 1) THEN
1831 4 : IF (unit_nr > 0) THEN
1832 : CALL cp_warn(__LOCATION__, &
1833 2 : "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 30 : CALL cp_fm_create(fm_mat_S_ia_bse, fm_mat_S(1)%matrix_struct)
1838 30 : 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 30 : IF (iter_sc_gw0 == 1) THEN
1841 : CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse, virtual(1), &
1842 18 : Eigenval_last(:, 1, 1), homo(1), omega)
1843 : ELSE
1844 : CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse, virtual(1), &
1845 12 : 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 30 : gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
1853 : ! Release BSE-copy of fm_mat_S
1854 30 : CALL cp_fm_release(fm_mat_S_ia_bse)
1855 : END IF
1856 :
1857 294 : 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 104 : kpoints, vec_Sigma_x_gw,.NOT. do_im_time)
1862 : END IF
1863 :
1864 294 : 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 134 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, mat_work, qs_env)
1876 :
1877 134 : 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 46 : mat_W, qs_env)
1884 : END IF
1885 :
1886 : END IF
1887 :
1888 294 : IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) CALL exchange_work%release()
1889 :
1890 294 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1891 236 : DEALLOCATE (tj)
1892 236 : DEALLOCATE (wj)
1893 236 : DEALLOCATE (trace_Qomega)
1894 : END IF
1895 :
1896 294 : IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
1897 162 : DEALLOCATE (tau_tj)
1898 162 : DEALLOCATE (tau_wj)
1899 : END IF
1900 :
1901 294 : IF (do_im_time .AND. calc_forces) THEN
1902 50 : CALL im_time_force_release(force_data)
1903 : END IF
1904 :
1905 294 : 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 40 : homo, virtual)
1908 :
1909 294 : CALL timestop(handle)
1910 :
1911 6002 : END SUBROUTINE rpa_num_int
1912 :
1913 : END MODULE rpa_main
|
