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