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: bse_solve,&
22 : prepare_bse_env_from_ri_rpa_gw
23 : USE bse_types, ONLY: bse_env_release
24 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
25 : cp_blacs_env_release,&
26 : cp_blacs_env_type
27 : USE cp_cfm_types, ONLY: cp_cfm_type
28 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
29 : dbcsr_clear,&
30 : dbcsr_get_info,&
31 : dbcsr_p_type,&
32 : dbcsr_type
33 : USE cp_fm_basic_linalg, ONLY: cp_fm_add_on_diag,&
34 : cp_fm_scale_and_add
35 : USE cp_fm_redistribution, ONLY: cp_array_to_fm
36 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37 : cp_fm_struct_release,&
38 : cp_fm_struct_type
39 : USE cp_fm_types, ONLY: cp_fm_create,&
40 : cp_fm_get_diag,&
41 : cp_fm_get_info,&
42 : cp_fm_release,&
43 : cp_fm_set_all,&
44 : cp_fm_to_fm,&
45 : cp_fm_type
46 : USE dbt_api, ONLY: dbt_type
47 : USE dgemm_counter_types, ONLY: dgemm_counter_init,&
48 : dgemm_counter_type,&
49 : dgemm_counter_write
50 : USE group_dist_types, ONLY: get_group_dist,&
51 : group_dist_d1_type,&
52 : maxsize,&
53 : release_group_dist
54 : USE hfx_types, ONLY: block_ind_type,&
55 : hfx_compression_type
56 : USE input_constants, ONLY: cohsex_qp_update,&
57 : rpa_exchange_axk,&
58 : rpa_exchange_none,&
59 : rpa_exchange_sosex,&
60 : sigma_none,&
61 : wfc_mm_style_gemm
62 : USE input_section_types, ONLY: section_vals_type,&
63 : section_vals_val_set
64 : USE kinds, ONLY: dp,&
65 : int_8
66 : USE kpoint_types, ONLY: get_kpoint_info,&
67 : kpoint_env_type,&
68 : kpoint_type
69 : USE machine, ONLY: m_flush,&
70 : m_memory
71 : USE mathconstants, ONLY: pi,&
72 : z_zero
73 : USE message_passing, ONLY: mp_comm_type,&
74 : mp_para_env_release,&
75 : mp_para_env_type
76 : USE minimax_exp, ONLY: check_exp_minimax_range
77 : USE mp2_laplace, ONLY: SOS_MP2_postprocessing
78 : USE mp2_types, ONLY: mp2_type,&
79 : three_dim_real_array,&
80 : two_dim_int_array,&
81 : two_dim_real_array
82 : USE qs_environment_types, ONLY: get_qs_env,&
83 : qs_environment_type
84 : USE qs_mo_types, ONLY: get_mo_set,&
85 : mo_set_type
86 : USE rpa_exchange, ONLY: rpa_exchange_needed_mem,&
87 : rpa_exchange_work_type
88 : USE rpa_grad, ONLY: rpa_grad_copy_Q,&
89 : rpa_grad_create,&
90 : rpa_grad_finalize,&
91 : rpa_grad_matrix_operations,&
92 : rpa_grad_needed_mem,&
93 : rpa_grad_type
94 : USE rpa_gw, ONLY: allocate_matrices_gw,&
95 : allocate_matrices_gw_im_time,&
96 : compute_GW_self_energy,&
97 : compute_QP_energies,&
98 : compute_W_cubic_GW,&
99 : deallocate_matrices_gw,&
100 : deallocate_matrices_gw_im_time,&
101 : get_fermi_level_offset
102 : USE rpa_gw_ic, ONLY: calculate_ic_correction
103 : USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_MOs,&
104 : invert_eps_compute_W_and_Erpa_kp
105 : USE rpa_gw_static_cohsex, ONLY: get_static_cohsex_mo_coeff,&
106 : rotate_bse_slab,&
107 : static_cohsex_qp_energies
108 : USE rpa_im_time, ONLY: compute_mat_P_omega
109 : USE rpa_im_time_force_methods, ONLY: calc_laplace_loop_forces,&
110 : calc_post_loop_forces,&
111 : calc_rpa_loop_forces,&
112 : init_im_time_forces,&
113 : keep_initial_quad
114 : USE rpa_im_time_force_types, ONLY: im_time_force_release,&
115 : im_time_force_type
116 : USE rpa_sigma_functional, ONLY: finalize_rpa_sigma,&
117 : rpa_sigma_create,&
118 : rpa_sigma_matrix_spectral,&
119 : rpa_sigma_type
120 : USE rpa_util, ONLY: alloc_im_time,&
121 : calc_mat_Q,&
122 : compute_Erpa_by_freq_int,&
123 : contract_P_omega_with_mat_L,&
124 : dealloc_im_time,&
125 : remove_scaling_factor_rpa
126 : USE time_frequency_grids, ONLY: build_clenshaw_grid,&
127 : build_minimax_time_frequency_grid,&
128 : time_frequency_grid_release,&
129 : time_frequency_grid_type
130 : USE util, ONLY: get_limit
131 : #include "./base/base_uses.f90"
132 :
133 : IMPLICIT NONE
134 :
135 : PRIVATE
136 :
137 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_main'
138 :
139 : PUBLIC :: rpa_ri_compute_en
140 :
141 : CONTAINS
142 :
143 : ! **************************************************************************************************
144 : !> \brief ...
145 : !> \param qs_env ...
146 : !> \param Erpa ...
147 : !> \param mp2_env ...
148 : !> \param BIb_C ...
149 : !> \param BIb_C_gw ...
150 : !> \param BIb_C_bse_ij ...
151 : !> \param BIb_C_bse_ab ...
152 : !> \param para_env ...
153 : !> \param para_env_sub ...
154 : !> \param color_sub ...
155 : !> \param gd_array ...
156 : !> \param gd_B_virtual ...
157 : !> \param gd_B_all ...
158 : !> \param gd_B_occ_bse ...
159 : !> \param gd_B_virt_bse ...
160 : !> \param mo_coeff ...
161 : !> \param fm_matrix_PQ ...
162 : !> \param cfm_matrix_L_kpoints ...
163 : !> \param cfm_matrix_Minv_L_kpoints ...
164 : !> \param fm_matrix_Minv_L_gamma ...
165 : !> \param fm_matrix_Minv ...
166 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
167 : !> \param kpoints ...
168 : !> \param Eigenval ...
169 : !> \param nmo ...
170 : !> \param homo ...
171 : !> \param dimen_RI ...
172 : !> \param dimen_RI_red ...
173 : !> \param gw_corr_lev_occ ...
174 : !> \param gw_corr_lev_virt ...
175 : !> \param bse_lev_virt ...
176 : !> \param unit_nr ...
177 : !> \param do_ri_sos_laplace_mp2 ...
178 : !> \param my_do_gw ...
179 : !> \param do_im_time ...
180 : !> \param do_bse ...
181 : !> \param matrix_s ...
182 : !> \param mat_munu ...
183 : !> \param mat_P_global ...
184 : !> \param t_3c_M ...
185 : !> \param t_3c_O ...
186 : !> \param t_3c_O_compressed ...
187 : !> \param t_3c_O_ind ...
188 : !> \param starts_array_mc ...
189 : !> \param ends_array_mc ...
190 : !> \param starts_array_mc_block ...
191 : !> \param ends_array_mc_block ...
192 : !> \param calc_forces ...
193 : ! **************************************************************************************************
194 342 : SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, &
195 : para_env, para_env_sub, color_sub, &
196 1026 : gd_array, gd_B_virtual, gd_B_all, gd_B_occ_bse, gd_B_virt_bse, &
197 342 : mo_coeff, fm_matrix_PQ, cfm_matrix_L_kpoints, cfm_matrix_Minv_L_kpoints, &
198 : fm_matrix_Minv_L_gamma, &
199 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, kpoints, &
200 684 : Eigenval, nmo, homo, dimen_RI, dimen_RI_red, gw_corr_lev_occ, gw_corr_lev_virt, &
201 342 : bse_lev_virt, &
202 : unit_nr, do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_bse, matrix_s, &
203 : mat_munu, mat_P_global, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
204 : starts_array_mc, ends_array_mc, &
205 : starts_array_mc_block, ends_array_mc_block, calc_forces)
206 :
207 : TYPE(qs_environment_type), POINTER :: qs_env
208 : REAL(KIND=dp), INTENT(OUT) :: Erpa
209 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
210 : TYPE(three_dim_real_array), DIMENSION(:), &
211 : INTENT(INOUT) :: BIb_C, BIb_C_gw, BIb_C_bse_ij, &
212 : BIb_C_bse_ab
213 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
214 : INTEGER, INTENT(INOUT) :: color_sub
215 : TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_array
216 : TYPE(group_dist_d1_type), DIMENSION(:), &
217 : INTENT(INOUT) :: gd_B_virtual
218 : TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_B_all
219 : TYPE(group_dist_d1_type), DIMENSION(:), &
220 : INTENT(INOUT) :: gd_B_occ_bse, gd_B_virt_bse
221 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
222 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_PQ
223 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_matrix_L_kpoints, &
224 : cfm_matrix_Minv_L_kpoints
225 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv_L_gamma, fm_matrix_Minv, &
226 : fm_matrix_Minv_Vtrunc_Minv
227 : TYPE(kpoint_type), POINTER :: kpoints
228 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
229 : INTENT(INOUT) :: Eigenval
230 : INTEGER, INTENT(IN) :: nmo
231 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
232 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
233 : INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, &
234 : bse_lev_virt
235 : INTEGER, INTENT(IN) :: unit_nr
236 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, my_do_gw, &
237 : do_im_time, do_bse
238 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
239 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
240 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_P_global
241 : TYPE(dbt_type) :: t_3c_M
242 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O
243 : TYPE(hfx_compression_type), ALLOCATABLE, &
244 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_compressed
245 : TYPE(block_ind_type), ALLOCATABLE, &
246 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_ind
247 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
248 : starts_array_mc_block, &
249 : ends_array_mc_block
250 : LOGICAL, INTENT(IN) :: calc_forces
251 :
252 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_ri_compute_en'
253 :
254 : INTEGER :: best_integ_group_size, best_num_integ_point, color_rpa_group, dimen_nm_gw, &
255 : dimen_virt_square, handle, handle2, handle3, ierr, iiB, input_num_integ_groups, &
256 : integ_group_size, ispin, jjB, min_integ_group_size, my_group_L_end, my_group_L_size, &
257 : my_group_L_start, my_nm_gw_end, my_nm_gw_size, my_nm_gw_start, ncol_block_mat, ngroup, &
258 : nrow_block_mat, nspins, num_integ_group, num_integ_points, pos_integ_group
259 : INTEGER(KIND=int_8) :: mem
260 684 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dimen_homo_square, dimen_ia, my_ab_comb_bse_end, &
261 342 : my_ab_comb_bse_size, my_ab_comb_bse_start, my_ia_end, my_ia_size, my_ia_start, &
262 342 : my_ij_comb_bse_end, my_ij_comb_bse_size, my_ij_comb_bse_start, virtual
263 : LOGICAL :: do_kpoints_from_Gamma, do_minimax_quad, &
264 : my_open_shell, skip_integ_group_opt
265 : REAL(KIND=dp) :: allowed_memory, avail_mem, E_Range, Emax, Emin, mem_for_iaK, mem_for_QK, &
266 : mem_min, mem_per_group, mem_per_rank, mem_per_repl, mem_real
267 342 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_kp
268 684 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_Q, fm_mat_Q_gemm, fm_mat_S, &
269 342 : fm_mat_S_ab_bse, fm_mat_S_gw, &
270 342 : fm_mat_S_ij_bse
271 684 : TYPE(cp_fm_type), DIMENSION(1) :: fm_mat_R_gw
272 : TYPE(mp_para_env_type), POINTER :: para_env_RPA
273 : TYPE(two_dim_real_array), ALLOCATABLE, &
274 342 : DIMENSION(:) :: BIb_C_2D, BIb_C_2D_bse_ab, &
275 342 : BIb_C_2D_bse_ij, BIb_C_2D_gw
276 :
277 342 : CALL timeset(routineN, handle)
278 :
279 342 : CALL cite_reference(DelBen2013)
280 342 : CALL cite_reference(DelBen2015)
281 :
282 342 : IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_axk) THEN
283 10 : CALL cite_reference(Bates2013)
284 332 : ELSE IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_sosex) THEN
285 2 : CALL cite_reference(Freeman1977)
286 2 : CALL cite_reference(Gruneis2009)
287 : END IF
288 342 : IF (mp2_env%ri_rpa%do_rse) THEN
289 8 : CALL cite_reference(Ren2011)
290 8 : CALL cite_reference(Ren2013)
291 : END IF
292 :
293 342 : IF (my_do_gw) THEN
294 134 : CALL cite_reference(Wilhelm2016a)
295 134 : CALL cite_reference(Wilhelm2017)
296 134 : CALL cite_reference(Wilhelm2018)
297 : END IF
298 :
299 342 : IF (do_im_time) THEN
300 144 : CALL cite_reference(Wilhelm2016b)
301 : END IF
302 :
303 342 : nspins = SIZE(homo)
304 342 : my_open_shell = (nspins == 2)
305 2394 : ALLOCATE (virtual(nspins), dimen_ia(nspins), my_ia_end(nspins), my_ia_start(nspins), my_ia_size(nspins))
306 758 : virtual(:) = nmo - homo(:)
307 758 : dimen_ia(:) = virtual(:)*homo(:)
308 :
309 1368 : ALLOCATE (Eigenval_kp(nmo, 1, nspins))
310 10988 : Eigenval_kp(:, 1, :) = Eigenval(:, :)
311 :
312 342 : IF (do_im_time) mp2_env%ri_rpa%minimax_quad = .TRUE.
313 342 : do_minimax_quad = mp2_env%ri_rpa%minimax_quad
314 :
315 342 : IF (do_ri_sos_laplace_mp2) THEN
316 64 : num_integ_points = mp2_env%ri_laplace%n_quadrature
317 64 : input_num_integ_groups = mp2_env%ri_laplace%num_integ_groups
318 :
319 : ! check the range for the minimax approximation
320 64 : E_Range = mp2_env%e_range
321 64 : IF (mp2_env%e_range <= 1.0_dp .OR. mp2_env%e_gap <= 0.0_dp) THEN
322 : Emin = HUGE(dp)
323 : Emax = 0.0_dp
324 104 : DO ispin = 1, nspins
325 104 : IF (homo(ispin) > 0) THEN
326 60 : Emin = MIN(Emin, 2.0_dp*(Eigenval(homo(ispin) + 1, ispin) - Eigenval(homo(ispin), ispin)))
327 3200 : Emax = MAX(Emax, 2.0_dp*(MAXVAL(Eigenval(:, ispin)) - MINVAL(Eigenval(:, ispin))))
328 : END IF
329 : END DO
330 44 : E_Range = Emax/Emin
331 : END IF
332 64 : IF (E_Range < 2.0_dp) E_Range = 2.0_dp
333 : ierr = 0
334 64 : CALL check_exp_minimax_range(num_integ_points, E_Range, ierr)
335 64 : IF (ierr /= 0) THEN
336 : jjB = num_integ_points - 1
337 0 : DO iiB = 1, jjB
338 0 : num_integ_points = num_integ_points - 1
339 : ierr = 0
340 0 : CALL check_exp_minimax_range(num_integ_points, E_Range, ierr)
341 0 : IF (ierr == 0) EXIT
342 : END DO
343 : END IF
344 64 : CPASSERT(num_integ_points >= 1)
345 : ELSE
346 278 : num_integ_points = mp2_env%ri_rpa%rpa_num_quad_points
347 278 : input_num_integ_groups = mp2_env%ri_rpa%rpa_num_integ_groups
348 278 : IF (my_do_gw .AND. do_minimax_quad) THEN
349 46 : IF (num_integ_points > 34) THEN
350 0 : IF (unit_nr > 0) THEN
351 : CALL cp_warn(__LOCATION__, &
352 : "The required number of quadrature point exceeds the maximum possible in the "// &
353 0 : "Minimax quadrature scheme. The number of quadrature point has been reset to 30.")
354 : END IF
355 0 : num_integ_points = 30
356 : END IF
357 : ELSE
358 232 : IF (do_minimax_quad .AND. num_integ_points > 20) THEN
359 0 : IF (unit_nr > 0) THEN
360 : CALL cp_warn(__LOCATION__, &
361 : "The required number of quadrature point exceeds the maximum possible in the "// &
362 0 : "Minimax quadrature scheme. The number of quadrature point has been reset to 20.")
363 : END IF
364 0 : num_integ_points = 20
365 : END IF
366 : END IF
367 : END IF
368 342 : allowed_memory = mp2_env%mp2_memory
369 :
370 342 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
371 :
372 342 : ngroup = para_env%num_pe/para_env_sub%num_pe
373 :
374 : ! for imaginary time or periodic GW or BSE, we use all processors for a single frequency/time point
375 342 : IF (do_im_time .OR. mp2_env%ri_g0w0%do_periodic .OR. do_bse) THEN
376 :
377 202 : integ_group_size = ngroup
378 202 : best_num_integ_point = num_integ_points
379 :
380 : ELSE
381 :
382 : ! Calculate available memory and create integral group according to that
383 : ! mem_for_iaK is the memory needed for storing the 3 centre integrals
384 310 : mem_for_iaK = REAL(SUM(dimen_ia), KIND=dp)*dimen_RI_red*8.0_dp/(1024_dp**2)
385 140 : mem_for_QK = REAL(dimen_RI_red, KIND=dp)*nspins*dimen_RI_red*8.0_dp/(1024_dp**2)
386 :
387 140 : CALL m_memory(mem)
388 140 : mem_real = (mem + 1024*1024 - 1)/(1024*1024)
389 140 : CALL para_env%min(mem_real)
390 :
391 140 : mem_per_rank = 0.0_dp
392 :
393 : ! B_ia_P
394 : mem_per_repl = mem_for_iaK
395 : ! Q (regular and for dgemm)
396 140 : mem_per_repl = mem_per_repl + 2.0_dp*mem_for_QK
397 :
398 140 : IF (calc_forces) CALL rpa_grad_needed_mem(homo, virtual, dimen_RI_red, mem_per_rank, mem_per_repl, do_ri_sos_laplace_mp2)
399 140 : CALL rpa_exchange_needed_mem(mp2_env, homo, virtual, dimen_RI_red, para_env, mem_per_rank, mem_per_repl)
400 :
401 140 : mem_min = mem_per_repl/para_env%num_pe + mem_per_rank
402 :
403 140 : IF (unit_nr > 0) THEN
404 70 : WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum required memory per MPI process:', mem_min, ' MiB'
405 70 : WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Available memory per MPI process:', mem_real, ' MiB'
406 : END IF
407 :
408 : ! Use only the allowed amount of memory
409 140 : mem_real = MIN(mem_real, allowed_memory)
410 : ! For the memory estimate, we require the amount of required memory per replication group and the available memory
411 140 : mem_real = mem_real - mem_per_rank
412 :
413 140 : mem_per_group = mem_real*para_env_sub%num_pe
414 :
415 : ! here we try to find the best rpa/laplace group size
416 140 : skip_integ_group_opt = .FALSE.
417 :
418 : ! Check the input number of integration groups
419 140 : IF (input_num_integ_groups > 0) THEN
420 2 : IF (num_integ_points < input_num_integ_groups) THEN
421 0 : IF (MOD(ngroup, input_num_integ_groups) == 0) THEN
422 0 : best_integ_group_size = ngroup/input_num_integ_groups
423 0 : best_num_integ_point = (num_integ_points + input_num_integ_groups - 1)/input_num_integ_groups
424 : skip_integ_group_opt = .TRUE.
425 : ELSE
426 0 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Total number of groups not multiple of NUM_INTEG_GROUPS'
427 : END IF
428 : ELSE
429 2 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Too many integration groups for the given number of quadrature points'
430 : END IF
431 : END IF
432 :
433 : IF (.NOT. skip_integ_group_opt) THEN
434 140 : best_integ_group_size = ngroup
435 140 : best_num_integ_point = num_integ_points
436 :
437 140 : min_integ_group_size = MAX(1, ngroup/num_integ_points)
438 :
439 140 : integ_group_size = min_integ_group_size - 1
440 140 : DO iiB = min_integ_group_size + 1, ngroup
441 118 : integ_group_size = integ_group_size + 1
442 :
443 : ! check that the ngroup is a multiple of integ_group_size
444 118 : IF (MOD(ngroup, integ_group_size) /= 0) CYCLE
445 :
446 : ! check for memory
447 118 : avail_mem = integ_group_size*mem_per_group
448 118 : IF (avail_mem < mem_per_repl) CYCLE
449 :
450 : ! check that the integration groups have the same size
451 118 : num_integ_group = ngroup/integ_group_size
452 :
453 118 : best_num_integ_point = (num_integ_points + num_integ_group - 1)/num_integ_group
454 118 : best_integ_group_size = integ_group_size
455 :
456 140 : EXIT
457 :
458 : END DO
459 : END IF
460 :
461 : integ_group_size = best_integ_group_size
462 :
463 : END IF
464 :
465 342 : IF (unit_nr > 0 .AND. .NOT. do_im_time) THEN
466 99 : IF (do_ri_sos_laplace_mp2) THEN
467 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
468 14 : "RI_INFO| Group size for laplace numerical integration:", integ_group_size*para_env_sub%num_pe
469 : WRITE (UNIT=unit_nr, FMT="(T3,A)") &
470 14 : "INTEG_INFO| MINIMAX approximation"
471 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
472 14 : "INTEG_INFO| Number of integration points:", num_integ_points
473 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
474 14 : "INTEG_INFO| Max. number of integration points per Laplace group:", best_num_integ_point
475 : ELSE
476 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
477 85 : "RI_INFO| Group size for frequency integration:", integ_group_size*para_env_sub%num_pe
478 85 : IF (do_minimax_quad) THEN
479 : WRITE (UNIT=unit_nr, FMT="(T3,A)") &
480 21 : "INTEG_INFO| MINIMAX quadrature"
481 : ELSE
482 : WRITE (UNIT=unit_nr, FMT="(T3,A)") &
483 64 : "INTEG_INFO| Clenshaw-Curtius quadrature"
484 : END IF
485 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
486 85 : "INTEG_INFO| Number of integration points:", num_integ_points
487 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
488 85 : "INTEG_INFO| Max. number of integration points per RPA group:", best_num_integ_point
489 : END IF
490 99 : CALL m_flush(unit_nr)
491 : END IF
492 :
493 342 : num_integ_group = ngroup/integ_group_size
494 :
495 342 : pos_integ_group = MOD(color_sub, integ_group_size)
496 342 : color_rpa_group = color_sub/integ_group_size
497 :
498 342 : CALL timeset(routineN//"_reorder", handle2)
499 :
500 : ! not necessary for imaginary time
501 :
502 1442 : ALLOCATE (BIb_C_2D(nspins))
503 :
504 342 : IF (.NOT. do_im_time) THEN
505 :
506 : ! reorder the local data in such a way to help the next stage of matrix creation
507 : ! now the data inside the group are divided into a ia x K matrix
508 438 : DO ispin = 1, nspins
509 : CALL calculate_BIb_C_2D(BIb_C_2D(ispin)%array, BIb_C(ispin)%array, para_env_sub, dimen_ia(ispin), &
510 : homo(ispin), virtual(ispin), gd_B_virtual(ispin), &
511 240 : my_ia_size(ispin), my_ia_start(ispin), my_ia_end(ispin), my_group_L_size)
512 :
513 240 : DEALLOCATE (BIb_C(ispin)%array)
514 438 : CALL release_group_dist(gd_B_virtual(ispin))
515 :
516 : END DO
517 :
518 : ! in the GW case, BIb_C_2D_gw is an nm x K matrix, with n: number of corr GW levels, m=nmo
519 198 : IF (my_do_gw) THEN
520 280 : ALLOCATE (BIb_C_2D_gw(nspins))
521 :
522 88 : CALL timeset(routineN//"_reorder_gw", handle3)
523 :
524 88 : dimen_nm_gw = nmo*(gw_corr_lev_occ(1) + gw_corr_lev_virt(1))
525 :
526 : ! The same for open shell
527 192 : DO ispin = 1, nspins
528 : CALL calculate_BIb_C_2D(BIb_C_2D_gw(ispin)%array, BIb_C_gw(ispin)%array, para_env_sub, dimen_nm_gw, &
529 : gw_corr_lev_occ(ispin) + gw_corr_lev_virt(ispin), nmo, gd_B_all, &
530 104 : my_nm_gw_size, my_nm_gw_start, my_nm_gw_end, my_group_L_size)
531 192 : DEALLOCATE (BIb_C_gw(ispin)%array)
532 : END DO
533 :
534 88 : CALL release_group_dist(gd_B_all)
535 :
536 176 : CALL timestop(handle3)
537 :
538 : END IF
539 : END IF
540 :
541 342 : IF (do_bse) THEN
542 :
543 56 : CALL timeset(routineN//"_reorder_bse1", handle3)
544 :
545 304 : ALLOCATE (BIb_C_2D_bse_ij(nspins), BIb_C_2D_bse_ab(nspins))
546 112 : ALLOCATE (dimen_homo_square(nspins))
547 224 : ALLOCATE (my_ij_comb_bse_size(nspins), my_ij_comb_bse_start(nspins), my_ij_comb_bse_end(nspins))
548 224 : ALLOCATE (my_ab_comb_bse_size(nspins), my_ab_comb_bse_start(nspins), my_ab_comb_bse_end(nspins))
549 :
550 : ! We do not implement an explicit bse_lev_occ different to homo here, because the small number of occupied levels
551 : ! does not critically influence the memory
552 124 : DO ispin = 1, nspins
553 68 : dimen_homo_square(ispin) = homo(ispin)**2
554 : CALL calculate_BIb_C_2D(BIb_C_2D_bse_ij(ispin)%array, BIb_C_bse_ij(ispin)%array, para_env_sub, &
555 : dimen_homo_square(ispin), homo(ispin), homo(ispin), gd_B_occ_bse(ispin), &
556 : my_ij_comb_bse_size(ispin), my_ij_comb_bse_start(ispin), &
557 68 : my_ij_comb_bse_end(ispin), my_group_L_size)
558 68 : DEALLOCATE (BIb_C_bse_ij(ispin)%array)
559 124 : CALL release_group_dist(gd_B_occ_bse(ispin))
560 : END DO
561 :
562 56 : CALL timestop(handle3)
563 :
564 56 : CALL timeset(routineN//"_reorder_bse2", handle3)
565 :
566 : ! bse_lev_virt(ispin) (hence dimen_virt_square) and gd_B_virt_bse(ispin) are per-spin
567 124 : DO ispin = 1, nspins
568 68 : dimen_virt_square = bse_lev_virt(ispin)**2
569 : CALL calculate_BIb_C_2D(BIb_C_2D_bse_ab(ispin)%array, BIb_C_bse_ab(ispin)%array, para_env_sub, &
570 : dimen_virt_square, bse_lev_virt(ispin), bse_lev_virt(ispin), gd_B_virt_bse(ispin), &
571 : my_ab_comb_bse_size(ispin), my_ab_comb_bse_start(ispin), &
572 68 : my_ab_comb_bse_end(ispin), my_group_L_size)
573 68 : DEALLOCATE (BIb_C_bse_ab(ispin)%array)
574 124 : CALL release_group_dist(gd_B_virt_bse(ispin))
575 : END DO
576 :
577 168 : CALL timestop(handle3)
578 :
579 : END IF
580 :
581 342 : CALL timestop(handle2)
582 :
583 342 : IF (num_integ_group > 1) THEN
584 118 : ALLOCATE (para_env_RPA)
585 118 : CALL para_env_RPA%from_split(para_env, color_rpa_group)
586 : ELSE
587 224 : para_env_RPA => para_env
588 : END IF
589 :
590 : ! now create the matrices needed for the calculation, Q, S and G
591 : ! Q and G will have omega dependence
592 :
593 342 : IF (do_im_time) THEN
594 896 : ALLOCATE (fm_mat_Q(nspins), fm_mat_Q_gemm(1), fm_mat_S(1))
595 : ELSE
596 1710 : ALLOCATE (fm_mat_Q(nspins), fm_mat_Q_gemm(nspins), fm_mat_S(nspins))
597 : END IF
598 :
599 : CALL create_integ_mat(BIb_C_2D, para_env, dimen_RI_red, dimen_ia, mp2_env%block_size_row, &
600 : mp2_env%block_size_col, unit_nr, my_ia_start, my_ia_end, my_group_L_start, &
601 : my_group_L_end, para_env_RPA, fm_mat_S, nrow_block_mat, ncol_block_mat, &
602 : dimen_ia_for_block_size=dimen_ia(1), do_im_time=do_im_time, &
603 342 : fm_mat_Q_gemm=fm_mat_Q_gemm, fm_mat_Q=fm_mat_Q, qs_env=qs_env)
604 :
605 758 : DEALLOCATE (BIb_C_2D, my_ia_end, my_ia_size, my_ia_start)
606 :
607 : ! for GW, we need other matrix fm_mat_S, always allocate the container to prevent crying compilers
608 1442 : ALLOCATE (fm_mat_S_gw(nspins))
609 342 : IF (my_do_gw .AND. .NOT. do_im_time) THEN
610 :
611 : CALL create_integ_mat(BIb_C_2D_gw, para_env, dimen_RI_red, [dimen_nm_gw, dimen_nm_gw], &
612 : mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
613 : [my_nm_gw_start, my_nm_gw_start], [my_nm_gw_end, my_nm_gw_end], &
614 : my_group_L_start, my_group_L_end, para_env_RPA, fm_mat_S_gw, nrow_block_mat, &
615 : ncol_block_mat, fm_mat_Q(1)%matrix_struct%context, &
616 616 : fm_mat_Q(1)%matrix_struct%context, fm_mat_Q=fm_mat_R_gw)
617 192 : DEALLOCATE (BIb_C_2D_gw)
618 :
619 : END IF
620 :
621 : ! for Bethe-Salpeter, we need other matrix fm_mat_S (per spin; the ab slab dimension is spin-independent)
622 342 : IF (do_bse) THEN
623 304 : ALLOCATE (fm_mat_S_ij_bse(nspins), fm_mat_S_ab_bse(nspins))
624 : CALL create_integ_mat(BIb_C_2D_bse_ij, para_env, dimen_RI_red, dimen_homo_square, &
625 : mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, my_ij_comb_bse_start, &
626 : my_ij_comb_bse_end, my_group_L_start, my_group_L_end, para_env_RPA, &
627 : fm_mat_S_ij_bse, nrow_block_mat, ncol_block_mat, &
628 56 : fm_mat_Q(1)%matrix_struct%context, fm_mat_Q(1)%matrix_struct%context)
629 :
630 : CALL create_integ_mat(BIb_C_2D_bse_ab, para_env, dimen_RI_red, &
631 : [(bse_lev_virt(ispin)**2, ispin=1, nspins)], mp2_env%block_size_row, &
632 : mp2_env%block_size_col, unit_nr, my_ab_comb_bse_start, my_ab_comb_bse_end, &
633 : my_group_L_start, my_group_L_end, para_env_RPA, fm_mat_S_ab_bse, &
634 : nrow_block_mat, ncol_block_mat, fm_mat_Q(1)%matrix_struct%context, &
635 248 : fm_mat_Q(1)%matrix_struct%context)
636 :
637 : END IF
638 :
639 342 : do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
640 342 : IF (do_kpoints_from_Gamma) THEN
641 16 : CALL get_bandstruc_and_k_dependent_MOs(qs_env, Eigenval_kp)
642 : END IF
643 :
644 : ! Now start the RPA calculation
645 : ! cfm_mo_coeff will be deallocated here
646 : CALL rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_sub, unit_nr, &
647 : homo, virtual, dimen_RI, dimen_RI_red, dimen_ia, dimen_nm_gw, &
648 : Eigenval_kp, num_integ_points, num_integ_group, color_rpa_group, &
649 : fm_matrix_PQ, fm_mat_S, fm_mat_Q_gemm, fm_mat_Q, fm_mat_S_gw, fm_mat_R_gw(1), &
650 : fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
651 : my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
652 : bse_lev_virt, &
653 : do_minimax_quad, &
654 : do_im_time, mo_coeff, &
655 : cfm_matrix_L_kpoints, cfm_matrix_Minv_L_kpoints, &
656 : fm_matrix_Minv_L_gamma, &
657 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, mat_munu, mat_P_global, &
658 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
659 : starts_array_mc, ends_array_mc, &
660 : starts_array_mc_block, ends_array_mc_block, &
661 : matrix_s, do_kpoints_from_Gamma, kpoints, gd_array, color_sub, &
662 342 : do_ri_sos_laplace_mp2=do_ri_sos_laplace_mp2, calc_forces=calc_forces)
663 :
664 342 : CALL release_group_dist(gd_array)
665 :
666 342 : IF (num_integ_group > 1) CALL mp_para_env_release(para_env_RPA)
667 :
668 342 : IF (.NOT. do_im_time) THEN
669 198 : CALL cp_fm_release(fm_mat_Q_gemm)
670 198 : CALL cp_fm_release(fm_mat_S)
671 : END IF
672 342 : CALL cp_fm_release(fm_mat_Q)
673 :
674 342 : IF (my_do_gw .AND. .NOT. do_im_time) THEN
675 88 : CALL cp_fm_release(fm_mat_S_gw)
676 88 : CALL cp_fm_release(fm_mat_R_gw(1))
677 : END IF
678 :
679 342 : IF (do_bse) THEN
680 124 : DO ispin = 1, nspins
681 68 : CALL cp_fm_release(fm_mat_S_ij_bse(ispin))
682 124 : CALL cp_fm_release(fm_mat_S_ab_bse(ispin))
683 : END DO
684 56 : DEALLOCATE (fm_mat_S_ij_bse, fm_mat_S_ab_bse)
685 : END IF
686 :
687 342 : CALL timestop(handle)
688 :
689 1504 : END SUBROUTINE rpa_ri_compute_en
690 :
691 : ! **************************************************************************************************
692 : !> \brief reorder the local data in such a way to help the next stage of matrix creation;
693 : !> now the data inside the group are divided into a ia x K matrix (BIb_C_2D);
694 : !> Subroutine created to avoid massive double coding
695 : !> \param BIb_C_2D ...
696 : !> \param BIb_C ...
697 : !> \param para_env_sub ...
698 : !> \param dimen_ia ...
699 : !> \param homo ...
700 : !> \param virtual ...
701 : !> \param gd_B_virtual ...
702 : !> \param my_ia_size ...
703 : !> \param my_ia_start ...
704 : !> \param my_ia_end ...
705 : !> \param my_group_L_size ...
706 : !> \author Jan Wilhelm, 03/2015
707 : ! **************************************************************************************************
708 480 : SUBROUTINE calculate_BIb_C_2D(BIb_C_2D, BIb_C, para_env_sub, dimen_ia, homo, virtual, &
709 : gd_B_virtual, &
710 : my_ia_size, my_ia_start, my_ia_end, my_group_L_size)
711 :
712 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
713 : INTENT(OUT) :: BIb_C_2D
714 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
715 : INTENT(IN) :: BIb_C
716 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
717 : INTEGER, INTENT(IN) :: dimen_ia, homo, virtual
718 : TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_B_virtual
719 : INTEGER :: my_ia_size, my_ia_start, my_ia_end, &
720 : my_group_L_size
721 :
722 : INTEGER, PARAMETER :: occ_chunk = 128
723 :
724 : INTEGER :: ia_global, iiB, itmp(2), jjB, my_B_size, my_B_virtual_start, occ_high, occ_low, &
725 : proc_receive, proc_send, proc_shift, rec_B_size, rec_B_virtual_end, rec_B_virtual_start
726 480 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: BIb_C_rec_1D
727 480 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: BIb_C_rec
728 :
729 480 : itmp = get_limit(dimen_ia, para_env_sub%num_pe, para_env_sub%mepos)
730 480 : my_ia_start = itmp(1)
731 480 : my_ia_end = itmp(2)
732 480 : my_ia_size = my_ia_end - my_ia_start + 1
733 :
734 480 : CALL get_group_dist(gd_B_virtual, para_env_sub%mepos, sizes=my_B_size, starts=my_B_virtual_start)
735 :
736 : ! reorder data
737 1914 : ALLOCATE (BIb_C_2D(my_group_L_size, my_ia_size))
738 :
739 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,ia_global) &
740 : !$OMP SHARED(homo,my_B_size,virtual,my_B_virtual_start,my_ia_start,my_ia_end,BIb_C,BIb_C_2D,&
741 480 : !$OMP my_group_L_size)
742 : DO iiB = 1, homo
743 : DO jjB = 1, my_B_size
744 : ia_global = (iiB - 1)*virtual + my_B_virtual_start + jjB - 1
745 : IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
746 : BIb_C_2D(1:my_group_L_size, ia_global - my_ia_start + 1) = BIb_C(1:my_group_L_size, jjB, iiB)
747 : END IF
748 : END DO
749 : END DO
750 :
751 480 : IF (para_env_sub%num_pe > 1) THEN
752 30 : ALLOCATE (BIb_C_rec_1D(INT(my_group_L_size, int_8)*maxsize(gd_B_virtual)*MIN(homo, occ_chunk)))
753 20 : DO proc_shift = 1, para_env_sub%num_pe - 1
754 10 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
755 10 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
756 :
757 10 : CALL get_group_dist(gd_B_virtual, proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
758 :
759 : ! do this in chunks to avoid high memory overhead
760 20 : DO occ_low = 1, homo, occ_chunk
761 10 : occ_high = MIN(homo, occ_low + occ_chunk - 1)
762 : BIb_C_rec(1:my_group_L_size, 1:rec_B_size, 1:occ_high - occ_low + 1) => &
763 10 : BIb_C_rec_1D(1:INT(my_group_L_size, int_8)*rec_B_size*(occ_high - occ_low + 1))
764 : CALL para_env_sub%sendrecv(BIb_C(:, :, occ_low:occ_high), proc_send, &
765 31970 : BIb_C_rec(:, :, 1:occ_high - occ_low + 1), proc_receive)
766 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,ia_global) &
767 : !$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,&
768 10 : !$OMP my_group_L_size)
769 : DO iiB = occ_low, occ_high
770 : DO jjB = 1, rec_B_size
771 : ia_global = (iiB - 1)*virtual + rec_B_virtual_start + jjB - 1
772 : IF (ia_global >= my_ia_start .AND. ia_global <= my_ia_end) THEN
773 : 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)
774 : END IF
775 : END DO
776 : END DO
777 : END DO
778 :
779 : END DO
780 10 : DEALLOCATE (BIb_C_rec_1D)
781 : END IF
782 :
783 480 : END SUBROUTINE calculate_BIb_C_2D
784 :
785 : ! **************************************************************************************************
786 : !> \brief ...
787 : !> \param BIb_C_2D ...
788 : !> \param para_env ...
789 : !> \param dimen_RI ...
790 : !> \param dimen_ia ...
791 : !> \param ext_row_block_size ...
792 : !> \param ext_col_block_size ...
793 : !> \param unit_nr ...
794 : !> \param my_ia_start ...
795 : !> \param my_ia_end ...
796 : !> \param my_group_L_start ...
797 : !> \param my_group_L_end ...
798 : !> \param para_env_RPA ...
799 : !> \param fm_mat_S ...
800 : !> \param nrow_block_mat ...
801 : !> \param ncol_block_mat ...
802 : !> \param blacs_env_ext ...
803 : !> \param blacs_env_ext_S ...
804 : !> \param dimen_ia_for_block_size ...
805 : !> \param do_im_time ...
806 : !> \param fm_mat_Q_gemm ...
807 : !> \param fm_mat_Q ...
808 : !> \param qs_env ...
809 : ! **************************************************************************************************
810 542 : SUBROUTINE create_integ_mat(BIb_C_2D, para_env, dimen_RI, dimen_ia, ext_row_block_size, ext_col_block_size, &
811 542 : unit_nr, my_ia_start, my_ia_end, my_group_L_start, my_group_L_end, para_env_RPA, &
812 542 : fm_mat_S, nrow_block_mat, ncol_block_mat, blacs_env_ext, blacs_env_ext_S, &
813 542 : dimen_ia_for_block_size, do_im_time, fm_mat_Q_gemm, fm_mat_Q, qs_env)
814 :
815 : TYPE(two_dim_real_array), DIMENSION(:), &
816 : INTENT(INOUT) :: BIb_C_2D
817 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
818 : INTEGER, INTENT(IN) :: dimen_RI
819 : INTEGER, DIMENSION(:), INTENT(IN) :: dimen_ia
820 : INTEGER, INTENT(IN) :: ext_row_block_size, ext_col_block_size, &
821 : unit_nr
822 : INTEGER, DIMENSION(:), INTENT(IN) :: my_ia_start, my_ia_end
823 : INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end
824 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_RPA
825 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_S
826 : INTEGER, INTENT(INOUT) :: nrow_block_mat, ncol_block_mat
827 : TYPE(cp_blacs_env_type), OPTIONAL, POINTER :: blacs_env_ext, blacs_env_ext_S
828 : INTEGER, INTENT(IN), OPTIONAL :: dimen_ia_for_block_size
829 : LOGICAL, INTENT(IN), OPTIONAL :: do_im_time
830 : TYPE(cp_fm_type), DIMENSION(:), OPTIONAL :: fm_mat_Q_gemm, fm_mat_Q
831 : TYPE(qs_environment_type), INTENT(IN), OPTIONAL, &
832 : POINTER :: qs_env
833 :
834 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_integ_mat'
835 :
836 : INTEGER :: col_row_proc_ratio, grid_2D(2), handle, &
837 : iproc, iproc_col, iproc_row, ispin
838 : LOGICAL :: my_blacs_ext, my_blacs_S_ext, &
839 : my_do_im_time
840 : TYPE(cp_blacs_env_type), POINTER :: blacs_env, blacs_env_Q
841 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
842 :
843 542 : CALL timeset(routineN, handle)
844 :
845 542 : CPASSERT(PRESENT(blacs_env_ext) .OR. PRESENT(dimen_ia_for_block_size))
846 :
847 542 : my_blacs_ext = .FALSE.
848 542 : IF (PRESENT(blacs_env_ext)) my_blacs_ext = .TRUE.
849 :
850 542 : my_blacs_S_ext = .FALSE.
851 542 : IF (PRESENT(blacs_env_ext_S)) my_blacs_S_ext = .TRUE.
852 :
853 542 : my_do_im_time = .FALSE.
854 542 : IF (PRESENT(do_im_time)) my_do_im_time = do_im_time
855 :
856 542 : NULLIFY (blacs_env)
857 : ! create the RPA blacs env
858 542 : IF (my_blacs_S_ext) THEN
859 200 : blacs_env => blacs_env_ext_S
860 : ELSE
861 342 : IF (para_env_RPA%num_pe > 1) THEN
862 224 : col_row_proc_ratio = MAX(1, dimen_ia_for_block_size/dimen_RI)
863 :
864 224 : 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
865 224 : DO iproc = 1, para_env_RPA%num_pe
866 224 : iproc_col = iproc_col - 1
867 224 : IF (MOD(para_env_RPA%num_pe, iproc_col) == 0) EXIT
868 : END DO
869 :
870 224 : iproc_row = para_env_RPA%num_pe/iproc_col
871 224 : grid_2D(1) = iproc_row
872 224 : grid_2D(2) = iproc_col
873 : ELSE
874 354 : grid_2D = 1
875 : END IF
876 342 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_RPA, grid_2d=grid_2D)
877 :
878 342 : IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
879 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
880 99 : "MATRIX_INFO| Number row processes:", grid_2D(1)
881 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
882 99 : "MATRIX_INFO| Number column processes:", grid_2D(2)
883 : END IF
884 :
885 : ! define the block_size for the row
886 342 : IF (ext_row_block_size > 0) THEN
887 0 : nrow_block_mat = ext_row_block_size
888 : ELSE
889 342 : nrow_block_mat = MAX(1, dimen_RI/grid_2D(1)/2)
890 : END IF
891 :
892 : ! define the block_size for the column
893 342 : IF (ext_col_block_size > 0) THEN
894 0 : ncol_block_mat = ext_col_block_size
895 : ELSE
896 342 : ncol_block_mat = MAX(1, dimen_ia_for_block_size/grid_2D(2)/2)
897 : END IF
898 :
899 342 : IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
900 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
901 99 : "MATRIX_INFO| Row block size:", nrow_block_mat
902 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
903 99 : "MATRIX_INFO| Column block size:", ncol_block_mat
904 : END IF
905 : END IF
906 :
907 470 : IF (.NOT. my_do_im_time) THEN
908 878 : DO ispin = 1, SIZE(BIb_C_2D)
909 480 : NULLIFY (fm_struct)
910 480 : IF (my_blacs_ext) THEN
911 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
912 240 : ncol_global=dimen_ia(ispin), para_env=para_env_RPA)
913 : ELSE
914 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
915 : ncol_global=dimen_ia(ispin), para_env=para_env_RPA, &
916 240 : nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.TRUE.)
917 :
918 : END IF ! external blacs_env
919 :
920 : CALL cp_array_to_fm(BIb_C_2D(ispin)%array, fm_struct, my_group_L_start, my_group_L_end, &
921 480 : my_ia_start(ispin), my_ia_end(ispin), fm_mat_S(ispin))
922 :
923 480 : CALL cp_fm_struct_release(fm_struct)
924 :
925 : ! sum the local data across processes belonging to different RPA group.
926 878 : IF (para_env_RPA%num_pe /= para_env%num_pe) THEN
927 : BLOCK
928 : TYPE(mp_comm_type) :: comm_exchange
929 180 : comm_exchange = fm_mat_S(ispin)%matrix_struct%context%interconnect(para_env)
930 180 : CALL comm_exchange%sum(fm_mat_S(ispin)%local_data)
931 360 : CALL comm_exchange%free()
932 : END BLOCK
933 : END IF
934 : END DO
935 : END IF
936 :
937 542 : IF (PRESENT(fm_mat_Q_gemm) .AND. .NOT. my_do_im_time) THEN
938 : ! create the Q matrix dimen_RIxdimen_RI where the result of the mat-mat-mult will be stored
939 198 : NULLIFY (fm_struct)
940 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
941 : ncol_global=dimen_RI, para_env=para_env_RPA, &
942 198 : nrow_block=nrow_block_mat, ncol_block=ncol_block_mat, force_block=.TRUE.)
943 438 : DO ispin = 1, SIZE(fm_mat_Q_gemm)
944 438 : CALL cp_fm_create(fm_mat_Q_gemm(ispin), fm_struct, name="fm_mat_Q_gemm")
945 : END DO
946 198 : CALL cp_fm_struct_release(fm_struct)
947 : END IF
948 :
949 542 : IF (PRESENT(fm_mat_Q)) THEN
950 430 : NULLIFY (blacs_env_Q)
951 430 : IF (my_blacs_ext) THEN
952 88 : blacs_env_Q => blacs_env_ext
953 342 : ELSE IF (para_env_RPA%num_pe == para_env%num_pe .AND. PRESENT(qs_env)) THEN
954 224 : CALL get_qs_env(qs_env, blacs_env=blacs_env_Q)
955 : ELSE
956 118 : CALL cp_blacs_env_create(blacs_env=blacs_env_Q, para_env=para_env_RPA)
957 : END IF
958 430 : NULLIFY (fm_struct)
959 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_Q, nrow_global=dimen_RI, &
960 430 : ncol_global=dimen_RI, para_env=para_env_RPA)
961 934 : DO ispin = 1, SIZE(fm_mat_Q)
962 934 : CALL cp_fm_create(fm_mat_Q(ispin), fm_struct, name="fm_mat_Q", set_zero=.TRUE.)
963 : END DO
964 :
965 430 : CALL cp_fm_struct_release(fm_struct)
966 :
967 430 : IF (.NOT. (my_blacs_ext .OR. (para_env_RPA%num_pe == para_env%num_pe .AND. PRESENT(qs_env)))) THEN
968 118 : CALL cp_blacs_env_release(blacs_env_Q)
969 : END IF
970 : END IF
971 :
972 : ! release blacs_env
973 542 : IF (.NOT. my_blacs_S_ext) THEN
974 342 : CALL cp_blacs_env_release(blacs_env)
975 : ELSE
976 200 : NULLIFY (blacs_env)
977 : END IF
978 :
979 542 : CALL timestop(handle)
980 :
981 542 : END SUBROUTINE create_integ_mat
982 :
983 : ! **************************************************************************************************
984 : !> \brief ...
985 : !> \param qs_env ...
986 : !> \param Erpa ...
987 : !> \param mp2_env ...
988 : !> \param para_env ...
989 : !> \param para_env_RPA ...
990 : !> \param para_env_sub ...
991 : !> \param unit_nr ...
992 : !> \param homo ...
993 : !> \param virtual ...
994 : !> \param dimen_RI ...
995 : !> \param dimen_RI_red ...
996 : !> \param dimen_ia ...
997 : !> \param dimen_nm_gw ...
998 : !> \param Eigenval ...
999 : !> \param num_integ_points ...
1000 : !> \param num_integ_group ...
1001 : !> \param color_rpa_group ...
1002 : !> \param fm_matrix_PQ ...
1003 : !> \param fm_mat_S ...
1004 : !> \param fm_mat_Q_gemm ...
1005 : !> \param fm_mat_Q ...
1006 : !> \param fm_mat_S_gw ...
1007 : !> \param fm_mat_R_gw ...
1008 : !> \param fm_mat_S_ij_bse ...
1009 : !> \param fm_mat_S_ab_bse ...
1010 : !> \param my_do_gw ...
1011 : !> \param do_bse ...
1012 : !> \param gw_corr_lev_occ ...
1013 : !> \param gw_corr_lev_virt ...
1014 : !> \param bse_lev_virt ...
1015 : !> \param do_minimax_quad ...
1016 : !> \param do_im_time ...
1017 : !> \param mo_coeff ...
1018 : !> \param cfm_matrix_L_kpoints ...
1019 : !> \param cfm_matrix_Minv_L_kpoints ...
1020 : !> \param fm_matrix_Minv_L_gamma ...
1021 : !> \param fm_matrix_Minv ...
1022 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
1023 : !> \param mat_munu ...
1024 : !> \param mat_P_global ...
1025 : !> \param t_3c_M ...
1026 : !> \param t_3c_O ...
1027 : !> \param t_3c_O_compressed ...
1028 : !> \param t_3c_O_ind ...
1029 : !> \param starts_array_mc ...
1030 : !> \param ends_array_mc ...
1031 : !> \param starts_array_mc_block ...
1032 : !> \param ends_array_mc_block ...
1033 : !> \param matrix_s ...
1034 : !> \param do_kpoints_from_Gamma ...
1035 : !> \param kpoints ...
1036 : !> \param gd_array ...
1037 : !> \param color_sub ...
1038 : !> \param do_ri_sos_laplace_mp2 ...
1039 : !> \param calc_forces ...
1040 : ! **************************************************************************************************
1041 342 : SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_sub, unit_nr, &
1042 342 : homo, virtual, dimen_RI, dimen_RI_red, dimen_ia, dimen_nm_gw, &
1043 : Eigenval, num_integ_points, num_integ_group, color_rpa_group, &
1044 684 : fm_matrix_PQ, fm_mat_S, fm_mat_Q_gemm, fm_mat_Q, fm_mat_S_gw, fm_mat_R_gw, &
1045 342 : fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
1046 342 : my_do_gw, do_bse, gw_corr_lev_occ, gw_corr_lev_virt, &
1047 342 : bse_lev_virt, &
1048 342 : do_minimax_quad, do_im_time, mo_coeff, &
1049 : cfm_matrix_L_kpoints, cfm_matrix_Minv_L_kpoints, &
1050 : fm_matrix_Minv_L_gamma, &
1051 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, mat_munu, mat_P_global, &
1052 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1053 : starts_array_mc, ends_array_mc, &
1054 : starts_array_mc_block, ends_array_mc_block, &
1055 : matrix_s, do_kpoints_from_Gamma, kpoints, gd_array, color_sub, &
1056 : do_ri_sos_laplace_mp2, calc_forces)
1057 :
1058 : TYPE(qs_environment_type), POINTER :: qs_env
1059 : REAL(KIND=dp), INTENT(OUT) :: Erpa
1060 : TYPE(mp2_type) :: mp2_env
1061 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_RPA, para_env_sub
1062 : INTEGER, INTENT(IN) :: unit_nr
1063 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
1064 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
1065 : INTEGER, DIMENSION(:), INTENT(IN) :: dimen_ia
1066 : INTEGER, INTENT(IN) :: dimen_nm_gw
1067 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1068 : INTENT(INOUT) :: Eigenval
1069 : INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, &
1070 : color_rpa_group
1071 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_PQ
1072 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_S
1073 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_Q_gemm
1074 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_Q
1075 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S_gw
1076 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_R_gw
1077 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S_ij_bse, fm_mat_S_ab_bse
1078 : LOGICAL, INTENT(IN) :: my_do_gw, do_bse
1079 : INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, &
1080 : bse_lev_virt
1081 : LOGICAL, INTENT(IN) :: do_minimax_quad, do_im_time
1082 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1083 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_matrix_L_kpoints, &
1084 : cfm_matrix_Minv_L_kpoints
1085 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv_L_gamma, fm_matrix_Minv, &
1086 : fm_matrix_Minv_Vtrunc_Minv
1087 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
1088 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_P_global
1089 : TYPE(dbt_type), INTENT(INOUT) :: t_3c_M
1090 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
1091 : INTENT(INOUT) :: t_3c_O
1092 : TYPE(hfx_compression_type), ALLOCATABLE, &
1093 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_compressed
1094 : TYPE(block_ind_type), ALLOCATABLE, &
1095 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_ind
1096 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
1097 : starts_array_mc_block, &
1098 : ends_array_mc_block
1099 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1100 : LOGICAL :: do_kpoints_from_Gamma
1101 : TYPE(kpoint_type), POINTER :: kpoints
1102 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1103 : INTEGER, INTENT(IN) :: color_sub
1104 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, calc_forces
1105 :
1106 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_num_int'
1107 :
1108 : COMPLEX(KIND=dp), ALLOCATABLE, &
1109 342 : DIMENSION(:, :, :, :) :: vec_Sigma_c_gw
1110 : INTEGER :: count_ev_sc_GW, cut_memory, group_size_P, gw_corr_lev_tot, handle, handle3, i, &
1111 : ikp_local, ispin, iter_evGW, iter_sc_GW0, j, jquad, min_bsize, mm_style, nkp, &
1112 : nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, Pspin, Qspin, &
1113 : size_P
1114 : INTEGER(int_8) :: dbcsr_nflop
1115 342 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c
1116 342 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c
1117 684 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, prim_blk_sizes, &
1118 342 : RI_blk_sizes
1119 : LOGICAL :: do_apply_ic_corr_to_gw, do_gw_im_time, do_ic_model, do_kpoints_cubic_RPA, &
1120 : do_periodic, do_print, do_ri_Sigma_x, do_static_cohsex, exit_ev_gw, first_cycle, &
1121 : first_cycle_periodic_correction, my_open_shell, print_ic_values
1122 342 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: has_mat_P_blocks
1123 : REAL(KIND=dp) :: alpha, dbcsr_time, e_exchange, e_exchange_corr, eps_filter, &
1124 : eps_filter_im_time, ext_scaling, fermi_level_offset, fermi_level_offset_input, &
1125 : my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old
1126 342 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, diag_Qomega, e_fermi, &
1127 342 : vec_omega_fit_gw, wkp_W
1128 342 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: static_cohsex_eigenval_cn, &
1129 684 : static_cohsex_eigenval_dg, vec_W_gw
1130 342 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_last, Eigenval_scf, &
1131 684 : static_cohsex_u, vec_Sigma_x_gw
1132 : TYPE(cp_cfm_type) :: cfm_mat_Q, cfm_mat_RI_global_work
1133 342 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_mat_L_kpoints, &
1134 342 : cfm_mat_Minv_L_kpoints, cfm_mo_coeff
1135 : TYPE(cp_fm_type) :: fm_mat_Minv_L_gamma, &
1136 : fm_mat_Q_static_bse_gemm, fm_mat_work
1137 342 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_S_gw_work, fm_mat_S_ia_bse, &
1138 684 : fm_mat_W, mo_coeff_bse
1139 : TYPE(dbcsr_p_type) :: mat_dm, mat_L, mat_M_P_munu_occ, &
1140 : mat_M_P_munu_virt, mat_MinvVMinv
1141 : TYPE(dbcsr_p_type), ALLOCATABLE, &
1142 342 : DIMENSION(:, :, :) :: mat_P_omega
1143 342 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_berry_im_mo_mo, &
1144 342 : matrix_berry_re_mo_mo
1145 342 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_kp
1146 : TYPE(dbcsr_type), POINTER :: mat_W, mat_work
1147 3078 : TYPE(dbt_type) :: t_3c_overl_int_ao_mo
1148 342 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_3c_overl_int_gw_AO, &
1149 342 : t_3c_overl_int_gw_RI, &
1150 342 : t_3c_overl_nnP_ic, &
1151 342 : t_3c_overl_nnP_ic_reflected
1152 : TYPE(dgemm_counter_type) :: dgemm_counter
1153 : TYPE(hfx_compression_type), ALLOCATABLE, &
1154 342 : DIMENSION(:) :: t_3c_O_mo_compressed
1155 28728 : TYPE(im_time_force_type) :: force_data
1156 342 : TYPE(rpa_exchange_work_type) :: exchange_work
1157 2052 : TYPE(rpa_grad_type) :: rpa_grad
1158 342 : TYPE(rpa_sigma_type) :: rpa_sigma
1159 342 : TYPE(time_frequency_grid_type) :: time_frequency_grid
1160 342 : TYPE(two_dim_int_array), ALLOCATABLE, DIMENSION(:) :: t_3c_O_mo_ind
1161 :
1162 342 : CALL timeset(routineN, handle)
1163 :
1164 342 : nspins = SIZE(homo)
1165 342 : nmo = homo(1) + virtual(1)
1166 :
1167 342 : my_open_shell = (nspins == 2)
1168 :
1169 342 : do_gw_im_time = my_do_gw .AND. do_im_time
1170 342 : do_ri_Sigma_x = mp2_env%ri_g0w0%do_ri_Sigma_x
1171 342 : do_ic_model = mp2_env%ri_g0w0%do_ic_model
1172 342 : print_ic_values = mp2_env%ri_g0w0%print_ic_values
1173 342 : do_periodic = mp2_env%ri_g0w0%do_periodic
1174 342 : do_kpoints_cubic_RPA = mp2_env%ri_rpa_im_time%do_im_time_kpoints
1175 342 : do_static_cohsex = my_do_gw .AND. mp2_env%ri_g0w0%do_static_cohsex
1176 :
1177 : ! For SOS-MP2 only gemm is implemented
1178 342 : mm_style = wfc_mm_style_gemm
1179 342 : IF (.NOT. do_ri_sos_laplace_mp2) mm_style = mp2_env%ri_rpa%mm_style
1180 :
1181 342 : IF (my_do_gw) THEN
1182 134 : ext_scaling = 0.2_dp
1183 134 : omega_max_fit = mp2_env%ri_g0w0%omega_max_fit
1184 134 : fermi_level_offset_input = mp2_env%ri_g0w0%fermi_level_offset
1185 134 : iter_evGW = mp2_env%ri_g0w0%iter_evGW
1186 134 : iter_sc_GW0 = mp2_env%ri_g0w0%iter_sc_GW0
1187 134 : IF ((.NOT. do_im_time)) THEN
1188 88 : IF (iter_sc_GW0 /= 1 .AND. iter_evGW /= 1) CPABORT("Mixed scGW0/evGW not implemented.")
1189 : ! in case of scGW0 with the N^4 algorithm, we use the evGW code but use the DFT eigenvalues for W
1190 88 : IF (iter_sc_GW0 /= 1) iter_evGW = iter_sc_GW0
1191 : END IF
1192 134 : IF (do_static_cohsex) THEN
1193 8 : IF (do_im_time) CPABORT("COHSEX is only implemented in the O(N^4) GW code.")
1194 8 : IF (do_periodic) CPABORT("COHSEX with PERIODIC_CORRECTION is not implemented.")
1195 : ! BSE after static COHSEX: W = v + α W^c(0) with one α for both A and B
1196 8 : IF (do_bse .AND. mp2_env%ri_g0w0%cohsex_scale_occ /= mp2_env%ri_g0w0%cohsex_scale_virt) THEN
1197 : CALL cp_abort(__LOCATION__, "BSE after COHSEX requires COHSEX_SCALE_OCC = COHSEX_SCALE_VIRT "// &
1198 0 : "(the scaling of W^c in the BSE).")
1199 : END IF
1200 8 : IF (do_ri_sos_laplace_mp2) CPABORT("COHSEX with SOS-MP2 is not implemented.")
1201 : ! static COHSEX is one-shot @DFT: Σ^COHSEX is static, no quasiparticle self-consistency
1202 8 : iter_evGW = 1
1203 8 : iter_sc_GW0 = 1
1204 : END IF
1205 : ELSE
1206 208 : ext_scaling = 0.0_dp
1207 208 : iter_evGW = 1
1208 208 : iter_sc_GW0 = 1
1209 : END IF
1210 :
1211 342 : IF (do_kpoints_cubic_RPA .AND. do_ri_sos_laplace_mp2) THEN
1212 0 : CPABORT("RI-SOS-Laplace-MP2 with k-point-sampling is not implemented.")
1213 : END IF
1214 :
1215 342 : do_apply_ic_corr_to_gw = .FALSE.
1216 342 : IF (mp2_env%ri_g0w0%ic_corr_list(1)%array(1) > 0.0_dp) do_apply_ic_corr_to_gw = .TRUE.
1217 :
1218 342 : IF (do_im_time) THEN
1219 144 : CPASSERT(do_minimax_quad .OR. do_ri_sos_laplace_mp2)
1220 :
1221 144 : group_size_P = mp2_env%ri_rpa_im_time%group_size_P
1222 144 : cut_memory = mp2_env%ri_rpa_im_time%cut_memory
1223 144 : eps_filter = mp2_env%ri_rpa_im_time%eps_filter
1224 : eps_filter_im_time = mp2_env%ri_rpa_im_time%eps_filter* &
1225 144 : mp2_env%ri_rpa_im_time%eps_filter_factor
1226 :
1227 144 : min_bsize = mp2_env%ri_rpa_im_time%min_bsize
1228 :
1229 : CALL alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, &
1230 : num_integ_points, nspins, calc_forces, fm_mat_Q(1), cfm_mo_coeff, &
1231 : cfm_matrix_Minv_L_kpoints, cfm_matrix_L_kpoints, fm_matrix_Minv_L_gamma, mat_P_global, &
1232 : t_3c_O, matrix_s, kpoints, eps_filter_im_time, &
1233 : cut_memory, nkp, num_cells_dm, num_3c_repl, &
1234 : size_P, ikp_local, &
1235 : index_to_cell_3c, &
1236 : cell_to_index_3c, &
1237 : col_blk_size, &
1238 : do_ic_model, do_kpoints_cubic_RPA, &
1239 : do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, &
1240 : has_mat_P_blocks, wkp_W, &
1241 : cfm_mat_Q, cfm_mat_Minv_L_kpoints, cfm_mat_L_kpoints, fm_mat_Minv_L_gamma, &
1242 : cfm_mat_RI_global_work, fm_mat_work, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, &
1243 144 : mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, mat_work, mo_coeff)
1244 :
1245 144 : IF (calc_forces) CALL init_im_time_forces(force_data, fm_matrix_PQ, t_3c_M, unit_nr, mp2_env, qs_env)
1246 :
1247 144 : IF (my_do_gw) THEN
1248 :
1249 : CALL dbcsr_get_info(mat_P_global%matrix, &
1250 46 : row_blk_size=RI_blk_sizes)
1251 :
1252 : CALL dbcsr_get_info(matrix_s(1)%matrix, &
1253 46 : row_blk_size=prim_blk_sizes)
1254 :
1255 46 : gw_corr_lev_tot = gw_corr_lev_occ(1) + gw_corr_lev_virt(1)
1256 :
1257 46 : IF (.NOT. do_kpoints_cubic_RPA) THEN
1258 : CALL allocate_matrices_gw_im_time(gw_corr_lev_occ, gw_corr_lev_virt, homo, nmo, &
1259 : num_integ_points, unit_nr, &
1260 : RI_blk_sizes, do_ic_model, &
1261 : para_env, fm_mat_W, fm_mat_Q(1), &
1262 : mo_coeff, &
1263 : t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
1264 : t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
1265 : starts_array_mc, ends_array_mc, &
1266 : t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
1267 : matrix_s, mat_W, t_3c_O, &
1268 : t_3c_O_compressed, t_3c_O_ind, &
1269 46 : qs_env)
1270 :
1271 : END IF
1272 : END IF
1273 :
1274 : END IF
1275 342 : IF (do_ic_model) THEN
1276 : ! image charge model only implemented for cubic scaling GW
1277 2 : CPASSERT(do_gw_im_time)
1278 2 : CPASSERT(.NOT. do_periodic)
1279 2 : IF (cut_memory /= 1) CPABORT("For IC, use MEMORY_CUT 1 in the LOW_SCALING section.")
1280 : END IF
1281 :
1282 1026 : ALLOCATE (e_fermi(nspins))
1283 342 : IF (do_minimax_quad .OR. do_ri_sos_laplace_mp2) THEN
1284 214 : do_print = .NOT. do_ic_model
1285 : CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, &
1286 : do_ri_sos_laplace_mp2, do_print, &
1287 214 : qs_env, do_gw_im_time, do_kpoints_cubic_RPA, e_fermi(1), time_frequency_grid)
1288 :
1289 : !For sos_laplace_mp2 and low-scaling RPA, potentially need to store/retrieve the initial weights
1290 214 : IF (qs_env%mp2_env%ri_rpa_im_time%keep_quad) THEN
1291 214 : CALL keep_initial_quad(time_frequency_grid, do_ri_sos_laplace_mp2, do_im_time, unit_nr, qs_env)
1292 : END IF
1293 : ELSE
1294 128 : IF (calc_forces) CPABORT("Forces with Clenshaw-Curtis grid not implemented.")
1295 : CALL get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, &
1296 : num_integ_group, color_rpa_group, fm_mat_S, my_do_gw, &
1297 128 : ext_scaling, time_frequency_grid)
1298 : END IF
1299 :
1300 : ! This array is needed for RPA
1301 342 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1302 834 : ALLOCATE (diag_Qomega(dimen_RI_red))
1303 : END IF
1304 :
1305 342 : IF (do_ri_sos_laplace_mp2 .AND. .NOT. do_im_time) THEN
1306 28 : alpha = 1.0_dp
1307 314 : ELSE IF (my_open_shell .OR. do_ri_sos_laplace_mp2) THEN
1308 90 : alpha = 2.0_dp
1309 : ELSE
1310 224 : alpha = 4.0_dp
1311 : END IF
1312 342 : IF (my_do_gw) THEN
1313 : CALL allocate_matrices_gw(vec_Sigma_c_gw, color_rpa_group, dimen_nm_gw, &
1314 : gw_corr_lev_occ, gw_corr_lev_virt, homo, &
1315 : nmo, num_integ_group, unit_nr, &
1316 : gw_corr_lev_tot, num_fit_points, omega_max_fit, &
1317 : do_minimax_quad, do_periodic, do_ri_Sigma_x,.NOT. do_im_time, &
1318 : first_cycle_periodic_correction, &
1319 : time_frequency_grid, Eigenval, vec_omega_fit_gw, vec_Sigma_x_gw, &
1320 : delta_corr, Eigenval_last, Eigenval_scf, vec_W_gw, &
1321 : fm_mat_S_gw, fm_mat_S_gw_work, &
1322 : para_env, mp2_env, kpoints, nkp, nkp_self_energy, &
1323 134 : do_kpoints_cubic_RPA, do_kpoints_from_Gamma)
1324 :
1325 134 : IF (do_bse) THEN
1326 :
1327 56 : CALL cp_fm_create(fm_mat_Q_static_bse_gemm, fm_mat_Q_gemm(1)%matrix_struct, set_zero=.TRUE.)
1328 :
1329 : END IF
1330 :
1331 : END IF
1332 :
1333 342 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_create(rpa_grad, fm_mat_Q(1), &
1334 : fm_mat_S, homo, virtual, mp2_env, Eigenval(:, 1, :), &
1335 44 : unit_nr, do_ri_sos_laplace_mp2)
1336 342 : IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) THEN
1337 : CALL exchange_work%create(qs_env, para_env_sub, mat_munu, dimen_RI_red, &
1338 170 : fm_mat_S, fm_mat_Q(1), fm_mat_Q_gemm(1), homo, virtual)
1339 : END IF
1340 342 : Erpa = 0.0_dp
1341 342 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) e_exchange = 0.0_dp
1342 342 : first_cycle = .TRUE.
1343 342 : omega_old = 0.0_dp
1344 342 : CALL dgemm_counter_init(dgemm_counter, unit_nr, mp2_env%ri_rpa%print_dgemm_info)
1345 :
1346 794 : DO count_ev_sc_GW = 1, iter_evGW
1347 472 : dbcsr_time = 0.0_dp
1348 472 : dbcsr_nflop = 0
1349 :
1350 472 : IF (do_ic_model) CYCLE
1351 :
1352 : ! reset some values, important when doing eigenvalue self-consistent GW
1353 470 : IF (my_do_gw) THEN
1354 262 : Erpa = 0.0_dp
1355 262 : vec_Sigma_c_gw = z_zero
1356 262 : first_cycle = .TRUE.
1357 : END IF
1358 :
1359 : ! calculate Q_PQ(it)
1360 470 : IF (do_im_time) THEN ! not using Imaginary time
1361 :
1362 156 : IF (.NOT. do_kpoints_cubic_RPA) THEN
1363 332 : DO ispin = 1, nspins
1364 332 : e_fermi(ispin) = (Eigenval(homo(ispin), 1, ispin) + Eigenval(homo(ispin) + 1, 1, ispin))*0.5_dp
1365 : END DO
1366 : END IF
1367 :
1368 156 : tau = 0.0_dp
1369 156 : tau_old = 0.0_dp
1370 :
1371 156 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(/T3,A,T66,i15)") &
1372 78 : "MEMORY_INFO| Memory cut:", cut_memory
1373 156 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
1374 78 : "SPARSITY_INFO| Eps filter for M virt/occ tensors:", eps_filter
1375 156 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
1376 78 : "SPARSITY_INFO| Eps filter for P matrix:", eps_filter_im_time
1377 156 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,i15)") &
1378 78 : "SPARSITY_INFO| Minimum tensor block size:", min_bsize
1379 :
1380 : ! compute the matrix Q(it) and Fourier transform it directly to mat_P_omega(iw)
1381 : CALL compute_mat_P_omega(mat_P_omega(:, :, 1), mo_coeff(1), homo(1), &
1382 : mat_P_global, matrix_s, 1, &
1383 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1384 : starts_array_mc, ends_array_mc, &
1385 : starts_array_mc_block, ends_array_mc_block, &
1386 : time_frequency_grid, e_fermi(1), eps_filter, alpha, &
1387 : eps_filter_im_time, Eigenval(:, 1, 1), nmo, &
1388 : cut_memory, &
1389 : unit_nr, mp2_env, para_env, &
1390 : qs_env, do_kpoints_from_Gamma, &
1391 : index_to_cell_3c, cell_to_index_3c, &
1392 : has_mat_P_blocks, do_ri_sos_laplace_mp2, &
1393 156 : dbcsr_time, dbcsr_nflop)
1394 :
1395 : ! the same for open shell, use the beta-spin MO coefficients
1396 156 : IF (my_open_shell) THEN
1397 : CALL compute_mat_P_omega(mat_P_omega(:, :, 2), mo_coeff(2), homo(2), &
1398 : mat_P_global, matrix_s, 2, &
1399 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1400 : starts_array_mc, ends_array_mc, &
1401 : starts_array_mc_block, ends_array_mc_block, &
1402 : time_frequency_grid, e_fermi(2), eps_filter, alpha, &
1403 : eps_filter_im_time, Eigenval(:, 1, 2), nmo, &
1404 : cut_memory, &
1405 : unit_nr, mp2_env, para_env, &
1406 : qs_env, do_kpoints_from_Gamma, &
1407 : index_to_cell_3c, cell_to_index_3c, &
1408 : has_mat_P_blocks, do_ri_sos_laplace_mp2, &
1409 32 : dbcsr_time, dbcsr_nflop)
1410 :
1411 : !For RPA, we sum up the P matrices. If no force needed, can clean-up the beta spin one
1412 32 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1413 90 : DO j = 1, SIZE(mat_P_omega, 2)
1414 598 : DO i = 1, SIZE(mat_P_omega, 1)
1415 508 : CALL dbcsr_add(mat_P_omega(i, j, 1)%matrix, mat_P_omega(i, j, 2)%matrix, 1.0_dp, 1.0_dp)
1416 578 : IF (.NOT. calc_forces) CALL dbcsr_clear(mat_P_omega(i, j, 2)%matrix)
1417 : END DO
1418 : END DO
1419 : END IF
1420 : END IF ! my_open_shell
1421 :
1422 : END IF ! do im time
1423 :
1424 470 : IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1425 10 : CALL rpa_sigma_create(rpa_sigma, mp2_env%ri_rpa%sigma_param, fm_mat_Q(1), unit_nr, para_env)
1426 : END IF
1427 :
1428 : ! ==============================================================================
1429 : ! static COHSEX: only the static screening at ω = 0, no frequency integration, no
1430 : ! analytic continuation and no RPA correlation energy (Erpa stays zero)
1431 : ! S^P_ia(0) = B^P_ia sqrt(1/(ϵ_a - ϵ_i))
1432 : ! Q_PQ(0) = 4 Σ_ia S^P_ia S^Q_ia (closed shell), 2 Σ_ia S^P_ia S^Q_ia per spin (open shell)
1433 : ! Q_PQ(0) = Σ_σ Q^σ_PQ(0) (spin sum, open shell)
1434 : ! ε_PQ(0) = δ_PQ + Q_PQ(0)
1435 : ! then W^c(0), Σ^c, Σ^x, H^QP and the diagonal/canonical static COHSEX energies
1436 : ! ==============================================================================
1437 470 : IF (do_static_cohsex) THEN
1438 8 : omega = 0.0_dp
1439 16 : DO ispin = 1, nspins
1440 : CALL calc_mat_Q(fm_mat_S(ispin), do_ri_sos_laplace_mp2, first_cycle, virtual(ispin), &
1441 : Eigenval(:, 1, ispin), homo(ispin), omega, omega_old, 1, mm_style, &
1442 : dimen_RI_red, dimen_ia(ispin), alpha, fm_mat_Q(ispin), &
1443 : fm_mat_Q_gemm(ispin), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
1444 16 : num_integ_points, count_ev_sc_GW)
1445 : END DO
1446 8 : DO ispin = 2, nspins
1447 8 : 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))
1448 : END DO
1449 : ! BSE: the static polarizability Q_PQ(0) = Σ_σ Q^σ_PQ(0);
1450 8 : IF (do_bse) THEN
1451 4 : CALL cp_fm_set_all(fm_mat_Q_static_bse_gemm, 0.0_dp)
1452 8 : DO ispin = 1, nspins
1453 8 : CALL cp_fm_scale_and_add(1.0_dp, fm_mat_Q_static_bse_gemm, 1.0_dp, fm_mat_Q_gemm(ispin))
1454 : END DO
1455 : END IF
1456 8 : CALL cp_fm_add_on_diag(fm_mat_Q(1), 1.0_dp, n_active=dimen_RI_red)
1457 8 : IF (do_bse) THEN
1458 : ! the BSE reads the quasiparticle energies from Eigenval: ϵ^CN (QP_UPDATE, with the
1459 : ! canonical MOs, see the BSE block) or ϵ^DG (DIAGONAL, with the KS MOs)
1460 4 : IF (mp2_env%ri_g0w0%cohsex_flavor == cohsex_qp_update) THEN
1461 16 : ALLOCATE (static_cohsex_eigenval_cn(nmo, nspins), static_cohsex_u(nmo, nmo, nspins))
1462 : CALL static_cohsex_qp_energies(fm_mat_Q(1), fm_mat_S_gw, fm_mat_S_gw_work, dimen_RI_red, nmo, &
1463 : homo, gw_corr_lev_occ, gw_corr_lev_virt, Eigenval(:, 1, :), &
1464 : mo_coeff, unit_nr, qs_env, mp2_env, &
1465 2 : eigenval_cn=static_cohsex_eigenval_cn, u_qp_all=static_cohsex_u)
1466 50 : Eigenval(:, 1, :) = static_cohsex_eigenval_cn(:, :)
1467 2 : DEALLOCATE (static_cohsex_eigenval_cn)
1468 : ELSE
1469 8 : ALLOCATE (static_cohsex_eigenval_dg(nmo, nspins))
1470 : CALL static_cohsex_qp_energies(fm_mat_Q(1), fm_mat_S_gw, fm_mat_S_gw_work, dimen_RI_red, nmo, &
1471 : homo, gw_corr_lev_occ, gw_corr_lev_virt, Eigenval(:, 1, :), &
1472 : mo_coeff, unit_nr, qs_env, mp2_env, &
1473 2 : eigenval_dg=static_cohsex_eigenval_dg)
1474 50 : Eigenval(:, 1, :) = static_cohsex_eigenval_dg(:, :)
1475 2 : DEALLOCATE (static_cohsex_eigenval_dg)
1476 : END IF
1477 : ! direct terms of the BSE with W = v + α W^c(0), α = COHSEX_SCALE_OCC (= COHSEX_SCALE_VIRT)
1478 4 : mp2_env%bse%wc_scaling = mp2_env%ri_g0w0%cohsex_scale_occ
1479 : ELSE
1480 : CALL static_cohsex_qp_energies(fm_mat_Q(1), fm_mat_S_gw, fm_mat_S_gw_work, dimen_RI_red, nmo, homo, &
1481 : gw_corr_lev_occ, gw_corr_lev_virt, Eigenval(:, 1, :), mo_coeff, &
1482 4 : unit_nr, qs_env, mp2_env)
1483 : END IF
1484 8 : first_cycle = .FALSE.
1485 : END IF
1486 :
1487 14302 : DO jquad = 1, num_integ_points
1488 : ! static COHSEX needs no frequency integral
1489 13840 : IF (do_static_cohsex) EXIT
1490 13832 : IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE
1491 :
1492 13059 : CALL timeset(routineN//"_RPA_matrix_operations", handle3)
1493 :
1494 13059 : IF (do_ri_sos_laplace_mp2) THEN
1495 206 : omega = time_frequency_grid%imaginary_time(jquad)
1496 : ELSE
1497 12853 : omega = time_frequency_grid%frequency(jquad)
1498 : END IF ! do_ri_sos_laplace_mp2
1499 :
1500 13059 : IF (do_im_time) THEN
1501 : ! in case we do imag time, we already calculated fm_mat_Q by a Fourier transform from im. time
1502 :
1503 1222 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
1504 :
1505 2428 : DO ispin = 1, SIZE(mat_P_omega, 3)
1506 : CALL contract_P_omega_with_mat_L(mat_P_omega(jquad, 1, ispin)%matrix, mat_L%matrix, mat_work, &
1507 : eps_filter_im_time, fm_mat_work, dimen_RI, dimen_RI_red, &
1508 2428 : fm_mat_Minv_L_gamma, fm_mat_Q(ispin))
1509 : END DO
1510 : END IF
1511 :
1512 : ELSE
1513 11837 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A, 1X, I3, 1X, A, 1X, I3)") &
1514 5915 : "INTEG_INFO| Started with Integration point", jquad, "of", num_integ_points
1515 :
1516 11837 : IF (first_cycle .AND. count_ev_sc_gw > 1) THEN
1517 116 : IF (iter_sc_gw0 == 1) THEN
1518 124 : DO ispin = 1, nspins
1519 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virtual(ispin), &
1520 124 : Eigenval_last(:, 1, ispin), homo(ispin), omega_old)
1521 : END DO
1522 : ELSE
1523 116 : DO ispin = 1, nspins
1524 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virtual(ispin), &
1525 116 : Eigenval_scf(:, 1, ispin), homo(ispin), omega_old)
1526 : END DO
1527 : END IF
1528 : END IF
1529 :
1530 11837 : IF (iter_sc_GW0 > 1) THEN
1531 12140 : DO ispin = 1, nspins
1532 : CALL calc_mat_Q(fm_mat_S(ispin), do_ri_sos_laplace_mp2, first_cycle, virtual(ispin), &
1533 : Eigenval_scf(:, 1, ispin), homo(ispin), omega, omega_old, jquad, mm_style, &
1534 : dimen_RI_red, dimen_ia(ispin), alpha, fm_mat_Q(ispin), &
1535 : fm_mat_Q_gemm(ispin), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
1536 12140 : num_integ_points, count_ev_sc_GW)
1537 : END DO
1538 :
1539 : ! For SOS-MP2 we need both matrices separately
1540 6070 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1541 6070 : DO ispin = 2, nspins
1542 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))
1543 : END DO
1544 : END IF
1545 : ELSE
1546 12208 : DO ispin = 1, nspins
1547 : CALL calc_mat_Q(fm_mat_S(ispin), do_ri_sos_laplace_mp2, first_cycle, virtual(ispin), &
1548 : Eigenval(:, 1, ispin), homo(ispin), omega, omega_old, jquad, mm_style, &
1549 : dimen_RI_red, dimen_ia(ispin), alpha, fm_mat_Q(ispin), &
1550 : fm_mat_Q_gemm(ispin), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
1551 12208 : num_integ_points, count_ev_sc_GW)
1552 : END DO
1553 : ! For open-shell BSE: the static screened-Coulomb polarizability is the
1554 : ! sum over both spin channels. calc_mat_Q overwrites fm_mat_Q_static_bse_gemm
1555 : ! per spin, so rebuild it here as the explicit spin sum at omega=0.
1556 5767 : IF (do_bse .AND. nspins > 1 .AND. jquad == num_integ_points .AND. &
1557 : count_ev_sc_GW == 1) THEN
1558 12 : CALL cp_fm_set_all(fm_mat_Q_static_bse_gemm, 0.0_dp)
1559 36 : DO ispin = 1, nspins
1560 : CALL cp_fm_scale_and_add(1.0_dp, fm_mat_Q_static_bse_gemm, &
1561 36 : 1.0_dp, fm_mat_Q_gemm(ispin))
1562 : END DO
1563 : END IF
1564 :
1565 : ! For SOS-MP2 we need both matrices separately
1566 5767 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1567 6353 : DO ispin = 2, nspins
1568 6353 : 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))
1569 : END DO
1570 : END IF
1571 :
1572 : END IF
1573 :
1574 : END IF ! im time
1575 :
1576 : ! Calculate RPA exchange energy correction
1577 13059 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
1578 12 : e_exchange_corr = 0.0_dp
1579 12 : CALL exchange_work%compute(fm_mat_Q(1), Eigenval(:, 1, :), fm_mat_S, omega, e_exchange_corr, mp2_env)
1580 :
1581 : ! Evaluate the final exchange energy correction
1582 12 : e_exchange = e_exchange + e_exchange_corr*time_frequency_grid%frequency_weights(jquad)
1583 : END IF
1584 :
1585 : ! for developing Sigma functional closed and open shell are taken cared for
1586 13059 : IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1587 30 : CALL rpa_sigma_matrix_spectral(rpa_sigma, fm_mat_Q(1), time_frequency_grid%frequency_weights(jquad), para_env_RPA)
1588 : END IF
1589 :
1590 13059 : IF (do_ri_sos_laplace_mp2) THEN
1591 :
1592 206 : CALL SOS_MP2_postprocessing(fm_mat_Q, Erpa, time_frequency_grid%time_weights_at_zero_frequency(jquad))
1593 :
1594 206 : IF (calc_forces .AND. .NOT. do_im_time) THEN
1595 : CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
1596 : fm_mat_Q, fm_mat_Q_gemm, dgemm_counter, fm_mat_S, omega, homo, virtual, &
1597 : Eigenval(:, 1, :), time_frequency_grid%time_weights_at_zero_frequency(jquad), &
1598 50 : unit_nr)
1599 : END IF
1600 : ELSE
1601 12853 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_copy_Q(fm_mat_Q(1), rpa_grad)
1602 :
1603 : ! epsilon = I + Q; retain the diagonal of Q for the RPA energy.
1604 12853 : CALL cp_fm_get_diag(fm_mat_Q(1), diag_Qomega, n_active=dimen_RI_red)
1605 12853 : CALL cp_fm_add_on_diag(fm_mat_Q(1), 1.0_dp, n_active=dimen_RI_red)
1606 :
1607 12853 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
1608 : CALL invert_eps_compute_W_and_Erpa_kp(dimen_RI, jquad, nkp, count_ev_sc_GW, para_env, &
1609 : Erpa, time_frequency_grid, &
1610 : wkp_W, do_gw_im_time, do_ri_Sigma_x, do_kpoints_from_Gamma, &
1611 : cfm_mat_Q, ikp_local, &
1612 : mat_P_omega(:, :, 1), mat_P_omega_kp, qs_env, eps_filter_im_time, unit_nr, &
1613 : kpoints, cfm_mat_Minv_L_kpoints, cfm_mat_L_kpoints, &
1614 : fm_mat_W, cfm_mat_RI_global_work, mat_MinvVMinv, &
1615 132 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv)
1616 : ELSE
1617 : CALL compute_Erpa_by_freq_int(dimen_RI_red, diag_Qomega, fm_mat_Q(1), para_env_RPA, Erpa, &
1618 12721 : time_frequency_grid%frequency_weights(jquad))
1619 : END IF
1620 :
1621 12853 : IF (calc_forces .AND. .NOT. do_im_time) THEN
1622 : CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
1623 : fm_mat_Q, fm_mat_Q_gemm, dgemm_counter, fm_mat_S, omega, homo, virtual, &
1624 56 : Eigenval(:, 1, :), time_frequency_grid%frequency_weights(jquad), unit_nr)
1625 : END IF
1626 : END IF ! do_ri_sos_laplace_mp2
1627 :
1628 : ! save omega and reset the first_cycle flag
1629 13059 : first_cycle = .FALSE.
1630 13059 : omega_old = omega
1631 :
1632 13059 : CALL timestop(handle3)
1633 :
1634 13059 : IF (my_do_gw) THEN
1635 :
1636 12188 : CALL get_fermi_level_offset(fermi_level_offset, fermi_level_offset_input, Eigenval(:, 1, :), homo)
1637 :
1638 : ! do_im_time = TRUE means low-scaling calculation
1639 12188 : IF (do_im_time) THEN
1640 : ! only for molecules
1641 818 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
1642 : CALL compute_W_cubic_GW(fm_mat_W, fm_mat_Q(1), fm_mat_work, dimen_RI, fm_mat_Minv_L_gamma, &
1643 722 : time_frequency_grid, jquad, omega)
1644 : END IF
1645 : ELSE
1646 : CALL compute_GW_self_energy(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI_red, gw_corr_lev_occ, &
1647 : gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, &
1648 : do_im_time, do_periodic, first_cycle_periodic_correction, &
1649 : fermi_level_offset, &
1650 : omega, Eigenval(:, 1, :), delta_corr, vec_omega_fit_gw, vec_W_gw, &
1651 : time_frequency_grid, &
1652 : fm_mat_Q(1), fm_mat_R_gw, fm_mat_S_gw, &
1653 : fm_mat_S_gw_work, mo_coeff(1), para_env, &
1654 : para_env_RPA, matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, &
1655 11370 : kpoints, qs_env, mp2_env)
1656 : END IF
1657 : END IF
1658 :
1659 13059 : IF (unit_nr > 0) CALL m_flush(unit_nr)
1660 27361 : CALL para_env_RPA%sync() ! sync to see output
1661 :
1662 : END DO ! jquad
1663 :
1664 470 : IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
1665 10 : CALL finalize_rpa_sigma(rpa_sigma, unit_nr, mp2_env%ri_rpa%e_sigma_corr, para_env, do_minimax_quad)
1666 10 : IF (do_minimax_quad) mp2_env%ri_rpa%e_sigma_corr = mp2_env%ri_rpa%e_sigma_corr/2.0_dp
1667 10 : CALL para_env%sum(mp2_env%ri_rpa%e_sigma_corr)
1668 : END IF
1669 :
1670 470 : CALL para_env%sum(Erpa)
1671 :
1672 470 : IF (.NOT. (do_ri_sos_laplace_mp2)) THEN
1673 406 : Erpa = Erpa/(pi*2.0_dp)
1674 406 : IF (do_minimax_quad) Erpa = Erpa/2.0_dp
1675 : END IF
1676 :
1677 470 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
1678 12 : CALL para_env%sum(E_exchange)
1679 12 : E_exchange = E_exchange/(pi*2.0_dp)
1680 12 : IF (do_minimax_quad) E_exchange = E_exchange/2.0_dp
1681 12 : mp2_env%ri_rpa%ener_exchange = E_exchange
1682 : END IF
1683 :
1684 470 : IF (calc_forces .AND. do_ri_sos_laplace_mp2 .AND. do_im_time) THEN
1685 22 : IF (my_open_shell) THEN
1686 4 : Pspin = 1
1687 4 : Qspin = 2
1688 : CALL calc_laplace_loop_forces(force_data, mat_P_omega(:, 1, :), t_3c_M, t_3c_O(1, 1), &
1689 : t_3c_O_compressed(1, 1, :), t_3c_O_ind(1, 1, :), mo_coeff, homo, &
1690 : starts_array_mc, ends_array_mc, starts_array_mc_block, &
1691 : ends_array_mc_block, nmo, Eigenval(:, 1, :), &
1692 : time_frequency_grid, &
1693 : cut_memory, Pspin, Qspin, my_open_shell, &
1694 4 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1695 4 : Pspin = 2
1696 4 : Qspin = 1
1697 : CALL calc_laplace_loop_forces(force_data, mat_P_omega(:, 1, :), t_3c_M, t_3c_O(1, 1), &
1698 : t_3c_O_compressed(1, 1, :), t_3c_O_ind(1, 1, :), mo_coeff, homo, &
1699 : starts_array_mc, ends_array_mc, starts_array_mc_block, &
1700 : ends_array_mc_block, nmo, Eigenval(:, 1, :), &
1701 : time_frequency_grid, &
1702 : cut_memory, Pspin, Qspin, my_open_shell, &
1703 4 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1704 :
1705 : ELSE
1706 18 : Pspin = 1
1707 18 : Qspin = 1
1708 : CALL calc_laplace_loop_forces(force_data, mat_P_omega(:, 1, :), t_3c_M, t_3c_O(1, 1), &
1709 : t_3c_O_compressed(1, 1, :), t_3c_O_ind(1, 1, :), mo_coeff, homo, &
1710 : starts_array_mc, ends_array_mc, starts_array_mc_block, &
1711 : ends_array_mc_block, nmo, Eigenval(:, 1, :), &
1712 : time_frequency_grid, &
1713 : cut_memory, Pspin, Qspin, my_open_shell, &
1714 18 : unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
1715 : END IF
1716 22 : CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1717 : END IF !laplace SOS-MP2
1718 :
1719 470 : IF (calc_forces .AND. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) THEN
1720 64 : DO ispin = 1, nspins
1721 : CALL calc_rpa_loop_forces(force_data, mat_P_omega(:, 1, :), t_3c_M, t_3c_O(1, 1), &
1722 : t_3c_O_compressed(1, 1, :), t_3c_O_ind(1, 1, :), mo_coeff, homo, &
1723 : starts_array_mc, ends_array_mc, starts_array_mc_block, &
1724 : ends_array_mc_block, nmo, Eigenval(:, 1, :), &
1725 : e_fermi(ispin), time_frequency_grid, cut_memory, ispin, my_open_shell, &
1726 : unit_nr, dbcsr_time, &
1727 64 : dbcsr_nflop, mp2_env, qs_env)
1728 : END DO
1729 28 : CALL calc_post_loop_forces(force_data, unit_nr, qs_env)
1730 : END IF
1731 :
1732 470 : IF (do_im_time) THEN
1733 :
1734 156 : my_flop_rate = REAL(dbcsr_nflop, dp)/(1.0E09_dp*dbcsr_time)
1735 156 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(/T3,A,T73,ES8.2)") &
1736 78 : "PERFORMANCE| DBCSR total number of flops:", REAL(dbcsr_nflop*para_env%num_pe, dp)
1737 156 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.2)") &
1738 78 : "PERFORMANCE| DBCSR total execution time:", dbcsr_time
1739 156 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.2)") &
1740 78 : "PERFORMANCE| DBCSR flop rate (Gflops / MPI rank):", my_flop_rate
1741 :
1742 : ELSE
1743 :
1744 314 : CALL dgemm_counter_write(dgemm_counter, para_env)
1745 :
1746 : END IF
1747 :
1748 : ! GW: for low-scaling calculation: Compute self-energy Sigma(i*tau), Sigma(i*omega)
1749 : ! for low-scaling and ordinary-scaling: analytic continuation from Sigma(iw) -> Sigma(w)
1750 : ! and correction of quasiparticle energies e_n^GW
1751 : ! (static COHSEX energies are already computed and printed in static_cohsex_qp_energies)
1752 812 : IF (my_do_gw .AND. .NOT. do_static_cohsex) THEN
1753 :
1754 : CALL compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, &
1755 : gw_corr_lev_tot, gw_corr_lev_virt, homo, &
1756 : nmo, num_fit_points, &
1757 : unit_nr, do_apply_ic_corr_to_gw, do_im_time, &
1758 : do_periodic, do_ri_Sigma_x, first_cycle_periodic_correction, &
1759 : e_fermi, eps_filter, fermi_level_offset, &
1760 : delta_corr, Eigenval, &
1761 : Eigenval_last, Eigenval_scf, iter_sc_GW0, exit_ev_gw, &
1762 : time_frequency_grid, vec_omega_fit_gw, vec_Sigma_x_gw, &
1763 : mp2_env%ri_g0w0%ic_corr_list, &
1764 : mo_coeff, fm_mat_W, para_env, &
1765 : para_env_RPA, mat_dm, mat_MinvVMinv, &
1766 : t_3c_O, t_3c_M, t_3c_overl_int_ao_mo, t_3c_O_compressed, t_3c_O_mo_compressed, &
1767 : t_3c_O_ind, t_3c_O_mo_ind, &
1768 : t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
1769 : matrix_berry_im_mo_mo, matrix_berry_re_mo_mo, mat_W, matrix_s, &
1770 : kpoints, mp2_env, qs_env, nkp_self_energy, do_kpoints_cubic_RPA, &
1771 254 : starts_array_mc, ends_array_mc)
1772 :
1773 : ! if HOMO-LUMO gap differs by less than mp2_env%ri_g0w0%eps_ev_sc_iter, exit ev sc GW loop
1774 254 : IF (exit_ev_gw) EXIT
1775 :
1776 : END IF ! my_do_gw if
1777 :
1778 : END DO ! evGW loop
1779 :
1780 342 : IF (do_ic_model) THEN
1781 :
1782 2 : IF (my_open_shell) THEN
1783 :
1784 : CALL calculate_ic_correction(Eigenval(:, 1, 1), mat_MinvVMinv%matrix, &
1785 : t_3c_overl_nnP_ic(1), t_3c_overl_nnP_ic_reflected(1), &
1786 : gw_corr_lev_tot, &
1787 : gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), unit_nr, &
1788 0 : print_ic_values, para_env, do_alpha=.TRUE.)
1789 :
1790 : CALL calculate_ic_correction(Eigenval(:, 1, 2), mat_MinvVMinv%matrix, &
1791 : t_3c_overl_nnP_ic(2), t_3c_overl_nnP_ic_reflected(2), &
1792 : gw_corr_lev_tot, &
1793 : gw_corr_lev_occ(2), gw_corr_lev_virt(2), homo(2), unit_nr, &
1794 0 : print_ic_values, para_env, do_beta=.TRUE.)
1795 :
1796 : ELSE
1797 :
1798 : CALL calculate_ic_correction(Eigenval(:, 1, 1), mat_MinvVMinv%matrix, &
1799 : t_3c_overl_nnP_ic(1), t_3c_overl_nnP_ic_reflected(1), &
1800 : gw_corr_lev_tot, &
1801 : gw_corr_lev_occ(1), gw_corr_lev_virt(1), homo(1), unit_nr, &
1802 2 : print_ic_values, para_env)
1803 :
1804 : END IF
1805 :
1806 : END IF
1807 :
1808 : ! postprocessing after GW for Bethe-Salpeter
1809 342 : IF (do_bse) THEN
1810 : ! Check used GW flavor; in Case of evGW we use W0 for BSE
1811 : ! Use environment variable, since local iter_evGW is overwritten if evGW0 is invoked
1812 56 : IF (mp2_env%ri_g0w0%iter_evGW > 1) THEN
1813 4 : IF (unit_nr > 0) THEN
1814 : CALL cp_warn(__LOCATION__, &
1815 2 : "BSE@evGW applies W0, i.e. screening with DFT energies to the BSE!")
1816 : END IF
1817 : END IF
1818 : ! Create a per-spin copy of fm_mat_S for usage in BSE
1819 236 : ALLOCATE (fm_mat_S_ia_bse(nspins))
1820 124 : DO ispin = 1, nspins
1821 68 : CALL cp_fm_create(fm_mat_S_ia_bse(ispin), fm_mat_S(ispin)%matrix_struct)
1822 68 : CALL cp_fm_to_fm(fm_mat_S(ispin), fm_mat_S_ia_bse(ispin))
1823 : ! Remove energy/frequency factor from 3c-Integral for BSE
1824 124 : IF (iter_sc_gw0 == 1) THEN
1825 : CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse(ispin), virtual(ispin), &
1826 56 : Eigenval_last(:, 1, ispin), homo(ispin), omega)
1827 : ELSE
1828 : CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse(ispin), virtual(ispin), &
1829 12 : Eigenval_scf(:, 1, ispin), homo(ispin), omega)
1830 : END IF
1831 : END DO
1832 : ! ==============================================================================
1833 : ! BSE with the canonical static COHSEX quasiparticles: the KS-MO slabs are rotated to the
1834 : ! canonical MOs C^QP = C U (U block diagonal in occ/virt, δ outside the window)
1835 : ! B̃^P_ia = Σ_jb U_ji U_ba B^P_jb, B̃^P_ij = Σ_kl U_ki U_lj B^P_kl, B̃^P_ab analogous
1836 : ! Q(0), ε^(-1)(0) and the kernel are invariant under these rotations, so the BSE with
1837 : ! ϵ^CN and C^QP equals the BSE with the full H^QP in the KS basis.
1838 : ! ==============================================================================
1839 56 : IF (do_static_cohsex .AND. mp2_env%ri_g0w0%cohsex_flavor == cohsex_qp_update) THEN
1840 6 : ALLOCATE (mo_coeff_bse(nspins))
1841 4 : DO ispin = 1, nspins
1842 : ! the ab slab holds the first bse_lev_virt virtuals: U must not mix beyond them
1843 2 : IF (gw_corr_lev_virt(ispin) > bse_lev_virt(ispin)) THEN
1844 0 : CPABORT("Canonical static COHSEX BSE: the corrected virtual MOs exceed the BSE (ab) slab.")
1845 : END IF
1846 : CALL rotate_bse_slab(fm_mat_S_ia_bse(ispin), dimen_RI_red, homo(ispin), virtual(ispin), &
1847 : static_cohsex_u(1:homo(ispin), 1:homo(ispin), ispin), &
1848 2 : static_cohsex_u(homo(ispin) + 1:nmo, homo(ispin) + 1:nmo, ispin))
1849 : CALL rotate_bse_slab(fm_mat_S_ij_bse(ispin), dimen_RI_red, homo(ispin), homo(ispin), &
1850 : static_cohsex_u(1:homo(ispin), 1:homo(ispin), ispin), &
1851 2 : static_cohsex_u(1:homo(ispin), 1:homo(ispin), ispin))
1852 : CALL rotate_bse_slab(fm_mat_S_ab_bse(ispin), dimen_RI_red, bse_lev_virt(ispin), bse_lev_virt(ispin), &
1853 : static_cohsex_u(homo(ispin) + 1:homo(ispin) + bse_lev_virt(ispin), &
1854 : homo(ispin) + 1:homo(ispin) + bse_lev_virt(ispin), ispin), &
1855 : static_cohsex_u(homo(ispin) + 1:homo(ispin) + bse_lev_virt(ispin), &
1856 2 : homo(ispin) + 1:homo(ispin) + bse_lev_virt(ispin), ispin))
1857 : ! transition dipoles, NTOs and exciton descriptors with C^QP
1858 4 : CALL get_static_cohsex_mo_coeff(mo_coeff(ispin), static_cohsex_u(:, :, ispin), mo_coeff_bse(ispin))
1859 : END DO
1860 : END IF
1861 56 : IF (ALLOCATED(static_cohsex_u)) DEALLOCATE (static_cohsex_u)
1862 :
1863 : ! the solver keeps [1+Q(0)]^-1 and the window-sized slabs; Q(0) and the full slabs stay the GW's
1864 : CALL prepare_bse_env_from_ri_rpa_gw(mp2_env%bse, fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
1865 : fm_mat_Q_static_bse_gemm, Eigenval, Eigenval_scf, &
1866 : homo, virtual, dimen_RI, dimen_RI_red, mp2_env%ri_g0w0%corr_mos_occ, &
1867 56 : bse_lev_virt, qs_env, unit_nr)
1868 56 : CALL cp_fm_release(fm_mat_Q_static_bse_gemm)
1869 56 : IF (ALLOCATED(mo_coeff_bse)) THEN
1870 2 : CALL bse_solve(mp2_env%bse, mo_coeff_bse, unit_nr)
1871 2 : CALL cp_fm_release(mo_coeff_bse)
1872 : ELSE
1873 54 : CALL bse_solve(mp2_env%bse, mo_coeff, unit_nr)
1874 : END IF
1875 56 : CALL bse_env_release(mp2_env%bse)
1876 : ! Release per-spin BSE-copy of fm_mat_S
1877 124 : DO ispin = 1, nspins
1878 124 : CALL cp_fm_release(fm_mat_S_ia_bse(ispin))
1879 : END DO
1880 56 : DEALLOCATE (fm_mat_S_ia_bse)
1881 : END IF
1882 :
1883 342 : IF (my_do_gw) THEN
1884 : CALL deallocate_matrices_gw(fm_mat_S_gw_work, vec_W_gw, vec_Sigma_c_gw, vec_omega_fit_gw, &
1885 : mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw, &
1886 : Eigenval_last, Eigenval_scf, do_periodic, matrix_berry_re_mo_mo, matrix_berry_im_mo_mo, &
1887 134 : kpoints, vec_Sigma_x_gw,.NOT. do_im_time)
1888 : END IF
1889 :
1890 342 : IF (do_im_time) THEN
1891 :
1892 : CALL dealloc_im_time(cfm_mo_coeff, index_to_cell_3c, &
1893 : cell_to_index_3c, do_ic_model, &
1894 : do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, &
1895 : has_mat_P_blocks, &
1896 : wkp_W, cfm_mat_Q, cfm_mat_Minv_L_kpoints, cfm_mat_L_kpoints, &
1897 : fm_mat_Minv_L_gamma, &
1898 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, cfm_mat_RI_global_work, fm_mat_work, &
1899 : mat_dm, mat_L, &
1900 : mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
1901 144 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, mat_work, qs_env)
1902 :
1903 144 : IF (my_do_gw) THEN
1904 : CALL deallocate_matrices_gw_im_time(do_ic_model, do_kpoints_cubic_RPA, fm_mat_W, &
1905 : t_3c_overl_int_ao_mo, t_3c_O_mo_compressed, t_3c_O_mo_ind, &
1906 : t_3c_overl_int_gw_RI, t_3c_overl_int_gw_AO, &
1907 : t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, &
1908 46 : mat_W, qs_env)
1909 : END IF
1910 :
1911 : END IF
1912 :
1913 342 : IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) CALL exchange_work%release()
1914 :
1915 342 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
1916 278 : DEALLOCATE (diag_Qomega)
1917 : END IF
1918 :
1919 342 : CALL time_frequency_grid_release(time_frequency_grid)
1920 :
1921 342 : IF (do_im_time .AND. calc_forces) THEN
1922 50 : CALL im_time_force_release(force_data)
1923 : END IF
1924 :
1925 342 : IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_finalize(rpa_grad, mp2_env, para_env_sub, para_env, &
1926 : qs_env, gd_array, color_sub, do_ri_sos_laplace_mp2, &
1927 44 : homo, virtual)
1928 :
1929 342 : CALL timestop(handle)
1930 :
1931 8330 : END SUBROUTINE rpa_num_int
1932 :
1933 : ! **************************************************************************************************
1934 : !> \brief ...
1935 : !> \param para_env ...
1936 : !> \param unit_nr ...
1937 : !> \param homo ...
1938 : !> \param Eigenval ...
1939 : !> \param num_integ_points ...
1940 : !> \param do_im_time ...
1941 : !> \param do_ri_sos_laplace_mp2 ...
1942 : !> \param do_print ...
1943 : !> \param qs_env ...
1944 : !> \param do_gw_im_time ...
1945 : !> \param do_kpoints_cubic_RPA ...
1946 : !> \param e_fermi ...
1947 : !> \param grid ...
1948 : ! **************************************************************************************************
1949 214 : SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, &
1950 : do_im_time, do_ri_sos_laplace_mp2, do_print, qs_env, do_gw_im_time, &
1951 : do_kpoints_cubic_RPA, e_fermi, grid)
1952 :
1953 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1954 : INTEGER, INTENT(IN) :: unit_nr
1955 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
1956 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval
1957 : INTEGER, INTENT(IN) :: num_integ_points
1958 : LOGICAL, INTENT(IN) :: do_im_time, do_ri_sos_laplace_mp2, &
1959 : do_print
1960 : TYPE(qs_environment_type), POINTER :: qs_env
1961 : LOGICAL, INTENT(IN) :: do_gw_im_time, do_kpoints_cubic_RPA
1962 : REAL(KIND=dp), INTENT(OUT) :: e_fermi
1963 : TYPE(time_frequency_grid_type), INTENT(OUT) :: grid
1964 :
1965 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_minimax_grid'
1966 : INTEGER, PARAMETER :: num_points_per_magnitude = 200
1967 :
1968 : INTEGER :: handle, jquad
1969 : LOGICAL :: used_external_backend
1970 : REAL(KIND=dp) :: E_Range, Emax, Emin, max_error_min
1971 :
1972 214 : CALL timeset(routineN, handle)
1973 :
1974 : CALL determine_energy_range(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, &
1975 214 : do_kpoints_cubic_RPA, Emin, Emax, e_range, e_fermi)
1976 :
1977 : ! Open-shell uses ONE minimax grid for the combined [min gap, max span] over both spins (the
1978 : ! superset covers each channel, so it is accurate; per-spin grids would only be more efficient).
1979 214 : IF (SIZE(homo) > 1) THEN
1980 : CALL cp_hint(__LOCATION__, &
1981 : "Open-shell RPA/GW uses one minimax grid spanning [min gap, max span] across "// &
1982 : "both spin channels; raise QUADRATURE_POINTS if QP convergence is marginal for "// &
1983 50 : "strongly spin-asymmetric systems.")
1984 : END IF
1985 :
1986 : CALL build_minimax_time_frequency_grid(num_integ_points, Emin, Emax, &
1987 : qs_env%mp2_env%ri_g0w0%regularization_minimax, &
1988 : num_points_per_magnitude, grid, &
1989 : build_frequency=.NOT. do_ri_sos_laplace_mp2, &
1990 : build_time=do_im_time .OR. do_ri_sos_laplace_mp2, &
1991 : build_transforms=do_im_time .AND. .NOT. do_ri_sos_laplace_mp2, &
1992 : build_sine=do_im_time .AND. (.NOT. do_ri_sos_laplace_mp2) .AND. do_gw_im_time, &
1993 : time_scaling=MERGE(1.0_dp, 2.0_dp, do_ri_sos_laplace_mp2), &
1994 : time_weight_scaling=MERGE(1.0_dp, 2.0_dp, do_ri_sos_laplace_mp2), &
1995 : max_fit_error=max_error_min, print_warning=.TRUE., unit_nr=unit_nr, &
1996 680 : prefer_external_backend=.TRUE., used_external_backend=used_external_backend)
1997 :
1998 : ! Keep the native diagnostics and warning behavior in this RPA policy wrapper. The external
1999 : ! backend reports its diagnostics from the grid builder.
2000 214 : IF (used_external_backend) THEN
2001 78 : CALL timestop(handle)
2002 78 : RETURN
2003 : END IF
2004 :
2005 136 : IF (num_integ_points > 20 .AND. e_range < 100.0_dp) THEN
2006 0 : IF (unit_nr > 0) THEN
2007 : CALL cp_warn(__LOCATION__, &
2008 : "You requested a large minimax grid (> 20 points) for a small minimax range R (R < 100). "// &
2009 : "That may lead to numerical "// &
2010 : "instabilities when computing minimax grid weights. You can prevent small ranges by choosing "// &
2011 0 : "a larger basis set with higher angular momenta or alternatively using all-electron calculations.")
2012 : END IF
2013 : END IF
2014 :
2015 136 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
2016 72 : IF (unit_nr > 0 .AND. do_print) THEN
2017 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
2018 35 : "MINIMAX_INFO| Number of integration points:", num_integ_points
2019 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
2020 35 : "MINIMAX_INFO| Gap for the minimax approximation:", Emin
2021 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
2022 35 : "MINIMAX_INFO| Range for the minimax approximation:", e_range
2023 35 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") "MINIMAX_INFO| Minimax parameters:", "Weights", "Abscissas"
2024 167 : DO jquad = 1, num_integ_points
2025 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") &
2026 167 : grid%frequency_weights(jquad)/Emin, grid%frequency(jquad)/Emin
2027 : END DO
2028 35 : CALL m_flush(unit_nr)
2029 : END IF
2030 : END IF
2031 :
2032 136 : IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
2033 106 : IF (unit_nr > 0 .AND. do_print) THEN
2034 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
2035 52 : "MINIMAX_INFO| Range for the minimax approximation:", e_range
2036 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
2037 52 : "MINIMAX_INFO| Gap:", Emin
2038 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
2039 52 : "MINIMAX_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas"
2040 246 : DO jquad = 1, num_integ_points
2041 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") &
2042 246 : grid%time_weights_at_zero_frequency(jquad)*Emin, grid%imaginary_time(jquad)*Emin
2043 : END DO
2044 52 : CALL m_flush(unit_nr)
2045 : END IF
2046 :
2047 106 : IF (unit_nr > 0 .AND. do_im_time .AND. do_gw_im_time .AND. .NOT. do_ri_sos_laplace_mp2) THEN
2048 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") &
2049 5 : "MINIMAX_INFO| Maximum deviation among requested minimax transform fits:", max_error_min
2050 : END IF
2051 : END IF
2052 :
2053 136 : CALL timestop(handle)
2054 :
2055 : END SUBROUTINE get_minimax_grid
2056 :
2057 : ! **************************************************************************************************
2058 : !> \brief Construct a Clenshaw-Curtis grid after determining its RPA scaling.
2059 : !> \param para_env ...
2060 : !> \param para_env_RPA ...
2061 : !> \param unit_nr ...
2062 : !> \param homo ...
2063 : !> \param virtual ...
2064 : !> \param Eigenval ...
2065 : !> \param num_integ_points ...
2066 : !> \param num_integ_group ...
2067 : !> \param color_rpa_group ...
2068 : !> \param fm_mat_S ...
2069 : !> \param my_do_gw ...
2070 : !> \param ext_scaling ...
2071 : !> \param grid ...
2072 : ! **************************************************************************************************
2073 256 : SUBROUTINE get_clenshaw_grid(para_env, para_env_RPA, unit_nr, homo, virtual, Eigenval, num_integ_points, &
2074 128 : num_integ_group, color_rpa_group, fm_mat_S, my_do_gw, &
2075 : ext_scaling, grid)
2076 :
2077 : TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_RPA
2078 : INTEGER, INTENT(IN) :: unit_nr
2079 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
2080 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval
2081 : INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, &
2082 : color_rpa_group
2083 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
2084 : LOGICAL, INTENT(IN) :: my_do_gw
2085 : REAL(KIND=dp), INTENT(IN) :: ext_scaling
2086 : TYPE(time_frequency_grid_type), INTENT(OUT) :: grid
2087 :
2088 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_clenshaw_grid'
2089 :
2090 : INTEGER :: handle
2091 : REAL(KIND=dp) :: a_scaling
2092 :
2093 128 : CALL timeset(routineN, handle)
2094 :
2095 128 : CALL build_clenshaw_grid(num_integ_points, grid)
2096 :
2097 128 : IF (my_do_gw .AND. ext_scaling > 0.0_dp) THEN
2098 88 : a_scaling = ext_scaling
2099 : ELSE
2100 : CALL calc_scaling_factor(a_scaling, para_env, para_env_RPA, homo, virtual, Eigenval, &
2101 : num_integ_points, num_integ_group, color_rpa_group, &
2102 40 : grid%frequency, grid%frequency_weights, fm_mat_S)
2103 : END IF
2104 :
2105 128 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.5)') 'INTEG_INFO| Scaling parameter:', a_scaling
2106 :
2107 5678 : grid%frequency_weights(:) = grid%frequency_weights(:)*a_scaling
2108 5678 : grid%frequency(:) = a_scaling/TAN(grid%frequency(:))
2109 :
2110 128 : CALL timestop(handle)
2111 :
2112 128 : END SUBROUTINE get_clenshaw_grid
2113 :
2114 : ! **************************************************************************************************
2115 : !> \brief ...
2116 : !> \param a_scaling_ext ...
2117 : !> \param para_env ...
2118 : !> \param para_env_RPA ...
2119 : !> \param homo ...
2120 : !> \param virtual ...
2121 : !> \param Eigenval ...
2122 : !> \param num_integ_points ...
2123 : !> \param num_integ_group ...
2124 : !> \param color_rpa_group ...
2125 : !> \param tj_ext ...
2126 : !> \param wj_ext ...
2127 : !> \param fm_mat_S ...
2128 : ! **************************************************************************************************
2129 40 : SUBROUTINE calc_scaling_factor(a_scaling_ext, para_env, para_env_RPA, homo, virtual, Eigenval, &
2130 : num_integ_points, num_integ_group, color_rpa_group, &
2131 40 : tj_ext, wj_ext, fm_mat_S)
2132 : REAL(KIND=dp), INTENT(OUT) :: a_scaling_ext
2133 : TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_RPA
2134 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
2135 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval
2136 : INTEGER, INTENT(IN) :: num_integ_points, num_integ_group, &
2137 : color_rpa_group
2138 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
2139 : INTENT(IN) :: tj_ext, wj_ext
2140 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
2141 :
2142 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_scaling_factor'
2143 :
2144 : INTEGER :: handle, icycle, jquad, ncol_local, &
2145 : ncol_local_beta, nspins
2146 : LOGICAL :: my_open_shell
2147 : REAL(KIND=dp) :: a_high, a_low, a_scaling, conv_param, eps, first_deriv, left_term, &
2148 : right_term, right_term_ref, right_term_ref_beta, step
2149 40 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cottj, D_ia, D_ia_beta, iaia_RI, &
2150 40 : iaia_RI_beta, M_ia, M_ia_beta
2151 : TYPE(mp_para_env_type), POINTER :: para_env_col, para_env_col_beta
2152 :
2153 40 : CALL timeset(routineN, handle)
2154 :
2155 40 : nspins = SIZE(homo)
2156 40 : my_open_shell = (nspins == 2)
2157 :
2158 40 : eps = 1.0E-10_dp
2159 :
2160 120 : ALLOCATE (cottj(num_integ_points))
2161 :
2162 : ! calculate the cotangent of the abscissa tj
2163 570 : DO jquad = 1, num_integ_points
2164 570 : cottj(jquad) = 1.0_dp/TAN(tj_ext(jquad))
2165 : END DO
2166 :
2167 : CALL calc_ia_ia_integrals(para_env_RPA, homo(1), virtual(1), ncol_local, right_term_ref, Eigenval(:, 1, 1), &
2168 40 : D_ia, iaia_RI, M_ia, fm_mat_S(1), para_env_col)
2169 :
2170 : ! In the open shell case do point 1-2-3 for the beta spin
2171 40 : IF (my_open_shell) THEN
2172 : CALL calc_ia_ia_integrals(para_env_RPA, homo(2), virtual(2), ncol_local_beta, right_term_ref_beta, Eigenval(:, 1, 2), &
2173 8 : D_ia_beta, iaia_RI_beta, M_ia_beta, fm_mat_S(2), para_env_col_beta)
2174 :
2175 8 : right_term_ref = right_term_ref + right_term_ref_beta
2176 : END IF
2177 :
2178 : ! bcast the result
2179 40 : IF (para_env%mepos == 0) THEN
2180 20 : CALL para_env%bcast(right_term_ref, 0)
2181 : ELSE
2182 20 : right_term_ref = 0.0_dp
2183 20 : CALL para_env%bcast(right_term_ref, 0)
2184 : END IF
2185 :
2186 : ! 5) start iteration for solving the non-linear equation by bisection
2187 : ! find limit, here step=0.5 seems a good compromise
2188 40 : conv_param = 100.0_dp*EPSILON(right_term_ref)
2189 40 : step = 0.5_dp
2190 40 : a_low = 0.0_dp
2191 40 : a_high = step
2192 40 : right_term = -right_term_ref
2193 104 : DO icycle = 1, num_integ_points*2
2194 96 : a_scaling = a_high
2195 :
2196 : CALL calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
2197 : M_ia, cottj, wj_ext, D_ia, D_ia_beta, M_ia_beta, &
2198 : ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
2199 96 : para_env, para_env_col, para_env_col_beta)
2200 96 : left_term = left_term/4.0_dp/pi*a_scaling
2201 :
2202 96 : IF (ABS(left_term) > ABS(right_term) .OR. ABS(left_term + right_term) <= conv_param) EXIT
2203 64 : a_low = a_high
2204 104 : a_high = a_high + step
2205 :
2206 : END DO
2207 :
2208 40 : IF (ABS(left_term + right_term) >= conv_param) THEN
2209 32 : IF (a_scaling >= 2*num_integ_points*step) THEN
2210 10 : a_scaling = 1.0_dp
2211 : ELSE
2212 :
2213 340 : DO icycle = 1, num_integ_points*2
2214 336 : a_scaling = (a_low + a_high)/2.0_dp
2215 :
2216 : CALL calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
2217 : M_ia, cottj, wj_ext, D_ia, D_ia_beta, M_ia_beta, &
2218 : ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
2219 336 : para_env, para_env_col, para_env_col_beta)
2220 336 : left_term = left_term/4.0_dp/pi*a_scaling
2221 :
2222 336 : IF (ABS(left_term) > ABS(right_term)) THEN
2223 : a_high = a_scaling
2224 : ELSE
2225 186 : a_low = a_scaling
2226 : END IF
2227 :
2228 340 : IF (ABS(a_high - a_low) < 1.0e-5_dp) EXIT
2229 :
2230 : END DO
2231 :
2232 : END IF
2233 : END IF
2234 :
2235 40 : a_scaling_ext = a_scaling
2236 40 : CALL para_env%bcast(a_scaling_ext, 0)
2237 :
2238 40 : DEALLOCATE (cottj)
2239 40 : DEALLOCATE (iaia_RI)
2240 40 : DEALLOCATE (D_ia)
2241 40 : DEALLOCATE (M_ia)
2242 40 : CALL mp_para_env_release(para_env_col)
2243 :
2244 40 : IF (my_open_shell) THEN
2245 8 : DEALLOCATE (iaia_RI_beta)
2246 8 : DEALLOCATE (D_ia_beta)
2247 8 : DEALLOCATE (M_ia_beta)
2248 8 : CALL mp_para_env_release(para_env_col_beta)
2249 : END IF
2250 :
2251 40 : CALL timestop(handle)
2252 :
2253 80 : END SUBROUTINE calc_scaling_factor
2254 :
2255 : ! **************************************************************************************************
2256 : !> \brief ...
2257 : !> \param para_env_RPA ...
2258 : !> \param homo ...
2259 : !> \param virtual ...
2260 : !> \param ncol_local ...
2261 : !> \param right_term_ref ...
2262 : !> \param Eigenval ...
2263 : !> \param D_ia ...
2264 : !> \param iaia_RI ...
2265 : !> \param M_ia ...
2266 : !> \param fm_mat_S ...
2267 : !> \param para_env_col ...
2268 : ! **************************************************************************************************
2269 48 : SUBROUTINE calc_ia_ia_integrals(para_env_RPA, homo, virtual, ncol_local, right_term_ref, Eigenval, &
2270 : D_ia, iaia_RI, M_ia, fm_mat_S, para_env_col)
2271 :
2272 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA
2273 : INTEGER, INTENT(IN) :: homo, virtual
2274 : INTEGER, INTENT(OUT) :: ncol_local
2275 : REAL(KIND=dp), INTENT(OUT) :: right_term_ref
2276 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
2277 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
2278 : INTENT(OUT) :: D_ia, iaia_RI, M_ia
2279 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
2280 : TYPE(mp_para_env_type), POINTER :: para_env_col
2281 :
2282 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_ia_ia_integrals'
2283 :
2284 : INTEGER :: avirt, color_col, color_row, handle, &
2285 : i_global, iiB, iocc, nrow_local
2286 48 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2287 : REAL(KIND=dp) :: eigen_diff
2288 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: iaia_RI_dp
2289 : TYPE(mp_para_env_type), POINTER :: para_env_row
2290 :
2291 48 : CALL timeset(routineN, handle)
2292 :
2293 : ! calculate the (ia|ia) RI integrals
2294 : ! ----------------------------------
2295 : ! 1) get info fm_mat_S
2296 : CALL cp_fm_get_info(matrix=fm_mat_S, &
2297 : nrow_local=nrow_local, &
2298 : ncol_local=ncol_local, &
2299 : row_indices=row_indices, &
2300 48 : col_indices=col_indices)
2301 :
2302 : ! allocate the local buffer of iaia_RI integrals (dp kind)
2303 142 : ALLOCATE (iaia_RI_dp(ncol_local))
2304 48 : iaia_RI_dp = 0.0_dp
2305 :
2306 : ! 2) perform the local multiplication SUM_K (ia|K)*(ia|K)
2307 2918 : DO iiB = 1, ncol_local
2308 200776 : iaia_RI_dp(iiB) = iaia_RI_dp(iiB) + DOT_PRODUCT(fm_mat_S%local_data(:, iiB), fm_mat_S%local_data(:, iiB))
2309 : END DO
2310 :
2311 : ! 3) sum the result with the processes of the RPA_group having the same columns
2312 : ! _______ia______ _
2313 : ! | | | | | | |
2314 : ! --> | 1 | 5 | 9 | 13| SUM --> | |
2315 : ! |___|__ |___|___| |_|
2316 : ! | | | | | | |
2317 : ! --> | 2 | 6 | 10| 14| SUM --> | |
2318 : ! K |___|___|___|___| |_| (ia|ia)_RI
2319 : ! | | | | | | |
2320 : ! --> | 3 | 7 | 11| 15| SUM --> | |
2321 : ! |___|___|___|___| |_|
2322 : ! | | | | | | |
2323 : ! --> | 4 | 8 | 12| 16| SUM --> | |
2324 : ! |___|___|___|___| |_|
2325 : !
2326 :
2327 48 : color_col = fm_mat_S%matrix_struct%context%mepos(2)
2328 48 : ALLOCATE (para_env_col)
2329 48 : CALL para_env_col%from_split(para_env_RPA, color_col)
2330 :
2331 48 : CALL para_env_col%sum(iaia_RI_dp)
2332 :
2333 : ! convert the iaia_RI_dp into double-double precision
2334 142 : ALLOCATE (iaia_RI(ncol_local))
2335 2918 : DO iiB = 1, ncol_local
2336 2918 : iaia_RI(iiB) = iaia_RI_dp(iiB)
2337 : END DO
2338 48 : DEALLOCATE (iaia_RI_dp)
2339 :
2340 : ! 4) calculate the right hand term, D_ia is the matrix containing the
2341 : ! orbital energy differences, M_ia is the diagonal of the full RPA 'excitation'
2342 : ! matrix
2343 142 : ALLOCATE (D_ia(ncol_local))
2344 :
2345 94 : ALLOCATE (M_ia(ncol_local))
2346 :
2347 2918 : DO iiB = 1, ncol_local
2348 2870 : i_global = col_indices(iiB)
2349 :
2350 2870 : iocc = MAX(1, i_global - 1)/virtual + 1
2351 2870 : avirt = i_global - (iocc - 1)*virtual
2352 2870 : eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
2353 :
2354 2918 : D_ia(iiB) = eigen_diff
2355 : END DO
2356 :
2357 2918 : DO iiB = 1, ncol_local
2358 2918 : M_ia(iiB) = D_ia(iiB)*D_ia(iiB) + 2.0_dp*D_ia(iiB)*iaia_RI(iiB)
2359 : END DO
2360 :
2361 48 : right_term_ref = 0.0_dp
2362 2918 : DO iiB = 1, ncol_local
2363 2918 : right_term_ref = right_term_ref + (SQRT(M_ia(iiB)) - D_ia(iiB) - iaia_RI(iiB))
2364 : END DO
2365 48 : right_term_ref = right_term_ref/2.0_dp
2366 :
2367 : ! sum the result with the processes of the RPA_group having the same row
2368 48 : color_row = fm_mat_S%matrix_struct%context%mepos(1)
2369 48 : ALLOCATE (para_env_row)
2370 48 : CALL para_env_row%from_split(para_env_RPA, color_row)
2371 :
2372 : ! allocate communication array for rows
2373 48 : CALL para_env_row%sum(right_term_ref)
2374 :
2375 48 : CALL mp_para_env_release(para_env_row)
2376 :
2377 48 : CALL timestop(handle)
2378 :
2379 48 : END SUBROUTINE calc_ia_ia_integrals
2380 :
2381 : ! **************************************************************************************************
2382 : !> \brief ...
2383 : !> \param a_scaling ...
2384 : !> \param left_term ...
2385 : !> \param first_deriv ...
2386 : !> \param num_integ_points ...
2387 : !> \param my_open_shell ...
2388 : !> \param M_ia ...
2389 : !> \param cottj ...
2390 : !> \param wj ...
2391 : !> \param D_ia ...
2392 : !> \param D_ia_beta ...
2393 : !> \param M_ia_beta ...
2394 : !> \param ncol_local ...
2395 : !> \param ncol_local_beta ...
2396 : !> \param num_integ_group ...
2397 : !> \param color_rpa_group ...
2398 : !> \param para_env ...
2399 : !> \param para_env_col ...
2400 : !> \param para_env_col_beta ...
2401 : ! **************************************************************************************************
2402 432 : SUBROUTINE calculate_objfunc(a_scaling, left_term, first_deriv, num_integ_points, my_open_shell, &
2403 : M_ia, cottj, wj, D_ia, D_ia_beta, M_ia_beta, &
2404 : ncol_local, ncol_local_beta, num_integ_group, color_rpa_group, &
2405 : para_env, para_env_col, para_env_col_beta)
2406 : REAL(KIND=dp), INTENT(IN) :: a_scaling
2407 : REAL(KIND=dp), INTENT(INOUT) :: left_term, first_deriv
2408 : INTEGER, INTENT(IN) :: num_integ_points
2409 : LOGICAL, INTENT(IN) :: my_open_shell
2410 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
2411 : INTENT(IN) :: M_ia, cottj, wj, D_ia, D_ia_beta, &
2412 : M_ia_beta
2413 : INTEGER, INTENT(IN) :: ncol_local, ncol_local_beta, &
2414 : num_integ_group, color_rpa_group
2415 : TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_col
2416 : TYPE(mp_para_env_type), POINTER :: para_env_col_beta
2417 :
2418 : INTEGER :: iiB, jquad
2419 : REAL(KIND=dp) :: first_deriv_beta, left_term_beta, omega
2420 :
2421 432 : left_term = 0.0_dp
2422 432 : first_deriv = 0.0_dp
2423 432 : left_term_beta = 0.0_dp
2424 432 : first_deriv_beta = 0.0_dp
2425 4616 : DO jquad = 1, num_integ_points
2426 : ! parallelize over integration points
2427 4184 : IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE
2428 2292 : omega = a_scaling*cottj(jquad)
2429 :
2430 168644 : DO iiB = 1, ncol_local
2431 : ! parallelize over ia elements in the para_env_row group
2432 166352 : IF (MODULO(iiB, para_env_col%num_pe) /= para_env_col%mepos) CYCLE
2433 : ! calculate left_term
2434 : left_term = left_term + wj(jquad)* &
2435 : (LOG(1.0_dp + (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2)) - &
2436 151152 : (M_ia(iiB) - D_ia(iiB)**2)/(omega**2 + D_ia(iiB)**2))
2437 : first_deriv = first_deriv + wj(jquad)*cottj(jquad)**2* &
2438 168644 : ((-M_ia(iiB) + D_ia(iiB)**2)**2/((omega**2 + D_ia(iiB)**2)**2*(omega**2 + M_ia(iiB))))
2439 : END DO
2440 :
2441 2724 : IF (my_open_shell) THEN
2442 14490 : DO iiB = 1, ncol_local_beta
2443 : ! parallelize over ia elements in the para_env_row group
2444 14140 : IF (MODULO(iiB, para_env_col_beta%num_pe) /= para_env_col_beta%mepos) CYCLE
2445 : ! calculate left_term
2446 : left_term_beta = left_term_beta + wj(jquad)* &
2447 : (LOG(1.0_dp + (M_ia_beta(iiB) - D_ia_beta(iiB)**2)/(omega**2 + D_ia_beta(iiB)**2)) - &
2448 14140 : (M_ia_beta(iiB) - D_ia_beta(iiB)**2)/(omega**2 + D_ia_beta(iiB)**2))
2449 : first_deriv_beta = &
2450 : first_deriv_beta + wj(jquad)*cottj(jquad)**2* &
2451 14490 : ((-M_ia_beta(iiB) + D_ia_beta(iiB)**2)**2/((omega**2 + D_ia_beta(iiB)**2)**2*(omega**2 + M_ia_beta(iiB))))
2452 : END DO
2453 : END IF
2454 :
2455 : END DO
2456 :
2457 : ! sum the contribution from all proc, starting form the row group
2458 432 : CALL para_env%sum(left_term)
2459 432 : CALL para_env%sum(first_deriv)
2460 :
2461 432 : IF (my_open_shell) THEN
2462 70 : CALL para_env%sum(left_term_beta)
2463 70 : CALL para_env%sum(first_deriv_beta)
2464 :
2465 70 : left_term = left_term + left_term_beta
2466 70 : first_deriv = first_deriv + first_deriv_beta
2467 : END IF
2468 :
2469 432 : END SUBROUTINE calculate_objfunc
2470 :
2471 : ! **************************************************************************************************
2472 : !> \brief ...
2473 : !> \param qs_env ...
2474 : !> \param para_env ...
2475 : !> \param gap ...
2476 : !> \param max_eig_diff ...
2477 : !> \param e_fermi ...
2478 : ! **************************************************************************************************
2479 12 : SUBROUTINE gap_and_max_eig_diff_kpoints(qs_env, para_env, gap, max_eig_diff, e_fermi)
2480 :
2481 : TYPE(qs_environment_type), POINTER :: qs_env
2482 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2483 : REAL(KIND=dp), INTENT(OUT) :: gap, max_eig_diff, e_fermi
2484 :
2485 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gap_and_max_eig_diff_kpoints'
2486 :
2487 : INTEGER :: handle, homo, ikpgr, ispin, kplocal, &
2488 : nmo, nspin
2489 : INTEGER, DIMENSION(2) :: kp_range
2490 : REAL(KIND=dp) :: e_homo, e_homo_temp, e_lumo, e_lumo_temp
2491 : REAL(KIND=dp), DIMENSION(3) :: tmp
2492 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
2493 : TYPE(kpoint_env_type), POINTER :: kp
2494 : TYPE(kpoint_type), POINTER :: kpoint
2495 : TYPE(mo_set_type), POINTER :: mo_set
2496 :
2497 6 : CALL timeset(routineN, handle)
2498 :
2499 : CALL get_qs_env(qs_env, &
2500 6 : kpoints=kpoint)
2501 :
2502 6 : mo_set => kpoint%kp_env(1)%kpoint_env%mos(1)
2503 6 : CALL get_mo_set(mo_set, nmo=nmo)
2504 :
2505 6 : CALL get_kpoint_info(kpoint, kp_range=kp_range)
2506 6 : kplocal = kp_range(2) - kp_range(1) + 1
2507 :
2508 6 : gap = 1000.0_dp
2509 6 : max_eig_diff = 0.0_dp
2510 6 : e_homo = -1000.0_dp
2511 6 : e_lumo = 1000.0_dp
2512 :
2513 18 : DO ikpgr = 1, kplocal
2514 12 : kp => kpoint%kp_env(ikpgr)%kpoint_env
2515 12 : nspin = SIZE(kp%mos)
2516 30 : DO ispin = 1, nspin
2517 12 : mo_set => kp%mos(ispin)
2518 12 : CALL get_mo_set(mo_set, eigenvalues=eigenvalues, homo=homo)
2519 12 : e_homo_temp = eigenvalues(homo)
2520 12 : e_lumo_temp = eigenvalues(homo + 1)
2521 :
2522 : IF (e_homo_temp > e_homo) e_homo = e_homo_temp
2523 : IF (e_lumo_temp < e_lumo) e_lumo = e_lumo_temp
2524 24 : IF (eigenvalues(nmo) - eigenvalues(1) > max_eig_diff) max_eig_diff = eigenvalues(nmo) - eigenvalues(1)
2525 :
2526 : END DO
2527 : END DO
2528 :
2529 : ! Collect all three numbers in an array
2530 : ! Reverse sign of lumo to reduce number of MPI calls
2531 6 : tmp(1) = e_homo
2532 6 : tmp(2) = -e_lumo
2533 6 : tmp(3) = max_eig_diff
2534 6 : CALL para_env%max(tmp)
2535 :
2536 6 : gap = -tmp(2) - tmp(1)
2537 6 : e_fermi = (tmp(1) - tmp(2))*0.5_dp
2538 6 : max_eig_diff = tmp(3)
2539 :
2540 6 : CALL timestop(handle)
2541 :
2542 6 : END SUBROUTINE gap_and_max_eig_diff_kpoints
2543 :
2544 : ! **************************************************************************************************
2545 : !> \brief returns minimal and maximal energy values for the E_range for the minimax grid selection
2546 : !> \param qs_env ...
2547 : !> \param para_env ...
2548 : !> \param homo index of the homo level for the respective spin channel
2549 : !> \param Eigenval eigenvalues
2550 : !> \param do_ri_sos_laplace_mp2 flag for SOS-MP2
2551 : !> \param do_kpoints_cubic_RPA flag for cubic-scaling RPA with k-points
2552 : !> \param Emin minimal eigenvalue difference (gap of the system)
2553 : !> \param Emax maximal eigenvalue difference
2554 : !> \param e_range ...
2555 : !> \param e_fermi Fermi level
2556 : ! **************************************************************************************************
2557 214 : SUBROUTINE determine_energy_range(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, &
2558 : do_kpoints_cubic_RPA, Emin, Emax, e_range, e_fermi)
2559 :
2560 : TYPE(qs_environment_type), POINTER :: qs_env
2561 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
2562 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
2563 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval
2564 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, &
2565 : do_kpoints_cubic_RPA
2566 : REAL(KIND=dp), INTENT(OUT) :: Emin, Emax, e_range, e_fermi
2567 :
2568 : CHARACTER(LEN=*), PARAMETER :: routineN = 'determine_energy_range'
2569 :
2570 : INTEGER :: handle, ispin, nspins
2571 : LOGICAL :: my_do_kpoints
2572 : TYPE(section_vals_type), POINTER :: input
2573 :
2574 214 : CALL timeset(routineN, handle)
2575 : ! Test for spin unrestricted
2576 214 : nspins = SIZE(homo)
2577 :
2578 : ! Test whether all necessary variables are available
2579 214 : my_do_kpoints = .FALSE.
2580 214 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
2581 150 : my_do_kpoints = do_kpoints_cubic_RPA
2582 : END IF
2583 :
2584 150 : IF (my_do_kpoints) THEN
2585 6 : CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi)
2586 6 : E_Range = Emax/Emin
2587 : ELSE
2588 208 : IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN
2589 160 : Emin = HUGE(dp)
2590 160 : Emax = 0.0_dp
2591 360 : DO ispin = 1, nspins
2592 360 : IF (homo(ispin) > 0) THEN
2593 196 : Emin = MIN(Emin, Eigenval(homo(ispin) + 1, 1, ispin) - Eigenval(homo(ispin), 1, ispin))
2594 15360 : Emax = MAX(Emax, MAXVAL(Eigenval(:, :, ispin)) - MINVAL(Eigenval(:, :, ispin)))
2595 : END IF
2596 : END DO
2597 160 : E_Range = Emax/Emin
2598 160 : qs_env%mp2_env%e_range = e_range
2599 160 : qs_env%mp2_env%e_gap = Emin
2600 :
2601 160 : CALL get_qs_env(qs_env, input=input)
2602 160 : CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=e_range)
2603 160 : CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_GAP", r_val=emin)
2604 : ELSE
2605 48 : E_range = qs_env%mp2_env%E_range
2606 48 : Emin = qs_env%mp2_env%E_gap
2607 48 : Emax = Emin*E_range
2608 : END IF
2609 : END IF
2610 :
2611 : ! When we perform SOS-MP2, we need an additional factor of 2 for the energies (compare with mp2_laplace.F)
2612 : ! We do not need weights etc. for the cosine transform
2613 : ! We do not scale Emax because it is not needed for SOS-MP2
2614 214 : IF (do_ri_sos_laplace_mp2) THEN
2615 64 : Emin = Emin*2.0_dp
2616 64 : Emax = Emax*2.0_dp
2617 : END IF
2618 :
2619 214 : CALL timestop(handle)
2620 214 : END SUBROUTINE determine_energy_range
2621 :
2622 : END MODULE rpa_main
|
