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