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 Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations
10 : !> \par History
11 : !> 11.2023 created [Maximilian Graml]
12 : ! **************************************************************************************************
13 : MODULE bse_util
14 : USE atomic_kind_types, ONLY: atomic_kind_type
15 : USE cell_types, ONLY: cell_type
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
19 : dbcsr_init_p,&
20 : dbcsr_p_type,&
21 : dbcsr_set,&
22 : dbcsr_type_symmetric
23 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
24 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
25 : dbcsr_allocate_matrix_set,&
26 : dbcsr_deallocate_matrix_set
27 : USE cp_fm_basic_linalg, ONLY: cp_fm_trace,&
28 : cp_fm_uplo_to_full
29 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
30 : cp_fm_cholesky_invert
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: cp_fm_create,&
35 : cp_fm_get_info,&
36 : cp_fm_release,&
37 : cp_fm_set_all,&
38 : cp_fm_to_fm_submat,&
39 : cp_fm_to_fm_submat_general,&
40 : cp_fm_type
41 : USE cp_log_handling, ONLY: cp_get_default_logger,&
42 : cp_logger_type
43 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
44 : cp_print_key_unit_nr
45 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
46 : USE input_constants, ONLY: bse_iterdiag,&
47 : bse_screening_alpha,&
48 : bse_screening_rpa,&
49 : bse_screening_tdhf,&
50 : use_mom_ref_coac
51 : USE input_section_types, ONLY: section_vals_type
52 : USE kinds, ONLY: default_path_length,&
53 : dp,&
54 : int_8
55 : USE message_passing, ONLY: mp_para_env_type,&
56 : mp_request_type
57 : USE mo_window, ONLY: mo_window_type
58 : USE moments_utils, ONLY: get_reference_point
59 : USE mp2_types, ONLY: integ_mat_buffer_type,&
60 : mp2_type
61 : USE parallel_gemm_api, ONLY: parallel_gemm
62 : USE particle_list_types, ONLY: particle_list_type
63 : USE particle_types, ONLY: particle_type
64 : USE physcon, ONLY: evolt
65 : USE pw_env_types, ONLY: pw_env_get,&
66 : pw_env_type
67 : USE pw_poisson_types, ONLY: pw_poisson_type
68 : USE pw_pool_types, ONLY: pw_pool_p_type,&
69 : pw_pool_type
70 : USE pw_types, ONLY: pw_c1d_gs_type,&
71 : pw_r3d_rs_type
72 : USE qs_collocate_density, ONLY: calculate_wavefunction
73 : USE qs_environment_types, ONLY: get_qs_env,&
74 : qs_environment_type
75 : USE qs_kind_types, ONLY: qs_kind_type
76 : USE qs_mo_types, ONLY: get_mo_set,&
77 : mo_set_type
78 : USE qs_moments, ONLY: build_local_moment_matrix
79 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
80 : USE qs_subsys_types, ONLY: qs_subsys_get,&
81 : qs_subsys_type
82 : USE rpa_communication, ONLY: communicate_buffer
83 : USE util, ONLY: sort,&
84 : sort_unique
85 : #include "./base/base_uses.f90"
86 :
87 : IMPLICIT NONE
88 :
89 : PRIVATE
90 :
91 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_util'
92 :
93 : PUBLIC :: mult_B_with_W, fm_general_add_bse, truncate_fm, &
94 : deallocate_matrices_bse, comp_eigvec_coeff_BSE, sort_excitations, &
95 : estimate_BSE_resources, filter_eigvec_contrib, truncate_BSE_matrices, &
96 : adapt_BSE_input_params, get_multipoles_mo, &
97 : reshuffle_eigvec, print_bse_nto_cubes, trace_exciton_descr, &
98 : get_bse_spin_block_layout, assemble_joint_ov_slab, &
99 : ia_of_occ_virt, occ_of_ia, virt_of_ia
100 :
101 : CONTAINS
102 :
103 : ! **************************************************************************************************
104 : !> \brief Compound transition index ia = (i-1)*virt + a of the pair (i, a): a is the fast index
105 : !> \param i_occ occupied level
106 : !> \param a_virt virtual level
107 : !> \param virt number of virtual levels
108 : !> \return ia
109 : ! **************************************************************************************************
110 12240 : PURE FUNCTION ia_of_occ_virt(i_occ, a_virt, virt) RESULT(ia)
111 :
112 : INTEGER, INTENT(IN) :: i_occ, a_virt, virt
113 : INTEGER :: ia
114 :
115 12240 : ia = (i_occ - 1)*virt + a_virt
116 :
117 12240 : END FUNCTION ia_of_occ_virt
118 :
119 : ! **************************************************************************************************
120 : !> \brief Occupied level of the compound transition index ia = (i-1)*virt + a
121 : !> \param ia compound transition index
122 : !> \param virt number of virtual levels
123 : !> \return i
124 : ! **************************************************************************************************
125 0 : PURE FUNCTION occ_of_ia(ia, virt) RESULT(i_occ)
126 :
127 : INTEGER, INTENT(IN) :: ia, virt
128 : INTEGER :: i_occ
129 :
130 0 : i_occ = (ia - 1)/virt + 1
131 :
132 0 : END FUNCTION occ_of_ia
133 :
134 : ! **************************************************************************************************
135 : !> \brief Virtual level of the compound transition index ia = (i-1)*virt + a
136 : !> \param ia compound transition index
137 : !> \param virt number of virtual levels
138 : !> \return a
139 : ! **************************************************************************************************
140 0 : PURE FUNCTION virt_of_ia(ia, virt) RESULT(a_virt)
141 :
142 : INTEGER, INTENT(IN) :: ia, virt
143 : INTEGER :: a_virt
144 :
145 0 : a_virt = MOD(ia - 1, virt) + 1
146 :
147 0 : END FUNCTION virt_of_ia
148 :
149 : ! **************************************************************************************************
150 : !> \brief Multiplies B-matrix (RI-3c-Integrals) with W (screening) to obtain \bar{B}
151 : !> \param fm_mat_S_ij_bse ...
152 : !> \param fm_mat_S_ia_bse ...
153 : !> \param fm_mat_S_bar_ia_bse ...
154 : !> \param fm_mat_S_bar_ij_bse ...
155 : !> \param fm_mat_Q_static_bse_gemm ...
156 : !> \param dimen_RI ...
157 : !> \param homo ...
158 : !> \param virtual ...
159 : ! **************************************************************************************************
160 336 : SUBROUTINE mult_B_with_W(fm_mat_S_ij_bse, fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, &
161 : fm_mat_S_bar_ij_bse, fm_mat_Q_static_bse_gemm, &
162 : dimen_RI, homo, virtual)
163 :
164 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ij_bse, fm_mat_S_ia_bse
165 : TYPE(cp_fm_type), INTENT(OUT) :: fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse
166 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q_static_bse_gemm
167 : INTEGER, INTENT(IN) :: dimen_RI, homo, virtual
168 :
169 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mult_B_with_W'
170 :
171 : INTEGER :: handle, i_global, iiB, info_chol, &
172 : j_global, jjB, ncol_local, nrow_local
173 56 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
174 : TYPE(cp_fm_type) :: fm_work
175 :
176 56 : CALL timeset(routineN, handle)
177 :
178 56 : CALL cp_fm_create(fm_mat_S_bar_ia_bse, fm_mat_S_ia_bse%matrix_struct)
179 56 : CALL cp_fm_set_all(fm_mat_S_bar_ia_bse, 0.0_dp)
180 :
181 56 : CALL cp_fm_create(fm_mat_S_bar_ij_bse, fm_mat_S_ij_bse%matrix_struct)
182 56 : CALL cp_fm_set_all(fm_mat_S_bar_ij_bse, 0.0_dp)
183 :
184 56 : CALL cp_fm_create(fm_work, fm_mat_Q_static_bse_gemm%matrix_struct)
185 56 : CALL cp_fm_set_all(fm_work, 0.0_dp)
186 :
187 : ! get info of fm_mat_Q_static_bse and compute ((1+Q(0))^-1-1)
188 : CALL cp_fm_get_info(matrix=fm_mat_Q_static_bse_gemm, &
189 : nrow_local=nrow_local, &
190 : ncol_local=ncol_local, &
191 : row_indices=row_indices, &
192 56 : col_indices=col_indices)
193 :
194 4278 : DO jjB = 1, ncol_local
195 4222 : j_global = col_indices(jjB)
196 238165 : DO iiB = 1, nrow_local
197 233887 : i_global = row_indices(iiB)
198 238109 : IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
199 2415 : fm_mat_Q_static_bse_gemm%local_data(iiB, jjB) = fm_mat_Q_static_bse_gemm%local_data(iiB, jjB) + 1.0_dp
200 : END IF
201 : END DO
202 : END DO
203 :
204 : ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition
205 56 : CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q_static_bse_gemm, n=dimen_RI, info_out=info_chol)
206 :
207 56 : IF (info_chol /= 0) THEN
208 0 : CALL cp_abort(__LOCATION__, 'Cholesky decomposition failed for static polarization in BSE')
209 : END IF
210 :
211 : ! calculate [1+Q(i0)]^-1
212 56 : CALL cp_fm_cholesky_invert(fm_mat_Q_static_bse_gemm)
213 :
214 : ! symmetrize the result
215 56 : CALL cp_fm_uplo_to_full(fm_mat_Q_static_bse_gemm, fm_work)
216 :
217 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=homo**2, k=dimen_RI, alpha=1.0_dp, &
218 : matrix_a=fm_mat_Q_static_bse_gemm, matrix_b=fm_mat_S_ij_bse, beta=0.0_dp, &
219 56 : matrix_c=fm_mat_S_bar_ij_bse)
220 :
221 : ! fm_mat_S_bar_ia_bse has a different blacs_env as fm_mat_S_ij_bse since we take
222 : ! fm_mat_S_ia_bse from RPA. Therefore, we also need a different fm_mat_Q_static_bse_gemm
223 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=homo*virtual, k=dimen_RI, alpha=1.0_dp, &
224 : matrix_a=fm_mat_Q_static_bse_gemm, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
225 56 : matrix_c=fm_mat_S_bar_ia_bse)
226 :
227 56 : CALL cp_fm_release(fm_work)
228 :
229 56 : CALL timestop(handle)
230 :
231 56 : END SUBROUTINE mult_B_with_W
232 :
233 : ! **************************************************************************************************
234 : !> \brief Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,ab
235 : !> to A_ia, jb which needs MPI communication.
236 : !> \param fm_out ...
237 : !> \param fm_in ...
238 : !> \param beta ...
239 : !> \param nrow_secidx_in ...
240 : !> \param ncol_secidx_in ...
241 : !> \param nrow_secidx_out ...
242 : !> \param ncol_secidx_out ...
243 : !> \param unit_nr ...
244 : !> \param reordering ...
245 : !> \param mp2_env ...
246 : !> \param row_offset ...
247 : !> \param col_offset ...
248 : ! **************************************************************************************************
249 536 : SUBROUTINE fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, &
250 : nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env, &
251 : row_offset, col_offset)
252 :
253 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_out
254 : TYPE(cp_fm_type), INTENT(IN) :: fm_in
255 : REAL(kind=dp) :: beta
256 : INTEGER, INTENT(IN) :: nrow_secidx_in, ncol_secidx_in, &
257 : nrow_secidx_out, ncol_secidx_out
258 : INTEGER :: unit_nr
259 : INTEGER, DIMENSION(4) :: reordering
260 : TYPE(mp2_type), INTENT(IN) :: mp2_env
261 : INTEGER, INTENT(IN), OPTIONAL :: row_offset, col_offset
262 :
263 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_general_add_bse'
264 :
265 : INTEGER :: col_idx_loc, dummy, handle, handle2, i_entry_rec, idx_col_out, idx_row_out, ii, &
266 : iproc, jj, my_col_offset, my_row_offset, ncol_block_in, ncol_block_out, ncol_local_in, &
267 : ncol_local_out, nprocs, nrow_block_in, nrow_block_out, nrow_local_in, nrow_local_out, &
268 : proc_send, row_idx_loc, send_pcol, send_prow
269 536 : INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, &
270 : num_entries_send
271 : INTEGER, DIMENSION(4) :: indices_in
272 536 : INTEGER, DIMENSION(:), POINTER :: col_indices_in, col_indices_out, &
273 536 : row_indices_in, row_indices_out
274 : TYPE(integ_mat_buffer_type), ALLOCATABLE, &
275 536 : DIMENSION(:) :: buffer_rec, buffer_send
276 : TYPE(mp_para_env_type), POINTER :: para_env_out
277 536 : TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
278 :
279 : ! Offsets place the reshuffled block into a sub-block of fm_out (open-shell joint matrix);
280 : ! both default 0, recovering the closed-shell single-block placement bit-identically.
281 :
282 536 : my_row_offset = 0
283 536 : my_col_offset = 0
284 72 : IF (PRESENT(row_offset)) my_row_offset = row_offset
285 536 : IF (PRESENT(col_offset)) my_col_offset = col_offset
286 :
287 536 : CALL timeset(routineN, handle)
288 536 : CALL timeset(routineN//"_1_setup", handle2)
289 :
290 536 : para_env_out => fm_out%matrix_struct%para_env
291 : ! A_iajb
292 : ! We start by moving data from local parts of W_ijab to the full matrix A_iajb using buffers
293 : CALL cp_fm_get_info(matrix=fm_out, &
294 : nrow_local=nrow_local_out, &
295 : ncol_local=ncol_local_out, &
296 : row_indices=row_indices_out, &
297 : col_indices=col_indices_out, &
298 : nrow_block=nrow_block_out, &
299 536 : ncol_block=ncol_block_out)
300 :
301 1608 : ALLOCATE (num_entries_rec(0:para_env_out%num_pe - 1))
302 1608 : ALLOCATE (num_entries_send(0:para_env_out%num_pe - 1))
303 :
304 536 : num_entries_rec(:) = 0
305 536 : num_entries_send(:) = 0
306 :
307 536 : dummy = 0
308 :
309 : CALL cp_fm_get_info(matrix=fm_in, &
310 : nrow_local=nrow_local_in, &
311 : ncol_local=ncol_local_in, &
312 : row_indices=row_indices_in, &
313 : col_indices=col_indices_in, &
314 : nrow_block=nrow_block_in, &
315 536 : ncol_block=ncol_block_in)
316 :
317 536 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
318 89 : WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_out%name, &
319 178 : fm_out%matrix_struct%nrow_global
320 89 : WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_out%name, &
321 178 : fm_out%matrix_struct%ncol_global
322 :
323 89 : WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_out%name, nrow_block_out
324 89 : WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_out%name, ncol_block_out
325 :
326 89 : WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_in%name, &
327 178 : fm_in%matrix_struct%nrow_global
328 89 : WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_in%name, &
329 178 : fm_in%matrix_struct%ncol_global
330 :
331 89 : WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_in%name, nrow_block_in
332 89 : WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_in%name, ncol_block_in
333 : END IF
334 :
335 : ! Use scalapack wrapper to find process index in fm_out
336 : ! To that end, we obtain the global index in fm_out from the level indices
337 536 : indices_in(:) = 0
338 12754 : DO row_idx_loc = 1, nrow_local_in
339 12218 : indices_in(1) = (row_indices_in(row_idx_loc) - 1)/nrow_secidx_in + 1
340 12218 : indices_in(2) = MOD(row_indices_in(row_idx_loc) - 1, nrow_secidx_in) + 1
341 1221584 : DO col_idx_loc = 1, ncol_local_in
342 1208830 : indices_in(3) = (col_indices_in(col_idx_loc) - 1)/ncol_secidx_in + 1
343 1208830 : indices_in(4) = MOD(col_indices_in(col_idx_loc) - 1, ncol_secidx_in) + 1
344 :
345 1208830 : idx_row_out = my_row_offset + indices_in(reordering(2)) + (indices_in(reordering(1)) - 1)*nrow_secidx_out
346 1208830 : idx_col_out = my_col_offset + indices_in(reordering(4)) + (indices_in(reordering(3)) - 1)*ncol_secidx_out
347 :
348 1208830 : send_prow = fm_out%matrix_struct%g2p_row(idx_row_out)
349 1208830 : send_pcol = fm_out%matrix_struct%g2p_col(idx_col_out)
350 :
351 1208830 : proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
352 :
353 1221048 : num_entries_send(proc_send) = num_entries_send(proc_send) + 1
354 :
355 : END DO
356 : END DO
357 :
358 536 : CALL timestop(handle2)
359 :
360 536 : CALL timeset(routineN//"_2_comm_entry_nums", handle2)
361 536 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
362 89 : WRITE (unit_nr, '(T2,A10,T13,A27)') 'BSE|DEBUG|', 'Communicating entry numbers'
363 : END IF
364 :
365 536 : CALL para_env_out%alltoall(num_entries_send, num_entries_rec, 1)
366 :
367 536 : CALL timestop(handle2)
368 :
369 536 : CALL timeset(routineN//"_3_alloc_buffer", handle2)
370 536 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
371 89 : WRITE (unit_nr, '(T2,A10,T13,A18)') 'BSE|DEBUG|', 'Allocating buffers'
372 : END IF
373 :
374 : ! Buffers for entries and their indices
375 2680 : ALLOCATE (buffer_rec(0:para_env_out%num_pe - 1))
376 2680 : ALLOCATE (buffer_send(0:para_env_out%num_pe - 1))
377 :
378 : ! allocate data message and corresponding indices
379 1608 : DO iproc = 0, para_env_out%num_pe - 1
380 :
381 3026 : ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc)))
382 1210438 : buffer_rec(iproc)%msg = 0.0_dp
383 :
384 : END DO
385 :
386 1608 : DO iproc = 0, para_env_out%num_pe - 1
387 :
388 3026 : ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc)))
389 1210438 : buffer_send(iproc)%msg = 0.0_dp
390 :
391 : END DO
392 :
393 1608 : DO iproc = 0, para_env_out%num_pe - 1
394 :
395 3026 : ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2))
396 2421412 : buffer_rec(iproc)%indx = 0
397 :
398 : END DO
399 :
400 1608 : DO iproc = 0, para_env_out%num_pe - 1
401 :
402 3026 : ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2))
403 2421412 : buffer_send(iproc)%indx = 0
404 :
405 : END DO
406 :
407 536 : CALL timestop(handle2)
408 :
409 536 : CALL timeset(routineN//"_4_buf_from_fmin_"//fm_out%name, handle2)
410 536 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
411 89 : WRITE (unit_nr, '(T2,A10,T13,A18,A10,A13)') 'BSE|DEBUG|', 'Writing data from ', fm_in%name, ' into buffers'
412 : END IF
413 :
414 1608 : ALLOCATE (entry_counter(0:para_env_out%num_pe - 1))
415 536 : entry_counter(:) = 0
416 :
417 : ! Now we can write the actual data and indices to the send-buffer
418 12754 : DO row_idx_loc = 1, nrow_local_in
419 12218 : indices_in(1) = (row_indices_in(row_idx_loc) - 1)/nrow_secidx_in + 1
420 12218 : indices_in(2) = MOD(row_indices_in(row_idx_loc) - 1, nrow_secidx_in) + 1
421 1221584 : DO col_idx_loc = 1, ncol_local_in
422 1208830 : indices_in(3) = (col_indices_in(col_idx_loc) - 1)/ncol_secidx_in + 1
423 1208830 : indices_in(4) = MOD(col_indices_in(col_idx_loc) - 1, ncol_secidx_in) + 1
424 :
425 1208830 : idx_row_out = my_row_offset + indices_in(reordering(2)) + (indices_in(reordering(1)) - 1)*nrow_secidx_out
426 1208830 : idx_col_out = my_col_offset + indices_in(reordering(4)) + (indices_in(reordering(3)) - 1)*ncol_secidx_out
427 :
428 1208830 : send_prow = fm_out%matrix_struct%g2p_row(idx_row_out)
429 1208830 : send_pcol = fm_out%matrix_struct%g2p_col(idx_col_out)
430 :
431 1208830 : proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
432 1208830 : entry_counter(proc_send) = entry_counter(proc_send) + 1
433 :
434 : buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
435 1208830 : fm_in%local_data(row_idx_loc, col_idx_loc)
436 :
437 1208830 : buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = idx_row_out
438 1221048 : buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = idx_col_out
439 :
440 : END DO
441 : END DO
442 :
443 8040 : ALLOCATE (req_array(1:para_env_out%num_pe, 4))
444 :
445 536 : CALL timestop(handle2)
446 :
447 536 : CALL timeset(routineN//"_5_comm_buffer", handle2)
448 536 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
449 89 : WRITE (unit_nr, '(T2,A10,T13,A21)') 'BSE|DEBUG|', 'Communicating buffers'
450 : END IF
451 :
452 : ! communicate the buffer
453 : CALL communicate_buffer(para_env_out, num_entries_rec, num_entries_send, buffer_rec, &
454 536 : buffer_send, req_array)
455 :
456 536 : CALL timestop(handle2)
457 :
458 536 : CALL timeset(routineN//"_6_buffer_to_fmout"//fm_out%name, handle2)
459 536 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
460 89 : WRITE (unit_nr, '(T2,A10,T13,A24,A10)') 'BSE|DEBUG|', 'Writing from buffers to ', fm_out%name
461 : END IF
462 :
463 : ! fill fm_out with the entries from buffer_rec, i.e. buffer_rec are parts of fm_in
464 536 : nprocs = para_env_out%num_pe
465 :
466 : !$OMP PARALLEL DO DEFAULT(NONE) &
467 : !$OMP SHARED(fm_out, nprocs, num_entries_rec, buffer_rec, beta) &
468 536 : !$OMP PRIVATE(iproc, i_entry_rec, ii, jj)
469 : DO iproc = 0, nprocs - 1
470 : DO i_entry_rec = 1, num_entries_rec(iproc)
471 : ii = fm_out%matrix_struct%g2l_row(buffer_rec(iproc)%indx(i_entry_rec, 1))
472 : jj = fm_out%matrix_struct%g2l_col(buffer_rec(iproc)%indx(i_entry_rec, 2))
473 :
474 : fm_out%local_data(ii, jj) = fm_out%local_data(ii, jj) + beta*buffer_rec(iproc)%msg(i_entry_rec)
475 : END DO
476 : END DO
477 : !$OMP END PARALLEL DO
478 :
479 536 : CALL timestop(handle2)
480 :
481 536 : CALL timeset(routineN//"_7_cleanup", handle2)
482 536 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
483 89 : WRITE (unit_nr, '(T2,A10,T13,A41)') 'BSE|DEBUG|', 'Starting cleanup of communication buffers'
484 : END IF
485 :
486 : !Clean up all the arrays from the communication process
487 1608 : DO iproc = 0, para_env_out%num_pe - 1
488 1072 : DEALLOCATE (buffer_rec(iproc)%msg)
489 1072 : DEALLOCATE (buffer_rec(iproc)%indx)
490 1072 : DEALLOCATE (buffer_send(iproc)%msg)
491 1608 : DEALLOCATE (buffer_send(iproc)%indx)
492 : END DO
493 2680 : DEALLOCATE (buffer_rec, buffer_send)
494 536 : DEALLOCATE (req_array)
495 536 : DEALLOCATE (entry_counter)
496 536 : DEALLOCATE (num_entries_rec, num_entries_send)
497 :
498 536 : CALL timestop(handle2)
499 536 : CALL timestop(handle)
500 :
501 4288 : END SUBROUTINE fm_general_add_bse
502 :
503 : ! **************************************************************************************************
504 : !> \brief Routine for truncating a full matrix as given by the energy cutoffs in the input file.
505 : !> Logic: Matrices have some dimension dimen_RI x nrow_in*ncol_in for the incoming (untruncated) matrix
506 : !> and dimen_RI x nrow_out*ncol_out for the truncated matrix. The truncation is done by resorting the indices
507 : !> via parallel communication.
508 : !> \param fm_out ...
509 : !> \param fm_in ...
510 : !> \param ncol_in ...
511 : !> \param nrow_out ...
512 : !> \param ncol_out ...
513 : !> \param unit_nr ...
514 : !> \param mp2_env ...
515 : !> \param nrow_offset ...
516 : !> \param ncol_offset ...
517 : ! **************************************************************************************************
518 132 : SUBROUTINE truncate_fm(fm_out, fm_in, ncol_in, &
519 : nrow_out, ncol_out, unit_nr, mp2_env, &
520 : nrow_offset, ncol_offset)
521 :
522 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_out
523 : TYPE(cp_fm_type), INTENT(IN) :: fm_in
524 : INTEGER :: ncol_in, nrow_out, ncol_out, unit_nr
525 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
526 : INTEGER, INTENT(IN), OPTIONAL :: nrow_offset, ncol_offset
527 :
528 : CHARACTER(LEN=*), PARAMETER :: routineN = 'truncate_fm'
529 :
530 : INTEGER :: col_idx_loc, dummy, handle, handle2, i_entry_rec, idx_col_first, idx_col_in, &
531 : idx_col_out, idx_col_sec, idx_row_in, ii, iproc, jj, ncol_block_in, ncol_block_out, &
532 : ncol_local_in, ncol_local_out, nprocs, nrow_block_in, nrow_block_out, nrow_local_in, &
533 : nrow_local_out, proc_send, row_idx_loc, send_pcol, send_prow
534 132 : INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, &
535 : num_entries_send
536 132 : INTEGER, DIMENSION(:), POINTER :: col_indices_in, col_indices_out, &
537 132 : row_indices_in, row_indices_out
538 : LOGICAL :: correct_ncol, correct_nrow
539 : TYPE(integ_mat_buffer_type), ALLOCATABLE, &
540 132 : DIMENSION(:) :: buffer_rec, buffer_send
541 : TYPE(mp_para_env_type), POINTER :: para_env_out
542 132 : TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
543 :
544 132 : CALL timeset(routineN, handle)
545 132 : CALL timeset(routineN//"_1_setup", handle2)
546 :
547 132 : correct_nrow = .FALSE.
548 132 : correct_ncol = .FALSE.
549 : !In case of truncation in the occupied space, we need to correct the interval of indices
550 132 : IF (PRESENT(nrow_offset)) THEN
551 88 : correct_nrow = .TRUE.
552 : END IF
553 132 : IF (PRESENT(ncol_offset)) THEN
554 44 : correct_ncol = .TRUE.
555 : END IF
556 :
557 132 : para_env_out => fm_out%matrix_struct%para_env
558 :
559 : CALL cp_fm_get_info(matrix=fm_out, &
560 : nrow_local=nrow_local_out, &
561 : ncol_local=ncol_local_out, &
562 : row_indices=row_indices_out, &
563 : col_indices=col_indices_out, &
564 : nrow_block=nrow_block_out, &
565 132 : ncol_block=ncol_block_out)
566 :
567 396 : ALLOCATE (num_entries_rec(0:para_env_out%num_pe - 1))
568 396 : ALLOCATE (num_entries_send(0:para_env_out%num_pe - 1))
569 :
570 132 : num_entries_rec(:) = 0
571 132 : num_entries_send(:) = 0
572 :
573 132 : dummy = 0
574 :
575 : CALL cp_fm_get_info(matrix=fm_in, &
576 : nrow_local=nrow_local_in, &
577 : ncol_local=ncol_local_in, &
578 : row_indices=row_indices_in, &
579 : col_indices=col_indices_in, &
580 : nrow_block=nrow_block_in, &
581 132 : ncol_block=ncol_block_in)
582 :
583 132 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
584 9 : WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_out%name, &
585 18 : fm_out%matrix_struct%nrow_global
586 9 : WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_out%name, &
587 18 : fm_out%matrix_struct%ncol_global
588 :
589 9 : WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_out%name, nrow_block_out
590 9 : WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_out%name, ncol_block_out
591 :
592 9 : WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_in%name, &
593 18 : fm_in%matrix_struct%nrow_global
594 9 : WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_in%name, &
595 18 : fm_in%matrix_struct%ncol_global
596 :
597 9 : WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_in%name, nrow_block_in
598 9 : WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_in%name, ncol_block_in
599 : END IF
600 :
601 : ! We find global indices in S with nrow_in and ncol_in for truncation
602 9240 : DO col_idx_loc = 1, ncol_local_in
603 9108 : idx_col_in = col_indices_in(col_idx_loc)
604 :
605 9108 : idx_col_first = (idx_col_in - 1)/ncol_in + 1
606 9108 : idx_col_sec = MOD(idx_col_in - 1, ncol_in) + 1
607 :
608 : ! If occupied orbitals are included, these have to be handled differently
609 : ! due to their reversed indexing
610 9108 : IF (correct_nrow) THEN
611 3284 : idx_col_first = idx_col_first - nrow_offset + 1
612 3284 : IF (idx_col_first <= 0) CYCLE
613 : ELSE
614 5824 : IF (idx_col_first > nrow_out) EXIT
615 : END IF
616 9108 : IF (correct_ncol) THEN
617 554 : idx_col_sec = idx_col_sec - ncol_offset + 1
618 554 : IF (idx_col_sec <= 0) CYCLE
619 : ELSE
620 8554 : IF (idx_col_sec > ncol_out) CYCLE
621 : END IF
622 :
623 8106 : idx_col_out = idx_col_sec + (idx_col_first - 1)*ncol_out
624 :
625 325508 : DO row_idx_loc = 1, nrow_local_in
626 317270 : idx_row_in = row_indices_in(row_idx_loc)
627 :
628 317270 : send_prow = fm_out%matrix_struct%g2p_row(idx_row_in)
629 317270 : send_pcol = fm_out%matrix_struct%g2p_col(idx_col_out)
630 :
631 317270 : proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
632 :
633 326378 : num_entries_send(proc_send) = num_entries_send(proc_send) + 1
634 :
635 : END DO
636 : END DO
637 :
638 132 : CALL timestop(handle2)
639 :
640 132 : CALL timeset(routineN//"_2_comm_entry_nums", handle2)
641 132 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
642 9 : WRITE (unit_nr, '(T2,A10,T13,A27)') 'BSE|DEBUG|', 'Communicating entry numbers'
643 : END IF
644 :
645 132 : CALL para_env_out%alltoall(num_entries_send, num_entries_rec, 1)
646 :
647 132 : CALL timestop(handle2)
648 :
649 132 : CALL timeset(routineN//"_3_alloc_buffer", handle2)
650 132 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
651 9 : WRITE (unit_nr, '(T2,A10,T13,A18)') 'BSE|DEBUG|', 'Allocating buffers'
652 : END IF
653 :
654 : ! Buffers for entries and their indices
655 660 : ALLOCATE (buffer_rec(0:para_env_out%num_pe - 1))
656 660 : ALLOCATE (buffer_send(0:para_env_out%num_pe - 1))
657 :
658 : ! allocate data message and corresponding indices
659 396 : DO iproc = 0, para_env_out%num_pe - 1
660 :
661 704 : ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc)))
662 317666 : buffer_rec(iproc)%msg = 0.0_dp
663 :
664 : END DO
665 :
666 396 : DO iproc = 0, para_env_out%num_pe - 1
667 :
668 704 : ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc)))
669 317666 : buffer_send(iproc)%msg = 0.0_dp
670 :
671 : END DO
672 :
673 396 : DO iproc = 0, para_env_out%num_pe - 1
674 :
675 704 : ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2))
676 635464 : buffer_rec(iproc)%indx = 0
677 :
678 : END DO
679 :
680 396 : DO iproc = 0, para_env_out%num_pe - 1
681 :
682 704 : ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2))
683 635464 : buffer_send(iproc)%indx = 0
684 :
685 : END DO
686 :
687 132 : CALL timestop(handle2)
688 :
689 132 : CALL timeset(routineN//"_4_buf_from_fmin_"//fm_out%name, handle2)
690 132 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
691 9 : WRITE (unit_nr, '(T2,A10,T13,A18,A10,A13)') 'BSE|DEBUG|', 'Writing data from ', fm_in%name, ' into buffers'
692 : END IF
693 :
694 396 : ALLOCATE (entry_counter(0:para_env_out%num_pe - 1))
695 132 : entry_counter(:) = 0
696 :
697 : ! Now we can write the actual data and indices to the send-buffer
698 9240 : DO col_idx_loc = 1, ncol_local_in
699 9108 : idx_col_in = col_indices_in(col_idx_loc)
700 :
701 9108 : idx_col_first = (idx_col_in - 1)/ncol_in + 1
702 9108 : idx_col_sec = MOD(idx_col_in - 1, ncol_in) + 1
703 :
704 : ! If occupied orbitals are included, these have to be handled differently
705 : ! due to their reversed indexing
706 9108 : IF (correct_nrow) THEN
707 3284 : idx_col_first = idx_col_first - nrow_offset + 1
708 3284 : IF (idx_col_first <= 0) CYCLE
709 : ELSE
710 5824 : IF (idx_col_first > nrow_out) EXIT
711 : END IF
712 9108 : IF (correct_ncol) THEN
713 554 : idx_col_sec = idx_col_sec - ncol_offset + 1
714 554 : IF (idx_col_sec <= 0) CYCLE
715 : ELSE
716 8554 : IF (idx_col_sec > ncol_out) CYCLE
717 : END IF
718 :
719 8106 : idx_col_out = idx_col_sec + (idx_col_first - 1)*ncol_out
720 :
721 325508 : DO row_idx_loc = 1, nrow_local_in
722 317270 : idx_row_in = row_indices_in(row_idx_loc)
723 :
724 317270 : send_prow = fm_out%matrix_struct%g2p_row(idx_row_in)
725 :
726 317270 : send_pcol = fm_out%matrix_struct%g2p_col(idx_col_out)
727 :
728 317270 : proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
729 317270 : entry_counter(proc_send) = entry_counter(proc_send) + 1
730 :
731 : buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
732 317270 : fm_in%local_data(row_idx_loc, col_idx_loc)
733 : !No need to create row_out, since it is identical to incoming
734 : !We dont change the RI index for any fm_mat_XX_BSE
735 317270 : buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = idx_row_in
736 326378 : buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = idx_col_out
737 :
738 : END DO
739 : END DO
740 :
741 1980 : ALLOCATE (req_array(1:para_env_out%num_pe, 4))
742 :
743 132 : CALL timestop(handle2)
744 :
745 132 : CALL timeset(routineN//"_5_comm_buffer", handle2)
746 132 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
747 9 : WRITE (unit_nr, '(T2,A10,T13,A21)') 'BSE|DEBUG|', 'Communicating buffers'
748 : END IF
749 :
750 : ! communicate the buffer
751 : CALL communicate_buffer(para_env_out, num_entries_rec, num_entries_send, buffer_rec, &
752 132 : buffer_send, req_array)
753 :
754 132 : CALL timestop(handle2)
755 :
756 132 : CALL timeset(routineN//"_6_buffer_to_fmout"//fm_out%name, handle2)
757 132 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
758 9 : WRITE (unit_nr, '(T2,A10,T13,A24,A10)') 'BSE|DEBUG|', 'Writing from buffers to ', fm_out%name
759 : END IF
760 :
761 : ! fill fm_out with the entries from buffer_rec, i.e. buffer_rec are parts of fm_in
762 132 : nprocs = para_env_out%num_pe
763 :
764 : !$OMP PARALLEL DO DEFAULT(NONE) &
765 : !$OMP SHARED(fm_out, nprocs, num_entries_rec, buffer_rec) &
766 132 : !$OMP PRIVATE(iproc, i_entry_rec, ii, jj)
767 : DO iproc = 0, nprocs - 1
768 : DO i_entry_rec = 1, num_entries_rec(iproc)
769 : ii = fm_out%matrix_struct%g2l_row(buffer_rec(iproc)%indx(i_entry_rec, 1))
770 : jj = fm_out%matrix_struct%g2l_col(buffer_rec(iproc)%indx(i_entry_rec, 2))
771 :
772 : fm_out%local_data(ii, jj) = fm_out%local_data(ii, jj) + buffer_rec(iproc)%msg(i_entry_rec)
773 : END DO
774 : END DO
775 : !$OMP END PARALLEL DO
776 :
777 132 : CALL timestop(handle2)
778 :
779 132 : CALL timeset(routineN//"_7_cleanup", handle2)
780 132 : IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
781 9 : WRITE (unit_nr, '(T2,A10,T13,A41)') 'BSE|DEBUG|', 'Starting cleanup of communication buffers'
782 : END IF
783 :
784 : !Clean up all the arrays from the communication process
785 396 : DO iproc = 0, para_env_out%num_pe - 1
786 264 : DEALLOCATE (buffer_rec(iproc)%msg)
787 264 : DEALLOCATE (buffer_rec(iproc)%indx)
788 264 : DEALLOCATE (buffer_send(iproc)%msg)
789 396 : DEALLOCATE (buffer_send(iproc)%indx)
790 : END DO
791 660 : DEALLOCATE (buffer_rec, buffer_send)
792 132 : DEALLOCATE (req_array)
793 132 : DEALLOCATE (entry_counter)
794 132 : DEALLOCATE (num_entries_rec, num_entries_send)
795 :
796 132 : CALL timestop(handle2)
797 132 : CALL timestop(handle)
798 :
799 1056 : END SUBROUTINE truncate_fm
800 :
801 : ! **************************************************************************************************
802 : !> \brief ...
803 : !> \param fm_mat_S_bar_ia_bse ...
804 : !> \param fm_mat_S_bar_ij_bse ...
805 : !> \param fm_mat_S_trunc ...
806 : !> \param fm_mat_S_ij_trunc ...
807 : !> \param fm_mat_S_ab_trunc ...
808 : !> \param fm_mat_Q_static_bse_gemm ...
809 : !> \param mp2_env ...
810 : ! **************************************************************************************************
811 40 : SUBROUTINE deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
812 : fm_mat_S_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
813 : fm_mat_Q_static_bse_gemm, mp2_env)
814 :
815 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_trunc, &
816 : fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, fm_mat_Q_static_bse_gemm
817 : TYPE(mp2_type) :: mp2_env
818 :
819 : CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_matrices_bse'
820 :
821 : INTEGER :: handle
822 :
823 40 : CALL timeset(routineN, handle)
824 :
825 40 : CALL cp_fm_release(fm_mat_S_bar_ia_bse)
826 40 : CALL cp_fm_release(fm_mat_S_bar_ij_bse)
827 40 : CALL cp_fm_release(fm_mat_S_trunc)
828 40 : CALL cp_fm_release(fm_mat_S_ij_trunc)
829 40 : CALL cp_fm_release(fm_mat_S_ab_trunc)
830 40 : CALL cp_fm_release(fm_mat_Q_static_bse_gemm)
831 40 : IF (mp2_env%bse%do_nto_analysis) THEN
832 4 : DEALLOCATE (mp2_env%bse%bse_nto_state_list_final)
833 : END IF
834 :
835 40 : CALL timestop(handle)
836 :
837 40 : END SUBROUTINE deallocate_matrices_bse
838 :
839 : ! **************************************************************************************************
840 : !> \brief Routine for computing the coefficients of the eigenvectors of the BSE matrix from a
841 : !> multiplication with the eigenvalues
842 : !> \param fm_work ...
843 : !> \param eig_vals ...
844 : !> \param beta ...
845 : !> \param gamma ...
846 : !> \param do_transpose ...
847 : ! **************************************************************************************************
848 88 : SUBROUTINE comp_eigvec_coeff_BSE(fm_work, eig_vals, beta, gamma, do_transpose)
849 :
850 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_work
851 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
852 : INTENT(IN) :: eig_vals
853 : REAL(KIND=dp), INTENT(IN) :: beta
854 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: gamma
855 : LOGICAL, INTENT(IN), OPTIONAL :: do_transpose
856 :
857 : CHARACTER(LEN=*), PARAMETER :: routineN = 'comp_eigvec_coeff_BSE'
858 :
859 : INTEGER :: handle, i_row_global, ii, j_col_global, &
860 : jj, ncol_local, nrow_local
861 44 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
862 : LOGICAL :: my_do_transpose
863 : REAL(KIND=dp) :: coeff, my_gamma
864 :
865 44 : CALL timeset(routineN, handle)
866 :
867 44 : IF (PRESENT(gamma)) THEN
868 44 : my_gamma = gamma
869 : ELSE
870 : my_gamma = 2.0_dp
871 : END IF
872 :
873 44 : IF (PRESENT(do_transpose)) THEN
874 44 : my_do_transpose = do_transpose
875 : ELSE
876 : my_do_transpose = .FALSE.
877 : END IF
878 :
879 : CALL cp_fm_get_info(matrix=fm_work, &
880 : nrow_local=nrow_local, &
881 : ncol_local=ncol_local, &
882 : row_indices=row_indices, &
883 44 : col_indices=col_indices)
884 :
885 44 : IF (my_do_transpose) THEN
886 3576 : DO jj = 1, ncol_local
887 3532 : j_col_global = col_indices(jj)
888 1531936 : DO ii = 1, nrow_local
889 1528360 : coeff = (eig_vals(j_col_global)**beta)/my_gamma
890 1531892 : fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff
891 : END DO
892 : END DO
893 : ELSE
894 0 : DO jj = 1, ncol_local
895 0 : DO ii = 1, nrow_local
896 0 : i_row_global = row_indices(ii)
897 0 : coeff = (eig_vals(i_row_global)**beta)/my_gamma
898 0 : fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff
899 : END DO
900 : END DO
901 : END IF
902 :
903 44 : CALL timestop(handle)
904 :
905 44 : END SUBROUTINE comp_eigvec_coeff_BSE
906 :
907 : ! **************************************************************************************************
908 : !> \brief Sorts excitation entries by ascending primary index, reordering the secondary index,
909 : !> the eigenvector coefficients and - open shell - the spin index alongside
910 : !> \param idx_prim Primary index of each entry; sorted in place and used as the sort key
911 : !> \param idx_sec Secondary index of each entry, reordered to follow idx_prim
912 : !> \param eigvec_entries Eigenvector coefficients of each entry, reordered to follow idx_prim
913 : !> \param idx_spin Optional spin index of each entry (open shell), reordered to follow idx_prim
914 : ! **************************************************************************************************
915 1414 : SUBROUTINE sort_excitations(idx_prim, idx_sec, eigvec_entries, idx_spin)
916 :
917 : INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim, idx_sec
918 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
919 : INTEGER, ALLOCATABLE, DIMENSION(:), OPTIONAL :: idx_spin
920 :
921 : CHARACTER(LEN=*), PARAMETER :: routineN = 'sort_excitations'
922 :
923 : INTEGER :: handle, ii, kk, num_entries, num_mults
924 1414 : INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim_work, idx_sec_work, &
925 1414 : idx_spin_work, tmp_index
926 : LOGICAL :: unique_entries
927 1414 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries_work
928 :
929 1414 : CALL timeset(routineN, handle)
930 :
931 1414 : num_entries = SIZE(idx_prim)
932 :
933 3726 : ALLOCATE (tmp_index(num_entries))
934 :
935 1414 : CALL sort(idx_prim, num_entries, tmp_index)
936 :
937 2312 : ALLOCATE (idx_sec_work(num_entries))
938 3726 : ALLOCATE (eigvec_entries_work(num_entries))
939 1528 : IF (PRESENT(idx_spin)) ALLOCATE (idx_spin_work(num_entries))
940 :
941 3970 : DO ii = 1, num_entries
942 2556 : idx_sec_work(ii) = idx_sec(tmp_index(ii))
943 2556 : eigvec_entries_work(ii) = eigvec_entries(tmp_index(ii))
944 3970 : IF (PRESENT(idx_spin)) idx_spin_work(ii) = idx_spin(tmp_index(ii))
945 : END DO
946 :
947 1414 : DEALLOCATE (tmp_index)
948 1414 : DEALLOCATE (idx_sec)
949 1414 : DEALLOCATE (eigvec_entries)
950 :
951 1414 : CALL MOVE_ALLOC(idx_sec_work, idx_sec)
952 1414 : CALL MOVE_ALLOC(eigvec_entries_work, eigvec_entries)
953 1414 : IF (PRESENT(idx_spin)) THEN
954 168 : DEALLOCATE (idx_spin)
955 168 : CALL MOVE_ALLOC(idx_spin_work, idx_spin)
956 : END IF
957 :
958 : !Now check for multiple entries in first idx to check necessity of sorting in second idx
959 1414 : CALL sort_unique(idx_prim, unique_entries)
960 1414 : IF (.NOT. unique_entries) THEN
961 768 : ALLOCATE (idx_prim_work(num_entries))
962 2048 : idx_prim_work(:) = idx_prim(:)
963 : ! Find duplicate entries in idx_prim
964 2048 : DO ii = 1, num_entries
965 1664 : IF (idx_prim_work(ii) == 0) CYCLE
966 6032 : num_mults = COUNT(idx_prim_work == idx_prim_work(ii))
967 990 : IF (num_mults > 1) THEN
968 : !Set all duplicate entries to 0
969 1714 : idx_prim_work(ii:ii + num_mults - 1) = 0
970 : !Start sorting in secondary index
971 1560 : ALLOCATE (idx_sec_work(num_mults))
972 1560 : ALLOCATE (eigvec_entries_work(num_mults))
973 1714 : idx_sec_work(:) = idx_sec(ii:ii + num_mults - 1)
974 1714 : eigvec_entries_work(:) = eigvec_entries(ii:ii + num_mults - 1)
975 520 : IF (PRESENT(idx_spin)) THEN
976 188 : ALLOCATE (idx_spin_work(num_mults))
977 332 : idx_spin_work(:) = idx_spin(ii:ii + num_mults - 1)
978 : END IF
979 1040 : ALLOCATE (tmp_index(num_mults))
980 520 : CALL sort(idx_sec_work, num_mults, tmp_index)
981 :
982 : !Now write newly sorted indices to original arrays
983 1714 : DO kk = ii, ii + num_mults - 1
984 1194 : idx_sec(kk) = idx_sec_work(kk - ii + 1)
985 1194 : eigvec_entries(kk) = eigvec_entries_work(tmp_index(kk - ii + 1))
986 1714 : IF (PRESENT(idx_spin)) idx_spin(kk) = idx_spin_work(tmp_index(kk - ii + 1))
987 : END DO
988 : !Deallocate work arrays
989 520 : DEALLOCATE (tmp_index)
990 520 : DEALLOCATE (idx_sec_work)
991 520 : DEALLOCATE (eigvec_entries_work)
992 520 : IF (PRESENT(idx_spin)) DEALLOCATE (idx_spin_work)
993 : END IF
994 2048 : idx_prim_work(ii) = idx_prim(ii)
995 : END DO
996 384 : DEALLOCATE (idx_prim_work)
997 : END IF
998 :
999 1414 : CALL timestop(handle)
1000 :
1001 4242 : END SUBROUTINE sort_excitations
1002 :
1003 : ! **************************************************************************************************
1004 : !> \brief Roughly estimates the needed runtime and memory during the BSE run
1005 : !> \param n_ov_joint ...
1006 : !> \param unit_nr ...
1007 : !> \param bse_abba ...
1008 : !> \param para_env ...
1009 : !> \param diag_runtime_est ...
1010 : ! **************************************************************************************************
1011 42 : SUBROUTINE estimate_BSE_resources(n_ov_joint, unit_nr, bse_abba, &
1012 : para_env, diag_runtime_est)
1013 :
1014 : INTEGER, INTENT(IN) :: n_ov_joint, unit_nr
1015 : LOGICAL :: bse_abba
1016 : TYPE(mp_para_env_type), POINTER :: para_env
1017 : REAL(KIND=dp) :: diag_runtime_est
1018 :
1019 : CHARACTER(LEN=*), PARAMETER :: routineN = 'estimate_BSE_resources'
1020 :
1021 : INTEGER :: handle, num_BSE_matrices
1022 : INTEGER(KIND=int_8) :: full_dim
1023 : REAL(KIND=dp) :: mem_est, mem_est_per_rank
1024 :
1025 42 : CALL timeset(routineN, handle)
1026 :
1027 : ! Number of matrices with size of A in TDA is 2 (A itself and W_ijab)
1028 42 : num_BSE_matrices = 2
1029 : ! With the full diagonalization of ABBA, we need several auxiliary matrices in the process
1030 : ! The maximum number is 2 + 2 + 6 (additional B and C matrix as well as 6 matrices to create C)
1031 42 : IF (bse_abba) THEN
1032 22 : num_BSE_matrices = 10
1033 : END IF
1034 :
1035 42 : full_dim = INT(n_ov_joint, KIND=int_8)**2*INT(num_BSE_matrices, KIND=int_8)
1036 42 : mem_est = REAL(8*full_dim, KIND=dp)/REAL(1024**3, KIND=dp)
1037 42 : mem_est_per_rank = REAL(mem_est/para_env%num_pe, KIND=dp)
1038 :
1039 42 : IF (unit_nr > 0) THEN
1040 : ! WRITE (unit_nr, '(T2,A4,T7,A40,T68,F13.3)') 'BSE|', 'Total peak memory estimate from BSE [GB]', &
1041 : ! mem_est
1042 21 : WRITE (unit_nr, '(T2,A4,T7,A40,T68,ES13.3)') 'BSE|', 'Total peak memory estimate from BSE [GB]', &
1043 42 : mem_est
1044 21 : WRITE (unit_nr, '(T2,A4,T7,A47,T68,F13.3)') 'BSE|', 'Peak memory estimate per MPI rank from BSE [GB]', &
1045 42 : mem_est_per_rank
1046 21 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1047 : END IF
1048 : ! Rough estimation of diagonalization runtimes. Baseline was a full BSE Naphthalene
1049 : ! run with 11000x11000 entries in A/B/C, which took 10s on 32 ranks
1050 : diag_runtime_est = REAL(INT(n_ov_joint, KIND=int_8)/11000_int_8, KIND=dp)**3* &
1051 42 : 10*32/REAL(para_env%num_pe, KIND=dp)
1052 :
1053 42 : CALL timestop(handle)
1054 :
1055 42 : END SUBROUTINE estimate_BSE_resources
1056 :
1057 : ! **************************************************************************************************
1058 : !> \brief Filters eigenvector entries above a given threshold to describe excitations in the
1059 : !> singleparticle basis
1060 : !> \param fm_eigvec ...
1061 : !> \param idx_homo ...
1062 : !> \param idx_virt ...
1063 : !> \param eigvec_entries ...
1064 : !> \param i_exc ...
1065 : !> \param virtual ...
1066 : !> \param num_entries ...
1067 : !> \param mp2_env ...
1068 : !> \param offsets ...
1069 : !> \param virtual_per_spin ...
1070 : !> \param idx_spin ...
1071 : ! **************************************************************************************************
1072 1414 : SUBROUTINE filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, &
1073 : i_exc, virtual, num_entries, mp2_env, &
1074 1414 : offsets, virtual_per_spin, idx_spin)
1075 :
1076 : TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec
1077 : INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt
1078 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
1079 : INTEGER :: i_exc, virtual, num_entries
1080 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1081 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: offsets, virtual_per_spin
1082 : INTEGER, ALLOCATABLE, DIMENSION(:), OPTIONAL :: idx_spin
1083 :
1084 : CHARACTER(LEN=*), PARAMETER :: routineN = 'filter_eigvec_contrib'
1085 :
1086 : INTEGER :: eigvec_idx, handle, ii, iproc, isp, jj, &
1087 : kk, ksp, ncol_local, nrow_local, &
1088 : num_entries_local, r_local, v_local
1089 : INTEGER, ALLOCATABLE, DIMENSION(:) :: num_entries_to_comm
1090 1414 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1091 : REAL(KIND=dp) :: eigvec_entry
1092 : TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1093 1414 : DIMENSION(:) :: buffer_entries
1094 : TYPE(mp_para_env_type), POINTER :: para_env
1095 :
1096 1414 : CALL timeset(routineN, handle)
1097 :
1098 1414 : para_env => fm_eigvec%matrix_struct%para_env
1099 :
1100 : CALL cp_fm_get_info(matrix=fm_eigvec, &
1101 : nrow_local=nrow_local, &
1102 : ncol_local=ncol_local, &
1103 : row_indices=row_indices, &
1104 1414 : col_indices=col_indices)
1105 :
1106 4242 : ALLOCATE (num_entries_to_comm(0:para_env%num_pe - 1))
1107 1414 : num_entries_to_comm(:) = 0
1108 :
1109 83720 : DO jj = 1, ncol_local
1110 : !First check if i is localized on this proc
1111 82306 : IF (col_indices(jj) /= i_exc) THEN
1112 : CYCLE
1113 : END IF
1114 54976 : DO ii = 1, nrow_local
1115 52172 : eigvec_idx = row_indices(ii)
1116 52172 : eigvec_entry = fm_eigvec%local_data(ii, jj)
1117 134478 : IF (ABS(eigvec_entry) > mp2_env%bse%eps_x) THEN
1118 1278 : num_entries_to_comm(para_env%mepos) = num_entries_to_comm(para_env%mepos) + 1
1119 : END IF
1120 : END DO
1121 : END DO
1122 :
1123 : !Gather number of entries of other processes
1124 1414 : CALL para_env%sum(num_entries_to_comm)
1125 :
1126 1414 : num_entries_local = num_entries_to_comm(para_env%mepos)
1127 :
1128 7070 : ALLOCATE (buffer_entries(0:para_env%num_pe - 1))
1129 :
1130 4242 : DO iproc = 0, para_env%num_pe - 1
1131 7106 : ALLOCATE (buffer_entries(iproc)%msg(num_entries_to_comm(iproc)))
1132 7106 : ALLOCATE (buffer_entries(iproc)%indx(num_entries_to_comm(iproc), 3))
1133 5384 : buffer_entries(iproc)%msg = 0.0_dp
1134 20394 : buffer_entries(iproc)%indx = 0
1135 : END DO
1136 :
1137 : kk = 1
1138 83720 : DO jj = 1, ncol_local
1139 : !First check if i is localized on this proc
1140 82306 : IF (col_indices(jj) /= i_exc) THEN
1141 : CYCLE
1142 : END IF
1143 54976 : DO ii = 1, nrow_local
1144 52172 : eigvec_idx = row_indices(ii)
1145 52172 : eigvec_entry = fm_eigvec%local_data(ii, jj)
1146 134478 : IF (ABS(eigvec_entry) > mp2_env%bse%eps_x) THEN
1147 : ! Decode spin block from the joint row index (blocks are contiguous; sigma is the
1148 : ! largest offset strictly below eigvec_idx). offsets absent -> closed shell, sigma=1.
1149 1278 : isp = 1
1150 1278 : r_local = eigvec_idx
1151 1278 : v_local = virtual
1152 1278 : IF (PRESENT(offsets)) THEN
1153 213 : DO ksp = SIZE(offsets), 1, -1
1154 213 : IF (eigvec_idx > offsets(ksp)) THEN
1155 : isp = ksp
1156 : EXIT
1157 : END IF
1158 : END DO
1159 146 : r_local = eigvec_idx - offsets(isp)
1160 146 : v_local = virtual_per_spin(isp)
1161 : END IF
1162 1278 : buffer_entries(para_env%mepos)%indx(kk, 1) = (r_local - 1)/v_local + 1
1163 1278 : buffer_entries(para_env%mepos)%indx(kk, 2) = MOD(r_local - 1, v_local) + 1
1164 1278 : buffer_entries(para_env%mepos)%indx(kk, 3) = isp
1165 1278 : buffer_entries(para_env%mepos)%msg(kk) = eigvec_entry
1166 1278 : kk = kk + 1
1167 : END IF
1168 : END DO
1169 : END DO
1170 :
1171 4242 : DO iproc = 0, para_env%num_pe - 1
1172 2828 : CALL para_env%sum(buffer_entries(iproc)%msg)
1173 4242 : CALL para_env%sum(buffer_entries(iproc)%indx)
1174 : END DO
1175 :
1176 : !Now sum up gathered information
1177 4242 : num_entries = SUM(num_entries_to_comm)
1178 3726 : ALLOCATE (idx_homo(num_entries))
1179 2312 : ALLOCATE (idx_virt(num_entries))
1180 3726 : ALLOCATE (eigvec_entries(num_entries))
1181 1528 : IF (PRESENT(idx_spin)) ALLOCATE (idx_spin(num_entries))
1182 :
1183 1414 : kk = 1
1184 4242 : DO iproc = 0, para_env%num_pe - 1
1185 4242 : IF (num_entries_to_comm(iproc) /= 0) THEN
1186 4006 : DO ii = 1, num_entries_to_comm(iproc)
1187 2556 : idx_homo(kk) = buffer_entries(iproc)%indx(ii, 1)
1188 2556 : idx_virt(kk) = buffer_entries(iproc)%indx(ii, 2)
1189 2556 : IF (PRESENT(idx_spin)) idx_spin(kk) = buffer_entries(iproc)%indx(ii, 3)
1190 2556 : eigvec_entries(kk) = buffer_entries(iproc)%msg(ii)
1191 4006 : kk = kk + 1
1192 : END DO
1193 : END IF
1194 : END DO
1195 :
1196 : !Deallocate all the used arrays
1197 4242 : DO iproc = 0, para_env%num_pe - 1
1198 2828 : DEALLOCATE (buffer_entries(iproc)%msg)
1199 4242 : DEALLOCATE (buffer_entries(iproc)%indx)
1200 : END DO
1201 5656 : DEALLOCATE (buffer_entries)
1202 1414 : DEALLOCATE (num_entries_to_comm)
1203 1414 : NULLIFY (row_indices)
1204 1414 : NULLIFY (col_indices)
1205 :
1206 : !Now sort the results according to the involved singleparticle orbitals
1207 : ! (homo first, then virtual). idx_spin is payload, permuted alongside the entries.
1208 1414 : IF (PRESENT(idx_spin)) THEN
1209 168 : CALL sort_excitations(idx_homo, idx_virt, eigvec_entries, idx_spin)
1210 : ELSE
1211 1246 : CALL sort_excitations(idx_homo, idx_virt, eigvec_entries)
1212 : END IF
1213 :
1214 1414 : CALL timestop(handle)
1215 :
1216 1414 : END SUBROUTINE filter_eigvec_contrib
1217 :
1218 : ! **************************************************************************************************
1219 : !> \brief Spin-block layout for the open-shell (joint) BSE matrix: per-spin OV-pair counts and the
1220 : !> block offsets into the joint matrix of dimension n_ov_joint = sum_sigma homo*virtual.
1221 : !> \param homo_red per-spin (reduced) number of occupied levels
1222 : !> \param virt_red per-spin (reduced) number of virtual levels
1223 : !> \param n_ov per-spin OV-pair count (OUT)
1224 : !> \param offsets per-spin block offset into the joint matrix (OUT)
1225 : !> \param n_ov_joint total joint dimension (OUT)
1226 : ! **************************************************************************************************
1227 126 : SUBROUTINE get_bse_spin_block_layout(homo_red, virt_red, n_ov, offsets, n_ov_joint)
1228 : INTEGER, DIMENSION(:), INTENT(IN) :: homo_red, virt_red
1229 : INTEGER, DIMENSION(:), INTENT(OUT) :: n_ov, offsets
1230 : INTEGER, INTENT(OUT) :: n_ov_joint
1231 :
1232 : INTEGER :: isp
1233 :
1234 126 : n_ov_joint = 0
1235 288 : DO isp = 1, SIZE(homo_red)
1236 162 : offsets(isp) = n_ov_joint
1237 162 : n_ov(isp) = homo_red(isp)*virt_red(isp)
1238 288 : n_ov_joint = n_ov_joint + n_ov(isp)
1239 : END DO
1240 :
1241 126 : END SUBROUTINE get_bse_spin_block_layout
1242 :
1243 : ! **************************************************************************************************
1244 : !> \brief Determines indices within the given energy cutoffs and truncates Eigenvalues and matrices
1245 : !> \param fm_mat_S_ia_bse ...
1246 : !> \param fm_mat_S_ij_bse ...
1247 : !> \param fm_mat_S_ab_bse ...
1248 : !> \param fm_mat_S_trunc ...
1249 : !> \param fm_mat_S_ij_trunc ...
1250 : !> \param fm_mat_S_ab_trunc ...
1251 : !> \param Eigenval_scf ...
1252 : !> \param Eigenval ...
1253 : !> \param Eigenval_reduced ...
1254 : !> \param homo ...
1255 : !> \param virtual ...
1256 : !> \param dimen_RI ...
1257 : !> \param unit_nr ...
1258 : !> \param bse_lev_virt ...
1259 : !> \param homo_red ...
1260 : !> \param virt_red ...
1261 : !> \param mp2_env ...
1262 : !> \param window Absolute MO window to retain.
1263 : !> \param print_window Whether to print the cutoff summary.
1264 : ! **************************************************************************************************
1265 112 : SUBROUTINE truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
1266 : fm_mat_S_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
1267 56 : Eigenval_scf, Eigenval, Eigenval_reduced, &
1268 : homo, virtual, dimen_RI, unit_nr, &
1269 : bse_lev_virt, &
1270 : homo_red, virt_red, &
1271 : mp2_env, &
1272 : window, print_window)
1273 :
1274 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, &
1275 : fm_mat_S_ab_bse
1276 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_S_trunc, fm_mat_S_ij_trunc, &
1277 : fm_mat_S_ab_trunc
1278 : REAL(KIND=dp), DIMENSION(:) :: Eigenval_scf, Eigenval
1279 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval_reduced
1280 : INTEGER, INTENT(IN) :: homo, virtual, dimen_RI, unit_nr, &
1281 : bse_lev_virt
1282 : INTEGER, INTENT(OUT) :: homo_red, virt_red
1283 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1284 : TYPE(mo_window_type), INTENT(IN) :: window
1285 : LOGICAL, INTENT(IN) :: print_window
1286 :
1287 : CHARACTER(LEN=*), PARAMETER :: routineN = 'truncate_BSE_matrices'
1288 :
1289 : INTEGER :: handle, homo_incl, virt_incl
1290 : TYPE(cp_blacs_env_type), POINTER :: context
1291 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_ab, fm_struct_ia, fm_struct_ij
1292 : TYPE(mp_para_env_type), POINTER :: para_env
1293 :
1294 56 : CALL timeset(routineN, handle)
1295 :
1296 56 : homo_incl = window%first_mo
1297 56 : virt_incl = window%last_mo - homo
1298 56 : homo_red = homo - homo_incl + 1
1299 56 : virt_red = virt_incl
1300 :
1301 56 : IF (print_window .AND. unit_nr > 0) THEN
1302 20 : IF (mp2_env%bse%bse_cutoff_occ > 0) THEN
1303 17 : WRITE (unit_nr, '(T2,A4,T7,A29,T71,F10.3)') 'BSE|', 'Cutoff occupied orbitals [eV]', &
1304 34 : mp2_env%bse%bse_cutoff_occ*evolt
1305 : ELSE
1306 3 : WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No cutoff given for occupied orbitals'
1307 : END IF
1308 20 : IF (mp2_env%bse%bse_cutoff_empty > 0) THEN
1309 18 : WRITE (unit_nr, '(T2,A4,T7,A26,T71,F10.3)') 'BSE|', 'Cutoff empty orbitals [eV]', &
1310 36 : mp2_env%bse%bse_cutoff_empty*evolt
1311 : ELSE
1312 2 : WRITE (unit_nr, '(T2,A4,T7,A34)') 'BSE|', 'No cutoff given for empty orbitals'
1313 : END IF
1314 20 : WRITE (unit_nr, '(T2,A4,T7,A20,T71,I10)') 'BSE|', 'First occupied index', homo_incl
1315 20 : WRITE (unit_nr, '(T2,A4,T7,A32,T71,I10)') 'BSE|', 'Last empty index (not MO index!)', virt_incl
1316 20 : WRITE (unit_nr, '(T2,A4,T7,A35,T71,F10.3)') 'BSE|', 'Energy of first occupied index [eV]', &
1317 40 : Eigenval(homo_incl)*evolt
1318 20 : WRITE (unit_nr, '(T2,A4,T7,A31,T71,F10.3)') 'BSE|', 'Energy of last empty index [eV]', &
1319 40 : Eigenval(homo + virt_incl)*evolt
1320 20 : WRITE (unit_nr, '(T2,A4,T7,A54,T71,F10.3)') 'BSE|', &
1321 20 : 'Energy difference of first occupied index to HOMO [eV]', &
1322 40 : -(Eigenval(homo_incl) - Eigenval(homo))*evolt
1323 20 : WRITE (unit_nr, '(T2,A4,T7,A50,T71,F10.3)') 'BSE|', &
1324 20 : 'Energy difference of last empty index to LUMO [eV]', &
1325 40 : (Eigenval(homo + virt_incl) - Eigenval(homo + 1))*evolt
1326 20 : WRITE (unit_nr, '(T2,A4,T7,A35,T71,I10)') 'BSE|', 'Number of GW-corrected occupied MOs', &
1327 40 : mp2_env%ri_g0w0%corr_mos_occ
1328 20 : WRITE (unit_nr, '(T2,A4,T7,A32,T71,I10)') 'BSE|', 'Number of GW-corrected empty MOs', &
1329 40 : bse_lev_virt
1330 20 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1331 : END IF
1332 56 : IF (unit_nr > 0) THEN
1333 28 : IF (homo - homo_incl + 1 > mp2_env%ri_g0w0%corr_mos_occ) THEN
1334 0 : CPABORT("Number of GW-corrected occupied MOs too small for chosen BSE cutoff")
1335 : END IF
1336 28 : IF (virt_incl > bse_lev_virt) THEN
1337 0 : CPABORT("Number of GW-corrected virtual MOs too small for chosen BSE cutoff")
1338 : END IF
1339 : END IF
1340 : !Truncate full fm_S matrices
1341 : !Allocate new truncated matrices of proper size
1342 56 : para_env => fm_mat_S_ia_bse%matrix_struct%para_env
1343 56 : context => fm_mat_S_ia_bse%matrix_struct%context
1344 :
1345 56 : CALL cp_fm_struct_create(fm_struct_ia, para_env, context, dimen_RI, homo_red*virt_red)
1346 56 : CALL cp_fm_struct_create(fm_struct_ij, para_env, context, dimen_RI, homo_red*homo_red)
1347 56 : CALL cp_fm_struct_create(fm_struct_ab, para_env, context, dimen_RI, virt_red*virt_red)
1348 :
1349 56 : CALL cp_fm_create(fm_mat_S_trunc, fm_struct_ia, name="fm_S_trunc", set_zero=.TRUE.)
1350 56 : CALL cp_fm_create(fm_mat_S_ij_trunc, fm_struct_ij, name="fm_S_ij_trunc", set_zero=.TRUE.)
1351 56 : CALL cp_fm_create(fm_mat_S_ab_trunc, fm_struct_ab, name="fm_S_ab_trunc", set_zero=.TRUE.)
1352 :
1353 : !Copy parts of original matrices to truncated ones
1354 56 : IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN
1355 : !Truncate eigenvals
1356 132 : ALLOCATE (Eigenval_reduced(homo_red + virt_red))
1357 : ! Include USE_KS_ENERGIES input
1358 44 : IF (mp2_env%bse%use_ks_energies) THEN
1359 48 : Eigenval_reduced(:) = Eigenval_scf(homo_incl:homo + virt_incl)
1360 : ELSE
1361 646 : Eigenval_reduced(:) = Eigenval(homo_incl:homo + virt_incl)
1362 : END IF
1363 :
1364 : CALL truncate_fm(fm_mat_S_trunc, fm_mat_S_ia_bse, virtual, &
1365 : homo_red, virt_red, unit_nr, mp2_env, &
1366 44 : nrow_offset=homo_incl)
1367 : CALL truncate_fm(fm_mat_S_ij_trunc, fm_mat_S_ij_bse, homo, &
1368 : homo_red, homo_red, unit_nr, mp2_env, &
1369 44 : homo_incl, homo_incl)
1370 : CALL truncate_fm(fm_mat_S_ab_trunc, fm_mat_S_ab_bse, bse_lev_virt, &
1371 44 : virt_red, virt_red, unit_nr, mp2_env)
1372 :
1373 : ELSE
1374 12 : IF (unit_nr > 0) THEN
1375 6 : WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No truncation of BSE matrices applied'
1376 6 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1377 : END IF
1378 36 : ALLOCATE (Eigenval_reduced(homo_red + virt_red))
1379 : ! Include USE_KS_ENERGIES input
1380 12 : IF (mp2_env%bse%use_ks_energies) THEN
1381 0 : Eigenval_reduced(:) = Eigenval_scf(:)
1382 : ELSE
1383 600 : Eigenval_reduced(:) = Eigenval(:)
1384 : END IF
1385 : CALL cp_fm_to_fm_submat_general(fm_mat_S_ia_bse, fm_mat_S_trunc, dimen_RI, homo_red*virt_red, &
1386 12 : 1, 1, 1, 1, context)
1387 : CALL cp_fm_to_fm_submat_general(fm_mat_S_ij_bse, fm_mat_S_ij_trunc, dimen_RI, homo_red*homo_red, &
1388 12 : 1, 1, 1, 1, context)
1389 : CALL cp_fm_to_fm_submat_general(fm_mat_S_ab_bse, fm_mat_S_ab_trunc, dimen_RI, virt_red*virt_red, &
1390 12 : 1, 1, 1, 1, context)
1391 : END IF
1392 :
1393 56 : CALL cp_fm_struct_release(fm_struct_ia)
1394 56 : CALL cp_fm_struct_release(fm_struct_ij)
1395 56 : CALL cp_fm_struct_release(fm_struct_ab)
1396 :
1397 56 : NULLIFY (para_env)
1398 56 : NULLIFY (context)
1399 :
1400 56 : CALL timestop(handle)
1401 :
1402 56 : END SUBROUTINE truncate_BSE_matrices
1403 :
1404 : ! **************************************************************************************************
1405 : !> \brief ...
1406 : !> \param fm_eigvec ...
1407 : !> \param fm_eigvec_reshuffled ...
1408 : !> \param homo ...
1409 : !> \param virtual ...
1410 : !> \param n_exc ...
1411 : !> \param do_transpose ...
1412 : !> \param unit_nr ...
1413 : !> \param mp2_env ...
1414 : ! **************************************************************************************************
1415 900 : SUBROUTINE reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, &
1416 : unit_nr, mp2_env)
1417 :
1418 : TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec
1419 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_eigvec_reshuffled
1420 : INTEGER, INTENT(IN) :: homo, virtual, n_exc
1421 : LOGICAL, INTENT(IN) :: do_transpose
1422 : INTEGER, INTENT(IN) :: unit_nr
1423 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1424 :
1425 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reshuffle_eigvec'
1426 :
1427 : INTEGER :: handle, my_m_col, my_n_row
1428 : INTEGER, DIMENSION(4) :: reordering
1429 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_eigvec_col, &
1430 : fm_struct_eigvec_reshuffled
1431 : TYPE(cp_fm_type) :: fm_eigvec_col
1432 :
1433 300 : CALL timeset(routineN, handle)
1434 :
1435 : ! Define reordering:
1436 : ! (ia,11) to (a1,i1) for transposition
1437 : ! (ia,11) to (i1,a1) for default
1438 300 : IF (do_transpose) THEN
1439 50 : reordering = [2, 3, 1, 4]
1440 50 : my_n_row = virtual
1441 50 : my_m_col = homo
1442 : ELSE
1443 250 : reordering = [1, 3, 2, 4]
1444 250 : my_n_row = homo
1445 250 : my_m_col = virtual
1446 : END IF
1447 :
1448 : CALL cp_fm_struct_create(fm_struct_eigvec_col, &
1449 : fm_eigvec%matrix_struct%para_env, fm_eigvec%matrix_struct%context, &
1450 300 : homo*virtual, 1)
1451 : CALL cp_fm_struct_create(fm_struct_eigvec_reshuffled, &
1452 : fm_eigvec%matrix_struct%para_env, fm_eigvec%matrix_struct%context, &
1453 300 : my_n_row, my_m_col)
1454 :
1455 : ! Resort indices
1456 300 : CALL cp_fm_create(fm_eigvec_col, fm_struct_eigvec_col, name="BSE_column_vector")
1457 300 : CALL cp_fm_set_all(fm_eigvec_col, 0.0_dp)
1458 300 : CALL cp_fm_create(fm_eigvec_reshuffled, fm_struct_eigvec_reshuffled, name="BSE_reshuffled_eigenvector")
1459 300 : CALL cp_fm_set_all(fm_eigvec_reshuffled, 0.0_dp)
1460 : ! Fill matrix
1461 : CALL cp_fm_to_fm_submat(fm_eigvec, fm_eigvec_col, &
1462 : homo*virtual, 1, &
1463 : 1, n_exc, &
1464 300 : 1, 1)
1465 : ! Reshuffle
1466 : CALL fm_general_add_bse(fm_eigvec_reshuffled, fm_eigvec_col, 1.0_dp, &
1467 : virtual, 1, &
1468 : 1, 1, &
1469 300 : unit_nr, reordering, mp2_env)
1470 :
1471 300 : CALL cp_fm_release(fm_eigvec_col)
1472 300 : CALL cp_fm_struct_release(fm_struct_eigvec_col)
1473 300 : CALL cp_fm_struct_release(fm_struct_eigvec_reshuffled)
1474 :
1475 300 : CALL timestop(handle)
1476 :
1477 300 : END SUBROUTINE reshuffle_eigvec
1478 :
1479 : ! **************************************************************************************************
1480 : !> \brief Borrowed from the tddfpt module with slight adaptions
1481 : !> \param qs_env ...
1482 : !> \param mos ...
1483 : !> \param istate ...
1484 : !> \param info_approximation ...
1485 : !> \param stride ...
1486 : !> \param append_cube ...
1487 : !> \param print_section ...
1488 : ! **************************************************************************************************
1489 0 : SUBROUTINE print_bse_nto_cubes(qs_env, mos, istate, info_approximation, &
1490 : stride, append_cube, print_section)
1491 :
1492 : TYPE(qs_environment_type), POINTER :: qs_env
1493 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1494 : INTEGER, INTENT(IN) :: istate
1495 : CHARACTER(LEN=10) :: info_approximation
1496 : INTEGER, DIMENSION(:), POINTER :: stride
1497 : LOGICAL, INTENT(IN) :: append_cube
1498 : TYPE(section_vals_type), POINTER :: print_section
1499 :
1500 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_bse_nto_cubes'
1501 :
1502 : CHARACTER(LEN=default_path_length) :: filename, info_approx_trunc, &
1503 : my_pos_cube, title
1504 : INTEGER :: handle, i, iset, nmo, unit_nr_cube
1505 : LOGICAL :: mpi_io
1506 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1507 : TYPE(cell_type), POINTER :: cell
1508 : TYPE(cp_fm_type), POINTER :: mo_coeff
1509 : TYPE(cp_logger_type), POINTER :: logger
1510 : TYPE(dft_control_type), POINTER :: dft_control
1511 : TYPE(particle_list_type), POINTER :: particles
1512 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1513 : TYPE(pw_c1d_gs_type) :: wf_g
1514 : TYPE(pw_env_type), POINTER :: pw_env
1515 0 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1516 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1517 : TYPE(pw_r3d_rs_type) :: wf_r
1518 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1519 : TYPE(qs_subsys_type), POINTER :: subsys
1520 :
1521 0 : logger => cp_get_default_logger()
1522 0 : CALL timeset(routineN, handle)
1523 :
1524 0 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
1525 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1526 0 : CALL auxbas_pw_pool%create_pw(wf_r)
1527 0 : CALL auxbas_pw_pool%create_pw(wf_g)
1528 :
1529 0 : CALL get_qs_env(qs_env, subsys=subsys)
1530 0 : CALL qs_subsys_get(subsys, particles=particles)
1531 :
1532 0 : my_pos_cube = "REWIND"
1533 0 : IF (append_cube) THEN
1534 0 : my_pos_cube = "APPEND"
1535 : END IF
1536 :
1537 : CALL get_qs_env(qs_env=qs_env, &
1538 : atomic_kind_set=atomic_kind_set, &
1539 : qs_kind_set=qs_kind_set, &
1540 : cell=cell, &
1541 0 : particle_set=particle_set)
1542 :
1543 0 : DO iset = 1, 2
1544 0 : CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
1545 0 : DO i = 1, nmo
1546 : CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1547 0 : cell, dft_control, particle_set, pw_env)
1548 0 : IF (iset == 1) THEN
1549 0 : WRITE (filename, '(A6,I3.3,A5,I2.2,a11)') "_NEXC_", istate, "_NTO_", i, "_Hole_State"
1550 0 : ELSE IF (iset == 2) THEN
1551 0 : WRITE (filename, '(A6,I3.3,A5,I2.2,a15)') "_NEXC_", istate, "_NTO_", i, "_Particle_State"
1552 : END IF
1553 0 : info_approx_trunc = TRIM(ADJUSTL(info_approximation))
1554 0 : info_approx_trunc = info_approx_trunc(2:LEN_TRIM(info_approx_trunc) - 1)
1555 0 : filename = TRIM(info_approx_trunc)//TRIM(filename)
1556 0 : mpi_io = .TRUE.
1557 : unit_nr_cube = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
1558 : middle_name=TRIM(filename), file_position=my_pos_cube, &
1559 0 : log_filename=.FALSE., ignore_should_output=.TRUE., mpi_io=mpi_io)
1560 0 : IF (iset == 1) THEN
1561 0 : WRITE (title, *) "Natural Transition Orbital Hole State", i
1562 0 : ELSE IF (iset == 2) THEN
1563 0 : WRITE (title, *) "Natural Transition Orbital Particle State", i
1564 : END IF
1565 0 : CALL cp_pw_to_cube(wf_r, unit_nr_cube, title, particles=particles, stride=stride, mpi_io=mpi_io)
1566 : CALL cp_print_key_finished_output(unit_nr_cube, logger, print_section, '', &
1567 0 : ignore_should_output=.TRUE., mpi_io=mpi_io)
1568 : END DO
1569 : END DO
1570 :
1571 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1572 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1573 :
1574 0 : CALL timestop(handle)
1575 0 : END SUBROUTINE print_bse_nto_cubes
1576 :
1577 : ! **************************************************************************************************
1578 : !> \brief Checks BSE input section and adapts them if necessary
1579 : !> \param homo ...
1580 : !> \param virtual ...
1581 : !> \param unit_nr ...
1582 : !> \param mp2_env ...
1583 : !> \param qs_env ...
1584 : ! **************************************************************************************************
1585 48 : SUBROUTINE adapt_BSE_input_params(homo, virtual, unit_nr, mp2_env, qs_env)
1586 :
1587 : INTEGER, INTENT(IN) :: homo, virtual, unit_nr
1588 : TYPE(mp2_type) :: mp2_env
1589 : TYPE(qs_environment_type), POINTER :: qs_env
1590 :
1591 : CHARACTER(LEN=*), PARAMETER :: routineN = 'adapt_BSE_input_params'
1592 :
1593 : INTEGER :: handle, i, j, n, n_exc_max, &
1594 : ndim_periodic_cell, &
1595 : ndim_periodic_poisson, &
1596 : num_state_list_exceptions
1597 : TYPE(cell_type), POINTER :: cell_ref
1598 : TYPE(pw_env_type), POINTER :: pw_env
1599 : TYPE(pw_poisson_type), POINTER :: poisson_env
1600 :
1601 48 : CALL timeset(routineN, handle)
1602 : ! Get environment infos for later usage
1603 48 : NULLIFY (pw_env, cell_ref, poisson_env)
1604 48 : CALL get_qs_env(qs_env, pw_env=pw_env, cell_ref=cell_ref)
1605 48 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
1606 192 : ndim_periodic_poisson = COUNT(poisson_env%parameters%periodic == 1)
1607 192 : ndim_periodic_cell = SUM(cell_ref%perd(1:3)) ! Borrowed from cell_methods.F/write_cell_low
1608 :
1609 : ! The iterative solver computes the NUM_EXC_EN lowest excitations only
1610 48 : n_exc_max = homo*virtual
1611 48 : IF (mp2_env%bse%bse_diag_method == bse_iterdiag) THEN
1612 6 : n_exc_max = MIN(n_exc_max, mp2_env%bse%num_exc_en)
1613 : END IF
1614 :
1615 : ! Handle negative NUM_PRINT_EXC
1616 48 : IF (mp2_env%bse%num_print_exc < 0 .OR. &
1617 : mp2_env%bse%num_print_exc > n_exc_max) THEN
1618 12 : mp2_env%bse%num_print_exc = n_exc_max
1619 12 : IF (unit_nr > 0) THEN
1620 : CALL cp_hint(__LOCATION__, &
1621 : "Keyword NUM_PRINT_EXC is either negative or too large. "// &
1622 6 : "Printing all computed excitations.")
1623 : END IF
1624 : END IF
1625 :
1626 : ! Default to NUM_PRINT_EXC if too large or negative,
1627 : ! but only if NTOs are called - would be confusing for the user otherwise
1628 : ! Prepare and adapt user inputs for NTO analysis
1629 : ! Logic: Explicit state list overrides NUM_PRINT_EXC_NTOS
1630 : ! If only NUM_PRINT_EXC_NTOS is given, we write the array 1,...,NUM_PRINT_EXC_NTOS to
1631 : ! bse_nto_state_list
1632 48 : IF (mp2_env%bse%do_nto_analysis) THEN
1633 4 : IF (mp2_env%bse%explicit_nto_list) THEN
1634 0 : IF (mp2_env%bse%num_print_exc_ntos > 0) THEN
1635 0 : IF (unit_nr > 0) THEN
1636 : CALL cp_hint(__LOCATION__, &
1637 : "Keywords NUM_PRINT_EXC_NTOS and STATE_LIST are both given in input. "// &
1638 0 : "Overriding NUM_PRINT_EXC_NTOS.")
1639 : END IF
1640 : END IF
1641 : ! Check if all states are within the range
1642 : ! Count them and initialize new array afterwards
1643 0 : num_state_list_exceptions = 0
1644 0 : DO i = 1, SIZE(mp2_env%bse%bse_nto_state_list)
1645 0 : IF (mp2_env%bse%bse_nto_state_list(i) < 1 .OR. &
1646 0 : mp2_env%bse%bse_nto_state_list(i) > mp2_env%bse%num_print_exc) THEN
1647 0 : num_state_list_exceptions = num_state_list_exceptions + 1
1648 : END IF
1649 : END DO
1650 0 : IF (num_state_list_exceptions > 0) THEN
1651 0 : IF (unit_nr > 0) THEN
1652 : CALL cp_hint(__LOCATION__, &
1653 : "STATE_LIST contains indices outside the range of included excitation levels. "// &
1654 0 : "Ignoring these states.")
1655 : END IF
1656 : END IF
1657 0 : n = SIZE(mp2_env%bse%bse_nto_state_list) - num_state_list_exceptions
1658 0 : ALLOCATE (mp2_env%bse%bse_nto_state_list_final(n))
1659 0 : mp2_env%bse%bse_nto_state_list_final(:) = 0
1660 : i = 1
1661 0 : DO j = 1, SIZE(mp2_env%bse%bse_nto_state_list)
1662 0 : IF (mp2_env%bse%bse_nto_state_list(j) >= 1 .AND. &
1663 0 : mp2_env%bse%bse_nto_state_list(j) <= mp2_env%bse%num_print_exc) THEN
1664 0 : mp2_env%bse%bse_nto_state_list_final(i) = mp2_env%bse%bse_nto_state_list(j)
1665 0 : i = i + 1
1666 : END IF
1667 : END DO
1668 :
1669 0 : mp2_env%bse%num_print_exc_ntos = SIZE(mp2_env%bse%bse_nto_state_list_final)
1670 : ELSE
1671 4 : IF (mp2_env%bse%num_print_exc_ntos > mp2_env%bse%num_print_exc .OR. &
1672 : mp2_env%bse%num_print_exc_ntos < 0) THEN
1673 4 : mp2_env%bse%num_print_exc_ntos = mp2_env%bse%num_print_exc
1674 : END IF
1675 12 : ALLOCATE (mp2_env%bse%bse_nto_state_list_final(mp2_env%bse%num_print_exc_ntos))
1676 104 : DO i = 1, mp2_env%bse%num_print_exc_ntos
1677 104 : mp2_env%bse%bse_nto_state_list_final(i) = i
1678 : END DO
1679 : END IF
1680 : END IF
1681 :
1682 : ! Takes care of triplet states, when oscillator strengths are 0
1683 48 : IF (mp2_env%bse%bse_spin_config /= 0 .AND. &
1684 : mp2_env%bse%eps_nto_osc_str > 0) THEN
1685 0 : IF (unit_nr > 0) THEN
1686 : CALL cp_warn(__LOCATION__, &
1687 : "Cannot apply EPS_OSC_STR for Triplet excitations. "// &
1688 0 : "Resetting EPS_OSC_STR to default.")
1689 : END IF
1690 0 : mp2_env%bse%eps_nto_osc_str = -1.0_dp
1691 : END IF
1692 :
1693 : ! Take care of number for computed exciton descriptors
1694 48 : IF (mp2_env%bse%num_print_exc_descr < 0 .OR. &
1695 : mp2_env%bse%num_print_exc_descr > mp2_env%bse%num_print_exc) THEN
1696 4 : IF (unit_nr > 0) THEN
1697 : CALL cp_hint(__LOCATION__, &
1698 : "Keyword NUM_PRINT_EXC_DESCR is either negative or too large. "// &
1699 2 : "Printing exciton descriptors up to NUM_PRINT_EXC.")
1700 : END IF
1701 4 : mp2_env%bse%num_print_exc_descr = mp2_env%bse%num_print_exc
1702 : END IF
1703 :
1704 : ! Handle screening factor options
1705 48 : IF (mp2_env%BSE%screening_factor > 0.0_dp) THEN
1706 2 : IF (mp2_env%BSE%screening_method /= bse_screening_alpha) THEN
1707 0 : IF (unit_nr > 0) THEN
1708 : CALL cp_warn(__LOCATION__, &
1709 : "Screening factor is only supported for &SCREENING_IN_W ALPHA. "// &
1710 0 : "Resetting SCREENING_IN_W to ALPHA.")
1711 : END IF
1712 0 : mp2_env%BSE%screening_method = bse_screening_alpha
1713 : END IF
1714 2 : IF (mp2_env%BSE%screening_factor > 1.0_dp) THEN
1715 0 : IF (unit_nr > 0) THEN
1716 : CALL cp_warn(__LOCATION__, &
1717 0 : "Screening factor is larger than 1.0. ")
1718 : END IF
1719 : END IF
1720 : END IF
1721 :
1722 48 : IF (mp2_env%BSE%screening_factor < 0.0_dp .AND. &
1723 : mp2_env%BSE%screening_method == bse_screening_alpha) THEN
1724 0 : IF (unit_nr > 0) THEN
1725 : CALL cp_warn(__LOCATION__, &
1726 0 : "Screening factor is negative. Defaulting to 0.25")
1727 : END IF
1728 0 : mp2_env%BSE%screening_factor = 0.25_dp
1729 : END IF
1730 :
1731 48 : IF (mp2_env%BSE%screening_factor == 0.0_dp) THEN
1732 : ! Use RPA internally in this case
1733 0 : mp2_env%BSE%screening_method = bse_screening_rpa
1734 : END IF
1735 48 : IF (mp2_env%BSE%screening_factor == 1.0_dp) THEN
1736 : ! Use TDHF internally in this case
1737 0 : mp2_env%BSE%screening_method = bse_screening_tdhf
1738 : END IF
1739 :
1740 : ! Add warning for usage of KS energies
1741 48 : IF (mp2_env%bse%use_ks_energies) THEN
1742 2 : IF (unit_nr > 0) THEN
1743 : CALL cp_warn(__LOCATION__, &
1744 : "Using KS energies for BSE calculations. Therefore, no quantities "// &
1745 1 : "of the preceeding GW calculation enter the BSE.")
1746 : END IF
1747 : END IF
1748 :
1749 : ! Add warning if periodic calculation is invoked
1750 48 : IF (ndim_periodic_poisson /= 0) THEN
1751 0 : IF (unit_nr > 0) THEN
1752 : CALL cp_warn(__LOCATION__, &
1753 : "Poisson solver should be invoked by PERIODIC NONE. "// &
1754 : "The applied length gauge might give misleading results for "// &
1755 0 : "oscillator strengths.")
1756 : END IF
1757 : END IF
1758 48 : IF (ndim_periodic_cell /= 0) THEN
1759 0 : IF (unit_nr > 0) THEN
1760 : CALL cp_warn(__LOCATION__, &
1761 : "CELL in SUBSYS should be invoked with PERIODIC NONE. "// &
1762 : "The applied length gauge might give misleading results for "// &
1763 0 : "oscillator strengths.")
1764 : END IF
1765 : END IF
1766 :
1767 48 : CALL timestop(handle)
1768 48 : END SUBROUTINE adapt_BSE_input_params
1769 :
1770 : ! **************************************************************************************************
1771 :
1772 : ! **************************************************************************************************
1773 : !> \brief ...
1774 : !> \param fm_multipole_ai_trunc ...
1775 : !> \param fm_multipole_ij_trunc ...
1776 : !> \param fm_multipole_ab_trunc ...
1777 : !> \param qs_env ...
1778 : !> \param mo_coeff ...
1779 : !> \param rpoint ...
1780 : !> \param n_moments ...
1781 : !> \param homo_red ...
1782 : !> \param virtual_red ...
1783 : !> \param context_BSE ...
1784 : !> \param ispin spin channel whose mo_set supplies homo/nao (default 1); open-shell beta needs 2
1785 : ! **************************************************************************************************
1786 60 : SUBROUTINE get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, &
1787 60 : qs_env, mo_coeff, rpoint, n_moments, &
1788 : homo_red, virtual_red, context_BSE, ispin)
1789 :
1790 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
1791 : INTENT(INOUT) :: fm_multipole_ai_trunc, &
1792 : fm_multipole_ij_trunc, &
1793 : fm_multipole_ab_trunc
1794 : TYPE(qs_environment_type), POINTER :: qs_env
1795 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1796 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: rpoint
1797 : INTEGER, INTENT(IN) :: n_moments, homo_red, virtual_red
1798 : TYPE(cp_blacs_env_type), POINTER :: context_BSE
1799 : INTEGER, INTENT(IN), OPTIONAL :: ispin
1800 :
1801 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_multipoles_mo'
1802 :
1803 : INTEGER :: handle, idir, my_ispin, n_multipole, &
1804 : n_occ, n_virt, nao, nmo_mp2
1805 60 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
1806 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_mp_ab_trunc, fm_struct_mp_ai_trunc, &
1807 : fm_struct_mp_ij_trunc, fm_struct_multipoles_ao, fm_struct_nao_nmo, fm_struct_nmo_nmo
1808 : TYPE(cp_fm_type) :: fm_work
1809 60 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_multipole_per_dir
1810 60 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_multipole, matrix_s
1811 60 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1812 : TYPE(mp_para_env_type), POINTER :: para_env_BSE
1813 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1814 60 : POINTER :: sab_orb
1815 :
1816 60 : CALL timeset(routineN, handle)
1817 :
1818 60 : my_ispin = 1
1819 60 : IF (PRESENT(ispin)) my_ispin = ispin
1820 :
1821 : !First, we calculate the AO dipoles
1822 60 : NULLIFY (sab_orb, matrix_s)
1823 : CALL get_qs_env(qs_env, &
1824 : mos=mos, &
1825 : matrix_s=matrix_s, &
1826 60 : sab_orb=sab_orb)
1827 :
1828 : ! Use the same blacs environment as for the MO coefficients to ensure correct multiplication dbcsr x fm later on
1829 60 : fm_struct_multipoles_ao => mos(my_ispin)%mo_coeff%matrix_struct
1830 : ! BSE has different contexts and blacsenvs
1831 60 : para_env_BSE => context_BSE%para_env
1832 : ! Get size of multipole tensor
1833 60 : n_multipole = (6 + 11*n_moments + 6*n_moments**2 + n_moments**3)/6 - 1
1834 60 : NULLIFY (matrix_multipole)
1835 60 : CALL dbcsr_allocate_matrix_set(matrix_multipole, n_multipole)
1836 396 : ALLOCATE (fm_multipole_per_dir(n_multipole))
1837 276 : DO idir = 1, n_multipole
1838 216 : CALL dbcsr_init_p(matrix_multipole(idir)%matrix)
1839 : CALL dbcsr_create(matrix_multipole(idir)%matrix, name="ao_multipole", &
1840 216 : template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
1841 216 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_multipole(idir)%matrix, sab_orb)
1842 276 : CALL dbcsr_set(matrix_multipole(idir)%matrix, 0._dp)
1843 : END DO
1844 :
1845 60 : CALL get_reference_point(rpoint, qs_env=qs_env, reference=use_mom_ref_coac, ref_point=ref_point)
1846 :
1847 60 : CALL build_local_moment_matrix(qs_env, matrix_multipole, n_moments, ref_point=rpoint)
1848 :
1849 60 : NULLIFY (sab_orb)
1850 :
1851 : ! Now we transform them to MO
1852 : ! n_occ is the number of occupied MOs, nao the number of all AOs
1853 : ! Writing homo to n_occ instead if nmo,
1854 : ! takes care of ADDED_MOS, which would overwrite nmo of qs_env-mos, if invoked
1855 60 : CALL get_mo_set(mo_set=mos(my_ispin), homo=n_occ, nao=nao)
1856 : ! Takes into account removed nullspace values from SVD
1857 60 : nmo_mp2 = mo_coeff(1)%matrix_struct%ncol_global
1858 60 : n_virt = nmo_mp2 - n_occ
1859 :
1860 : ! At the end, we need four different layouts of matrices in this multiplication, e.g. for a dipole:
1861 : ! D_pq = full multipole matrix for occupied and unoccupied
1862 : ! Final result:D_pq= C_{mu p} <\mu|\vec{r}|\nu> C_{\nu q} EQ.I
1863 : ! \_______/ \___________/ \______/
1864 : ! fm_coeff matrix_multipole fm_coeff
1865 : ! (EQ.Ia) (EQ.Ib) (EQ.Ia)
1866 : ! Intermediate work matrices:
1867 : ! fm_work = <\mu|\vec{r}|\nu> C_{\nu q} EQ.II
1868 :
1869 : ! Struct for the full multipole matrix
1870 : CALL cp_fm_struct_create(fm_struct_nao_nmo, &
1871 : fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, &
1872 60 : nao, nmo_mp2)
1873 : CALL cp_fm_struct_create(fm_struct_nmo_nmo, &
1874 : fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, &
1875 60 : nmo_mp2, nmo_mp2)
1876 :
1877 : ! At the very end, we copy the multipoles corresponding to truncated BSE indices in i and a
1878 : CALL cp_fm_struct_create(fm_struct_mp_ai_trunc, para_env_BSE, &
1879 60 : context_BSE, virtual_red, homo_red)
1880 : CALL cp_fm_struct_create(fm_struct_mp_ij_trunc, para_env_BSE, &
1881 60 : context_BSE, homo_red, homo_red)
1882 : CALL cp_fm_struct_create(fm_struct_mp_ab_trunc, para_env_BSE, &
1883 60 : context_BSE, virtual_red, virtual_red)
1884 276 : DO idir = 1, n_multipole
1885 : CALL cp_fm_create(fm_multipole_ai_trunc(idir), matrix_struct=fm_struct_mp_ai_trunc, &
1886 216 : name="dipoles_mo_ai_trunc")
1887 216 : CALL cp_fm_set_all(fm_multipole_ai_trunc(idir), 0.0_dp)
1888 : CALL cp_fm_create(fm_multipole_ij_trunc(idir), matrix_struct=fm_struct_mp_ij_trunc, &
1889 216 : name="dipoles_mo_ij_trunc")
1890 216 : CALL cp_fm_set_all(fm_multipole_ij_trunc(idir), 0.0_dp)
1891 : CALL cp_fm_create(fm_multipole_ab_trunc(idir), matrix_struct=fm_struct_mp_ab_trunc, &
1892 216 : name="dipoles_mo_ab_trunc")
1893 276 : CALL cp_fm_set_all(fm_multipole_ab_trunc(idir), 0.0_dp)
1894 : END DO
1895 :
1896 : ! Need another temporary matrix to store intermediate result from right multiplication
1897 : ! D = C_{mu a} <\mu|\vec{r}|\nu> C_{\nu i}
1898 60 : CALL cp_fm_create(fm_work, matrix_struct=fm_struct_nao_nmo, name="multipole_work")
1899 60 : CALL cp_fm_set_all(fm_work, 0.0_dp)
1900 :
1901 276 : DO idir = 1, n_multipole
1902 : ! Create the full multipole matrix per direction
1903 216 : CALL cp_fm_create(fm_multipole_per_dir(idir), matrix_struct=fm_struct_nmo_nmo, name="multipoles_mo")
1904 216 : CALL cp_fm_set_all(fm_multipole_per_dir(idir), 0.0_dp)
1905 : ! Fill final (MO) multipole matrix
1906 : CALL cp_dbcsr_sm_fm_multiply(matrix_multipole(idir)%matrix, mo_coeff(1), &
1907 216 : fm_work, ncol=nmo_mp2)
1908 : ! Now obtain the multipoles by the final multiplication;
1909 : ! We do that inside the loop to obtain multipoles per axis for print
1910 216 : CALL parallel_gemm('T', 'N', nmo_mp2, nmo_mp2, nao, 1.0_dp, mo_coeff(1), fm_work, 0.0_dp, fm_multipole_per_dir(idir))
1911 :
1912 : ! Truncate full matrix to the BSE indices
1913 : ! D_ai
1914 : CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1915 : fm_multipole_ai_trunc(idir), &
1916 : virtual_red, &
1917 : homo_red, &
1918 : n_occ + 1, &
1919 : n_occ - homo_red + 1, &
1920 : 1, &
1921 : 1, &
1922 216 : fm_multipole_per_dir(idir)%matrix_struct%context)
1923 : ! D_ij
1924 : CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1925 : fm_multipole_ij_trunc(idir), &
1926 : homo_red, &
1927 : homo_red, &
1928 : n_occ - homo_red + 1, &
1929 : n_occ - homo_red + 1, &
1930 : 1, &
1931 : 1, &
1932 216 : fm_multipole_per_dir(idir)%matrix_struct%context)
1933 : ! D_ab
1934 : CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1935 : fm_multipole_ab_trunc(idir), &
1936 : virtual_red, &
1937 : virtual_red, &
1938 : n_occ + 1, &
1939 : n_occ + 1, &
1940 : 1, &
1941 : 1, &
1942 276 : fm_multipole_per_dir(idir)%matrix_struct%context)
1943 : END DO
1944 :
1945 : !Release matrices and structs
1946 60 : NULLIFY (fm_struct_multipoles_ao)
1947 60 : CALL cp_fm_struct_release(fm_struct_mp_ai_trunc)
1948 60 : CALL cp_fm_struct_release(fm_struct_mp_ij_trunc)
1949 60 : CALL cp_fm_struct_release(fm_struct_mp_ab_trunc)
1950 60 : CALL cp_fm_struct_release(fm_struct_nao_nmo)
1951 60 : CALL cp_fm_struct_release(fm_struct_nmo_nmo)
1952 276 : DO idir = 1, n_multipole
1953 276 : CALL cp_fm_release(fm_multipole_per_dir(idir))
1954 : END DO
1955 60 : DEALLOCATE (fm_multipole_per_dir)
1956 60 : CALL cp_fm_release(fm_work)
1957 60 : CALL dbcsr_deallocate_matrix_set(matrix_multipole)
1958 :
1959 60 : CALL timestop(handle)
1960 :
1961 180 : END SUBROUTINE get_multipoles_mo
1962 :
1963 : ! **************************************************************************************************
1964 : !> \brief Computes trace of form Tr{A^T B C} for exciton descriptors
1965 : !> \param fm_A Full Matrix, typically X or Y, in format homo x virtual
1966 : !> \param fm_B ...
1967 : !> \param fm_C ...
1968 : !> \param alpha ...
1969 : ! **************************************************************************************************
1970 11520 : SUBROUTINE trace_exciton_descr(fm_A, fm_B, fm_C, alpha)
1971 :
1972 : TYPE(cp_fm_type), INTENT(IN) :: fm_A, fm_B, fm_C
1973 : REAL(KIND=dp), INTENT(OUT) :: alpha
1974 :
1975 : CHARACTER(LEN=*), PARAMETER :: routineN = 'trace_exciton_descr'
1976 :
1977 : INTEGER :: handle, ncol_A, ncol_B, ncol_C, nrow_A, &
1978 : nrow_B, nrow_C
1979 : TYPE(cp_fm_type) :: fm_work_ia
1980 :
1981 1920 : CALL timeset(routineN, handle)
1982 :
1983 1920 : CALL cp_fm_create(fm_work_ia, fm_A%matrix_struct)
1984 1920 : CALL cp_fm_get_info(fm_A, nrow_global=nrow_A, ncol_global=ncol_A)
1985 1920 : CALL cp_fm_get_info(fm_B, nrow_global=nrow_B, ncol_global=ncol_B)
1986 1920 : CALL cp_fm_get_info(fm_C, nrow_global=nrow_C, ncol_global=ncol_C)
1987 :
1988 : ! Check matrix sizes
1989 1920 : CPASSERT(nrow_A == nrow_B .AND. ncol_A == ncol_C .AND. ncol_B == nrow_C)
1990 :
1991 1920 : CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
1992 :
1993 : CALL parallel_gemm("N", "N", nrow_A, ncol_A, nrow_C, 1.0_dp, &
1994 1920 : fm_B, fm_C, 0.0_dp, fm_work_ia)
1995 :
1996 1920 : CALL cp_fm_trace(fm_A, fm_work_ia, alpha)
1997 :
1998 1920 : CALL cp_fm_release(fm_work_ia)
1999 :
2000 1920 : CALL timestop(handle)
2001 :
2002 1920 : END SUBROUTINE trace_exciton_descr
2003 :
2004 : ! **************************************************************************************************
2005 : !> \brief Column-concatenate per-spin ia-slabs into the joint dimen_RI x n_ov_joint slab.
2006 : !> Sigma-block of spin isp occupies columns offsets(isp)+1 .. offsets(isp)+n_ov(isp).
2007 : !> fm_S_joint must be pre-created and zeroed by the caller.
2008 : !> \param fm_S_ia per-spin ia-slabs, shape (dimen_RI, n_ov(isp)) per spin
2009 : !> \param offsets per-spin column offsets into fm_S_joint (0-based)
2010 : !> \param n_ov per-spin OV-pair counts
2011 : !> \param dimen_RI RI auxiliary basis dimension (row count)
2012 : !> \param fm_S_joint pre-created output slab (dimen_RI x n_ov_joint)
2013 : ! **************************************************************************************************
2014 12 : SUBROUTINE assemble_joint_ov_slab(fm_S_ia, offsets, n_ov, dimen_RI, fm_S_joint)
2015 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_S_ia
2016 : INTEGER, DIMENSION(:), INTENT(IN) :: offsets, n_ov
2017 : INTEGER, INTENT(IN) :: dimen_RI
2018 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_S_joint
2019 :
2020 : CHARACTER(LEN=*), PARAMETER :: routineN = 'assemble_joint_ov_slab'
2021 :
2022 : INTEGER :: handle, isp
2023 :
2024 12 : CALL timeset(routineN, handle)
2025 36 : DO isp = 1, SIZE(fm_S_ia)
2026 36 : CALL cp_fm_to_fm_submat(fm_S_ia(isp), fm_S_joint, dimen_RI, n_ov(isp), 1, 1, 1, offsets(isp) + 1)
2027 : END DO
2028 12 : CALL timestop(handle)
2029 :
2030 12 : END SUBROUTINE assemble_joint_ov_slab
2031 :
2032 : END MODULE bse_util
|
