Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Utility functions for RPA calculations
10 : !> \par History
11 : !> 06.2019 Moved from rpa_ri_gpw.F [Frederick Stein]
12 : ! **************************************************************************************************
13 : MODULE rpa_util
14 :
15 : USE cell_types, ONLY: cell_type,&
16 : get_cell
17 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
18 : cp_blacs_env_release,&
19 : cp_blacs_env_type
20 : USE cp_cfm_types, ONLY: cp_cfm_create,&
21 : cp_cfm_release,&
22 : cp_cfm_set_all,&
23 : cp_cfm_type
24 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
25 : copy_fm_to_dbcsr,&
26 : dbcsr_allocate_matrix_set,&
27 : dbcsr_deallocate_matrix_set
28 : USE cp_fm_basic_linalg, ONLY: cp_fm_syrk,&
29 : cp_fm_transpose
30 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
31 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32 : cp_fm_struct_release,&
33 : cp_fm_struct_type
34 : USE cp_fm_types, ONLY: &
35 : cp_fm_copy_general, cp_fm_create, cp_fm_get_info, cp_fm_release, cp_fm_set_all, &
36 : cp_fm_set_element, cp_fm_to_fm, cp_fm_to_fm_submat_general, cp_fm_type
37 : USE dbcsr_api, ONLY: &
38 : dbcsr_create, dbcsr_deallocate_matrix, dbcsr_filter, dbcsr_get_info, dbcsr_multiply, &
39 : dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
40 : USE dbt_api, ONLY: dbt_destroy,&
41 : dbt_type
42 : USE dgemm_counter_types, ONLY: dgemm_counter_start,&
43 : dgemm_counter_stop,&
44 : dgemm_counter_type
45 : USE hfx_types, ONLY: block_ind_type,&
46 : dealloc_containers,&
47 : hfx_compression_type
48 : USE input_constants, ONLY: wfc_mm_style_gemm,&
49 : wfc_mm_style_syrk
50 : USE kinds, ONLY: dp
51 : USE kpoint_types, ONLY: get_kpoint_info,&
52 : kpoint_release,&
53 : kpoint_type
54 : USE mathconstants, ONLY: z_zero
55 : USE message_passing, ONLY: mp_para_env_release,&
56 : mp_para_env_type
57 : USE mp2_laplace, ONLY: calc_fm_mat_S_laplace
58 : USE parallel_gemm_api, ONLY: parallel_gemm
59 : USE qs_environment_types, ONLY: get_qs_env,&
60 : qs_environment_type
61 : USE rpa_gw_kpoints_util, ONLY: compute_wkp_W
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_util'
69 :
70 : PUBLIC :: compute_Erpa_by_freq_int, alloc_im_time, calc_mat_Q, Q_trace_and_add_unit_matrix, &
71 : dealloc_im_time, contract_P_omega_with_mat_L, calc_fm_mat_S_rpa, remove_scaling_factor_rpa
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief ...
77 : !> \param qs_env ...
78 : !> \param para_env ...
79 : !> \param dimen_RI ...
80 : !> \param dimen_RI_red ...
81 : !> \param num_integ_points ...
82 : !> \param nspins ...
83 : !> \param fm_mat_Q ...
84 : !> \param fm_mo_coeff_occ ...
85 : !> \param fm_mo_coeff_virt ...
86 : !> \param fm_matrix_Minv_L_kpoints ...
87 : !> \param fm_matrix_L_kpoints ...
88 : !> \param mat_P_global ...
89 : !> \param t_3c_O ...
90 : !> \param matrix_s ...
91 : !> \param kpoints ...
92 : !> \param eps_filter_im_time ...
93 : !> \param cut_memory ...
94 : !> \param nkp ...
95 : !> \param num_cells_dm ...
96 : !> \param num_3c_repl ...
97 : !> \param size_P ...
98 : !> \param ikp_local ...
99 : !> \param index_to_cell_3c ...
100 : !> \param cell_to_index_3c ...
101 : !> \param col_blk_size ...
102 : !> \param do_ic_model ...
103 : !> \param do_kpoints_cubic_RPA ...
104 : !> \param do_kpoints_from_Gamma ...
105 : !> \param do_ri_Sigma_x ...
106 : !> \param my_open_shell ...
107 : !> \param has_mat_P_blocks ...
108 : !> \param wkp_W ...
109 : !> \param cfm_mat_Q ...
110 : !> \param fm_mat_Minv_L_kpoints ...
111 : !> \param fm_mat_L_kpoints ...
112 : !> \param fm_mat_RI_global_work ...
113 : !> \param fm_mat_work ...
114 : !> \param fm_mo_coeff_occ_scaled ...
115 : !> \param fm_mo_coeff_virt_scaled ...
116 : !> \param mat_dm ...
117 : !> \param mat_L ...
118 : !> \param mat_M_P_munu_occ ...
119 : !> \param mat_M_P_munu_virt ...
120 : !> \param mat_MinvVMinv ...
121 : !> \param mat_P_omega ...
122 : !> \param mat_P_omega_kp ...
123 : !> \param mat_work ...
124 : !> \param mo_coeff ...
125 : !> \param fm_scaled_dm_occ_tau ...
126 : !> \param fm_scaled_dm_virt_tau ...
127 : !> \param homo ...
128 : !> \param nmo ...
129 : ! **************************************************************************************************
130 268 : SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_points, nspins, &
131 : fm_mat_Q, fm_mo_coeff_occ, fm_mo_coeff_virt, &
132 : fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, mat_P_global, &
133 : t_3c_O, matrix_s, kpoints, eps_filter_im_time, &
134 : cut_memory, nkp, num_cells_dm, num_3c_repl, &
135 : size_P, ikp_local, &
136 : index_to_cell_3c, &
137 : cell_to_index_3c, &
138 : col_blk_size, &
139 : do_ic_model, do_kpoints_cubic_RPA, &
140 : do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, &
141 : has_mat_P_blocks, wkp_W, &
142 : cfm_mat_Q, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
143 : fm_mat_RI_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, &
144 : fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, &
145 : mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
146 134 : mat_work, mo_coeff, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, homo, nmo)
147 :
148 : TYPE(qs_environment_type), POINTER :: qs_env
149 : TYPE(mp_para_env_type), POINTER :: para_env
150 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, &
151 : num_integ_points, nspins
152 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
153 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mo_coeff_occ, fm_mo_coeff_virt
154 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv_L_kpoints, &
155 : fm_matrix_L_kpoints
156 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_P_global
157 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
158 : INTENT(INOUT) :: t_3c_O
159 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
160 : TYPE(kpoint_type), POINTER :: kpoints
161 : REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time
162 : INTEGER, INTENT(IN) :: cut_memory
163 : INTEGER, INTENT(OUT) :: nkp, num_cells_dm, num_3c_repl, size_P, &
164 : ikp_local
165 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell_3c
166 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
167 : INTENT(OUT) :: cell_to_index_3c
168 : INTEGER, DIMENSION(:), POINTER :: col_blk_size
169 : LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA, &
170 : do_kpoints_from_Gamma, do_ri_Sigma_x, &
171 : my_open_shell
172 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), &
173 : INTENT(OUT) :: has_mat_P_blocks
174 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
175 : INTENT(OUT) :: wkp_W
176 : TYPE(cp_cfm_type), INTENT(OUT) :: cfm_mat_Q
177 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_Minv_L_kpoints, fm_mat_L_kpoints
178 : TYPE(cp_fm_type), INTENT(OUT) :: fm_mat_RI_global_work, fm_mat_work, &
179 : fm_mo_coeff_occ_scaled, &
180 : fm_mo_coeff_virt_scaled
181 : TYPE(dbcsr_p_type), INTENT(OUT) :: mat_dm, mat_L, mat_M_P_munu_occ, &
182 : mat_M_P_munu_virt, mat_MinvVMinv
183 : TYPE(dbcsr_p_type), ALLOCATABLE, &
184 : DIMENSION(:, :, :) :: mat_P_omega
185 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_kp
186 : TYPE(dbcsr_type), POINTER :: mat_work
187 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
188 : TYPE(cp_fm_type), INTENT(OUT) :: fm_scaled_dm_occ_tau, &
189 : fm_scaled_dm_virt_tau
190 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
191 : INTEGER, INTENT(IN) :: nmo
192 :
193 : CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_im_time'
194 :
195 : INTEGER :: cell_grid_dm(3), first_ikp_local, &
196 : handle, i_dim, i_kp, ispin, jquad, &
197 : nspins_P_omega, periodic(3)
198 134 : INTEGER, DIMENSION(:), POINTER :: row_blk_size
199 134 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wkp_V
200 : TYPE(cell_type), POINTER :: cell
201 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_sub_kp
202 :
203 134 : CALL timeset(routineN, handle)
204 :
205 998 : ALLOCATE (fm_mo_coeff_occ(nspins), fm_mo_coeff_virt(nspins))
206 :
207 134 : CALL cp_fm_create(fm_scaled_dm_occ_tau, mo_coeff(1)%matrix_struct)
208 134 : CALL cp_fm_set_all(fm_scaled_dm_occ_tau, 0.0_dp)
209 :
210 134 : CALL cp_fm_create(fm_scaled_dm_virt_tau, mo_coeff(1)%matrix_struct)
211 134 : CALL cp_fm_set_all(fm_scaled_dm_virt_tau, 0.0_dp)
212 :
213 298 : DO ispin = 1, SIZE(mo_coeff)
214 : CALL create_occ_virt_mo_coeffs(fm_mo_coeff_occ(ispin), fm_mo_coeff_virt(ispin), mo_coeff(ispin), &
215 298 : nmo, homo(ispin))
216 : END DO
217 :
218 134 : num_3c_repl = SIZE(t_3c_O, 2)
219 :
220 134 : IF (do_kpoints_cubic_RPA) THEN
221 : ! we always use an odd number of image cells
222 : ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical
223 16 : DO i_dim = 1, 3
224 16 : cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
225 : END DO
226 4 : num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
227 12 : ALLOCATE (index_to_cell_3c(3, SIZE(kpoints%index_to_cell, 2)))
228 4 : CPASSERT(SIZE(kpoints%index_to_cell, 1) == 3)
229 84 : index_to_cell_3c(:, :) = kpoints%index_to_cell(:, :)
230 0 : ALLOCATE (cell_to_index_3c(LBOUND(kpoints%cell_to_index, 1):UBOUND(kpoints%cell_to_index, 1), &
231 : LBOUND(kpoints%cell_to_index, 2):UBOUND(kpoints%cell_to_index, 2), &
232 20 : LBOUND(kpoints%cell_to_index, 3):UBOUND(kpoints%cell_to_index, 3)))
233 64 : cell_to_index_3c(:, :, :) = kpoints%cell_to_index(:, :, :)
234 :
235 : ELSE
236 130 : ALLOCATE (index_to_cell_3c(3, 1))
237 520 : index_to_cell_3c(:, 1) = 0
238 130 : ALLOCATE (cell_to_index_3c(0:0, 0:0, 0:0))
239 130 : cell_to_index_3c(0, 0, 0) = 1
240 130 : num_cells_dm = 1
241 : END IF
242 :
243 134 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
244 :
245 22 : CALL get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_RI, ikp_local, first_ikp_local)
246 :
247 22 : CALL cp_cfm_create(cfm_mat_Q, fm_struct_sub_kp)
248 22 : CALL cp_cfm_set_all(cfm_mat_Q, z_zero)
249 : ELSE
250 112 : first_ikp_local = 1
251 : END IF
252 :
253 : ! if we do kpoints, mat_P has a kpoint and mat_P_omega has the inted
254 : ! mat_P(tau, kpoint)
255 134 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
256 :
257 22 : NULLIFY (cell)
258 22 : CALL get_qs_env(qs_env, cell=cell)
259 22 : CALL get_cell(cell=cell, periodic=periodic)
260 :
261 22 : CALL get_kpoint_info(kpoints, nkp=nkp)
262 : ! compute k-point weights such that functions 1/k^2, 1/k and const function are
263 : ! integrated correctly
264 22 : CALL compute_wkp_W(qs_env, wkp_W, wkp_V, kpoints, cell%h_inv, periodic)
265 22 : DEALLOCATE (wkp_V)
266 :
267 : ELSE
268 112 : nkp = 1
269 : END IF
270 :
271 134 : IF (do_kpoints_cubic_RPA) THEN
272 4 : size_P = MAX(num_cells_dm/2 + 1, nkp)
273 130 : ELSE IF (do_kpoints_from_Gamma) THEN
274 18 : size_P = MAX(3**(periodic(1) + periodic(2) + periodic(3)), nkp)
275 : ELSE
276 112 : size_P = 1
277 : END IF
278 :
279 134 : nspins_P_omega = 1
280 134 : IF (my_open_shell) nspins_P_omega = 2
281 :
282 5920 : ALLOCATE (mat_P_omega(num_integ_points, size_P, nspins_P_omega))
283 298 : DO ispin = 1, nspins_P_omega
284 1016 : DO i_kp = 1, size_P
285 5250 : DO jquad = 1, num_integ_points
286 4368 : NULLIFY (mat_P_omega(jquad, i_kp, ispin)%matrix)
287 4368 : ALLOCATE (mat_P_omega(jquad, i_kp, ispin)%matrix)
288 : CALL dbcsr_create(matrix=mat_P_omega(jquad, i_kp, ispin)%matrix, &
289 4368 : template=mat_P_global%matrix)
290 5086 : CALL dbcsr_set(mat_P_omega(jquad, i_kp, ispin)%matrix, 0.0_dp)
291 : END DO
292 : END DO
293 : END DO
294 :
295 134 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
296 22 : CALL alloc_mat_P_omega(mat_P_omega_kp, 2, size_P, mat_P_global%matrix)
297 : END IF
298 :
299 134 : CALL cp_fm_create(fm_mo_coeff_occ_scaled, fm_mo_coeff_occ(1)%matrix_struct)
300 134 : CALL cp_fm_to_fm(fm_mo_coeff_occ(1), fm_mo_coeff_occ_scaled)
301 134 : CALL cp_fm_set_all(matrix=fm_mo_coeff_occ_scaled, alpha=0.0_dp)
302 :
303 134 : CALL cp_fm_create(fm_mo_coeff_virt_scaled, fm_mo_coeff_virt(1)%matrix_struct)
304 134 : CALL cp_fm_to_fm(fm_mo_coeff_virt(1), fm_mo_coeff_virt_scaled)
305 134 : CALL cp_fm_set_all(matrix=fm_mo_coeff_virt_scaled, alpha=0.0_dp)
306 :
307 134 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
308 22 : CALL cp_fm_create(fm_mat_RI_global_work, fm_matrix_Minv_L_kpoints(1, 1)%matrix_struct)
309 22 : CALL cp_fm_to_fm(fm_matrix_Minv_L_kpoints(1, 1), fm_mat_RI_global_work)
310 22 : CALL cp_fm_set_all(fm_mat_RI_global_work, 0.0_dp)
311 : END IF
312 :
313 938 : ALLOCATE (has_mat_P_blocks(num_cells_dm/2 + 1, cut_memory, cut_memory, num_3c_repl, num_3c_repl))
314 3036 : has_mat_P_blocks = .TRUE.
315 :
316 134 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
317 : CALL reorder_mat_L(fm_mat_Minv_L_kpoints, fm_matrix_Minv_L_kpoints, fm_mat_Q%matrix_struct, para_env, mat_L, &
318 : mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp, &
319 22 : allocate_mat_L=.FALSE.)
320 :
321 : CALL reorder_mat_L(fm_mat_L_kpoints, fm_matrix_L_kpoints, fm_mat_Q%matrix_struct, para_env, mat_L, &
322 22 : mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp)
323 :
324 22 : CALL cp_fm_struct_release(fm_struct_sub_kp)
325 :
326 : ELSE
327 : CALL reorder_mat_L(fm_mat_Minv_L_kpoints, fm_matrix_Minv_L_kpoints, fm_mat_Q%matrix_struct, para_env, mat_L, &
328 112 : mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local)
329 : END IF
330 :
331 : ! Create Scalapack working matrix for the contraction with the metric
332 134 : IF (dimen_RI == dimen_RI_red) THEN
333 130 : CALL cp_fm_create(fm_mat_work, fm_mat_Q%matrix_struct)
334 130 : CALL cp_fm_to_fm(fm_mat_Q, fm_mat_work)
335 130 : CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp)
336 :
337 : ELSE
338 4 : NULLIFY (fm_struct)
339 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_mat_Q%matrix_struct, &
340 4 : ncol_global=dimen_RI_red, nrow_global=dimen_RI)
341 :
342 4 : CALL cp_fm_create(fm_mat_work, fm_struct)
343 4 : CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp)
344 :
345 4 : CALL cp_fm_struct_release(fm_struct)
346 :
347 : END IF
348 :
349 : ! Then its DBCSR counter part
350 134 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
351 112 : CALL dbcsr_get_info(mat_L%matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size)
352 :
353 : ! Create mat_work having the shape of the transposed of mat_L (compare with contract_P_omega_with_mat_L)
354 : NULLIFY (mat_work)
355 112 : ALLOCATE (mat_work)
356 112 : CALL dbcsr_create(mat_work, template=mat_L%matrix, row_blk_size=col_blk_size, col_blk_size=row_blk_size)
357 : END IF
358 :
359 134 : IF (do_ri_Sigma_x .OR. do_ic_model) THEN
360 :
361 : NULLIFY (mat_MinvVMinv%matrix)
362 106 : ALLOCATE (mat_MinvVMinv%matrix)
363 106 : CALL dbcsr_create(mat_MinvVMinv%matrix, template=mat_P_global%matrix)
364 106 : CALL dbcsr_set(mat_MinvVMinv%matrix, 0.0_dp)
365 :
366 : ! for kpoints we compute SinvVSinv later with kpoints
367 106 : IF (.NOT. do_kpoints_from_Gamma) THEN
368 :
369 : ! get the Coulomb matrix for Sigma_x = G*V
370 : CALL dbcsr_multiply("T", "N", 1.0_dp, mat_L%matrix, mat_L%matrix, &
371 92 : 0.0_dp, mat_MinvVMinv%matrix, filter_eps=eps_filter_im_time)
372 :
373 : END IF
374 :
375 : END IF
376 :
377 134 : IF (do_ri_Sigma_x) THEN
378 :
379 : NULLIFY (mat_dm%matrix)
380 106 : ALLOCATE (mat_dm%matrix)
381 106 : CALL dbcsr_create(mat_dm%matrix, template=matrix_s(1)%matrix)
382 :
383 : END IF
384 :
385 134 : CALL timestop(handle)
386 :
387 268 : END SUBROUTINE alloc_im_time
388 :
389 : ! **************************************************************************************************
390 : !> \brief ...
391 : !> \param fm_mo_coeff_occ ...
392 : !> \param fm_mo_coeff_virt ...
393 : !> \param mo_coeff ...
394 : !> \param nmo ...
395 : !> \param homo ...
396 : ! **************************************************************************************************
397 164 : SUBROUTINE create_occ_virt_mo_coeffs(fm_mo_coeff_occ, fm_mo_coeff_virt, mo_coeff, &
398 : nmo, homo)
399 :
400 : TYPE(cp_fm_type), INTENT(OUT) :: fm_mo_coeff_occ, fm_mo_coeff_virt
401 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
402 : INTEGER, INTENT(IN) :: nmo, homo
403 :
404 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_occ_virt_mo_coeffs'
405 :
406 : INTEGER :: handle, icol_global, irow_global
407 :
408 164 : CALL timeset(routineN, handle)
409 :
410 164 : CALL cp_fm_create(fm_mo_coeff_occ, mo_coeff%matrix_struct)
411 164 : CALL cp_fm_set_all(fm_mo_coeff_occ, 0.0_dp)
412 164 : CALL cp_fm_to_fm(mo_coeff, fm_mo_coeff_occ)
413 :
414 : ! set all virtual MO coeffs to zero
415 4176 : DO irow_global = 1, nmo
416 98924 : DO icol_global = homo + 1, nmo
417 98760 : CALL cp_fm_set_element(fm_mo_coeff_occ, irow_global, icol_global, 0.0_dp)
418 : END DO
419 : END DO
420 :
421 164 : CALL cp_fm_create(fm_mo_coeff_virt, mo_coeff%matrix_struct)
422 164 : CALL cp_fm_set_all(fm_mo_coeff_virt, 0.0_dp)
423 164 : CALL cp_fm_to_fm(mo_coeff, fm_mo_coeff_virt)
424 :
425 : ! set all occupied MO coeffs to zero
426 4176 : DO irow_global = 1, nmo
427 19260 : DO icol_global = 1, homo
428 19096 : CALL cp_fm_set_element(fm_mo_coeff_virt, irow_global, icol_global, 0.0_dp)
429 : END DO
430 : END DO
431 :
432 164 : CALL timestop(handle)
433 :
434 164 : END SUBROUTINE create_occ_virt_mo_coeffs
435 :
436 : ! **************************************************************************************************
437 : !> \brief ...
438 : !> \param mat_P_omega ...
439 : !> \param num_integ_points ...
440 : !> \param size_P ...
441 : !> \param template ...
442 : ! **************************************************************************************************
443 22 : SUBROUTINE alloc_mat_P_omega(mat_P_omega, num_integ_points, size_P, template)
444 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega
445 : INTEGER, INTENT(IN) :: num_integ_points, size_P
446 : TYPE(dbcsr_type), POINTER :: template
447 :
448 : CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_mat_P_omega'
449 :
450 : INTEGER :: handle, i_kp, jquad
451 :
452 22 : CALL timeset(routineN, handle)
453 :
454 22 : NULLIFY (mat_P_omega)
455 22 : CALL dbcsr_allocate_matrix_set(mat_P_omega, num_integ_points, size_P)
456 498 : DO i_kp = 1, size_P
457 1450 : DO jquad = 1, num_integ_points
458 952 : ALLOCATE (mat_P_omega(jquad, i_kp)%matrix)
459 : CALL dbcsr_create(matrix=mat_P_omega(jquad, i_kp)%matrix, &
460 952 : template=template)
461 1428 : CALL dbcsr_set(mat_P_omega(jquad, i_kp)%matrix, 0.0_dp)
462 : END DO
463 : END DO
464 :
465 22 : CALL timestop(handle)
466 :
467 22 : END SUBROUTINE alloc_mat_P_omega
468 :
469 : ! **************************************************************************************************
470 : !> \brief ...
471 : !> \param fm_mat_L ...
472 : !> \param fm_matrix_Minv_L_kpoints ...
473 : !> \param fm_struct_template ...
474 : !> \param para_env ...
475 : !> \param mat_L ...
476 : !> \param mat_template ...
477 : !> \param dimen_RI ...
478 : !> \param dimen_RI_red ...
479 : !> \param first_ikp_local ...
480 : !> \param ikp_local ...
481 : !> \param fm_struct_sub_kp ...
482 : !> \param allocate_mat_L ...
483 : ! **************************************************************************************************
484 156 : SUBROUTINE reorder_mat_L(fm_mat_L, fm_matrix_Minv_L_kpoints, fm_struct_template, para_env, mat_L, mat_template, &
485 : dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp, allocate_mat_L)
486 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_L, fm_matrix_Minv_L_kpoints
487 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_template
488 : TYPE(mp_para_env_type), POINTER :: para_env
489 : TYPE(dbcsr_p_type), INTENT(OUT) :: mat_L
490 : TYPE(dbcsr_type), INTENT(IN) :: mat_template
491 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red, first_ikp_local
492 : INTEGER, OPTIONAL :: ikp_local
493 : TYPE(cp_fm_struct_type), OPTIONAL, POINTER :: fm_struct_sub_kp
494 : LOGICAL, INTENT(IN), OPTIONAL :: allocate_mat_L
495 :
496 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reorder_mat_L'
497 :
498 : INTEGER :: handle, ikp, j_size, nblk
499 156 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
500 : LOGICAL :: do_kpoints, my_allocate_mat_L
501 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
502 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
503 : TYPE(cp_fm_type) :: fm_mat_L_transposed, fmdummy
504 :
505 156 : CALL timeset(routineN, handle)
506 :
507 156 : do_kpoints = .FALSE.
508 156 : IF (PRESENT(ikp_local) .AND. PRESENT(fm_struct_sub_kp)) THEN
509 44 : do_kpoints = .TRUE.
510 : END IF
511 :
512 : ! Get the fm_struct for fm_mat_L
513 156 : NULLIFY (fm_struct)
514 156 : IF (dimen_RI == dimen_RI_red) THEN
515 152 : fm_struct => fm_struct_template
516 : ELSE
517 : ! The template is assumed to be square such that we need a new fm_struct if dimensions are not equal
518 4 : CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI_red, ncol_global=dimen_RI, template_fmstruct=fm_struct_template)
519 : END IF
520 :
521 : ! Start to allocate the new full matrix
522 2840 : ALLOCATE (fm_mat_L(SIZE(fm_matrix_Minv_L_kpoints, 1), SIZE(fm_matrix_Minv_L_kpoints, 2)))
523 1220 : DO ikp = 1, SIZE(fm_matrix_Minv_L_kpoints, 1)
524 3236 : DO j_size = 1, SIZE(fm_matrix_Minv_L_kpoints, 2)
525 3080 : IF (do_kpoints) THEN
526 1904 : IF (ikp == first_ikp_local .OR. ikp_local == -1) THEN
527 1904 : CALL cp_fm_create(fm_mat_L(ikp, j_size), fm_struct_sub_kp)
528 1904 : CALL cp_fm_set_all(fm_mat_L(ikp, j_size), 0.0_dp)
529 : END IF
530 : ELSE
531 112 : CALL cp_fm_create(fm_mat_L(ikp, j_size), fm_struct)
532 112 : CALL cp_fm_set_all(fm_mat_L(ikp, j_size), 0.0_dp)
533 : END IF
534 : END DO
535 : END DO
536 :
537 : ! For the transposed matric we need a different fm_struct
538 156 : IF (dimen_RI == dimen_RI_red) THEN
539 152 : fm_struct => fm_mat_L(first_ikp_local, 1)%matrix_struct
540 : ELSE
541 4 : CALL cp_fm_struct_release(fm_struct)
542 :
543 : ! Create a fm_struct with transposed sizes
544 4 : NULLIFY (fm_struct)
545 : CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI, ncol_global=dimen_RI_red, &
546 4 : template_fmstruct=fm_mat_L(first_ikp_local, 1)%matrix_struct) !, force_block=.TRUE.)
547 : END IF
548 :
549 : ! Allocate buffer matrix
550 156 : CALL cp_fm_create(fm_mat_L_transposed, fm_struct)
551 156 : CALL cp_fm_set_all(matrix=fm_mat_L_transposed, alpha=0.0_dp)
552 :
553 156 : IF (dimen_RI /= dimen_RI_red) CALL cp_fm_struct_release(fm_struct)
554 :
555 156 : CALL cp_fm_get_info(fm_mat_L_transposed, context=blacs_env)
556 :
557 : ! For k-points copy matrices of your group
558 : ! Without kpoints, transpose matrix
559 : ! without kpoints, the size of fm_mat_L is 1x1. with kpoints, the size is N_kpoints x 2 (2 for real/complex)
560 1220 : DO ikp = 1, SIZE(fm_matrix_Minv_L_kpoints, 1)
561 3236 : DO j_size = 1, SIZE(fm_matrix_Minv_L_kpoints, 2)
562 3080 : IF (do_kpoints) THEN
563 1904 : IF (ikp_local == ikp .OR. ikp_local == -1) THEN
564 1904 : CALL cp_fm_copy_general(fm_matrix_Minv_L_kpoints(ikp, j_size), fm_mat_L_transposed, para_env)
565 1904 : CALL cp_fm_to_fm(fm_mat_L_transposed, fm_mat_L(ikp, j_size))
566 : ELSE
567 0 : CALL cp_fm_copy_general(fm_matrix_Minv_L_kpoints(ikp, j_size), fmdummy, para_env)
568 : END IF
569 : ELSE
570 112 : CALL cp_fm_copy_general(fm_matrix_Minv_L_kpoints(ikp, j_size), fm_mat_L_transposed, blacs_env%para_env)
571 112 : CALL cp_fm_transpose(fm_mat_L_transposed, fm_mat_L(ikp, j_size))
572 : END IF
573 : END DO
574 : END DO
575 :
576 : ! Release old matrix
577 156 : CALL cp_fm_release(fm_matrix_Minv_L_kpoints)
578 : ! Release buffer
579 156 : CALL cp_fm_release(fm_mat_L_transposed)
580 :
581 156 : my_allocate_mat_L = .TRUE.
582 156 : IF (PRESENT(allocate_mat_L)) my_allocate_mat_L = allocate_mat_L
583 :
584 22 : IF (my_allocate_mat_L) THEN
585 : ! Create sparse variant of L
586 : NULLIFY (mat_L%matrix)
587 134 : ALLOCATE (mat_L%matrix)
588 134 : IF (dimen_RI == dimen_RI_red) THEN
589 130 : CALL dbcsr_create(mat_L%matrix, template=mat_template)
590 : ELSE
591 4 : CALL dbcsr_get_info(mat_template, nblkrows_total=nblk, col_blk_size=col_blk_size)
592 :
593 4 : CALL calculate_equal_blk_size(row_blk_size, dimen_RI_red, nblk)
594 :
595 4 : CALL dbcsr_create(mat_L%matrix, template=mat_template, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
596 :
597 4 : DEALLOCATE (row_blk_size)
598 : END IF
599 :
600 134 : IF (.NOT. (do_kpoints)) THEN
601 112 : CALL copy_fm_to_dbcsr(fm_mat_L(1, 1), mat_L%matrix)
602 : END IF
603 :
604 : END IF
605 :
606 156 : CALL timestop(handle)
607 :
608 156 : END SUBROUTINE reorder_mat_L
609 :
610 : ! **************************************************************************************************
611 : !> \brief ...
612 : !> \param blk_size_new ...
613 : !> \param dimen_RI_red ...
614 : !> \param nblk ...
615 : ! **************************************************************************************************
616 4 : SUBROUTINE calculate_equal_blk_size(blk_size_new, dimen_RI_red, nblk)
617 : INTEGER, DIMENSION(:), POINTER :: blk_size_new
618 : INTEGER, INTENT(IN) :: dimen_RI_red, nblk
619 :
620 : INTEGER :: col_per_blk, remainder
621 :
622 4 : NULLIFY (blk_size_new)
623 12 : ALLOCATE (blk_size_new(nblk))
624 :
625 4 : remainder = MOD(dimen_RI_red, nblk)
626 4 : col_per_blk = dimen_RI_red/nblk
627 :
628 : ! Determine a new distribution for the columns (corresponding to the number of columns)
629 10 : IF (remainder > 0) blk_size_new(1:remainder) = col_per_blk + 1
630 10 : blk_size_new(remainder + 1:nblk) = col_per_blk
631 :
632 4 : END SUBROUTINE calculate_equal_blk_size
633 :
634 : ! **************************************************************************************************
635 : !> \brief ...
636 : !> \param fm_mat_S ...
637 : !> \param do_ri_sos_laplace_mp2 ...
638 : !> \param first_cycle ...
639 : !> \param iter_sc_GW0 ...
640 : !> \param virtual ...
641 : !> \param Eigenval ...
642 : !> \param Eigenval_scf ...
643 : !> \param homo ...
644 : !> \param omega ...
645 : !> \param omega_old ...
646 : !> \param jquad ...
647 : !> \param mm_style ...
648 : !> \param dimen_RI ...
649 : !> \param dimen_ia ...
650 : !> \param alpha ...
651 : !> \param fm_mat_Q ...
652 : !> \param fm_mat_Q_gemm ...
653 : !> \param do_bse ...
654 : !> \param fm_mat_Q_static_bse_gemm ...
655 : !> \param dgemm_counter ...
656 : !> \param num_integ_points ...
657 : ! **************************************************************************************************
658 2812 : SUBROUTINE calc_mat_Q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, iter_sc_GW0, virtual, &
659 1406 : Eigenval, Eigenval_scf, &
660 : homo, omega, omega_old, jquad, mm_style, dimen_RI, dimen_ia, alpha, fm_mat_Q, fm_mat_Q_gemm, &
661 : do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, num_integ_points)
662 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
663 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, first_cycle
664 : INTEGER, INTENT(IN) :: iter_sc_GW0, virtual
665 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval, Eigenval_scf
666 : INTEGER, INTENT(IN) :: homo
667 : REAL(KIND=dp), INTENT(IN) :: omega, omega_old
668 : INTEGER, INTENT(IN) :: jquad, mm_style, dimen_RI, dimen_ia
669 : REAL(KIND=dp), INTENT(IN) :: alpha
670 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q, fm_mat_Q_gemm
671 : LOGICAL, INTENT(IN) :: do_bse
672 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q_static_bse_gemm
673 : TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter
674 : INTEGER, INTENT(IN) :: num_integ_points
675 :
676 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mat_Q'
677 :
678 : INTEGER :: handle
679 :
680 1406 : CALL timeset(routineN, handle)
681 :
682 1406 : IF (do_ri_sos_laplace_mp2) THEN
683 : ! the first index of tau_tj starts with 0 (see mp2_weights)
684 128 : CALL calc_fm_mat_S_laplace(fm_mat_S, homo, virtual, Eigenval, omega - omega_old)
685 : ELSE
686 1278 : IF (iter_sc_GW0 == 1) THEN
687 : CALL calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval, &
688 1248 : homo, omega, omega_old)
689 : ELSE
690 : CALL calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval_scf, &
691 30 : homo, omega, omega_old)
692 : END IF
693 : END IF
694 :
695 : CALL contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, &
696 1406 : fm_mat_Q, dgemm_counter)
697 : ! fm_mat_Q_static_bse_gemm does not enter W_ijab (A matrix in TDA), but only full ABBA
698 : ! (since only B_ij_bar enters W_ijab)
699 : ! Changing jquad, since omega=0 is at last idx
700 1406 : IF (do_bse .AND. jquad == num_integ_points) THEN
701 4 : CALL cp_fm_to_fm(fm_mat_Q_gemm, fm_mat_Q_static_bse_gemm)
702 : END IF
703 1406 : CALL timestop(handle)
704 :
705 1406 : END SUBROUTINE calc_mat_Q
706 :
707 : ! **************************************************************************************************
708 : !> \brief ...
709 : !> \param fm_mat_S ...
710 : !> \param virtual ...
711 : !> \param Eigenval_last ...
712 : !> \param homo ...
713 : !> \param omega_old ...
714 : ! **************************************************************************************************
715 110 : SUBROUTINE remove_scaling_factor_rpa(fm_mat_S, virtual, Eigenval_last, homo, omega_old)
716 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
717 : INTEGER, INTENT(IN) :: virtual
718 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval_last
719 : INTEGER, INTENT(IN) :: homo
720 : REAL(KIND=dp), INTENT(IN) :: omega_old
721 :
722 : CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_scaling_factor_rpa'
723 :
724 : INTEGER :: avirt, handle, i_global, iiB, iocc, &
725 : ncol_local
726 110 : INTEGER, DIMENSION(:), POINTER :: col_indices
727 : REAL(KIND=dp) :: eigen_diff
728 :
729 110 : CALL timeset(routineN, handle)
730 :
731 : ! get info of fm_mat_S
732 : CALL cp_fm_get_info(matrix=fm_mat_S, &
733 : ncol_local=ncol_local, &
734 110 : col_indices=col_indices)
735 :
736 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
737 110 : !$OMP SHARED(ncol_local,col_indices,Eigenval_last,fm_mat_S,virtual,homo,omega_old)
738 : DO iiB = 1, ncol_local
739 : i_global = col_indices(iiB)
740 :
741 : iocc = MAX(1, i_global - 1)/virtual + 1
742 : avirt = i_global - (iocc - 1)*virtual
743 : eigen_diff = Eigenval_last(avirt + homo) - Eigenval_last(iocc)
744 :
745 : fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)/ &
746 : SQRT(eigen_diff/(eigen_diff**2 + omega_old**2))
747 :
748 : END DO
749 :
750 110 : CALL timestop(handle)
751 :
752 110 : END SUBROUTINE remove_scaling_factor_rpa
753 :
754 : ! **************************************************************************************************
755 : !> \brief ...
756 : !> \param fm_mat_S ...
757 : !> \param first_cycle ...
758 : !> \param virtual ...
759 : !> \param Eigenval ...
760 : !> \param homo ...
761 : !> \param omega ...
762 : !> \param omega_old ...
763 : ! **************************************************************************************************
764 1348 : SUBROUTINE calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval, homo, &
765 : omega, omega_old)
766 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
767 : LOGICAL, INTENT(IN) :: first_cycle
768 : INTEGER, INTENT(IN) :: virtual
769 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
770 : INTEGER, INTENT(IN) :: homo
771 : REAL(KIND=dp), INTENT(IN) :: omega, omega_old
772 :
773 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_fm_mat_S_rpa'
774 :
775 : INTEGER :: avirt, handle, i_global, iiB, iocc, &
776 : ncol_local
777 1348 : INTEGER, DIMENSION(:), POINTER :: col_indices
778 : REAL(KIND=dp) :: eigen_diff
779 :
780 1348 : CALL timeset(routineN, handle)
781 :
782 : ! get info of fm_mat_S
783 : CALL cp_fm_get_info(matrix=fm_mat_S, &
784 : ncol_local=ncol_local, &
785 1348 : col_indices=col_indices)
786 :
787 : ! update G matrix with the new value of omega
788 1348 : IF (first_cycle) THEN
789 : ! In this case just update the matrix (symmetric form) with
790 : ! SQRT((epsi_a-epsi_i)/((epsi_a-epsi_i)**2+omega**2))
791 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
792 204 : !$OMP SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega)
793 : DO iiB = 1, ncol_local
794 : i_global = col_indices(iiB)
795 :
796 : iocc = MAX(1, i_global - 1)/virtual + 1
797 : avirt = i_global - (iocc - 1)*virtual
798 : eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
799 :
800 : fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)* &
801 : SQRT(eigen_diff/(eigen_diff**2 + omega**2))
802 :
803 : END DO
804 : ELSE
805 : ! In this case the update has to remove the old omega component thus
806 : ! SQRT(((epsi_a-epsi_i)**2+omega_old**2)/((epsi_a-epsi_i)**2+omega**2))
807 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
808 1144 : !$OMP SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega,omega_old)
809 : DO iiB = 1, ncol_local
810 : i_global = col_indices(iiB)
811 :
812 : iocc = MAX(1, i_global - 1)/virtual + 1
813 : avirt = i_global - (iocc - 1)*virtual
814 : eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
815 :
816 : fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)* &
817 : SQRT((eigen_diff**2 + omega_old**2)/(eigen_diff**2 + omega**2))
818 :
819 : END DO
820 : END IF
821 :
822 1348 : CALL timestop(handle)
823 :
824 1348 : END SUBROUTINE calc_fm_mat_S_rpa
825 :
826 : ! **************************************************************************************************
827 : !> \brief ...
828 : !> \param mm_style ...
829 : !> \param dimen_RI ...
830 : !> \param dimen_ia ...
831 : !> \param alpha ...
832 : !> \param fm_mat_S ...
833 : !> \param fm_mat_Q_gemm ...
834 : !> \param fm_mat_Q ...
835 : !> \param dgemm_counter ...
836 : ! **************************************************************************************************
837 1406 : SUBROUTINE contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, &
838 : fm_mat_Q, dgemm_counter)
839 :
840 : INTEGER, INTENT(IN) :: mm_style, dimen_RI, dimen_ia
841 : REAL(KIND=dp), INTENT(IN) :: alpha
842 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S, fm_mat_Q_gemm, fm_mat_Q
843 : TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter
844 :
845 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_S_to_Q'
846 :
847 : INTEGER :: handle
848 :
849 1406 : CALL timeset(routineN, handle)
850 :
851 1406 : CALL dgemm_counter_start(dgemm_counter)
852 2802 : SELECT CASE (mm_style)
853 : CASE (wfc_mm_style_gemm)
854 : ! waste-fully computes the full symmetrix matrix, but maybe faster than cp_fm_syrk for optimized cp_fm_gemm
855 : CALL parallel_gemm(transa="N", transb="T", m=dimen_RI, n=dimen_RI, k=dimen_ia, alpha=alpha, &
856 : matrix_a=fm_mat_S, matrix_b=fm_mat_S, beta=0.0_dp, &
857 1396 : matrix_c=fm_mat_Q_gemm)
858 : CASE (wfc_mm_style_syrk)
859 : ! will only compute the upper half of the matrix, which is fine, since we only use it for cholesky later
860 : CALL cp_fm_syrk(uplo='U', trans='N', k=dimen_ia, alpha=alpha, matrix_a=fm_mat_S, &
861 10 : ia=1, ja=1, beta=0.0_dp, matrix_c=fm_mat_Q_gemm)
862 : CASE DEFAULT
863 1406 : CPABORT("")
864 : END SELECT
865 1406 : CALL dgemm_counter_stop(dgemm_counter, dimen_RI, dimen_RI, dimen_ia)
866 :
867 : ! copy/redistribute fm_mat_Q_gemm to fm_mat_Q
868 1406 : CALL cp_fm_set_all(matrix=fm_mat_Q, alpha=0.0_dp)
869 : CALL cp_fm_to_fm_submat_general(fm_mat_Q_gemm, fm_mat_Q, dimen_RI, dimen_RI, 1, 1, 1, 1, &
870 1406 : fm_mat_Q_gemm%matrix_struct%context)
871 :
872 1406 : CALL timestop(handle)
873 :
874 1406 : END SUBROUTINE contract_S_to_Q
875 :
876 : ! **************************************************************************************************
877 : !> \brief ...
878 : !> \param dimen_RI ...
879 : !> \param trace_Qomega ...
880 : !> \param fm_mat_Q ...
881 : !> \param para_env_RPA ...
882 : ! **************************************************************************************************
883 3936 : SUBROUTINE Q_trace_and_add_unit_matrix(dimen_RI, trace_Qomega, fm_mat_Q, para_env_RPA)
884 :
885 : INTEGER, INTENT(IN) :: dimen_RI
886 : REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(OUT) :: trace_Qomega
887 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
888 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA
889 :
890 : CHARACTER(LEN=*), PARAMETER :: routineN = 'Q_trace_and_add_unit_matrix'
891 :
892 : INTEGER :: handle, i_global, iiB, j_global, jjB, &
893 : ncol_local, nrow_local
894 1968 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
895 :
896 1968 : CALL timeset(routineN, handle)
897 :
898 : CALL cp_fm_get_info(matrix=fm_mat_Q, &
899 : nrow_local=nrow_local, &
900 : ncol_local=ncol_local, &
901 : row_indices=row_indices, &
902 1968 : col_indices=col_indices)
903 :
904 : ! calculate the trace of Q and add 1 on the diagonal
905 157874 : trace_Qomega = 0.0_dp
906 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
907 1968 : !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,trace_Qomega,fm_mat_Q,dimen_RI)
908 : DO jjB = 1, ncol_local
909 : j_global = col_indices(jjB)
910 : DO iiB = 1, nrow_local
911 : i_global = row_indices(iiB)
912 : IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
913 : trace_Qomega(i_global) = fm_mat_Q%local_data(iiB, jjB)
914 : fm_mat_Q%local_data(iiB, jjB) = fm_mat_Q%local_data(iiB, jjB) + 1.0_dp
915 : END IF
916 : END DO
917 : END DO
918 1968 : CALL para_env_RPA%sum(trace_Qomega)
919 :
920 1968 : CALL timestop(handle)
921 :
922 1968 : END SUBROUTINE Q_trace_and_add_unit_matrix
923 :
924 : ! **************************************************************************************************
925 : !> \brief ...
926 : !> \param dimen_RI ...
927 : !> \param trace_Qomega ...
928 : !> \param fm_mat_Q ...
929 : !> \param para_env_RPA ...
930 : !> \param Erpa ...
931 : !> \param wjquad ...
932 : ! **************************************************************************************************
933 1836 : SUBROUTINE compute_Erpa_by_freq_int(dimen_RI, trace_Qomega, fm_mat_Q, para_env_RPA, Erpa, wjquad)
934 :
935 : INTEGER, INTENT(IN) :: dimen_RI
936 : REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(IN) :: trace_Qomega
937 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
938 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA
939 : REAL(KIND=dp), INTENT(INOUT) :: Erpa
940 : REAL(KIND=dp), INTENT(IN) :: wjquad
941 :
942 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Erpa_by_freq_int'
943 :
944 : INTEGER :: handle, i_global, iiB, info_chol, &
945 : j_global, jjB, ncol_local, nrow_local
946 1836 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
947 : REAL(KIND=dp) :: FComega
948 1836 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Q_log
949 :
950 1836 : CALL timeset(routineN, handle)
951 :
952 : CALL cp_fm_get_info(matrix=fm_mat_Q, &
953 : nrow_local=nrow_local, &
954 : ncol_local=ncol_local, &
955 : row_indices=row_indices, &
956 1836 : col_indices=col_indices)
957 :
958 : ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition
959 1836 : CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q, n=dimen_RI, info_out=info_chol)
960 1836 : IF (info_chol .NE. 0) THEN
961 : CALL cp_warn(__LOCATION__, &
962 : "The Cholesky decomposition before inverting the RPA matrix / dielectric "// &
963 : "function failed. "// &
964 : "In case of low-scaling RPA/GW, decreasing EPS_FILTER in the &LOW_SCALING "// &
965 : "section might "// &
966 : "increase the overall accuracy making the matrix positive definite. "// &
967 0 : "Code will abort.")
968 : END IF
969 :
970 1836 : CPASSERT(info_chol == 0)
971 :
972 5508 : ALLOCATE (Q_log(dimen_RI))
973 148550 : Q_log = 0.0_dp
974 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
975 1836 : !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,Q_log,fm_mat_Q,dimen_RI)
976 : DO jjB = 1, ncol_local
977 : j_global = col_indices(jjB)
978 : DO iiB = 1, nrow_local
979 : i_global = row_indices(iiB)
980 : IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
981 : Q_log(i_global) = 2.0_dp*LOG(fm_mat_Q%local_data(iiB, jjB))
982 : END IF
983 : END DO
984 : END DO
985 1836 : CALL para_env_RPA%sum(Q_log)
986 :
987 : ! the following frequency integration is Eq. (27) in M. Del Ben et al., JCTC 9, 2654 (2013)
988 : ! (https://doi.org/10.1021/ct4002202)
989 1836 : FComega = 0.0_dp
990 148550 : DO iiB = 1, dimen_RI
991 146714 : IF (MODULO(iiB, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE
992 148550 : FComega = FComega + (Q_log(iiB) - trace_Qomega(iiB))/2.0_dp
993 : END DO
994 1836 : Erpa = Erpa + FComega*wjquad
995 :
996 1836 : DEALLOCATE (Q_log)
997 :
998 1836 : CALL timestop(handle)
999 :
1000 3672 : END SUBROUTINE compute_Erpa_by_freq_int
1001 :
1002 : ! **************************************************************************************************
1003 : !> \brief ...
1004 : !> \param fm_struct_sub_kp ...
1005 : !> \param para_env ...
1006 : !> \param dimen_RI ...
1007 : !> \param ikp_local ...
1008 : !> \param first_ikp_local ...
1009 : ! **************************************************************************************************
1010 22 : SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_RI, &
1011 : ikp_local, first_ikp_local)
1012 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_sub_kp
1013 : TYPE(mp_para_env_type), POINTER :: para_env
1014 : INTEGER, INTENT(IN) :: dimen_RI
1015 : INTEGER, INTENT(OUT) :: ikp_local, first_ikp_local
1016 :
1017 : CHARACTER(len=*), PARAMETER :: routineN = 'get_sub_para_kp'
1018 :
1019 : INTEGER :: color_sub_kp, handle, num_proc_per_kp
1020 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub_kp
1021 : TYPE(mp_para_env_type), POINTER :: para_env_sub_kp
1022 :
1023 22 : CALL timeset(routineN, handle)
1024 :
1025 : ! we use all processors for every k-point, subgroups for cp_cfm_heevd only seems to work for
1026 : ! very small subgroups with 1, 2, or 3 MPI ranks. For more MPI-ranks, eigenvalues and
1027 : ! eigenvectors coming out of cp_cfm_heevd are totally wrong unfortunately.
1028 22 : num_proc_per_kp = para_env%num_pe
1029 :
1030 : ! IF(nkp > para_env%num_pe) THEN
1031 : ! num_proc_per_kp = para_env%num_pe
1032 : ! ELSE
1033 : ! num_proc_per_kp = para_env%num_pe/nkp
1034 : ! END IF
1035 :
1036 22 : color_sub_kp = para_env%mepos/num_proc_per_kp
1037 22 : ALLOCATE (para_env_sub_kp)
1038 22 : CALL para_env_sub_kp%from_split(para_env, color_sub_kp)
1039 :
1040 : ! grid_2d(1) = 1
1041 : ! grid_2d(2) = para_env_sub_kp%num_pe
1042 :
1043 22 : NULLIFY (blacs_env_sub_kp)
1044 : ! CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp, grid_2d=grid_2d)
1045 22 : CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp)
1046 :
1047 22 : NULLIFY (fm_struct_sub_kp)
1048 : CALL cp_fm_struct_create(fm_struct_sub_kp, context=blacs_env_sub_kp, nrow_global=dimen_RI, &
1049 22 : ncol_global=dimen_RI, para_env=para_env_sub_kp)
1050 :
1051 22 : CALL cp_blacs_env_release(blacs_env_sub_kp)
1052 :
1053 : ! IF(nkp > para_env%num_pe) THEN
1054 : ! every processor has all ikp's
1055 22 : ikp_local = -1
1056 22 : first_ikp_local = 1
1057 : ! ELSE
1058 : ! ikp_local = 0
1059 : ! first_ikp_local = 1
1060 : ! DO ikp = 1, nkp
1061 : ! IF(MOD(ikp-1, para_env%num_pe/num_proc_per_kp) == color_sub_kp) THEN
1062 : ! ikp_local = ikp
1063 : ! first_ikp_local = ikp
1064 : ! END IF
1065 : ! END DO
1066 : ! END IF
1067 :
1068 22 : CALL mp_para_env_release(para_env_sub_kp)
1069 :
1070 22 : CALL timestop(handle)
1071 :
1072 22 : END SUBROUTINE get_sub_para_kp
1073 :
1074 : ! **************************************************************************************************
1075 : !> \brief ...
1076 : !> \param fm_mo_coeff_occ ...
1077 : !> \param fm_mo_coeff_virt ...
1078 : !> \param fm_scaled_dm_occ_tau ...
1079 : !> \param fm_scaled_dm_virt_tau ...
1080 : !> \param index_to_cell_3c ...
1081 : !> \param cell_to_index_3c ...
1082 : !> \param do_ic_model ...
1083 : !> \param do_kpoints_cubic_RPA ...
1084 : !> \param do_kpoints_from_Gamma ...
1085 : !> \param do_ri_Sigma_x ...
1086 : !> \param has_mat_P_blocks ...
1087 : !> \param wkp_W ...
1088 : !> \param cfm_mat_Q ...
1089 : !> \param fm_mat_Minv_L_kpoints ...
1090 : !> \param fm_mat_L_kpoints ...
1091 : !> \param fm_matrix_Minv ...
1092 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
1093 : !> \param fm_mat_RI_global_work ...
1094 : !> \param fm_mat_work ...
1095 : !> \param fm_mo_coeff_occ_scaled ...
1096 : !> \param fm_mo_coeff_virt_scaled ...
1097 : !> \param mat_dm ...
1098 : !> \param mat_L ...
1099 : !> \param mat_MinvVMinv ...
1100 : !> \param mat_P_omega ...
1101 : !> \param mat_P_omega_kp ...
1102 : !> \param t_3c_M ...
1103 : !> \param t_3c_O ...
1104 : !> \param t_3c_O_compressed ...
1105 : !> \param t_3c_O_ind ...
1106 : !> \param mat_work ...
1107 : !> \param qs_env ...
1108 : ! **************************************************************************************************
1109 134 : SUBROUTINE dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, &
1110 : fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, index_to_cell_3c, &
1111 : cell_to_index_3c, do_ic_model, &
1112 : do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, &
1113 : has_mat_P_blocks, &
1114 : wkp_W, cfm_mat_Q, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
1115 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
1116 : fm_mat_RI_global_work, fm_mat_work, &
1117 : fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_L, &
1118 : mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
1119 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
1120 : mat_work, qs_env)
1121 :
1122 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mo_coeff_occ, fm_mo_coeff_virt
1123 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_scaled_dm_occ_tau, &
1124 : fm_scaled_dm_virt_tau
1125 : INTEGER, ALLOCATABLE, DIMENSION(:, :), &
1126 : INTENT(INOUT) :: index_to_cell_3c
1127 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
1128 : INTENT(INOUT) :: cell_to_index_3c
1129 : LOGICAL, INTENT(IN) :: do_ic_model, do_kpoints_cubic_RPA, &
1130 : do_kpoints_from_Gamma, do_ri_Sigma_x
1131 : LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), &
1132 : INTENT(INOUT) :: has_mat_P_blocks
1133 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1134 : INTENT(INOUT) :: wkp_W
1135 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm_mat_Q
1136 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
1137 : fm_matrix_Minv, &
1138 : fm_matrix_Minv_Vtrunc_Minv
1139 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_RI_global_work, fm_mat_work, &
1140 : fm_mo_coeff_occ_scaled, &
1141 : fm_mo_coeff_virt_scaled
1142 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_dm, mat_L, mat_MinvVMinv
1143 : TYPE(dbcsr_p_type), ALLOCATABLE, &
1144 : DIMENSION(:, :, :), INTENT(INOUT) :: mat_P_omega
1145 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_P_omega_kp
1146 : TYPE(dbt_type) :: t_3c_M
1147 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O
1148 : TYPE(hfx_compression_type), ALLOCATABLE, &
1149 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_compressed
1150 : TYPE(block_ind_type), ALLOCATABLE, &
1151 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_ind
1152 : TYPE(dbcsr_type), POINTER :: mat_work
1153 : TYPE(qs_environment_type), POINTER :: qs_env
1154 :
1155 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dealloc_im_time'
1156 :
1157 : INTEGER :: cut_memory, handle, i_kp, i_mem, i_size, &
1158 : ispin, j_size, jquad, nspins, unused
1159 : LOGICAL :: my_open_shell
1160 :
1161 134 : CALL timeset(routineN, handle)
1162 :
1163 134 : nspins = SIZE(fm_mo_coeff_occ)
1164 134 : my_open_shell = (nspins == 2)
1165 :
1166 134 : CALL cp_fm_release(fm_scaled_dm_occ_tau)
1167 134 : CALL cp_fm_release(fm_scaled_dm_virt_tau)
1168 298 : DO ispin = 1, SIZE(fm_mo_coeff_occ)
1169 164 : CALL cp_fm_release(fm_mo_coeff_occ(ispin))
1170 298 : CALL cp_fm_release(fm_mo_coeff_virt(ispin))
1171 : END DO
1172 134 : CALL cp_fm_release(fm_mo_coeff_occ_scaled)
1173 134 : CALL cp_fm_release(fm_mo_coeff_virt_scaled)
1174 :
1175 134 : CALL cp_fm_release(fm_mat_Minv_L_kpoints)
1176 134 : CALL cp_fm_release(fm_mat_L_kpoints)
1177 :
1178 134 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
1179 22 : CALL cp_fm_release(fm_matrix_Minv_Vtrunc_Minv)
1180 22 : CALL cp_fm_release(fm_matrix_Minv)
1181 : END IF
1182 :
1183 134 : CALL cp_fm_release(fm_mat_work)
1184 :
1185 134 : IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
1186 112 : CALL dbcsr_release(mat_work)
1187 112 : DEALLOCATE (mat_work)
1188 : END IF
1189 :
1190 134 : CALL dbcsr_release(mat_L%matrix)
1191 134 : DEALLOCATE (mat_L%matrix)
1192 :
1193 134 : IF (do_ri_Sigma_x .OR. do_ic_model) THEN
1194 106 : CALL dbcsr_release(mat_MinvVMinv%matrix)
1195 106 : DEALLOCATE (mat_MinvVMinv%matrix)
1196 : END IF
1197 134 : IF (do_ri_Sigma_x) THEN
1198 106 : CALL dbcsr_release(mat_dm%matrix)
1199 106 : DEALLOCATE (mat_dm%matrix)
1200 : END IF
1201 :
1202 134 : DEALLOCATE (index_to_cell_3c, cell_to_index_3c)
1203 :
1204 134 : IF (ALLOCATED(mat_P_omega)) THEN
1205 298 : DO ispin = 1, SIZE(mat_P_omega, 3)
1206 1016 : DO i_kp = 1, SIZE(mat_P_omega, 2)
1207 5250 : DO jquad = 1, SIZE(mat_P_omega, 1)
1208 5086 : CALL dbcsr_deallocate_matrix(mat_P_omega(jquad, i_kp, ispin)%matrix)
1209 : END DO
1210 : END DO
1211 : END DO
1212 134 : DEALLOCATE (mat_P_omega)
1213 : END IF
1214 :
1215 284 : DO i_size = 1, SIZE(t_3c_O, 1)
1216 514 : DO j_size = 1, SIZE(t_3c_O, 2)
1217 380 : CALL dbt_destroy(t_3c_O(i_size, j_size))
1218 : END DO
1219 : END DO
1220 :
1221 364 : DEALLOCATE (t_3c_O)
1222 134 : CALL dbt_destroy(t_3c_M)
1223 :
1224 134 : DEALLOCATE (has_mat_P_blocks)
1225 :
1226 134 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
1227 22 : CALL cp_cfm_release(cfm_mat_Q)
1228 22 : CALL cp_fm_release(fm_mat_RI_global_work)
1229 22 : CALL dbcsr_deallocate_matrix_set(mat_P_omega_kp)
1230 22 : DEALLOCATE (wkp_W)
1231 : END IF
1232 :
1233 134 : cut_memory = SIZE(t_3c_O_compressed, 3)
1234 :
1235 550 : DEALLOCATE (t_3c_O_ind)
1236 284 : DO i_size = 1, SIZE(t_3c_O_compressed, 1)
1237 514 : DO j_size = 1, SIZE(t_3c_O_compressed, 2)
1238 796 : DO i_mem = 1, cut_memory
1239 646 : CALL dealloc_containers(t_3c_O_compressed(i_size, j_size, i_mem), unused)
1240 : END DO
1241 : END DO
1242 : END DO
1243 134 : DEALLOCATE (t_3c_O_compressed)
1244 :
1245 134 : IF (do_kpoints_from_Gamma) THEN
1246 18 : CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_G)
1247 18 : IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
1248 18 : CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma)
1249 18 : CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc)
1250 : END IF
1251 : END IF
1252 :
1253 134 : CALL timestop(handle)
1254 :
1255 134 : END SUBROUTINE dealloc_im_time
1256 :
1257 : ! **************************************************************************************************
1258 : !> \brief ...
1259 : !> \param mat_P_omega ...
1260 : !> \param mat_L ...
1261 : !> \param mat_work ...
1262 : !> \param eps_filter_im_time ...
1263 : !> \param fm_mat_work ...
1264 : !> \param dimen_RI ...
1265 : !> \param dimen_RI_red ...
1266 : !> \param fm_mat_L ...
1267 : !> \param fm_mat_Q ...
1268 : ! **************************************************************************************************
1269 978 : SUBROUTINE contract_P_omega_with_mat_L(mat_P_omega, mat_L, mat_work, eps_filter_im_time, fm_mat_work, dimen_RI, &
1270 : dimen_RI_red, fm_mat_L, fm_mat_Q)
1271 :
1272 : TYPE(dbcsr_type), POINTER :: mat_P_omega, mat_L, mat_work
1273 : REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time
1274 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_work
1275 : INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
1276 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_L, fm_mat_Q
1277 :
1278 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_P_omega_with_mat_L'
1279 :
1280 : INTEGER :: handle
1281 :
1282 978 : CALL timeset(routineN, handle)
1283 :
1284 : ! multiplication with RI metric/Coulomb operator
1285 : CALL dbcsr_multiply("N", "T", 1.0_dp, mat_P_omega, mat_L, &
1286 978 : 0.0_dp, mat_work, filter_eps=eps_filter_im_time)
1287 :
1288 978 : CALL copy_dbcsr_to_fm(mat_work, fm_mat_work)
1289 :
1290 : CALL parallel_gemm('N', 'N', dimen_RI_red, dimen_RI_red, dimen_RI, 1.0_dp, fm_mat_L, fm_mat_work, &
1291 978 : 0.0_dp, fm_mat_Q)
1292 :
1293 : ! Reset mat_work to save memory
1294 978 : CALL dbcsr_set(mat_work, 0.0_dp)
1295 978 : CALL dbcsr_filter(mat_work, 1.0_dp)
1296 :
1297 978 : CALL timestop(handle)
1298 :
1299 978 : END SUBROUTINE contract_P_omega_with_mat_L
1300 :
1301 : END MODULE rpa_util
|
