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