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