Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to calculate 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 : wfc_mm_style_gemm
55 : USE kinds, ONLY: dp,&
56 : int_8
57 : USE kpoint_types, ONLY: kpoint_type
58 : USE machine, ONLY: m_flush,&
59 : m_memory
60 : USE mathconstants, ONLY: pi,&
61 : z_zero
62 : USE message_passing, ONLY: mp_comm_type,&
63 : mp_para_env_release,&
64 : mp_para_env_type
65 : USE minimax_exp, ONLY: check_exp_minimax_range
66 : USE mp2_grids, ONLY: get_clenshaw_grid,&
67 : get_minimax_grid
68 : USE mp2_laplace, ONLY: SOS_MP2_postprocessing
69 : USE mp2_ri_grad_util, ONLY: array2fm
70 : USE mp2_types, ONLY: mp2_type,&
71 : three_dim_real_array,&
72 : two_dim_int_array,&
73 : two_dim_real_array
74 : USE qs_environment_types, ONLY: get_qs_env,&
75 : qs_environment_type
76 : USE rpa_exchange, ONLY: rpa_exchange_needed_mem,&
77 : rpa_exchange_work_type
78 : USE rpa_grad, ONLY: rpa_grad_copy_Q,&
79 : rpa_grad_create,&
80 : rpa_grad_finalize,&
81 : rpa_grad_matrix_operations,&
82 : rpa_grad_needed_mem,&
83 : rpa_grad_type
84 : USE rpa_gw, ONLY: allocate_matrices_gw,&
85 : allocate_matrices_gw_im_time,&
86 : compute_GW_self_energy,&
87 : compute_QP_energies,&
88 : compute_W_cubic_GW,&
89 : deallocate_matrices_gw,&
90 : deallocate_matrices_gw_im_time,&
91 : get_fermi_level_offset
92 : USE rpa_gw_ic, ONLY: calculate_ic_correction
93 : USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_MOs,&
94 : invert_eps_compute_W_and_Erpa_kp
95 : USE rpa_im_time, ONLY: compute_mat_P_omega,&
96 : zero_mat_P_omega
97 : USE rpa_im_time_force_methods, ONLY: calc_laplace_loop_forces,&
98 : calc_post_loop_forces,&
99 : calc_rpa_loop_forces,&
100 : init_im_time_forces,&
101 : keep_initial_quad
102 : USE rpa_im_time_force_types, ONLY: im_time_force_release,&
103 : im_time_force_type
104 : USE rpa_util, ONLY: Q_trace_and_add_unit_matrix,&
105 : alloc_im_time,&
106 : calc_mat_Q,&
107 : compute_Erpa_by_freq_int,&
108 : contract_P_omega_with_mat_L,&
109 : dealloc_im_time,&
110 : remove_scaling_factor_rpa
111 : USE util, ONLY: get_limit
112 : #include "./base/base_uses.f90"
113 :
114 : IMPLICIT NONE
115 :
116 : PRIVATE
117 :
118 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_main'
119 :
120 : PUBLIC :: rpa_ri_compute_en
121 :
122 : CONTAINS
123 :
124 : ! **************************************************************************************************
125 : !> \brief ...
126 : !> \param qs_env ...
127 : !> \param Erpa ...
128 : !> \param mp2_env ...
129 : !> \param BIb_C ...
130 : !> \param BIb_C_gw ...
131 : !> \param BIb_C_bse_ij ...
132 : !> \param BIb_C_bse_ab ...
133 : !> \param para_env ...
134 : !> \param para_env_sub ...
135 : !> \param color_sub ...
136 : !> \param gd_array ...
137 : !> \param gd_B_virtual ...
138 : !> \param gd_B_all ...
139 : !> \param gd_B_occ_bse ...
140 : !> \param gd_B_virt_bse ...
141 : !> \param mo_coeff ...
142 : !> \param fm_matrix_PQ ...
143 : !> \param fm_matrix_L_kpoints ...
144 : !> \param fm_matrix_Minv_L_kpoints ...
145 : !> \param fm_matrix_Minv ...
146 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
147 : !> \param kpoints ...
148 : !> \param Eigenval ...
149 : !> \param nmo ...
150 : !> \param homo ...
151 : !> \param dimen_RI ...
152 : !> \param dimen_RI_red ...
153 : !> \param gw_corr_lev_occ ...
154 : !> \param gw_corr_lev_virt ...
155 : !> \param bse_lev_virt ...
156 : !> \param unit_nr ...
157 : !> \param do_ri_sos_laplace_mp2 ...
158 : !> \param my_do_gw ...
159 : !> \param do_im_time ...
160 : !> \param do_bse ...
161 : !> \param matrix_s ...
162 : !> \param mat_munu ...
163 : !> \param mat_P_global ...
164 : !> \param t_3c_M ...
165 : !> \param t_3c_O ...
166 : !> \param t_3c_O_compressed ...
167 : !> \param t_3c_O_ind ...
168 : !> \param starts_array_mc ...
169 : !> \param ends_array_mc ...
170 : !> \param starts_array_mc_block ...
171 : !> \param ends_array_mc_block ...
172 : !> \param calc_forces ...
173 : ! **************************************************************************************************
174 1360 : SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, &
175 : para_env, para_env_sub, color_sub, &
176 272 : gd_array, gd_B_virtual, gd_B_all, gd_B_occ_bse, gd_B_virt_bse, &
177 272 : mo_coeff, fm_matrix_PQ, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
178 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, kpoints, &
179 544 : Eigenval, nmo, homo, dimen_RI, dimen_RI_red, gw_corr_lev_occ, gw_corr_lev_virt, &
180 : bse_lev_virt, &
181 : unit_nr, do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_bse, matrix_s, &
182 : mat_munu, mat_P_global, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
183 : starts_array_mc, ends_array_mc, &
184 : starts_array_mc_block, ends_array_mc_block, calc_forces)
185 :
186 : TYPE(qs_environment_type), POINTER :: qs_env
187 : REAL(KIND=dp), INTENT(OUT) :: Erpa
188 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
189 : TYPE(three_dim_real_array), DIMENSION(:), &
190 : INTENT(INOUT) :: BIb_C, BIb_C_gw
191 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
192 : INTENT(INOUT) :: BIb_C_bse_ij, BIb_C_bse_ab
193 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
194 : INTEGER, INTENT(INOUT) :: color_sub
195 : TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_array
196 : TYPE(group_dist_d1_type), DIMENSION(:), &
197 : INTENT(INOUT) :: gd_B_virtual
198 : TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_B_all, gd_B_occ_bse, gd_B_virt_bse
199 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
200 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_PQ
201 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_L_kpoints, &
202 : fm_matrix_Minv_L_kpoints, &
203 : fm_matrix_Minv, &
204 : fm_matrix_Minv_Vtrunc_Minv
205 : TYPE(kpoint_type), POINTER :: kpoints
206 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
207 : INTENT(INOUT) :: Eigenval
208 : INTEGER, INTENT(IN) :: nmo
209 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
210 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
211 : INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt
212 : INTEGER, INTENT(IN) :: bse_lev_virt, unit_nr
213 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, my_do_gw, &
214 : do_im_time, do_bse
215 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
216 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
217 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_P_global
218 : TYPE(dbt_type) :: t_3c_M
219 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O
220 : TYPE(hfx_compression_type), ALLOCATABLE, &
221 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_compressed
222 : TYPE(block_ind_type), ALLOCATABLE, &
223 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_ind
224 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
225 : starts_array_mc_block, &
226 : ends_array_mc_block
227 : LOGICAL, INTENT(IN) :: calc_forces
228 :
229 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_ri_compute_en'
230 :
231 : INTEGER :: best_integ_group_size, best_num_integ_point, color_rpa_group, dimen_homo_square, &
232 : dimen_nm_gw, dimen_virt_square, handle, handle2, handle3, ierr, iiB, &
233 : input_num_integ_groups, integ_group_size, ispin, jjB, min_integ_group_size, &
234 : my_ab_comb_bse_end, my_ab_comb_bse_size, my_ab_comb_bse_start, my_group_L_end, &
235 : my_group_L_size, my_group_L_start, my_ij_comb_bse_end, my_ij_comb_bse_size, &
236 : my_ij_comb_bse_start, my_nm_gw_end, my_nm_gw_size, my_nm_gw_start, ncol_block_mat, &
237 : ngroup, nrow_block_mat, nspins, num_integ_group, num_integ_points, pos_integ_group
238 : INTEGER(KIND=int_8) :: mem
239 544 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dimen_ia, my_ia_end, my_ia_size, &
240 272 : my_ia_start, virtual
241 : LOGICAL :: do_kpoints_from_Gamma, do_minimax_quad, &
242 : my_open_shell, skip_integ_group_opt
243 : REAL(KIND=dp) :: allowed_memory, avail_mem, E_Range, Emax, Emin, mem_for_iaK, mem_for_QK, &
244 : mem_min, mem_per_group, mem_per_rank, mem_per_repl, mem_real
245 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_kp
246 544 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_Q, fm_mat_Q_gemm, fm_mat_S, &
247 272 : fm_mat_S_gw
248 816 : TYPE(cp_fm_type), DIMENSION(1) :: fm_mat_R_gw, fm_mat_S_ab_bse, &
249 544 : fm_mat_S_ij_bse
250 : TYPE(mp_para_env_type), POINTER :: para_env_RPA
251 : TYPE(two_dim_real_array), ALLOCATABLE, &
252 544 : DIMENSION(:) :: BIb_C_2D, BIb_C_2D_gw
253 2720 : TYPE(two_dim_real_array), DIMENSION(1) :: BIb_C_2D_bse_ab, BIb_C_2D_bse_ij
254 :
255 272 : CALL timeset(routineN, handle)
256 :
257 272 : CALL cite_reference(DelBen2013)
258 272 : CALL cite_reference(DelBen2015)
259 :
260 272 : IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_axk) THEN
261 10 : CALL cite_reference(Bates2013)
262 262 : ELSE IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_sosex) THEN
263 2 : CALL cite_reference(Freeman1977)
264 2 : CALL cite_reference(Gruneis2009)
265 : END IF
266 272 : IF (mp2_env%ri_rpa%do_rse) THEN
267 4 : CALL cite_reference(Ren2011)
268 4 : CALL cite_reference(Ren2013)
269 : END IF
270 :
271 272 : IF (my_do_gw) THEN
272 92 : CALL cite_reference(Wilhelm2016a)
273 92 : CALL cite_reference(Wilhelm2017)
274 92 : CALL cite_reference(Wilhelm2018)
275 : END IF
276 :
277 272 : IF (do_im_time) THEN
278 134 : CALL cite_reference(Wilhelm2016b)
279 : END IF
280 :
281 272 : nspins = SIZE(homo)
282 272 : my_open_shell = (nspins == 2)
283 1904 : ALLOCATE (virtual(nspins), dimen_ia(nspins), my_ia_end(nspins), my_ia_start(nspins), my_ia_size(nspins))
284 598 : virtual(:) = nmo - homo(:)
285 598 : dimen_ia(:) = virtual(:)*homo(:)
286 :
287 1088 : ALLOCATE (Eigenval_kp(nmo, 1, nspins))
288 8378 : Eigenval_kp(:, 1, :) = Eigenval(:, :)
289 :
290 272 : IF (do_im_time) mp2_env%ri_rpa%minimax_quad = .TRUE.
291 272 : do_minimax_quad = mp2_env%ri_rpa%minimax_quad
292 :
293 272 : IF (do_ri_sos_laplace_mp2) THEN
294 58 : num_integ_points = mp2_env%ri_laplace%n_quadrature
295 58 : input_num_integ_groups = mp2_env%ri_laplace%num_integ_groups
296 :
297 : ! check the range for the minimax approximation
298 58 : E_Range = mp2_env%e_range
299 58 : IF (mp2_env%e_range <= 1.0_dp .OR. mp2_env%e_gap <= 0.0_dp) THEN
300 : Emin = HUGE(dp)
301 : Emax = 0.0_dp
302 88 : DO ispin = 1, nspins
303 88 : IF (homo(ispin) > 0) THEN
304 50 : Emin = MIN(Emin, 2.0_dp*(Eigenval(homo(ispin) + 1, ispin) - Eigenval(homo(ispin), ispin)))
305 2740 : Emax = MAX(Emax, 2.0_dp*(MAXVAL(Eigenval(:, ispin)) - MINVAL(Eigenval(:, ispin))))
306 : END IF
307 : END DO
308 38 : E_Range = Emax/Emin
309 : END IF
310 58 : IF (E_Range < 2.0_dp) E_Range = 2.0_dp
311 : ierr = 0
312 58 : CALL check_exp_minimax_range(num_integ_points, E_Range, ierr)
313 58 : IF (ierr /= 0) THEN
314 : jjB = num_integ_points - 1
315 0 : DO iiB = 1, jjB
316 0 : num_integ_points = num_integ_points - 1
317 : ierr = 0
318 0 : CALL check_exp_minimax_range(num_integ_points, E_Range, ierr)
319 0 : IF (ierr == 0) EXIT
320 : END DO
321 : END IF
322 58 : CPASSERT(num_integ_points >= 1)
323 : ELSE
324 214 : num_integ_points = mp2_env%ri_rpa%rpa_num_quad_points
325 214 : input_num_integ_groups = mp2_env%ri_rpa%rpa_num_integ_groups
326 :
327 214 : IF (my_do_gw .AND. do_minimax_quad) THEN
328 46 : IF (num_integ_points > 34) THEN
329 0 : IF (unit_nr > 0) &
330 : CALL cp_warn(__LOCATION__, &
331 : "The required number of quadrature point exceeds the maximum possible in the "// &
332 0 : "Minimax quadrature scheme. The number of quadrature point has been reset to 30.")
333 0 : num_integ_points = 30
334 : END IF
335 : ELSE
336 168 : IF (do_minimax_quad .AND. num_integ_points > 20) THEN
337 0 : IF (unit_nr > 0) &
338 : CALL cp_warn(__LOCATION__, &
339 : "The required number of quadrature point exceeds the maximum possible in the "// &
340 0 : "Minimax quadrature scheme. The number of quadrature point has been reset to 20.")
341 0 : num_integ_points = 20
342 : END IF
343 : END IF
344 : END IF
345 272 : allowed_memory = mp2_env%mp2_memory
346 :
347 272 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
348 :
349 272 : ngroup = para_env%num_pe/para_env_sub%num_pe
350 :
351 : ! for imaginary time or periodic GW or BSE, we use all processors for a single frequency/time point
352 272 : IF (do_im_time .OR. mp2_env%ri_g0w0%do_periodic .OR. do_bse) THEN
353 :
354 154 : integ_group_size = ngroup
355 154 : best_num_integ_point = num_integ_points
356 :
357 : ELSE
358 :
359 : ! Calculate available memory and create integral group according to that
360 : ! mem_for_iaK is the memory needed for storing the 3 centre integrals
361 260 : mem_for_iaK = REAL(SUM(dimen_ia), KIND=dp)*dimen_RI_red*8.0_dp/(1024_dp**2)
362 118 : mem_for_QK = REAL(dimen_RI_red, KIND=dp)*nspins*dimen_RI_red*8.0_dp/(1024_dp**2)
363 :
364 118 : CALL m_memory(mem)
365 118 : mem_real = (mem + 1024*1024 - 1)/(1024*1024)
366 118 : CALL para_env%min(mem_real)
367 :
368 118 : mem_per_rank = 0.0_dp
369 :
370 : ! B_ia_P
371 : mem_per_repl = mem_for_iaK
372 : ! Q (regular and for dgemm)
373 118 : mem_per_repl = mem_per_repl + 2.0_dp*mem_for_QK
374 :
375 118 : IF (calc_forces) CALL rpa_grad_needed_mem(homo, virtual, dimen_RI_red, mem_per_rank, mem_per_repl, do_ri_sos_laplace_mp2)
376 118 : CALL rpa_exchange_needed_mem(mp2_env, homo, virtual, dimen_RI_red, para_env, mem_per_rank, mem_per_repl)
377 :
378 118 : mem_min = mem_per_repl/para_env%num_pe + mem_per_rank
379 :
380 118 : IF (unit_nr > 0) THEN
381 59 : WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum required memory per MPI process:', mem_min, ' MiB'
382 59 : WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Available memory per MPI process:', mem_real, ' MiB'
383 : END IF
384 :
385 : ! Use only the allowed amount of memory
386 118 : mem_real = MIN(mem_real, allowed_memory)
387 : ! For the memory estimate, we require the amount of required memory per replication group and the available memory
388 118 : mem_real = mem_real - mem_per_rank
389 :
390 118 : mem_per_group = mem_real*para_env_sub%num_pe
391 :
392 : ! here we try to find the best rpa/laplace group size
393 118 : skip_integ_group_opt = .FALSE.
394 :
395 : ! Check the input number of integration groups
396 118 : IF (input_num_integ_groups > 0) THEN
397 0 : IF (MOD(num_integ_points, input_num_integ_groups) == 0) THEN
398 0 : IF (MOD(ngroup, input_num_integ_groups) == 0) THEN
399 0 : best_integ_group_size = ngroup/input_num_integ_groups
400 0 : best_num_integ_point = num_integ_points/input_num_integ_groups
401 0 : skip_integ_group_opt = .TRUE.
402 : ELSE
403 0 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Total number of groups not multiple of NUM_INTEG_GROUPS'
404 : END IF
405 : ELSE
406 0 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'NUM_QUAD_POINTS not multiple of NUM_INTEG_GROUPS'
407 : END IF
408 : END IF
409 :
410 118 : IF (.NOT. skip_integ_group_opt) THEN
411 118 : best_integ_group_size = ngroup
412 118 : best_num_integ_point = num_integ_points
413 :
414 118 : min_integ_group_size = MAX(1, ngroup/num_integ_points)
415 :
416 118 : integ_group_size = min_integ_group_size - 1
417 132 : DO iiB = min_integ_group_size + 1, ngroup
418 96 : integ_group_size = integ_group_size + 1
419 :
420 : ! check that the ngroup is a multiple of integ_group_size
421 96 : IF (MOD(ngroup, integ_group_size) /= 0) CYCLE
422 :
423 : ! check for memory
424 96 : avail_mem = integ_group_size*mem_per_group
425 96 : IF (avail_mem < mem_per_repl) CYCLE
426 :
427 : ! check that the integration groups have the same size
428 96 : num_integ_group = ngroup/integ_group_size
429 96 : IF (MOD(num_integ_points, num_integ_group) /= 0) CYCLE
430 :
431 82 : best_num_integ_point = num_integ_points/num_integ_group
432 82 : best_integ_group_size = integ_group_size
433 :
434 132 : EXIT
435 :
436 : END DO
437 : END IF
438 :
439 118 : integ_group_size = best_integ_group_size
440 :
441 : END IF
442 :
443 272 : IF (unit_nr > 0 .AND. .NOT. do_im_time) THEN
444 69 : IF (do_ri_sos_laplace_mp2) THEN
445 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
446 14 : "RI_INFO| Group size for laplace numerical integration:", integ_group_size*para_env_sub%num_pe
447 : WRITE (UNIT=unit_nr, FMT="(T3,A)") &
448 14 : "INTEG_INFO| MINIMAX approximation"
449 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
450 14 : "INTEG_INFO| Number of integration points:", num_integ_points
451 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
452 14 : "INTEG_INFO| Number of integration points per Laplace group:", best_num_integ_point
453 : ELSE
454 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
455 55 : "RI_INFO| Group size for frequency integration:", integ_group_size*para_env_sub%num_pe
456 55 : IF (do_minimax_quad) THEN
457 : WRITE (UNIT=unit_nr, FMT="(T3,A)") &
458 16 : "INTEG_INFO| MINIMAX quadrature"
459 : ELSE
460 : WRITE (UNIT=unit_nr, FMT="(T3,A)") &
461 39 : "INTEG_INFO| Clenshaw-Curtius quadrature"
462 : END IF
463 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
464 55 : "INTEG_INFO| Number of integration points:", num_integ_points
465 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
466 55 : "INTEG_INFO| Number of integration points per RPA group:", best_num_integ_point
467 : END IF
468 69 : CALL m_flush(unit_nr)
469 : END IF
470 :
471 272 : num_integ_group = ngroup/integ_group_size
472 :
473 272 : pos_integ_group = MOD(color_sub, integ_group_size)
474 272 : color_rpa_group = color_sub/integ_group_size
475 :
476 272 : CALL timeset(routineN//"_reorder", handle2)
477 :
478 : ! not necessary for imaginary time
479 :
480 1142 : ALLOCATE (BIb_C_2D(nspins))
481 272 : IF (.NOT. do_im_time) THEN
482 :
483 : ! reorder the local data in such a way to help the next stage of matrix creation
484 : ! now the data inside the group are divided into a ia x K matrix
485 300 : DO ispin = 1, nspins
486 : CALL calculate_BIb_C_2D(BIb_C_2D(ispin)%array, BIb_C(ispin)%array, para_env_sub, dimen_ia(ispin), &
487 : homo(ispin), virtual(ispin), gd_B_virtual(ispin), &
488 162 : my_ia_size(ispin), my_ia_start(ispin), my_ia_end(ispin), my_group_L_size)
489 :
490 162 : DEALLOCATE (BIb_C(ispin)%array)
491 462 : CALL release_group_dist(gd_B_virtual(ispin))
492 :
493 : END DO
494 :
495 : ! in the GW case, BIb_C_2D_gw is an nm x K matrix, with n: number of corr GW levels, m=nmo
496 138 : IF (my_do_gw) THEN
497 142 : ALLOCATE (BIb_C_2D_gw(nspins))
498 :
499 46 : CALL timeset(routineN//"_reorder_gw", handle3)
500 :
501 46 : dimen_nm_gw = nmo*(gw_corr_lev_occ(1) + gw_corr_lev_virt(1))
502 :
503 : ! The same for open shell
504 96 : DO ispin = 1, nspins
505 : CALL calculate_BIb_C_2D(BIb_C_2D_gw(ispin)%array, BIb_C_gw(ispin)%array, para_env_sub, dimen_nm_gw, &
506 : gw_corr_lev_occ(ispin) + gw_corr_lev_virt(ispin), nmo, gd_B_all, &
507 50 : my_nm_gw_size, my_nm_gw_start, my_nm_gw_end, my_group_L_size)
508 96 : DEALLOCATE (BIb_C_gw(ispin)%array)
509 : END DO
510 :
511 46 : CALL release_group_dist(gd_B_all)
512 :
513 92 : CALL timestop(handle3)
514 :
515 : END IF
516 : END IF
517 :
518 308 : IF (do_bse) THEN
519 :
520 18 : CALL timeset(routineN//"_reorder_bse1", handle3)
521 :
522 18 : dimen_homo_square = homo(1)**2
523 : ! We do not implement an explicit bse_lev_occ different to homo here, because the small number of occupied levels
524 : ! does not critically influence the memory
525 : CALL calculate_BIb_C_2D(BIb_C_2D_bse_ij(1)%array, BIb_C_bse_ij, para_env_sub, dimen_homo_square, &
526 : homo(1), homo(1), gd_B_occ_bse, &
527 18 : my_ij_comb_bse_size, my_ij_comb_bse_start, my_ij_comb_bse_end, my_group_L_size)
528 :
529 18 : DEALLOCATE (BIb_C_bse_ij)
530 18 : CALL release_group_dist(gd_B_occ_bse)
531 :
532 18 : CALL timestop(handle3)
533 :
534 18 : CALL timeset(routineN//"_reorder_bse2", handle3)
535 :
536 18 : dimen_virt_square = bse_lev_virt**2
537 :
538 : CALL calculate_BIb_C_2D(BIb_C_2D_bse_ab(1)%array, BIb_C_bse_ab, para_env_sub, dimen_virt_square, &
539 : bse_lev_virt, bse_lev_virt, gd_B_virt_bse, &
540 18 : my_ab_comb_bse_size, my_ab_comb_bse_start, my_ab_comb_bse_end, my_group_L_size)
541 :
542 18 : DEALLOCATE (BIb_C_bse_ab)
543 18 : CALL release_group_dist(gd_B_virt_bse)
544 :
545 18 : CALL timestop(handle3)
546 :
547 : END IF
548 :
549 272 : CALL timestop(handle2)
550 :
551 272 : IF (num_integ_group > 1) THEN
552 82 : ALLOCATE (para_env_RPA)
553 82 : CALL para_env_RPA%from_split(para_env, color_rpa_group)
554 : ELSE
555 190 : para_env_RPA => para_env
556 : END IF
557 :
558 : ! now create the matrices needed for the calculation, Q, S and G
559 : ! Q and G will have omega dependence
560 :
561 272 : IF (do_im_time) THEN
562 834 : ALLOCATE (fm_mat_Q(nspins), fm_mat_Q_gemm(1), fm_mat_S(1))
563 : ELSE
564 1176 : ALLOCATE (fm_mat_Q(nspins), fm_mat_Q_gemm(nspins), fm_mat_S(nspins))
565 : END IF
566 :
567 : CALL create_integ_mat(BIb_C_2D, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
568 : dimen_RI_red, dimen_ia, color_rpa_group, &
569 : mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
570 : my_ia_size, my_ia_start, my_ia_end, &
571 : my_group_L_size, my_group_L_start, my_group_L_end, &
572 : para_env_RPA, fm_mat_S, nrow_block_mat, ncol_block_mat, &
573 : dimen_ia_for_block_size=dimen_ia(1), &
574 272 : do_im_time=do_im_time, fm_mat_Q_gemm=fm_mat_Q_gemm, fm_mat_Q=fm_mat_Q, qs_env=qs_env)
575 :
576 598 : DEALLOCATE (BIb_C_2D, my_ia_end, my_ia_size, my_ia_start)
577 :
578 : ! for GW, we need other matrix fm_mat_S, always allocate the container to prevent crying compilers
579 1142 : ALLOCATE (fm_mat_S_gw(nspins))
580 272 : IF (my_do_gw .AND. .NOT. do_im_time) THEN
581 :
582 : CALL create_integ_mat(BIb_C_2D_gw, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
583 : dimen_RI_red, [dimen_nm_gw, dimen_nm_gw], color_rpa_group, &
584 : mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
585 : [my_nm_gw_size, my_nm_gw_size], [my_nm_gw_start, my_nm_gw_start], [my_nm_gw_end, my_nm_gw_end], &
586 : my_group_L_size, my_group_L_start, my_group_L_end, &
587 : para_env_RPA, fm_mat_S_gw, nrow_block_mat, ncol_block_mat, &
588 : fm_mat_Q(1)%matrix_struct%context, fm_mat_Q(1)%matrix_struct%context, &
589 414 : fm_mat_Q=fm_mat_R_gw)
590 96 : DEALLOCATE (BIb_C_2D_gw)
591 :
592 : END IF
593 :
594 : ! for Bethe-Salpeter, we need other matrix fm_mat_S
595 272 : IF (do_bse) THEN
596 : CALL create_integ_mat(BIb_C_2D_bse_ij, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
597 : dimen_RI_red, [dimen_homo_square], color_rpa_group, &
598 : mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
599 : [my_ij_comb_bse_size], [my_ij_comb_bse_start], [my_ij_comb_bse_end], &
600 : my_group_L_size, my_group_L_start, my_group_L_end, &
601 : para_env_RPA, fm_mat_S_ij_bse, nrow_block_mat, ncol_block_mat, &
602 90 : fm_mat_Q(1)%matrix_struct%context, fm_mat_Q(1)%matrix_struct%context)
603 :
604 : CALL create_integ_mat(BIb_C_2D_bse_ab, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
605 : dimen_RI_red, [dimen_virt_square], color_rpa_group, &
606 : mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
607 : [my_ab_comb_bse_size], [my_ab_comb_bse_start], [my_ab_comb_bse_end], &
608 : my_group_L_size, my_group_L_start, my_group_L_end, &
609 : para_env_RPA, fm_mat_S_ab_bse, nrow_block_mat, ncol_block_mat, &
610 90 : fm_mat_Q(1)%matrix_struct%context, fm_mat_Q(1)%matrix_struct%context)
611 :
612 : END IF
613 :
614 272 : do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
615 272 : IF (do_kpoints_from_Gamma) THEN
616 18 : CALL get_bandstruc_and_k_dependent_MOs(qs_env, Eigenval_kp)
617 : END IF
618 :
619 : ! Now start the RPA calculation
620 : ! fm_mo_coeff_occ, fm_mo_coeff_virt will be deallocated here
621 : CALL rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_sub, unit_nr, &
622 : homo, virtual, dimen_RI, dimen_RI_red, dimen_ia, dimen_nm_gw, &
623 : Eigenval_kp, num_integ_points, num_integ_group, color_rpa_group, &
624 : fm_matrix_PQ, fm_mat_S, fm_mat_Q_gemm, fm_mat_Q, fm_mat_S_gw, fm_mat_R_gw(1), &
625 : fm_mat_S_ij_bse(1), fm_mat_S_ab_bse(1), &
626 : my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
627 : bse_lev_virt, &
628 : do_minimax_quad, &
629 : do_im_time, mo_coeff, &
630 : fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
631 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, mat_munu, mat_P_global, &
632 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
633 : starts_array_mc, ends_array_mc, &
634 : starts_array_mc_block, ends_array_mc_block, &
635 : matrix_s, do_kpoints_from_Gamma, kpoints, gd_array, color_sub, &
636 272 : do_ri_sos_laplace_mp2=do_ri_sos_laplace_mp2, calc_forces=calc_forces)
637 :
638 272 : CALL release_group_dist(gd_array)
639 :
640 272 : IF (num_integ_group > 1) CALL mp_para_env_release(para_env_RPA)
641 :
642 272 : IF (.NOT. do_im_time) THEN
643 138 : CALL cp_fm_release(fm_mat_Q_gemm)
644 138 : CALL cp_fm_release(fm_mat_S)
645 : END IF
646 272 : CALL cp_fm_release(fm_mat_Q)
647 :
648 272 : IF (my_do_gw .AND. .NOT. do_im_time) THEN
649 46 : CALL cp_fm_release(fm_mat_S_gw)
650 46 : CALL cp_fm_release(fm_mat_R_gw(1))
651 : END IF
652 :
653 272 : IF (do_bse) THEN
654 18 : CALL cp_fm_release(fm_mat_S_ij_bse(1))
655 18 : CALL cp_fm_release(fm_mat_S_ab_bse(1))
656 : END IF
657 :
658 272 : CALL timestop(handle)
659 :
660 1088 : END SUBROUTINE rpa_ri_compute_en
661 :
662 : ! **************************************************************************************************
663 : !> \brief reorder the local data in such a way to help the next stage of matrix creation;
664 : !> now the data inside the group are divided into a ia x K matrix (BIb_C_2D);
665 : !> Subroutine created to avoid massive double coding
666 : !> \param BIb_C_2D ...
667 : !> \param BIb_C ...
668 : !> \param para_env_sub ...
669 : !> \param dimen_ia ...
670 : !> \param homo ...
671 : !> \param virtual ...
672 : !> \param gd_B_virtual ...
673 : !> \param my_ia_size ...
674 : !> \param my_ia_start ...
675 : !> \param my_ia_end ...
676 : !> \param my_group_L_size ...
677 : !> \author Jan Wilhelm, 03/2015
678 : ! **************************************************************************************************
679 248 : SUBROUTINE calculate_BIb_C_2D(BIb_C_2D, BIb_C, para_env_sub, dimen_ia, homo, virtual, &
680 : gd_B_virtual, &
681 : my_ia_size, my_ia_start, my_ia_end, my_group_L_size)
682 :
683 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
684 : INTENT(OUT) :: BIb_C_2D
685 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
686 : INTENT(IN) :: BIb_C
687 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
688 : INTEGER, INTENT(IN) :: dimen_ia, homo, virtual
689 : TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_B_virtual
690 : INTEGER :: my_ia_size, my_ia_start, my_ia_end, &
691 : my_group_L_size
692 :
693 : INTEGER, PARAMETER :: occ_chunk = 128
694 :
695 : INTEGER :: ia_global, iiB, itmp(2), jjB, my_B_size, my_B_virtual_start, occ_high, occ_low, &
696 : proc_receive, proc_send, proc_shift, rec_B_size, rec_B_virtual_end, rec_B_virtual_start
697 248 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: BIb_C_rec_1D
698 248 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: BIb_C_rec
699 :
700 248 : itmp = get_limit(dimen_ia, para_env_sub%num_pe, para_env_sub%mepos)
701 248 : my_ia_start = itmp(1)
702 248 : my_ia_end = itmp(2)
703 248 : my_ia_size = my_ia_end - my_ia_start + 1
704 :
705 248 : CALL get_group_dist(gd_B_virtual, para_env_sub%mepos, sizes=my_B_size, starts=my_B_virtual_start)
706 :
707 : ! reorder data
708 990 : ALLOCATE (BIb_C_2D(my_group_L_size, my_ia_size))
709 :
710 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,ia_global) &
711 : !$OMP SHARED(homo,my_B_size,virtual,my_B_virtual_start,my_ia_start,my_ia_end,BIb_C,BIb_C_2D,&
712 248 : !$OMP my_group_L_size)
713 : DO iiB = 1, homo
714 : DO jjB = 1, my_B_size
715 : ia_global = (iiB - 1)*virtual + my_B_virtual_start + jjB - 1
716 : IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
717 : BIb_C_2D(1:my_group_L_size, ia_global - my_ia_start + 1) = BIb_C(1:my_group_L_size, jjB, iiB)
718 : END IF
719 : END DO
720 : END DO
721 :
722 248 : IF (para_env_sub%num_pe > 1) THEN
723 30 : ALLOCATE (BIb_C_rec_1D(INT(my_group_L_size, int_8)*maxsize(gd_B_virtual)*MIN(homo, occ_chunk)))
724 20 : DO proc_shift = 1, para_env_sub%num_pe - 1
725 10 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
726 10 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
727 :
728 10 : CALL get_group_dist(gd_B_virtual, proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
729 :
730 : ! do this in chunks to avoid high memory overhead
731 20 : DO occ_low = 1, homo, occ_chunk
732 10 : occ_high = MIN(homo, occ_low + occ_chunk - 1)
733 : BIb_C_rec(1:my_group_L_size, 1:rec_B_size, 1:occ_high - occ_low + 1) => &
734 10 : BIb_C_rec_1D(1:INT(my_group_L_size, int_8)*rec_B_size*(occ_high - occ_low + 1))
735 : CALL para_env_sub%sendrecv(BIb_C(:, :, occ_low:occ_high), proc_send, &
736 31970 : BIb_C_rec(:, :, 1:occ_high - occ_low + 1), proc_receive)
737 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,ia_global) &
738 : !$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,&
739 10 : !$OMP my_group_L_size)
740 : DO iiB = occ_low, occ_high
741 : DO jjB = 1, rec_B_size
742 : ia_global = (iiB - 1)*virtual + rec_B_virtual_start + jjB - 1
743 : IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
744 : 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)
745 : END IF
746 : END DO
747 : END DO
748 : END DO
749 :
750 : END DO
751 10 : DEALLOCATE (BIb_C_rec_1D)
752 : END IF
753 :
754 248 : END SUBROUTINE calculate_BIb_C_2D
755 :
756 : ! **************************************************************************************************
757 : !> \brief ...
758 : !> \param BIb_C_2D ...
759 : !> \param para_env ...
760 : !> \param para_env_sub ...
761 : !> \param color_sub ...
762 : !> \param ngroup ...
763 : !> \param integ_group_size ...
764 : !> \param dimen_RI ...
765 : !> \param dimen_ia ...
766 : !> \param color_rpa_group ...
767 : !> \param ext_row_block_size ...
768 : !> \param ext_col_block_size ...
769 : !> \param unit_nr ...
770 : !> \param my_ia_size ...
771 : !> \param my_ia_start ...
772 : !> \param my_ia_end ...
773 : !> \param my_group_L_size ...
774 : !> \param my_group_L_start ...
775 : !> \param my_group_L_end ...
776 : !> \param para_env_RPA ...
777 : !> \param fm_mat_S ...
778 : !> \param nrow_block_mat ...
779 : !> \param ncol_block_mat ...
780 : !> \param blacs_env_ext ...
781 : !> \param blacs_env_ext_S ...
782 : !> \param dimen_ia_for_block_size ...
783 : !> \param do_im_time ...
784 : !> \param fm_mat_Q_gemm ...
785 : !> \param fm_mat_Q ...
786 : !> \param qs_env ...
787 : ! **************************************************************************************************
788 354 : SUBROUTINE create_integ_mat(BIb_C_2D, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
789 354 : dimen_RI, dimen_ia, color_rpa_group, &
790 : ext_row_block_size, ext_col_block_size, unit_nr, &
791 354 : my_ia_size, my_ia_start, my_ia_end, &
792 : my_group_L_size, my_group_L_start, my_group_L_end, &
793 354 : para_env_RPA, fm_mat_S, nrow_block_mat, ncol_block_mat, &
794 : blacs_env_ext, blacs_env_ext_S, dimen_ia_for_block_size, &
795 354 : do_im_time, fm_mat_Q_gemm, fm_mat_Q, qs_env)
796 :
797 : TYPE(two_dim_real_array), DIMENSION(:), &
798 : INTENT(INOUT) :: BIb_C_2D
799 : TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_sub
800 : INTEGER, INTENT(IN) :: color_sub, ngroup, integ_group_size, &
801 : dimen_RI
802 : INTEGER, DIMENSION(:), INTENT(IN) :: dimen_ia
803 : INTEGER, INTENT(IN) :: color_rpa_group, ext_row_block_size, &
804 : ext_col_block_size, unit_nr
805 : INTEGER, DIMENSION(:), INTENT(IN) :: my_ia_size, my_ia_start, my_ia_end
806 : INTEGER, INTENT(IN) :: my_group_L_size, my_group_L_start, &
807 : my_group_L_end
808 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_RPA
809 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_S
810 : INTEGER, INTENT(INOUT) :: nrow_block_mat, ncol_block_mat
811 : TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: blacs_env_ext, blacs_env_ext_S
812 : INTEGER, INTENT(IN), OPTIONAL :: dimen_ia_for_block_size
813 : LOGICAL, INTENT(IN), OPTIONAL :: do_im_time
814 : TYPE(cp_fm_type), DIMENSION(:), OPTIONAL :: fm_mat_Q_gemm, fm_mat_Q
815 : TYPE(qs_environment_type), INTENT(IN), OPTIONAL, &
816 : POINTER :: qs_env
817 :
818 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_integ_mat'
819 :
820 : INTEGER :: col_row_proc_ratio, grid_2D(2), handle, &
821 : iproc, iproc_col, iproc_row, ispin, &
822 : mepos_in_RPA_group
823 354 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: group_grid_2_mepos
824 : LOGICAL :: my_blacs_ext, my_blacs_S_ext, &
825 : my_do_im_time
826 : TYPE(cp_blacs_env_type), POINTER :: blacs_env, blacs_env_Q
827 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
828 354 : TYPE(group_dist_d1_type) :: gd_ia, gd_L
829 :
830 354 : CALL timeset(routineN, handle)
831 :
832 354 : CPASSERT(PRESENT(blacs_env_ext) .OR. PRESENT(dimen_ia_for_block_size))
833 :
834 354 : my_blacs_ext = .FALSE.
835 354 : IF (PRESENT(blacs_env_ext)) my_blacs_ext = .TRUE.
836 :
837 354 : my_blacs_S_ext = .FALSE.
838 354 : IF (PRESENT(blacs_env_ext_S)) my_blacs_S_ext = .TRUE.
839 :
840 354 : my_do_im_time = .FALSE.
841 354 : IF (PRESENT(do_im_time)) my_do_im_time = do_im_time
842 :
843 354 : NULLIFY (blacs_env)
844 : ! create the RPA blacs env
845 354 : IF (my_blacs_S_ext) THEN
846 82 : blacs_env => blacs_env_ext_S
847 : ELSE
848 272 : IF (para_env_RPA%num_pe > 1) THEN
849 190 : col_row_proc_ratio = MAX(1, dimen_ia_for_block_size/dimen_RI)
850 :
851 190 : 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
852 190 : DO iproc = 1, para_env_RPA%num_pe
853 190 : iproc_col = iproc_col - 1
854 190 : IF (MOD(para_env_RPA%num_pe, iproc_col) == 0) EXIT
855 : END DO
856 :
857 190 : iproc_row = para_env_RPA%num_pe/iproc_col
858 190 : grid_2D(1) = iproc_row
859 190 : grid_2D(2) = iproc_col
860 : ELSE
861 246 : grid_2D = 1
862 : END IF
863 272 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_RPA, grid_2d=grid_2D)
864 :
865 272 : IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
866 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
867 69 : "MATRIX_INFO| Number row processes:", grid_2D(1)
868 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
869 69 : "MATRIX_INFO| Number column processes:", grid_2D(2)
870 : END IF
871 :
872 : ! define the block_size for the row
873 272 : IF (ext_row_block_size > 0) THEN
874 0 : nrow_block_mat = ext_row_block_size
875 : ELSE
876 272 : nrow_block_mat = MAX(1, dimen_RI/grid_2D(1)/2)
877 : END IF
878 :
879 : ! define the block_size for the column
880 272 : IF (ext_col_block_size > 0) THEN
881 0 : ncol_block_mat = ext_col_block_size
882 : ELSE
883 272 : ncol_block_mat = MAX(1, dimen_ia_for_block_size/grid_2D(2)/2)
884 : END IF
885 :
886 272 : IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
887 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
888 69 : "MATRIX_INFO| Row block size:", nrow_block_mat
889 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
890 69 : "MATRIX_INFO| Column block size:", ncol_block_mat
891 : END IF
892 : END IF
893 :
894 287 : IF (.NOT. my_do_im_time) THEN
895 468 : DO ispin = 1, SIZE(BIb_C_2D)
896 248 : NULLIFY (fm_struct)
897 248 : IF (my_blacs_ext) THEN
898 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
899 86 : ncol_global=dimen_ia(ispin), para_env=para_env_RPA)
900 : ELSE
901 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
902 : ncol_global=dimen_ia(ispin), para_env=para_env_RPA, &
903 162 : nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.TRUE.)
904 :
905 : END IF ! external blacs_env
906 :
907 248 : CALL create_group_dist(gd_ia, my_ia_start(ispin), my_ia_end(ispin), my_ia_size(ispin), para_env_RPA)
908 248 : CALL create_group_dist(gd_L, my_group_L_start, my_group_L_end, my_group_L_size, para_env_RPA)
909 :
910 : ! create the info array
911 :
912 248 : mepos_in_RPA_group = MOD(color_sub, integ_group_size)
913 992 : ALLOCATE (group_grid_2_mepos(0:integ_group_size - 1, 0:para_env_sub%num_pe - 1))
914 872 : group_grid_2_mepos = 0
915 248 : group_grid_2_mepos(mepos_in_RPA_group, para_env_sub%mepos) = para_env_RPA%mepos
916 248 : CALL para_env_RPA%sum(group_grid_2_mepos)
917 :
918 : CALL array2fm(BIb_C_2D(ispin)%array, fm_struct, my_group_L_start, my_group_L_end, &
919 : my_ia_start(ispin), my_ia_end(ispin), gd_L, gd_ia, &
920 : group_grid_2_mepos, ngroup, para_env_sub%num_pe, fm_mat_S(ispin), &
921 248 : integ_group_size, color_rpa_group)
922 :
923 248 : DEALLOCATE (group_grid_2_mepos)
924 248 : CALL cp_fm_struct_release(fm_struct)
925 :
926 : ! deallocate the info array
927 248 : CALL release_group_dist(gd_L)
928 248 : CALL release_group_dist(gd_ia)
929 :
930 : ! sum the local data across processes belonging to different RPA group.
931 468 : IF (para_env_RPA%num_pe /= para_env%num_pe) THEN
932 : BLOCK
933 : TYPE(mp_comm_type) :: comm_exchange
934 130 : comm_exchange = fm_mat_S(ispin)%matrix_struct%context%interconnect(para_env)
935 130 : CALL comm_exchange%sum(fm_mat_S(ispin)%local_data)
936 260 : CALL comm_exchange%free()
937 : END BLOCK
938 : END IF
939 : END DO
940 : END IF
941 :
942 354 : IF (PRESENT(fm_mat_Q_gemm) .AND. .NOT. my_do_im_time) THEN
943 : ! create the Q matrix dimen_RIxdimen_RI where the result of the mat-mat-mult will be stored
944 138 : NULLIFY (fm_struct)
945 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
946 : ncol_global=dimen_RI, para_env=para_env_RPA, &
947 138 : nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.TRUE.)
948 300 : DO ispin = 1, SIZE(fm_mat_Q_gemm)
949 300 : CALL cp_fm_create(fm_mat_Q_gemm(ispin), fm_struct, name="fm_mat_Q_gemm")
950 : END DO
951 138 : CALL cp_fm_struct_release(fm_struct)
952 : END IF
953 :
954 354 : IF (PRESENT(fm_mat_Q)) THEN
955 318 : NULLIFY (blacs_env_Q)
956 318 : IF (my_blacs_ext) THEN
957 46 : blacs_env_Q => blacs_env_ext
958 272 : ELSE IF (para_env_RPA%num_pe == para_env%num_pe .AND. PRESENT(qs_env)) THEN
959 190 : CALL get_qs_env(qs_env, blacs_env=blacs_env_Q)
960 : ELSE
961 82 : CALL cp_blacs_env_create(blacs_env=blacs_env_Q, para_env=para_env_RPA)
962 : END IF
963 318 : NULLIFY (fm_struct)
964 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_Q, nrow_global=dimen_RI, &
965 318 : ncol_global=dimen_RI, para_env=para_env_RPA)
966 690 : DO ispin = 1, SIZE(fm_mat_Q)
967 690 : CALL cp_fm_create(fm_mat_Q(ispin), fm_struct, name="fm_mat_Q")
968 : END DO
969 :
970 318 : CALL cp_fm_struct_release(fm_struct)
971 :
972 318 : IF (.NOT. (my_blacs_ext .OR. (para_env_RPA%num_pe == para_env%num_pe .AND. PRESENT(qs_env)))) &
973 82 : CALL cp_blacs_env_release(blacs_env_Q)
974 : END IF
975 :
976 : ! release blacs_env
977 354 : IF (.NOT. my_blacs_S_ext) CALL cp_blacs_env_release(blacs_env)
978 :
979 354 : CALL timestop(handle)
980 :
981 354 : END SUBROUTINE create_integ_mat
982 :
983 : ! **************************************************************************************************
984 : !> \brief ...
985 : !> \param qs_env ...
986 : !> \param Erpa ...
987 : !> \param mp2_env ...
988 : !> \param para_env ...
989 : !> \param para_env_RPA ...
990 : !> \param para_env_sub ...
991 : !> \param unit_nr ...
992 : !> \param homo ...
993 : !> \param virtual ...
994 : !> \param dimen_RI ...
995 : !> \param dimen_RI_red ...
996 : !> \param dimen_ia ...
997 : !> \param dimen_nm_gw ...
998 : !> \param Eigenval ...
999 : !> \param num_integ_points ...
1000 : !> \param num_integ_group ...
1001 : !> \param color_rpa_group ...
1002 : !> \param fm_matrix_PQ ...
1003 : !> \param fm_mat_S ...
1004 : !> \param fm_mat_Q_gemm ...
1005 : !> \param fm_mat_Q ...
1006 : !> \param fm_mat_S_gw ...
1007 : !> \param fm_mat_R_gw ...
1008 : !> \param fm_mat_S_ij_bse ...
1009 : !> \param fm_mat_S_ab_bse ...
1010 : !> \param my_do_gw ...
1011 : !> \param do_bse ...
1012 : !> \param gw_corr_lev_occ ...
1013 : !> \param gw_corr_lev_virt ...
1014 : !> \param bse_lev_virt ...
1015 : !> \param do_minimax_quad ...
1016 : !> \param do_im_time ...
1017 : !> \param mo_coeff ...
1018 : !> \param fm_matrix_L_kpoints ...
1019 : !> \param fm_matrix_Minv_L_kpoints ...
1020 : !> \param fm_matrix_Minv ...
1021 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
1022 : !> \param mat_munu ...
1023 : !> \param mat_P_global ...
1024 : !> \param t_3c_M ...
1025 : !> \param t_3c_O ...
1026 : !> \param t_3c_O_compressed ...
1027 : !> \param t_3c_O_ind ...
1028 : !> \param starts_array_mc ...
1029 : !> \param ends_array_mc ...
1030 : !> \param starts_array_mc_block ...
1031 : !> \param ends_array_mc_block ...
1032 : !> \param matrix_s ...
1033 : !> \param do_kpoints_from_Gamma ...
1034 : !> \param kpoints ...
1035 : !> \param gd_array ...
1036 : !> \param color_sub ...
1037 : !> \param do_ri_sos_laplace_mp2 ...
1038 : !> \param calc_forces ...
1039 : ! **************************************************************************************************
1040 272 : SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_sub, unit_nr, &
1041 272 : homo, virtual, dimen_RI, dimen_RI_red, dimen_ia, dimen_nm_gw, &
1042 : Eigenval, num_integ_points, num_integ_group, color_rpa_group, &
1043 544 : fm_matrix_PQ, fm_mat_S, fm_mat_Q_gemm, fm_mat_Q, fm_mat_S_gw, fm_mat_R_gw, &
1044 : fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
1045 272 : my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
1046 : bse_lev_virt, &
1047 272 : do_minimax_quad, do_im_time, mo_coeff, &
1048 : fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
1049 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, mat_munu, mat_P_global, &
1050 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1051 : starts_array_mc, ends_array_mc, &
1052 : starts_array_mc_block, ends_array_mc_block, &
1053 : matrix_s, do_kpoints_from_Gamma, kpoints, gd_array, color_sub, &
1054 : do_ri_sos_laplace_mp2, calc_forces)
1055 :
1056 : TYPE(qs_environment_type), POINTER :: qs_env
1057 : REAL(KIND=dp), INTENT(OUT) :: Erpa
1058 : TYPE(mp2_type) :: mp2_env
1059 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_RPA, para_env_sub
1060 : INTEGER, INTENT(IN) :: unit_nr
1061 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
1062 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
1063 : INTEGER, DIMENSION(:), INTENT(IN) :: dimen_ia
1064 : INTEGER, INTENT(IN) :: dimen_nm_gw
1065 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1066 : INTENT(INOUT) :: Eigenval
1067 : INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, &
1068 : color_rpa_group
1069 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_PQ
1070 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_S
1071 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_Q_gemm, fm_mat_Q, fm_mat_S_gw
1072 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_R_gw, fm_mat_S_ij_bse, &
1073 : fm_mat_S_ab_bse
1074 : LOGICAL, INTENT(IN) :: my_do_gw, do_bse
1075 : INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt
1076 : INTEGER, INTENT(IN) :: bse_lev_virt
1077 : LOGICAL, INTENT(IN) :: do_minimax_quad, do_im_time
1078 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1079 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_L_kpoints, &
1080 : fm_matrix_Minv_L_kpoints, &
1081 : fm_matrix_Minv, &
1082 : fm_matrix_Minv_Vtrunc_Minv
1083 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
1084 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_P_global
1085 : TYPE(dbt_type), INTENT(INOUT) :: t_3c_M
1086 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
1087 : INTENT(INOUT) :: t_3c_O
1088 : TYPE(hfx_compression_type), ALLOCATABLE, &
1089 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_compressed
1090 : TYPE(block_ind_type), ALLOCATABLE, &
1091 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_ind
1092 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
1093 : starts_array_mc_block, &
1094 : ends_array_mc_block
1095 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1096 : LOGICAL :: do_kpoints_from_Gamma
1097 : TYPE(kpoint_type), POINTER :: kpoints
1098 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1099 : INTEGER, INTENT(IN) :: color_sub
1100 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, calc_forces
1101 :
1102 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_num_int'
1103 :
1104 : COMPLEX(KIND=dp), ALLOCATABLE, &
1105 272 : DIMENSION(:, :, :, :) :: vec_Sigma_c_gw
1106 : INTEGER :: count_ev_sc_GW, cut_memory, group_size_P, gw_corr_lev_tot, handle, handle3, i, &
1107 : ikp_local, ispin, iter_evGW, iter_sc_GW0, j, jquad, min_bsize, mm_style, nkp, &
1108 : nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, Pspin, Qspin, &
1109 : size_P
1110 : INTEGER(int_8) :: dbcsr_nflop
1111 272 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c
1112 272 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c
1113 544 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, prim_blk_sizes, &
1114 272 : RI_blk_sizes
1115 : LOGICAL :: do_apply_ic_corr_to_gw, do_gw_im_time, do_ic_model, do_kpoints_cubic_RPA, &
1116 : do_periodic, do_print, do_ri_Sigma_x, exit_ev_gw, first_cycle, &
1117 : first_cycle_periodic_correction, my_open_shell, print_ic_values
1118 272 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: has_mat_P_blocks
1119 : REAL(KIND=dp) :: a_scaling, alpha, dbcsr_time, e_exchange, e_exchange_corr, eps_filter, &
1120 : eps_filter_im_time, ext_scaling, fermi_level_offset, fermi_level_offset_input, &
1121 : my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old
1122 544 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, e_fermi, tau_tj, tau_wj, tj, &
1123 272 : trace_Qomega, vec_omega_fit_gw, wj, &
1124 272 : wkp_W
1125 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_W_gw, weights_cos_tf_t_to_w, &
1126 272 : weights_cos_tf_w_to_t, &
1127 272 : weights_sin_tf_t_to_w
1128 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_last, Eigenval_scf, &
1129 272 : vec_Sigma_x_gw
1130 : TYPE(cp_cfm_type) :: cfm_mat_Q
1131 : TYPE(cp_fm_type) :: fm_mat_Q_static_bse_gemm, fm_mat_RI_global_work, fm_mat_S_ia_bse, &
1132 : fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, &
1133 : fm_scaled_dm_virt_tau
1134 272 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_S_gw_work, fm_mat_W, &
1135 272 : fm_mo_coeff_occ, fm_mo_coeff_virt
1136 272 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_L_kpoints, fm_mat_Minv_L_kpoints
1137 : TYPE(dbcsr_p_type) :: mat_dm, mat_L, mat_M_P_munu_occ, &
1138 : mat_M_P_munu_virt, mat_MinvVMinv
1139 : TYPE(dbcsr_p_type), ALLOCATABLE, &
1140 272 : DIMENSION(:, :, :) :: mat_P_omega
1141 272 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_im_mo_mo, &
1142 272 : matrix_berry_re_mo_mo
1143 272 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_kp
1144 : TYPE(dbcsr_type), POINTER :: mat_W, mat_work
1145 1904 : TYPE(dbt_type) :: t_3c_overl_int_ao_mo
1146 272 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_3c_overl_int_gw_AO, &
1147 272 : t_3c_overl_int_gw_RI, &
1148 272 : t_3c_overl_nnP_ic, &
1149 272 : t_3c_overl_nnP_ic_reflected
1150 : TYPE(dgemm_counter_type) :: dgemm_counter
1151 : TYPE(hfx_compression_type), ALLOCATABLE, &
1152 272 : DIMENSION(:) :: t_3c_O_mo_compressed
1153 20128 : TYPE(im_time_force_type) :: force_data
1154 272 : TYPE(rpa_exchange_work_type) :: exchange_work
1155 1360 : TYPE(rpa_grad_type) :: rpa_grad
1156 272 : TYPE(two_dim_int_array), ALLOCATABLE, DIMENSION(:) :: t_3c_O_mo_ind
1157 :
1158 272 : CALL timeset(routineN, handle)
1159 :
1160 272 : nspins = SIZE(homo)
1161 272 : nmo = homo(1) + virtual(1)
1162 272 : my_open_shell = (nspins == 2)
1163 :
1164 272 : do_gw_im_time = my_do_gw .AND. do_im_time
1165 272 : do_ri_Sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x
1166 272 : do_ic_model = mp2_env%ri_g0w0%do_ic_model
1167 272 : print_ic_values = mp2_env%ri_g0w0%print_ic_values
1168 272 : do_periodic = mp2_env%ri_g0w0%do_periodic
1169 272 : do_kpoints_cubic_RPA = mp2_env%ri_rpa_im_time%do_im_time_kpoints
1170 :
1171 : ! For SOS-MP2 only gemm is implemented
1172 272 : mm_style = wfc_mm_style_gemm
1173 272 : IF (.NOT. do_ri_sos_laplace_mp2) mm_style = mp2_env%ri_rpa%mm_style
1174 :
1175 272 : IF (my_do_gw) THEN
1176 92 : ext_scaling = 0.2_dp
1177 92 : omega_max_fit = mp2_env%ri_g0w0%omega_max_fit
1178 92 : fermi_level_offset_input = mp2_env%ri_g0w0%fermi_level_offset
1179 92 : iter_evGW = mp2_env%ri_g0w0%iter_evGW
1180 92 : iter_sc_GW0 = mp2_env%ri_g0w0%iter_sc_GW0
1181 92 : IF ((.NOT. do_im_time)) THEN
1182 46 : IF (iter_sc_GW0 .NE. 1 .AND. iter_evGW .NE. 1) CPABORT("Mixed scGW0/evGW not implemented.")
1183 : ! in case of scGW0 with the N^4 algorithm, we use the evGW code but use the DFT eigenvalues for W
1184 46 : IF (iter_sc_GW0 .NE. 1) iter_evGW = iter_sc_GW0
1185 : END IF
1186 : ELSE
1187 180 : ext_scaling = 0.0_dp
1188 180 : iter_evGW = 1
1189 180 : iter_sc_GW0 = 1
1190 : END IF
1191 :
1192 272 : IF (do_kpoints_cubic_RPA .AND. do_ri_sos_laplace_mp2) THEN
1193 0 : CPABORT("RI-SOS-Laplace-MP2 with k-point-sampling is not implemented.")
1194 : END IF
1195 :
1196 272 : do_apply_ic_corr_to_gw = .FALSE.
1197 272 : IF (mp2_env%ri_g0w0%ic_corr_list(1)%array(1) > 0.0_dp) do_apply_ic_corr_to_gw = .TRUE.
1198 :
1199 272 : IF (do_im_time) THEN
1200 134 : CPASSERT(do_minimax_quad .OR. do_ri_sos_laplace_mp2)
1201 :
1202 134 : group_size_P = mp2_env%ri_rpa_im_time%group_size_P
1203 134 : cut_memory = mp2_env%ri_rpa_im_time%cut_memory
1204 134 : eps_filter = mp2_env%ri_rpa_im_time%eps_filter
1205 : eps_filter_im_time = mp2_env%ri_rpa_im_time%eps_filter* &
1206 134 : mp2_env%ri_rpa_im_time%eps_filter_factor
1207 :
1208 134 : min_bsize = mp2_env%ri_rpa_im_time%min_bsize
1209 :
1210 : CALL alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, &
1211 : num_integ_points, nspins, fm_mat_Q(1), fm_mo_coeff_occ, fm_mo_coeff_virt, &
1212 : fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, mat_P_global, &
1213 : t_3c_O, matrix_s, kpoints, eps_filter_im_time, &
1214 : cut_memory, nkp, num_cells_dm, num_3c_repl, &
1215 : size_P, ikp_local, &
1216 : index_to_cell_3c, &
1217 : cell_to_index_3c, &
1218 : col_blk_size, &
1219 : do_ic_model, do_kpoints_cubic_RPA, &
1220 : do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, &
1221 : has_mat_P_blocks, wkp_W, &
1222 : cfm_mat_Q, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
1223 : fm_mat_RI_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, &
1224 : fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, &
1225 : mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, mat_work, mo_coeff, &
1226 134 : fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, homo, nmo)
1227 :
1228 134 : IF (calc_forces) CALL init_im_time_forces(force_data, fm_matrix_PQ, t_3c_M, unit_nr, mp2_env, qs_env)
1229 :
1230 134 : IF (my_do_gw) THEN
1231 :
1232 : CALL dbcsr_get_info(mat_P_global%matrix, &
1233 46 : row_blk_size=RI_blk_sizes)
1234 :
1235 : CALL dbcsr_get_info(matrix_s(1)%matrix, &
1236 46 : row_blk_size=prim_blk_sizes)
1237 :
1238 46 : gw_corr_lev_tot = gw_corr_lev_occ(1) + gw_corr_lev_virt(1)
1239 :
1240 46 : IF (.NOT. do_kpoints_cubic_RPA) THEN
1241 : CALL allocate_matrices_gw_im_time(gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, &
1242 : num_integ_points, unit_nr, &
1243 : RI_blk_sizes, do_ic_model, &
1244 : para_env, fm_mat_W, fm_mat_Q(1), &
1245 : mo_coeff, &
1246 : t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
1247 : t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
1248 : starts_array_mc, ends_array_mc, &
1249 : t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
1250 : matrix_s, mat_W, t_3c_O, &
1251 : t_3c_O_compressed, t_3c_O_ind, &
1252 46 : qs_env)
1253 :
1254 : END IF
1255 : END IF
1256 :
1257 : END IF
1258 :
1259 272 : IF (do_ic_model) THEN
1260 : ! image charge model only implemented for cubic scaling GW
1261 2 : CPASSERT(do_gw_im_time)
1262 2 : CPASSERT(.NOT. do_periodic)
1263 2 : IF (cut_memory .NE. 1) CPABORT("For IC, use MEMORY_CUT 1 in the LOW_SCALING section.")
1264 : END IF
1265 :
1266 816 : ALLOCATE (e_fermi(nspins))
1267 272 : IF (do_minimax_quad .OR. do_ri_sos_laplace_mp2) THEN
1268 194 : do_print = .NOT. do_ic_model
1269 : CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, &
1270 : do_ri_sos_laplace_mp2, do_print, &
1271 : tau_tj, tau_wj, qs_env, do_gw_im_time, &
1272 : do_kpoints_cubic_RPA, e_fermi(1), tj, wj, &
1273 : weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
1274 194 : qs_env%mp2_env%ri_g0w0%regularization_minimax)
1275 :
1276 : !For sos_laplace_mp2 and low-scaling RPA, potentially need to store/retrieve the initial weights
1277 194 : IF (qs_env%mp2_env%ri_rpa_im_time%keep_quad) THEN
1278 : CALL keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, &
1279 : weights_cos_tf_w_to_t, do_ri_sos_laplace_mp2, do_im_time, &
1280 194 : num_integ_points, unit_nr, qs_env)
1281 : END IF
1282 : ELSE
1283 78 : IF (calc_forces) CPABORT("Forces with Clenshaw-Curtis grid not implemented.")
1284 : CALL get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, &
1285 : num_integ_group, color_rpa_group, fm_mat_S, my_do_gw, &
1286 78 : ext_scaling, a_scaling, tj, wj)
1287 : END IF
1288 :
1289 : ! This array is needed for RPA
1290 272 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1291 642 : ALLOCATE (trace_Qomega(dimen_RI_red))
1292 : END IF
1293 :
1294 272 : IF (do_ri_sos_laplace_mp2 .AND. .NOT. do_im_time) THEN
1295 28 : alpha = 1.0_dp
1296 244 : ELSE IF (my_open_shell .OR. do_ri_sos_laplace_mp2) THEN
1297 68 : alpha = 2.0_dp
1298 : ELSE
1299 176 : alpha = 4.0_dp
1300 : END IF
1301 :
1302 272 : IF (my_do_gw) THEN
1303 : CALL allocate_matrices_gw(vec_Sigma_c_gw, color_rpa_group, dimen_nm_gw, &
1304 : gw_corr_lev_occ, gw_corr_lev_virt, homo, &
1305 : nmo, num_integ_group, num_integ_points, unit_nr, &
1306 : gw_corr_lev_tot, num_fit_points, omega_max_fit, &
1307 : do_minimax_quad, do_periodic, do_ri_Sigma_x,.NOT. do_im_time, &
1308 : first_cycle_periodic_correction, &
1309 : a_scaling, Eigenval, tj, vec_omega_fit_gw, vec_Sigma_x_gw, &
1310 : delta_corr, Eigenval_last, Eigenval_scf, vec_W_gw, &
1311 : fm_mat_S_gw, fm_mat_S_gw_work, &
1312 : para_env, mp2_env, kpoints, nkp, nkp_self_energy, &
1313 92 : do_kpoints_cubic_RPA, do_kpoints_from_Gamma)
1314 :
1315 92 : IF (do_bse) THEN
1316 :
1317 18 : CALL cp_fm_create(fm_mat_Q_static_bse_gemm, fm_mat_Q_gemm(1)%matrix_struct)
1318 18 : CALL cp_fm_to_fm(fm_mat_Q_gemm(1), fm_mat_Q_static_bse_gemm)
1319 18 : CALL cp_fm_set_all(fm_mat_Q_static_bse_gemm, 0.0_dp)
1320 :
1321 : END IF
1322 :
1323 : END IF
1324 :
1325 272 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_create(rpa_grad, fm_mat_Q(1), &
1326 : fm_mat_S, homo, virtual, mp2_env, Eigenval(:, 1, :), &
1327 40 : unit_nr, do_ri_sos_laplace_mp2)
1328 :
1329 272 : IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) &
1330 : CALL exchange_work%create(qs_env, para_env_sub, mat_munu, dimen_RI_red, &
1331 110 : fm_mat_S, fm_mat_Q(1), fm_mat_Q_gemm(1), homo, virtual)
1332 :
1333 272 : Erpa = 0.0_dp
1334 272 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) e_exchange = 0.0_dp
1335 272 : first_cycle = .TRUE.
1336 272 : omega_old = 0.0_dp
1337 272 : CALL dgemm_counter_init(dgemm_counter, unit_nr, mp2_env%ri_rpa%print_dgemm_info)
1338 :
1339 642 : DO count_ev_sc_GW = 1, iter_evGW
1340 :
1341 386 : dbcsr_time = 0.0_dp
1342 386 : dbcsr_nflop = 0
1343 :
1344 386 : IF (do_ic_model) CYCLE
1345 :
1346 : ! reset some values, important when doing eigenvalue self-consistent GW
1347 384 : IF (my_do_gw) THEN
1348 204 : Erpa = 0.0_dp
1349 115390 : vec_Sigma_c_gw = z_zero
1350 204 : first_cycle = .TRUE.
1351 : END IF
1352 :
1353 : ! calculate Q_PQ(it)
1354 384 : IF (do_im_time) THEN
1355 :
1356 146 : IF (.NOT. do_kpoints_cubic_RPA) THEN
1357 314 : DO ispin = 1, nspins
1358 314 : e_fermi(ispin) = (Eigenval(homo(ispin), 1, ispin) + Eigenval(homo(ispin) + 1, 1, ispin))*0.5_dp
1359 : END DO
1360 : END IF
1361 :
1362 146 : tau = 0.0_dp
1363 146 : tau_old = 0.0_dp
1364 :
1365 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(/T3,A,T66,i15)") &
1366 73 : "MEMORY_INFO| Memory cut:", cut_memory
1367 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
1368 73 : "SPARSITY_INFO| Eps filter for M virt/occ tensors:", eps_filter
1369 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
1370 73 : "SPARSITY_INFO| Eps filter for P matrix:", eps_filter_im_time
1371 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,i15)") &
1372 73 : "SPARSITY_INFO| Minimum tensor block size:", min_bsize
1373 :
1374 : ! for evGW, we have to ensure that mat_P_omega is zero
1375 146 : CALL zero_mat_P_omega(mat_P_omega(:, :, 1))
1376 :
1377 : ! compute the matrix Q(it) and Fourier transform it directly to mat_P_omega(iw)
1378 : CALL compute_mat_P_omega(mat_P_omega(:, :, 1), fm_scaled_dm_occ_tau, &
1379 : fm_scaled_dm_virt_tau, fm_mo_coeff_occ(1), fm_mo_coeff_virt(1), &
1380 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1381 : mat_P_global, matrix_s, 1, &
1382 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1383 : starts_array_mc, ends_array_mc, &
1384 : starts_array_mc_block, ends_array_mc_block, &
1385 : weights_cos_tf_t_to_w, tj, tau_tj, e_fermi(1), eps_filter, alpha, &
1386 : eps_filter_im_time, Eigenval(:, 1, 1), nmo, &
1387 : num_integ_points, cut_memory, &
1388 : unit_nr, mp2_env, para_env, &
1389 : qs_env, do_kpoints_from_Gamma, &
1390 : index_to_cell_3c, cell_to_index_3c, &
1391 : has_mat_P_blocks, do_ri_sos_laplace_mp2, &
1392 146 : dbcsr_time, dbcsr_nflop)
1393 :
1394 : ! the same for open shell, use fm_mo_coeff_occ_beta and fm_mo_coeff_virt_beta
1395 146 : IF (my_open_shell) THEN
1396 30 : CALL zero_mat_P_omega(mat_P_omega(:, :, 2))
1397 : CALL compute_mat_P_omega(mat_P_omega(:, :, 2), fm_scaled_dm_occ_tau, &
1398 : fm_scaled_dm_virt_tau, fm_mo_coeff_occ(2), &
1399 : fm_mo_coeff_virt(2), &
1400 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1401 : mat_P_global, matrix_s, 2, &
1402 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1403 : starts_array_mc, ends_array_mc, &
1404 : starts_array_mc_block, ends_array_mc_block, &
1405 : weights_cos_tf_t_to_w, tj, tau_tj, e_fermi(2), eps_filter, alpha, &
1406 : eps_filter_im_time, Eigenval(:, 1, 2), nmo, &
1407 : num_integ_points, cut_memory, &
1408 : unit_nr, mp2_env, para_env, &
1409 : qs_env, do_kpoints_from_Gamma, &
1410 : index_to_cell_3c, cell_to_index_3c, &
1411 : has_mat_P_blocks, do_ri_sos_laplace_mp2, &
1412 30 : dbcsr_time, dbcsr_nflop)
1413 :
1414 : !For RPA, we sum up the P matrices. If no force needed, can clean-up the beta spin one
1415 30 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1416 144 : DO j = 1, SIZE(mat_P_omega, 2)
1417 928 : DO i = 1, SIZE(mat_P_omega, 1)
1418 784 : CALL dbcsr_add(mat_P_omega(i, j, 1)%matrix, mat_P_omega(i, j, 2)%matrix, 1.0_dp, 1.0_dp)
1419 906 : IF (.NOT. calc_forces) CALL dbcsr_clear(mat_P_omega(i, j, 2)%matrix)
1420 : END DO
1421 : END DO
1422 : END IF
1423 : END IF ! my_open_shell
1424 :
1425 : END IF ! do im time
1426 :
1427 9526 : DO jquad = 1, num_integ_points
1428 :
1429 9142 : IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE
1430 :
1431 8532 : CALL timeset(routineN//"_RPA_matrix_operations", handle3)
1432 :
1433 8532 : IF (do_ri_sos_laplace_mp2) THEN
1434 206 : omega = tau_tj(jquad)
1435 : ELSE
1436 8326 : IF (do_minimax_quad) THEN
1437 928 : omega = tj(jquad)
1438 : ELSE
1439 7398 : omega = a_scaling/TAN(tj(jquad))
1440 : END IF
1441 : END IF ! do_ri_sos_laplace_mp2
1442 :
1443 8532 : IF (do_im_time) THEN
1444 : ! in case we do imag time, we already calculated fm_mat_Q by a Fourier transform from im. time
1445 :
1446 928 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
1447 :
1448 1784 : DO ispin = 1, SIZE(mat_P_omega, 3)
1449 : CALL contract_P_omega_with_mat_L(mat_P_omega(jquad, 1, ispin)%matrix, mat_L%matrix, mat_work, &
1450 : eps_filter_im_time, fm_mat_work, dimen_RI, dimen_RI_red, &
1451 1784 : fm_mat_Minv_L_kpoints(1, 1), fm_mat_Q(ispin))
1452 : END DO
1453 : END IF
1454 :
1455 : ELSE
1456 :
1457 7604 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A, 1X, I3, 1X, A, 1X, I3)") &
1458 3802 : "INTEG_INFO| Started with Integration point", jquad, "of", num_integ_points
1459 :
1460 7604 : IF (first_cycle .AND. count_ev_sc_gw > 1) THEN
1461 100 : IF (iter_sc_gw0 == 1) THEN
1462 140 : DO ispin = 1, nspins
1463 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virtual(ispin), &
1464 140 : Eigenval_last(:, 1, ispin), homo(ispin), omega_old)
1465 : END DO
1466 : ELSE
1467 68 : DO ispin = 1, nspins
1468 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virtual(ispin), &
1469 68 : Eigenval_scf(:, 1, ispin), homo(ispin), omega_old)
1470 : END DO
1471 : END IF
1472 : END IF
1473 :
1474 7604 : IF (iter_sc_GW0 > 1) THEN
1475 6140 : DO ispin = 1, nspins
1476 : CALL calc_mat_Q(fm_mat_S(ispin), do_ri_sos_laplace_mp2, first_cycle, virtual(ispin), &
1477 : Eigenval_scf(:, 1, ispin), homo(ispin), omega, omega_old, jquad, mm_style, &
1478 : dimen_RI_red, dimen_ia(ispin), alpha, fm_mat_Q(ispin), &
1479 : fm_mat_Q_gemm(ispin), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
1480 6140 : num_integ_points, count_ev_sc_GW)
1481 : END DO
1482 :
1483 : ! For SOS-MP2 we need both matrices separately
1484 3070 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1485 3070 : DO ispin = 2, nspins
1486 3070 : 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))
1487 : END DO
1488 : END IF
1489 : ELSE
1490 9220 : DO ispin = 1, nspins
1491 : CALL calc_mat_Q(fm_mat_S(ispin), do_ri_sos_laplace_mp2, first_cycle, virtual(ispin), &
1492 : Eigenval(:, 1, ispin), homo(ispin), omega, omega_old, jquad, mm_style, &
1493 : dimen_RI_red, dimen_ia(ispin), alpha, fm_mat_Q(ispin), &
1494 : fm_mat_Q_gemm(ispin), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
1495 9220 : num_integ_points, count_ev_sc_GW)
1496 : END DO
1497 :
1498 : ! For SOS-MP2 we need both matrices separately
1499 4534 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1500 4558 : DO ispin = 2, nspins
1501 4558 : 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))
1502 : END DO
1503 : END IF
1504 :
1505 : END IF
1506 :
1507 : END IF ! im time
1508 :
1509 : ! Calculate RPA exchange energy correction
1510 8532 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
1511 12 : e_exchange_corr = 0.0_dp
1512 12 : CALL exchange_work%compute(fm_mat_Q(1), Eigenval(:, 1, :), fm_mat_S, omega, e_exchange_corr, mp2_env)
1513 :
1514 : ! Evaluate the final exchange energy correction
1515 12 : e_exchange = e_exchange + e_exchange_corr*wj(jquad)
1516 : END IF
1517 :
1518 8532 : IF (do_ri_sos_laplace_mp2) THEN
1519 206 : CALL SOS_MP2_postprocessing(fm_mat_Q, Erpa, tau_wj(jquad))
1520 :
1521 206 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
1522 : fm_mat_Q, fm_mat_Q_gemm, dgemm_counter, fm_mat_S, omega, homo, virtual, &
1523 60 : Eigenval(:, 1, :), tau_wj(jquad), unit_nr)
1524 : ELSE
1525 8326 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_copy_Q(fm_mat_Q(1), rpa_grad)
1526 :
1527 8326 : CALL Q_trace_and_add_unit_matrix(dimen_RI_red, trace_Qomega, fm_mat_Q(1))
1528 :
1529 8326 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
1530 : CALL invert_eps_compute_W_and_Erpa_kp(dimen_RI, num_integ_points, jquad, nkp, count_ev_sc_GW, para_env, &
1531 : Erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, &
1532 : wkp_W, do_gw_im_time, do_ri_Sigma_x, do_kpoints_from_Gamma, &
1533 : cfm_mat_Q, ikp_local, &
1534 : mat_P_omega(:, :, 1), mat_P_omega_kp, qs_env, eps_filter_im_time, unit_nr, &
1535 : kpoints, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
1536 : fm_mat_W, fm_mat_RI_global_work, mat_MinvVMinv, &
1537 132 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv)
1538 : ELSE
1539 8194 : CALL compute_Erpa_by_freq_int(dimen_RI_red, trace_Qomega, fm_mat_Q(1), para_env_RPA, Erpa, wj(jquad))
1540 : END IF
1541 :
1542 8326 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
1543 : fm_mat_Q, fm_mat_Q_gemm, dgemm_counter, fm_mat_S, omega, homo, virtual, &
1544 48 : Eigenval(:, 1, :), wj(jquad), unit_nr)
1545 : END IF ! do_ri_sos_laplace_mp2
1546 :
1547 : ! save omega and reset the first_cycle flag
1548 8532 : first_cycle = .FALSE.
1549 8532 : omega_old = omega
1550 :
1551 8532 : CALL timestop(handle3)
1552 :
1553 8532 : IF (my_do_gw) THEN
1554 :
1555 7788 : CALL get_fermi_level_offset(fermi_level_offset, fermi_level_offset_input, Eigenval(:, 1, :), homo)
1556 :
1557 : ! do_im_time = TRUE means low-scaling calculation
1558 7788 : IF (do_im_time) THEN
1559 : ! only for molecules
1560 578 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
1561 : 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, &
1562 470 : tj, tau_tj, weights_cos_tf_w_to_t, jquad, omega)
1563 : END IF
1564 : ELSE
1565 : CALL compute_GW_self_energy(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI_red, gw_corr_lev_occ, &
1566 : gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, &
1567 : do_im_time, do_periodic, first_cycle_periodic_correction, &
1568 : fermi_level_offset, &
1569 : omega, Eigenval(:, 1, :), delta_corr, vec_omega_fit_gw, vec_W_gw, wj, &
1570 : fm_mat_Q(1), fm_mat_R_gw, fm_mat_S_gw, &
1571 : fm_mat_S_gw_work, mo_coeff(1), para_env, &
1572 : para_env_RPA, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, &
1573 7210 : kpoints, qs_env, mp2_env)
1574 : END IF
1575 : END IF
1576 :
1577 8532 : IF (unit_nr > 0) CALL m_flush(unit_nr)
1578 18058 : CALL para_env%sync() ! sync to see output
1579 :
1580 : END DO ! jquad
1581 :
1582 384 : CALL para_env%sum(Erpa)
1583 :
1584 384 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1585 326 : Erpa = Erpa/(pi*2.0_dp)
1586 326 : IF (do_minimax_quad) Erpa = Erpa/2.0_dp
1587 : END IF
1588 :
1589 384 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
1590 12 : CALL para_env%sum(E_exchange)
1591 12 : E_exchange = E_exchange/(pi*2.0_dp)
1592 12 : IF (do_minimax_quad) E_exchange = E_exchange/2.0_dp
1593 12 : mp2_env%ri_rpa%ener_exchange = E_exchange
1594 : END IF
1595 :
1596 384 : IF (calc_forces .AND. do_ri_sos_laplace_mp2 .AND. do_im_time) THEN
1597 22 : IF (my_open_shell) THEN
1598 4 : Pspin = 1
1599 4 : Qspin = 2
1600 : CALL calc_laplace_loop_forces(force_data, mat_P_omega(:, 1, :), t_3c_M, t_3c_O(1, 1), &
1601 : t_3c_O_compressed(1, 1, :), t_3c_O_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1602 : fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1603 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1604 : starts_array_mc, ends_array_mc, starts_array_mc_block, &
1605 : ends_array_mc_block, num_integ_points, nmo, Eigenval(:, 1, :), &
1606 : tau_tj, tau_wj, cut_memory, Pspin, Qspin, my_open_shell, &
1607 4 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1608 4 : Pspin = 2
1609 4 : Qspin = 1
1610 : CALL calc_laplace_loop_forces(force_data, mat_P_omega(:, 1, :), t_3c_M, t_3c_O(1, 1), &
1611 : t_3c_O_compressed(1, 1, :), t_3c_O_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1612 : fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1613 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1614 : starts_array_mc, ends_array_mc, starts_array_mc_block, &
1615 : ends_array_mc_block, num_integ_points, nmo, Eigenval(:, 1, :), &
1616 : tau_tj, tau_wj, cut_memory, Pspin, Qspin, my_open_shell, &
1617 4 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1618 :
1619 : ELSE
1620 18 : Pspin = 1
1621 18 : Qspin = 1
1622 : CALL calc_laplace_loop_forces(force_data, mat_P_omega(:, 1, :), t_3c_M, t_3c_O(1, 1), &
1623 : t_3c_O_compressed(1, 1, :), t_3c_O_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1624 : fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1625 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1626 : starts_array_mc, ends_array_mc, starts_array_mc_block, &
1627 : ends_array_mc_block, num_integ_points, nmo, Eigenval(:, 1, :), &
1628 : tau_tj, tau_wj, cut_memory, Pspin, Qspin, my_open_shell, &
1629 18 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1630 : END IF
1631 22 : CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1632 : END IF !laplace SOS-MP2
1633 :
1634 384 : IF (calc_forces .AND. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) THEN
1635 64 : DO ispin = 1, nspins
1636 : CALL calc_rpa_loop_forces(force_data, mat_P_omega(:, 1, :), t_3c_M, t_3c_O(1, 1), &
1637 : t_3c_O_compressed(1, 1, :), t_3c_O_ind(1, 1, :), fm_scaled_dm_occ_tau, &
1638 : fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
1639 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
1640 : starts_array_mc, ends_array_mc, starts_array_mc_block, &
1641 : ends_array_mc_block, num_integ_points, nmo, Eigenval(:, 1, :), &
1642 : e_fermi(ispin), weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, tj, &
1643 : wj, tau_tj, cut_memory, ispin, my_open_shell, unit_nr, dbcsr_time, &
1644 64 : dbcsr_nflop, mp2_env, qs_env)
1645 : END DO
1646 28 : CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1647 : END IF
1648 :
1649 384 : IF (do_im_time) THEN
1650 :
1651 146 : my_flop_rate = REAL(dbcsr_nflop, dp)/(1.0E09_dp*dbcsr_time)
1652 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(/T3,A,T73,ES8.2)") &
1653 73 : "PERFORMANCE| DBCSR total number of flops:", REAL(dbcsr_nflop*para_env%num_pe, dp)
1654 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.2)") &
1655 73 : "PERFORMANCE| DBCSR total execution time:", dbcsr_time
1656 146 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.2)") &
1657 73 : "PERFORMANCE| DBCSR flop rate (Gflops / MPI rank):", my_flop_rate
1658 :
1659 : ELSE
1660 :
1661 238 : CALL dgemm_counter_write(dgemm_counter, para_env)
1662 :
1663 : END IF
1664 :
1665 : ! GW: for low-scaling calculation: Compute self-energy Sigma(i*tau), Sigma(i*omega)
1666 : ! for low-scaling and ordinary-scaling: analytic continuation from Sigma(iw) -> Sigma(w)
1667 : ! and correction of quasiparticle energies e_n^GW
1668 656 : IF (my_do_gw) THEN
1669 :
1670 : CALL compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, &
1671 : gw_corr_lev_tot, gw_corr_lev_virt, homo, &
1672 : nmo, num_fit_points, num_integ_points, &
1673 : unit_nr, do_apply_ic_corr_to_gw, do_im_time, &
1674 : do_periodic, do_ri_Sigma_x, first_cycle_periodic_correction, &
1675 : e_fermi, eps_filter, fermi_level_offset, &
1676 : delta_corr, Eigenval, &
1677 : Eigenval_last, Eigenval_scf, iter_sc_GW0, exit_ev_gw, tau_tj, tj, &
1678 : vec_omega_fit_gw, vec_Sigma_x_gw, mp2_env%ri_g0w0%ic_corr_list, &
1679 : weights_cos_tf_t_to_w, weights_sin_tf_t_to_w, &
1680 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_mo_coeff_occ, &
1681 : fm_mo_coeff_virt, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
1682 : mo_coeff(1), fm_mat_W, para_env, para_env_RPA, mat_dm, mat_MinvVMinv, &
1683 : t_3c_O, t_3c_M, t_3c_overl_int_ao_mo, t_3c_O_compressed, t_3c_O_mo_compressed, &
1684 : t_3c_O_ind, t_3c_O_mo_ind, &
1685 : t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
1686 : matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, mat_W, matrix_s, &
1687 : kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_RPA, &
1688 204 : starts_array_mc, ends_array_mc)
1689 :
1690 : ! if HOMO-LUMO gap differs by less than mp2_env%ri_g0w0%eps_ev_sc_iter, exit ev sc GW loop
1691 204 : IF (exit_ev_gw) EXIT
1692 :
1693 : END IF ! my_do_gw if
1694 :
1695 : END DO ! evGW loop
1696 :
1697 272 : IF (do_ic_model) THEN
1698 :
1699 2 : IF (my_open_shell) THEN
1700 :
1701 : CALL calculate_ic_correction(Eigenval(:, 1, 1), mat_MinvVMinv%matrix, &
1702 : t_3c_overl_nnP_ic(1), t_3c_overl_nnP_ic_reflected(1), &
1703 : gw_corr_lev_tot, &
1704 : gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), unit_nr, &
1705 0 : print_ic_values, para_env, do_alpha=.TRUE.)
1706 :
1707 : CALL calculate_ic_correction(Eigenval(:, 1, 2), mat_MinvVMinv%matrix, &
1708 : t_3c_overl_nnP_ic(2), t_3c_overl_nnP_ic_reflected(2), &
1709 : gw_corr_lev_tot, &
1710 : gw_corr_lev_occ(2), gw_corr_lev_virt(2), homo(2), unit_nr, &
1711 0 : print_ic_values, para_env, do_beta=.TRUE.)
1712 :
1713 : ELSE
1714 :
1715 : CALL calculate_ic_correction(Eigenval(:, 1, 1), mat_MinvVMinv%matrix, &
1716 : t_3c_overl_nnP_ic(1), t_3c_overl_nnP_ic_reflected(1), &
1717 : gw_corr_lev_tot, &
1718 : gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), unit_nr, &
1719 2 : print_ic_values, para_env)
1720 :
1721 : END IF
1722 :
1723 : END IF
1724 :
1725 : ! postprocessing after GW for Bethe-Salpeter
1726 272 : IF (do_bse) THEN
1727 : ! Check used GW flavor; in Case of evGW we use W0 for BSE
1728 : ! Use environment variable, since local iter_evGW is overwritten if evGW0 is invoked
1729 18 : IF (mp2_env%ri_g0w0%iter_evGW > 1) THEN
1730 6 : IF (unit_nr > 0) THEN
1731 : CALL cp_warn(__LOCATION__, &
1732 : "BSE@evGW applies W0, i.e. screening with DFT energies! "// &
1733 3 : "Use BSE@evGW0 to avoid this warning!")
1734 : END IF
1735 : END IF
1736 : ! Create a copy of fm_mat_S for usage in BSE
1737 18 : CALL cp_fm_create(fm_mat_S_ia_bse, fm_mat_S(1)%matrix_struct)
1738 18 : CALL cp_fm_to_fm(fm_mat_S(1), fm_mat_S_ia_bse)
1739 : ! Remove energy/frequency factor from 3c-Integral for BSE
1740 18 : IF (iter_sc_gw0 == 1) THEN
1741 : CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse, virtual(1), &
1742 12 : Eigenval_last(:, 1, 1), homo(1), omega)
1743 : ELSE
1744 : CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse, virtual(1), &
1745 6 : Eigenval_scf(:, 1, 1), homo(1), omega)
1746 : END IF
1747 : ! Main routine for all BSE postprocessing
1748 : CALL start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
1749 : fm_mat_Q_static_bse_gemm, &
1750 : Eigenval, Eigenval_scf, &
1751 : homo, virtual, dimen_RI, dimen_RI_red, bse_lev_virt, &
1752 18 : gd_array, color_sub, mp2_env, qs_env, mo_coeff, unit_nr)
1753 : ! Release BSE-copy of fm_mat_S
1754 18 : CALL cp_fm_release(fm_mat_S_ia_bse)
1755 : END IF
1756 :
1757 272 : IF (my_do_gw) THEN
1758 : CALL deallocate_matrices_gw(fm_mat_S_gw_work, vec_W_gw, vec_Sigma_c_gw, vec_omega_fit_gw, &
1759 : mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw, &
1760 : Eigenval_last, Eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
1761 92 : kpoints, vec_Sigma_x_gw,.NOT. do_im_time)
1762 : END IF
1763 :
1764 272 : IF (do_im_time) THEN
1765 :
1766 : CALL dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, &
1767 : fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, index_to_cell_3c, &
1768 : cell_to_index_3c, do_ic_model, &
1769 : do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, &
1770 : has_mat_P_blocks, &
1771 : wkp_W, cfm_mat_Q, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
1772 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, fm_mat_RI_global_work, fm_mat_work, &
1773 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_L, &
1774 : mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
1775 134 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, mat_work, qs_env)
1776 :
1777 134 : IF (my_do_gw) THEN
1778 : CALL deallocate_matrices_gw_im_time(weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, do_ic_model, &
1779 : do_kpoints_cubic_RPA, fm_mat_W, &
1780 : t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
1781 : t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
1782 : t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
1783 46 : mat_W, qs_env)
1784 : END IF
1785 :
1786 : END IF
1787 :
1788 272 : IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) CALL exchange_work%release()
1789 :
1790 272 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1791 214 : DEALLOCATE (tj)
1792 214 : DEALLOCATE (wj)
1793 214 : DEALLOCATE (trace_Qomega)
1794 : END IF
1795 :
1796 272 : IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
1797 162 : DEALLOCATE (tau_tj)
1798 162 : DEALLOCATE (tau_wj)
1799 : END IF
1800 :
1801 272 : IF (do_im_time .AND. calc_forces) THEN
1802 50 : CALL im_time_force_release(force_data)
1803 : END IF
1804 :
1805 272 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_finalize(rpa_grad, mp2_env, para_env_sub, para_env, &
1806 : qs_env, gd_array, color_sub, do_ri_sos_laplace_mp2, &
1807 40 : homo, virtual)
1808 :
1809 272 : CALL timestop(handle)
1810 :
1811 5570 : END SUBROUTINE rpa_num_int
1812 :
1813 : END MODULE rpa_main
|
