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