(git:ed6f26b)
Loading...
Searching...
No Matches
bse_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 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! **************************************************************************************************
13MODULE 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
97CONTAINS
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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
124 TYPE(cp_fm_type) :: fm_work
125
126 CALL timeset(routinen, handle)
127
128 CALL cp_fm_create(fm_mat_s_bar_ia_bse, fm_mat_s_ia_bse%matrix_struct)
129 CALL cp_fm_set_all(fm_mat_s_bar_ia_bse, 0.0_dp)
130
131 CALL cp_fm_create(fm_mat_s_bar_ij_bse, fm_mat_s_ij_bse%matrix_struct)
132 CALL cp_fm_set_all(fm_mat_s_bar_ij_bse, 0.0_dp)
133
134 CALL cp_fm_create(fm_work, fm_mat_q_static_bse_gemm%matrix_struct)
135 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 col_indices=col_indices)
143
144 DO jjb = 1, ncol_local
145 j_global = col_indices(jjb)
146 DO iib = 1, nrow_local
147 i_global = row_indices(iib)
148 IF (j_global == i_global .AND. i_global <= dimen_ri) THEN
149 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 CALL cp_fm_cholesky_decompose(matrix=fm_mat_q_static_bse_gemm, n=dimen_ri, info_out=info_chol)
156
157 IF (info_chol /= 0) THEN
158 CALL cp_abort(__location__, 'Cholesky decomposition failed for static polarization in BSE')
159 END IF
160
161 ! calculate [1+Q(i0)]^-1
162 CALL cp_fm_cholesky_invert(fm_mat_q_static_bse_gemm)
163
164 ! symmetrize the result
165 CALL cp_fm_uplo_to_full(fm_mat_q_static_bse_gemm, fm_work)
166
167 CALL parallel_gemm(transa="N", transb="N", m=dimen_ri, n=homo**2, k=dimen_ri, alpha=1.0_dp, &
168 matrix_a=fm_mat_q_static_bse_gemm, matrix_b=fm_mat_s_ij_bse, beta=0.0_dp, &
169 matrix_c=fm_mat_s_bar_ij_bse)
170
171 ! fm_mat_S_bar_ia_bse has a different blacs_env as fm_mat_S_ij_bse since we take
172 ! fm_mat_S_ia_bse from RPA. Therefore, we also need a different fm_mat_Q_static_bse_gemm
173 CALL parallel_gemm(transa="N", transb="N", m=dimen_ri, n=homo*virtual, k=dimen_ri, alpha=1.0_dp, &
174 matrix_a=fm_mat_q_static_bse_gemm, matrix_b=fm_mat_s_ia_bse, beta=0.0_dp, &
175 matrix_c=fm_mat_s_bar_ia_bse)
176
177 CALL cp_fm_release(fm_work)
178
179 CALL timestop(handle)
180
181 END SUBROUTINE
182
183! **************************************************************************************************
184!> \brief Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,ab
185!> to A_ia, jb which needs MPI communication.
186!> \param fm_out ...
187!> \param fm_in ...
188!> \param beta ...
189!> \param nrow_secidx_in ...
190!> \param ncol_secidx_in ...
191!> \param nrow_secidx_out ...
192!> \param ncol_secidx_out ...
193!> \param unit_nr ...
194!> \param reordering ...
195!> \param mp2_env ...
196! **************************************************************************************************
197 SUBROUTINE fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, &
198 nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env)
199
200 TYPE(cp_fm_type), INTENT(INOUT) :: fm_out
201 TYPE(cp_fm_type), INTENT(IN) :: fm_in
202 REAL(kind=dp) :: beta
203 INTEGER, INTENT(IN) :: nrow_secidx_in, ncol_secidx_in, &
204 nrow_secidx_out, ncol_secidx_out
205 INTEGER :: unit_nr
206 INTEGER, DIMENSION(4) :: reordering
207 TYPE(mp2_type), INTENT(IN) :: mp2_env
208
209 CHARACTER(LEN=*), PARAMETER :: routinen = 'fm_general_add_bse'
210
211 INTEGER :: col_idx_loc, dummy, handle, handle2, i_entry_rec, idx_col_out, idx_row_out, ii, &
212 iproc, jj, ncol_block_in, ncol_block_out, ncol_local_in, ncol_local_out, nprocs, &
213 nrow_block_in, nrow_block_out, nrow_local_in, nrow_local_out, proc_send, row_idx_loc, &
214 send_pcol, send_prow
215 INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, &
216 num_entries_send
217 INTEGER, DIMENSION(4) :: indices_in
218 INTEGER, DIMENSION(:), POINTER :: col_indices_in, col_indices_out, &
219 row_indices_in, row_indices_out
220 TYPE(integ_mat_buffer_type), ALLOCATABLE, &
221 DIMENSION(:) :: buffer_rec, buffer_send
222 TYPE(mp_para_env_type), POINTER :: para_env_out
223 TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
224
225 CALL timeset(routinen, handle)
226 CALL timeset(routinen//"_1_setup", handle2)
227
228 para_env_out => fm_out%matrix_struct%para_env
229 ! A_iajb
230 ! We start by moving data from local parts of W_ijab to the full matrix A_iajb using buffers
231 CALL cp_fm_get_info(matrix=fm_out, &
232 nrow_local=nrow_local_out, &
233 ncol_local=ncol_local_out, &
234 row_indices=row_indices_out, &
235 col_indices=col_indices_out, &
236 nrow_block=nrow_block_out, &
237 ncol_block=ncol_block_out)
238
239 ALLOCATE (num_entries_rec(0:para_env_out%num_pe - 1))
240 ALLOCATE (num_entries_send(0:para_env_out%num_pe - 1))
241
242 num_entries_rec(:) = 0
243 num_entries_send(:) = 0
244
245 dummy = 0
246
247 CALL cp_fm_get_info(matrix=fm_in, &
248 nrow_local=nrow_local_in, &
249 ncol_local=ncol_local_in, &
250 row_indices=row_indices_in, &
251 col_indices=col_indices_in, &
252 nrow_block=nrow_block_in, &
253 ncol_block=ncol_block_in)
254
255 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
256 WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_out%name, &
257 fm_out%matrix_struct%nrow_global
258 WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_out%name, &
259 fm_out%matrix_struct%ncol_global
260
261 WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_out%name, nrow_block_out
262 WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_out%name, ncol_block_out
263
264 WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_in%name, &
265 fm_in%matrix_struct%nrow_global
266 WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_in%name, &
267 fm_in%matrix_struct%ncol_global
268
269 WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_in%name, nrow_block_in
270 WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_in%name, ncol_block_in
271 END IF
272
273 ! Use scalapack wrapper to find process index in fm_out
274 ! To that end, we obtain the global index in fm_out from the level indices
275 indices_in(:) = 0
276 DO row_idx_loc = 1, nrow_local_in
277 indices_in(1) = (row_indices_in(row_idx_loc) - 1)/nrow_secidx_in + 1
278 indices_in(2) = mod(row_indices_in(row_idx_loc) - 1, nrow_secidx_in) + 1
279 DO col_idx_loc = 1, ncol_local_in
280 indices_in(3) = (col_indices_in(col_idx_loc) - 1)/ncol_secidx_in + 1
281 indices_in(4) = mod(col_indices_in(col_idx_loc) - 1, ncol_secidx_in) + 1
282
283 idx_row_out = indices_in(reordering(2)) + (indices_in(reordering(1)) - 1)*nrow_secidx_out
284 idx_col_out = indices_in(reordering(4)) + (indices_in(reordering(3)) - 1)*ncol_secidx_out
285
286 send_prow = fm_out%matrix_struct%g2p_row(idx_row_out)
287 send_pcol = fm_out%matrix_struct%g2p_col(idx_col_out)
288
289 proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
290
291 num_entries_send(proc_send) = num_entries_send(proc_send) + 1
292
293 END DO
294 END DO
295
296 CALL timestop(handle2)
297
298 CALL timeset(routinen//"_2_comm_entry_nums", handle2)
299 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
300 WRITE (unit_nr, '(T2,A10,T13,A27)') 'BSE|DEBUG|', 'Communicating entry numbers'
301 END IF
302
303 CALL para_env_out%alltoall(num_entries_send, num_entries_rec, 1)
304
305 CALL timestop(handle2)
306
307 CALL timeset(routinen//"_3_alloc_buffer", handle2)
308 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
309 WRITE (unit_nr, '(T2,A10,T13,A18)') 'BSE|DEBUG|', 'Allocating buffers'
310 END IF
311
312 ! Buffers for entries and their indices
313 ALLOCATE (buffer_rec(0:para_env_out%num_pe - 1))
314 ALLOCATE (buffer_send(0:para_env_out%num_pe - 1))
315
316 ! allocate data message and corresponding indices
317 DO iproc = 0, para_env_out%num_pe - 1
318
319 ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc)))
320 buffer_rec(iproc)%msg = 0.0_dp
321
322 END DO
323
324 DO iproc = 0, para_env_out%num_pe - 1
325
326 ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc)))
327 buffer_send(iproc)%msg = 0.0_dp
328
329 END DO
330
331 DO iproc = 0, para_env_out%num_pe - 1
332
333 ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2))
334 buffer_rec(iproc)%indx = 0
335
336 END DO
337
338 DO iproc = 0, para_env_out%num_pe - 1
339
340 ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2))
341 buffer_send(iproc)%indx = 0
342
343 END DO
344
345 CALL timestop(handle2)
346
347 CALL timeset(routinen//"_4_buf_from_fmin_"//fm_out%name, handle2)
348 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
349 WRITE (unit_nr, '(T2,A10,T13,A18,A10,A13)') 'BSE|DEBUG|', 'Writing data from ', fm_in%name, ' into buffers'
350 END IF
351
352 ALLOCATE (entry_counter(0:para_env_out%num_pe - 1))
353 entry_counter(:) = 0
354
355 ! Now we can write the actual data and indices to the send-buffer
356 DO row_idx_loc = 1, nrow_local_in
357 indices_in(1) = (row_indices_in(row_idx_loc) - 1)/nrow_secidx_in + 1
358 indices_in(2) = mod(row_indices_in(row_idx_loc) - 1, nrow_secidx_in) + 1
359 DO col_idx_loc = 1, ncol_local_in
360 indices_in(3) = (col_indices_in(col_idx_loc) - 1)/ncol_secidx_in + 1
361 indices_in(4) = mod(col_indices_in(col_idx_loc) - 1, ncol_secidx_in) + 1
362
363 idx_row_out = indices_in(reordering(2)) + (indices_in(reordering(1)) - 1)*nrow_secidx_out
364 idx_col_out = indices_in(reordering(4)) + (indices_in(reordering(3)) - 1)*ncol_secidx_out
365
366 send_prow = fm_out%matrix_struct%g2p_row(idx_row_out)
367 send_pcol = fm_out%matrix_struct%g2p_col(idx_col_out)
368
369 proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
370 entry_counter(proc_send) = entry_counter(proc_send) + 1
371
372 buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
373 fm_in%local_data(row_idx_loc, col_idx_loc)
374
375 buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = idx_row_out
376 buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = idx_col_out
377
378 END DO
379 END DO
380
381 ALLOCATE (req_array(1:para_env_out%num_pe, 4))
382
383 CALL timestop(handle2)
384
385 CALL timeset(routinen//"_5_comm_buffer", handle2)
386 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
387 WRITE (unit_nr, '(T2,A10,T13,A21)') 'BSE|DEBUG|', 'Communicating buffers'
388 END IF
389
390 ! communicate the buffer
391 CALL communicate_buffer(para_env_out, num_entries_rec, num_entries_send, buffer_rec, &
392 buffer_send, req_array)
393
394 CALL timestop(handle2)
395
396 CALL timeset(routinen//"_6_buffer_to_fmout"//fm_out%name, handle2)
397 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
398 WRITE (unit_nr, '(T2,A10,T13,A24,A10)') 'BSE|DEBUG|', 'Writing from buffers to ', fm_out%name
399 END IF
400
401 ! fill fm_out with the entries from buffer_rec, i.e. buffer_rec are parts of fm_in
402 nprocs = para_env_out%num_pe
403
404!$OMP PARALLEL DO DEFAULT(NONE) &
405!$OMP SHARED(fm_out, nprocs, num_entries_rec, buffer_rec, beta) &
406!$OMP PRIVATE(iproc, i_entry_rec, ii, jj)
407 DO iproc = 0, nprocs - 1
408 DO i_entry_rec = 1, num_entries_rec(iproc)
409 ii = fm_out%matrix_struct%g2l_row(buffer_rec(iproc)%indx(i_entry_rec, 1))
410 jj = fm_out%matrix_struct%g2l_col(buffer_rec(iproc)%indx(i_entry_rec, 2))
411
412 fm_out%local_data(ii, jj) = fm_out%local_data(ii, jj) + beta*buffer_rec(iproc)%msg(i_entry_rec)
413 END DO
414 END DO
415!$OMP END PARALLEL DO
416
417 CALL timestop(handle2)
418
419 CALL timeset(routinen//"_7_cleanup", handle2)
420 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
421 WRITE (unit_nr, '(T2,A10,T13,A41)') 'BSE|DEBUG|', 'Starting cleanup of communication buffers'
422 END IF
423
424 !Clean up all the arrays from the communication process
425 DO iproc = 0, para_env_out%num_pe - 1
426 DEALLOCATE (buffer_rec(iproc)%msg)
427 DEALLOCATE (buffer_rec(iproc)%indx)
428 DEALLOCATE (buffer_send(iproc)%msg)
429 DEALLOCATE (buffer_send(iproc)%indx)
430 END DO
431 DEALLOCATE (buffer_rec, buffer_send)
432 DEALLOCATE (req_array)
433 DEALLOCATE (entry_counter)
434 DEALLOCATE (num_entries_rec, num_entries_send)
435
436 CALL timestop(handle2)
437 CALL timestop(handle)
438
439 END SUBROUTINE fm_general_add_bse
440
441! **************************************************************************************************
442!> \brief Routine for truncating a full matrix as given by the energy cutoffs in the input file.
443!> Logic: Matrices have some dimension dimen_RI x nrow_in*ncol_in for the incoming (untruncated) matrix
444!> and dimen_RI x nrow_out*ncol_out for the truncated matrix. The truncation is done by resorting the indices
445!> via parallel communication.
446!> \param fm_out ...
447!> \param fm_in ...
448!> \param ncol_in ...
449!> \param nrow_out ...
450!> \param ncol_out ...
451!> \param unit_nr ...
452!> \param mp2_env ...
453!> \param nrow_offset ...
454!> \param ncol_offset ...
455! **************************************************************************************************
456 SUBROUTINE truncate_fm(fm_out, fm_in, ncol_in, &
457 nrow_out, ncol_out, unit_nr, mp2_env, &
458 nrow_offset, ncol_offset)
459
460 TYPE(cp_fm_type), INTENT(INOUT) :: fm_out
461 TYPE(cp_fm_type), INTENT(IN) :: fm_in
462 INTEGER :: ncol_in, nrow_out, ncol_out, unit_nr
463 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
464 INTEGER, INTENT(IN), OPTIONAL :: nrow_offset, ncol_offset
465
466 CHARACTER(LEN=*), PARAMETER :: routinen = 'truncate_fm'
467
468 INTEGER :: col_idx_loc, dummy, handle, handle2, i_entry_rec, idx_col_first, idx_col_in, &
469 idx_col_out, idx_col_sec, idx_row_in, ii, iproc, jj, ncol_block_in, ncol_block_out, &
470 ncol_local_in, ncol_local_out, nprocs, nrow_block_in, nrow_block_out, nrow_local_in, &
471 nrow_local_out, proc_send, row_idx_loc, send_pcol, send_prow
472 INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, &
473 num_entries_send
474 INTEGER, DIMENSION(:), POINTER :: col_indices_in, col_indices_out, &
475 row_indices_in, row_indices_out
476 LOGICAL :: correct_ncol, correct_nrow
477 TYPE(integ_mat_buffer_type), ALLOCATABLE, &
478 DIMENSION(:) :: buffer_rec, buffer_send
479 TYPE(mp_para_env_type), POINTER :: para_env_out
480 TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
481
482 CALL timeset(routinen, handle)
483 CALL timeset(routinen//"_1_setup", handle2)
484
485 correct_nrow = .false.
486 correct_ncol = .false.
487 !In case of truncation in the occupied space, we need to correct the interval of indices
488 IF (PRESENT(nrow_offset)) THEN
489 correct_nrow = .true.
490 END IF
491 IF (PRESENT(ncol_offset)) THEN
492 correct_ncol = .true.
493 END IF
494
495 para_env_out => fm_out%matrix_struct%para_env
496
497 CALL cp_fm_get_info(matrix=fm_out, &
498 nrow_local=nrow_local_out, &
499 ncol_local=ncol_local_out, &
500 row_indices=row_indices_out, &
501 col_indices=col_indices_out, &
502 nrow_block=nrow_block_out, &
503 ncol_block=ncol_block_out)
504
505 ALLOCATE (num_entries_rec(0:para_env_out%num_pe - 1))
506 ALLOCATE (num_entries_send(0:para_env_out%num_pe - 1))
507
508 num_entries_rec(:) = 0
509 num_entries_send(:) = 0
510
511 dummy = 0
512
513 CALL cp_fm_get_info(matrix=fm_in, &
514 nrow_local=nrow_local_in, &
515 ncol_local=ncol_local_in, &
516 row_indices=row_indices_in, &
517 col_indices=col_indices_in, &
518 nrow_block=nrow_block_in, &
519 ncol_block=ncol_block_in)
520
521 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
522 WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_out%name, &
523 fm_out%matrix_struct%nrow_global
524 WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_out%name, &
525 fm_out%matrix_struct%ncol_global
526
527 WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_out%name, nrow_block_out
528 WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_out%name, ncol_block_out
529
530 WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_in%name, &
531 fm_in%matrix_struct%nrow_global
532 WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_in%name, &
533 fm_in%matrix_struct%ncol_global
534
535 WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_in%name, nrow_block_in
536 WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_in%name, ncol_block_in
537 END IF
538
539 ! We find global indices in S with nrow_in and ncol_in for truncation
540 DO col_idx_loc = 1, ncol_local_in
541 idx_col_in = col_indices_in(col_idx_loc)
542
543 idx_col_first = (idx_col_in - 1)/ncol_in + 1
544 idx_col_sec = mod(idx_col_in - 1, ncol_in) + 1
545
546 ! If occupied orbitals are included, these have to be handled differently
547 ! due to their reversed indexing
548 IF (correct_nrow) THEN
549 idx_col_first = idx_col_first - nrow_offset + 1
550 IF (idx_col_first .LE. 0) cycle
551 ELSE
552 IF (idx_col_first > nrow_out) EXIT
553 END IF
554 IF (correct_ncol) THEN
555 idx_col_sec = idx_col_sec - ncol_offset + 1
556 IF (idx_col_sec .LE. 0) cycle
557 ELSE
558 IF (idx_col_sec > ncol_out) cycle
559 END IF
560
561 idx_col_out = idx_col_sec + (idx_col_first - 1)*ncol_out
562
563 DO row_idx_loc = 1, nrow_local_in
564 idx_row_in = row_indices_in(row_idx_loc)
565
566 send_prow = fm_out%matrix_struct%g2p_row(idx_row_in)
567 send_pcol = fm_out%matrix_struct%g2p_col(idx_col_out)
568
569 proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
570
571 num_entries_send(proc_send) = num_entries_send(proc_send) + 1
572
573 END DO
574 END DO
575
576 CALL timestop(handle2)
577
578 CALL timeset(routinen//"_2_comm_entry_nums", handle2)
579 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
580 WRITE (unit_nr, '(T2,A10,T13,A27)') 'BSE|DEBUG|', 'Communicating entry numbers'
581 END IF
582
583 CALL para_env_out%alltoall(num_entries_send, num_entries_rec, 1)
584
585 CALL timestop(handle2)
586
587 CALL timeset(routinen//"_3_alloc_buffer", handle2)
588 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
589 WRITE (unit_nr, '(T2,A10,T13,A18)') 'BSE|DEBUG|', 'Allocating buffers'
590 END IF
591
592 ! Buffers for entries and their indices
593 ALLOCATE (buffer_rec(0:para_env_out%num_pe - 1))
594 ALLOCATE (buffer_send(0:para_env_out%num_pe - 1))
595
596 ! allocate data message and corresponding indices
597 DO iproc = 0, para_env_out%num_pe - 1
598
599 ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc)))
600 buffer_rec(iproc)%msg = 0.0_dp
601
602 END DO
603
604 DO iproc = 0, para_env_out%num_pe - 1
605
606 ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc)))
607 buffer_send(iproc)%msg = 0.0_dp
608
609 END DO
610
611 DO iproc = 0, para_env_out%num_pe - 1
612
613 ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2))
614 buffer_rec(iproc)%indx = 0
615
616 END DO
617
618 DO iproc = 0, para_env_out%num_pe - 1
619
620 ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2))
621 buffer_send(iproc)%indx = 0
622
623 END DO
624
625 CALL timestop(handle2)
626
627 CALL timeset(routinen//"_4_buf_from_fmin_"//fm_out%name, handle2)
628 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
629 WRITE (unit_nr, '(T2,A10,T13,A18,A10,A13)') 'BSE|DEBUG|', 'Writing data from ', fm_in%name, ' into buffers'
630 END IF
631
632 ALLOCATE (entry_counter(0:para_env_out%num_pe - 1))
633 entry_counter(:) = 0
634
635 ! Now we can write the actual data and indices to the send-buffer
636 DO col_idx_loc = 1, ncol_local_in
637 idx_col_in = col_indices_in(col_idx_loc)
638
639 idx_col_first = (idx_col_in - 1)/ncol_in + 1
640 idx_col_sec = mod(idx_col_in - 1, ncol_in) + 1
641
642 ! If occupied orbitals are included, these have to be handled differently
643 ! due to their reversed indexing
644 IF (correct_nrow) THEN
645 idx_col_first = idx_col_first - nrow_offset + 1
646 IF (idx_col_first .LE. 0) cycle
647 ELSE
648 IF (idx_col_first > nrow_out) EXIT
649 END IF
650 IF (correct_ncol) THEN
651 idx_col_sec = idx_col_sec - ncol_offset + 1
652 IF (idx_col_sec .LE. 0) cycle
653 ELSE
654 IF (idx_col_sec > ncol_out) cycle
655 END IF
656
657 idx_col_out = idx_col_sec + (idx_col_first - 1)*ncol_out
658
659 DO row_idx_loc = 1, nrow_local_in
660 idx_row_in = row_indices_in(row_idx_loc)
661
662 send_prow = fm_out%matrix_struct%g2p_row(idx_row_in)
663
664 send_pcol = fm_out%matrix_struct%g2p_col(idx_col_out)
665
666 proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
667 entry_counter(proc_send) = entry_counter(proc_send) + 1
668
669 buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
670 fm_in%local_data(row_idx_loc, col_idx_loc)
671 !No need to create row_out, since it is identical to incoming
672 !We dont change the RI index for any fm_mat_XX_BSE
673 buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = idx_row_in
674 buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = idx_col_out
675
676 END DO
677 END DO
678
679 ALLOCATE (req_array(1:para_env_out%num_pe, 4))
680
681 CALL timestop(handle2)
682
683 CALL timeset(routinen//"_5_comm_buffer", handle2)
684 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
685 WRITE (unit_nr, '(T2,A10,T13,A21)') 'BSE|DEBUG|', 'Communicating buffers'
686 END IF
687
688 ! communicate the buffer
689 CALL communicate_buffer(para_env_out, num_entries_rec, num_entries_send, buffer_rec, &
690 buffer_send, req_array)
691
692 CALL timestop(handle2)
693
694 CALL timeset(routinen//"_6_buffer_to_fmout"//fm_out%name, handle2)
695 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
696 WRITE (unit_nr, '(T2,A10,T13,A24,A10)') 'BSE|DEBUG|', 'Writing from buffers to ', fm_out%name
697 END IF
698
699 ! fill fm_out with the entries from buffer_rec, i.e. buffer_rec are parts of fm_in
700 nprocs = para_env_out%num_pe
701
702!$OMP PARALLEL DO DEFAULT(NONE) &
703!$OMP SHARED(fm_out, nprocs, num_entries_rec, buffer_rec) &
704!$OMP PRIVATE(iproc, i_entry_rec, ii, jj)
705 DO iproc = 0, nprocs - 1
706 DO i_entry_rec = 1, num_entries_rec(iproc)
707 ii = fm_out%matrix_struct%g2l_row(buffer_rec(iproc)%indx(i_entry_rec, 1))
708 jj = fm_out%matrix_struct%g2l_col(buffer_rec(iproc)%indx(i_entry_rec, 2))
709
710 fm_out%local_data(ii, jj) = fm_out%local_data(ii, jj) + buffer_rec(iproc)%msg(i_entry_rec)
711 END DO
712 END DO
713!$OMP END PARALLEL DO
714
715 CALL timestop(handle2)
716
717 CALL timeset(routinen//"_7_cleanup", handle2)
718 IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
719 WRITE (unit_nr, '(T2,A10,T13,A41)') 'BSE|DEBUG|', 'Starting cleanup of communication buffers'
720 END IF
721
722 !Clean up all the arrays from the communication process
723 DO iproc = 0, para_env_out%num_pe - 1
724 DEALLOCATE (buffer_rec(iproc)%msg)
725 DEALLOCATE (buffer_rec(iproc)%indx)
726 DEALLOCATE (buffer_send(iproc)%msg)
727 DEALLOCATE (buffer_send(iproc)%indx)
728 END DO
729 DEALLOCATE (buffer_rec, buffer_send)
730 DEALLOCATE (req_array)
731 DEALLOCATE (entry_counter)
732 DEALLOCATE (num_entries_rec, num_entries_send)
733
734 CALL timestop(handle2)
735 CALL timestop(handle)
736
737 END SUBROUTINE truncate_fm
738
739! **************************************************************************************************
740!> \brief ...
741!> \param fm_mat_S_bar_ia_bse ...
742!> \param fm_mat_S_bar_ij_bse ...
743!> \param fm_mat_S_trunc ...
744!> \param fm_mat_S_ij_trunc ...
745!> \param fm_mat_S_ab_trunc ...
746!> \param fm_mat_Q_static_bse_gemm ...
747!> \param mp2_env ...
748! **************************************************************************************************
749 SUBROUTINE deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
750 fm_mat_S_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
751 fm_mat_Q_static_bse_gemm, mp2_env)
752
753 TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_s_bar_ia_bse, fm_mat_s_bar_ij_bse, fm_mat_s_trunc, &
754 fm_mat_s_ij_trunc, fm_mat_s_ab_trunc, fm_mat_q_static_bse_gemm
755 TYPE(mp2_type) :: mp2_env
756
757 CHARACTER(LEN=*), PARAMETER :: routinen = 'deallocate_matrices_bse'
758
759 INTEGER :: handle
760
761 CALL timeset(routinen, handle)
762
763 CALL cp_fm_release(fm_mat_s_bar_ia_bse)
764 CALL cp_fm_release(fm_mat_s_bar_ij_bse)
765 CALL cp_fm_release(fm_mat_s_trunc)
766 CALL cp_fm_release(fm_mat_s_ij_trunc)
767 CALL cp_fm_release(fm_mat_s_ab_trunc)
768 CALL cp_fm_release(fm_mat_q_static_bse_gemm)
769 IF (mp2_env%bse%do_nto_analysis) THEN
770 DEALLOCATE (mp2_env%bse%bse_nto_state_list_final)
771 END IF
772
773 CALL timestop(handle)
774
775 END SUBROUTINE deallocate_matrices_bse
776
777! **************************************************************************************************
778!> \brief Routine for computing the coefficients of the eigenvectors of the BSE matrix from a
779!> multiplication with the eigenvalues
780!> \param fm_work ...
781!> \param eig_vals ...
782!> \param beta ...
783!> \param gamma ...
784!> \param do_transpose ...
785! **************************************************************************************************
786 SUBROUTINE comp_eigvec_coeff_bse(fm_work, eig_vals, beta, gamma, do_transpose)
787
788 TYPE(cp_fm_type), INTENT(INOUT) :: fm_work
789 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
790 INTENT(IN) :: eig_vals
791 REAL(kind=dp), INTENT(IN) :: beta
792 REAL(kind=dp), INTENT(IN), OPTIONAL :: gamma
793 LOGICAL, INTENT(IN), OPTIONAL :: do_transpose
794
795 CHARACTER(LEN=*), PARAMETER :: routinen = 'comp_eigvec_coeff_BSE'
796
797 INTEGER :: handle, i_row_global, ii, j_col_global, &
798 jj, ncol_local, nrow_local
799 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
800 LOGICAL :: my_do_transpose
801 REAL(kind=dp) :: coeff, my_gamma
802
803 CALL timeset(routinen, handle)
804
805 IF (PRESENT(gamma)) THEN
806 my_gamma = gamma
807 ELSE
808 my_gamma = 2.0_dp
809 END IF
810
811 IF (PRESENT(do_transpose)) THEN
812 my_do_transpose = do_transpose
813 ELSE
814 my_do_transpose = .false.
815 END IF
816
817 CALL cp_fm_get_info(matrix=fm_work, &
818 nrow_local=nrow_local, &
819 ncol_local=ncol_local, &
820 row_indices=row_indices, &
821 col_indices=col_indices)
822
823 IF (my_do_transpose) THEN
824 DO jj = 1, ncol_local
825 j_col_global = col_indices(jj)
826 DO ii = 1, nrow_local
827 coeff = (eig_vals(j_col_global)**beta)/my_gamma
828 fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff
829 END DO
830 END DO
831 ELSE
832 DO jj = 1, ncol_local
833 DO ii = 1, nrow_local
834 i_row_global = row_indices(ii)
835 coeff = (eig_vals(i_row_global)**beta)/my_gamma
836 fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff
837 END DO
838 END DO
839 END IF
840
841 CALL timestop(handle)
842
843 END SUBROUTINE
844
845! **************************************************************************************************
846!> \brief ...
847!> \param idx_prim ...
848!> \param idx_sec ...
849!> \param eigvec_entries ...
850! **************************************************************************************************
851 SUBROUTINE sort_excitations(idx_prim, idx_sec, eigvec_entries)
852
853 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim, idx_sec
854 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
855
856 CHARACTER(LEN=*), PARAMETER :: routinen = 'sort_excitations'
857
858 INTEGER :: handle, ii, kk, num_entries, num_mults
859 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim_work, idx_sec_work, tmp_index
860 LOGICAL :: unique_entries
861 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries_work
862
863 CALL timeset(routinen, handle)
864
865 num_entries = SIZE(idx_prim)
866
867 ALLOCATE (tmp_index(num_entries))
868
869 CALL sort(idx_prim, num_entries, tmp_index)
870
871 ALLOCATE (idx_sec_work(num_entries))
872 ALLOCATE (eigvec_entries_work(num_entries))
873
874 DO ii = 1, num_entries
875 idx_sec_work(ii) = idx_sec(tmp_index(ii))
876 eigvec_entries_work(ii) = eigvec_entries(tmp_index(ii))
877 END DO
878
879 DEALLOCATE (tmp_index)
880 DEALLOCATE (idx_sec)
881 DEALLOCATE (eigvec_entries)
882
883 CALL move_alloc(idx_sec_work, idx_sec)
884 CALL move_alloc(eigvec_entries_work, eigvec_entries)
885
886 !Now check for multiple entries in first idx to check necessity of sorting in second idx
887 CALL sort_unique(idx_prim, unique_entries)
888 IF (.NOT. unique_entries) THEN
889 ALLOCATE (idx_prim_work(num_entries))
890 idx_prim_work(:) = idx_prim(:)
891 ! Find duplicate entries in idx_prim
892 DO ii = 1, num_entries
893 IF (idx_prim_work(ii) == 0) cycle
894 num_mults = count(idx_prim_work == idx_prim_work(ii))
895 IF (num_mults > 1) THEN
896 !Set all duplicate entries to 0
897 idx_prim_work(ii:ii + num_mults - 1) = 0
898 !Start sorting in secondary index
899 ALLOCATE (idx_sec_work(num_mults))
900 ALLOCATE (eigvec_entries_work(num_mults))
901 idx_sec_work(:) = idx_sec(ii:ii + num_mults - 1)
902 eigvec_entries_work(:) = eigvec_entries(ii:ii + num_mults - 1)
903 ALLOCATE (tmp_index(num_mults))
904 CALL sort(idx_sec_work, num_mults, tmp_index)
905
906 !Now write newly sorted indices to original arrays
907 DO kk = ii, ii + num_mults - 1
908 idx_sec(kk) = idx_sec_work(kk - ii + 1)
909 eigvec_entries(kk) = eigvec_entries_work(tmp_index(kk - ii + 1))
910 END DO
911 !Deallocate work arrays
912 DEALLOCATE (tmp_index)
913 DEALLOCATE (idx_sec_work)
914 DEALLOCATE (eigvec_entries_work)
915 END IF
916 idx_prim_work(ii) = idx_prim(ii)
917 END DO
918 DEALLOCATE (idx_prim_work)
919 END IF
920
921 CALL timestop(handle)
922
923 END SUBROUTINE sort_excitations
924
925! **************************************************************************************************
926!> \brief Roughly estimates the needed runtime and memory during the BSE run
927!> \param homo_red ...
928!> \param virtual_red ...
929!> \param unit_nr ...
930!> \param bse_abba ...
931!> \param para_env ...
932!> \param diag_runtime_est ...
933! **************************************************************************************************
934 SUBROUTINE estimate_bse_resources(homo_red, virtual_red, unit_nr, bse_abba, &
935 para_env, diag_runtime_est)
936
937 INTEGER :: homo_red, virtual_red, unit_nr
938 LOGICAL :: bse_abba
939 TYPE(mp_para_env_type), POINTER :: para_env
940 REAL(kind=dp) :: diag_runtime_est
941
942 CHARACTER(LEN=*), PARAMETER :: routinen = 'estimate_BSE_resources'
943
944 INTEGER :: handle, num_bse_matrices
945 INTEGER(KIND=int_8) :: full_dim
946 REAL(kind=dp) :: mem_est, mem_est_per_rank
947
948 CALL timeset(routinen, handle)
949
950 ! Number of matrices with size of A in TDA is 2 (A itself and W_ijab)
951 num_bse_matrices = 2
952 ! With the full diagonalization of ABBA, we need several auxiliary matrices in the process
953 ! The maximum number is 2 + 2 + 6 (additional B and C matrix as well as 6 matrices to create C)
954 IF (bse_abba) THEN
955 num_bse_matrices = 10
956 END IF
957
958 full_dim = (int(homo_red, kind=int_8)**2*int(virtual_red, kind=int_8)**2)*int(num_bse_matrices, kind=int_8)
959 mem_est = real(8*full_dim, kind=dp)/real(1024**3, kind=dp)
960 mem_est_per_rank = real(mem_est/para_env%num_pe, kind=dp)
961
962 IF (unit_nr > 0) THEN
963 ! WRITE (unit_nr, '(T2,A4,T7,A40,T68,F13.3)') 'BSE|', 'Total peak memory estimate from BSE [GB]', &
964 ! mem_est
965 WRITE (unit_nr, '(T2,A4,T7,A40,T68,ES13.3)') 'BSE|', 'Total peak memory estimate from BSE [GB]', &
966 mem_est
967 WRITE (unit_nr, '(T2,A4,T7,A47,T68,F13.3)') 'BSE|', 'Peak memory estimate per MPI rank from BSE [GB]', &
968 mem_est_per_rank
969 WRITE (unit_nr, '(T2,A4)') 'BSE|'
970 END IF
971 ! Rough estimation of diagonalization runtimes. Baseline was a full BSE Naphthalene
972 ! run with 11000x11000 entries in A/B/C, which took 10s on 32 ranks
973 diag_runtime_est = real(int(homo_red, kind=int_8)*int(virtual_red, kind=int_8)/11000_int_8, kind=dp)**3* &
974 10*32/real(para_env%num_pe, kind=dp)
975
976 CALL timestop(handle)
977
978 END SUBROUTINE estimate_bse_resources
979
980! **************************************************************************************************
981!> \brief Filters eigenvector entries above a given threshold to describe excitations in the
982!> singleparticle basis
983!> \param fm_eigvec ...
984!> \param idx_homo ...
985!> \param idx_virt ...
986!> \param eigvec_entries ...
987!> \param i_exc ...
988!> \param virtual ...
989!> \param num_entries ...
990!> \param mp2_env ...
991! **************************************************************************************************
992 SUBROUTINE filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, &
993 i_exc, virtual, num_entries, mp2_env)
994
995 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec
996 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt
997 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
998 INTEGER :: i_exc, virtual, num_entries
999 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1000
1001 CHARACTER(LEN=*), PARAMETER :: routinen = 'filter_eigvec_contrib'
1002
1003 INTEGER :: eigvec_idx, handle, ii, iproc, jj, kk, &
1004 ncol_local, nrow_local, &
1005 num_entries_local
1006 INTEGER, ALLOCATABLE, DIMENSION(:) :: num_entries_to_comm
1007 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1008 REAL(kind=dp) :: eigvec_entry
1009 TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1010 DIMENSION(:) :: buffer_entries
1011 TYPE(mp_para_env_type), POINTER :: para_env
1012
1013 CALL timeset(routinen, handle)
1014
1015 para_env => fm_eigvec%matrix_struct%para_env
1016
1017 CALL cp_fm_get_info(matrix=fm_eigvec, &
1018 nrow_local=nrow_local, &
1019 ncol_local=ncol_local, &
1020 row_indices=row_indices, &
1021 col_indices=col_indices)
1022
1023 ALLOCATE (num_entries_to_comm(0:para_env%num_pe - 1))
1024 num_entries_to_comm(:) = 0
1025
1026 DO jj = 1, ncol_local
1027 !First check if i is localized on this proc
1028 IF (col_indices(jj) /= i_exc) THEN
1029 cycle
1030 END IF
1031 DO ii = 1, nrow_local
1032 eigvec_idx = row_indices(ii)
1033 eigvec_entry = fm_eigvec%local_data(ii, jj)/sqrt(2.0_dp)
1034 IF (abs(eigvec_entry) > mp2_env%bse%eps_x) THEN
1035 num_entries_to_comm(para_env%mepos) = num_entries_to_comm(para_env%mepos) + 1
1036 END IF
1037 END DO
1038 END DO
1039
1040 !Gather number of entries of other processes
1041 CALL para_env%sum(num_entries_to_comm)
1042
1043 num_entries_local = num_entries_to_comm(para_env%mepos)
1044
1045 ALLOCATE (buffer_entries(0:para_env%num_pe - 1))
1046
1047 DO iproc = 0, para_env%num_pe - 1
1048 ALLOCATE (buffer_entries(iproc)%msg(num_entries_to_comm(iproc)))
1049 ALLOCATE (buffer_entries(iproc)%indx(num_entries_to_comm(iproc), 2))
1050 buffer_entries(iproc)%msg = 0.0_dp
1051 buffer_entries(iproc)%indx = 0
1052 END DO
1053
1054 kk = 1
1055 DO jj = 1, ncol_local
1056 !First check if i is localized on this proc
1057 IF (col_indices(jj) /= i_exc) THEN
1058 cycle
1059 END IF
1060 DO ii = 1, nrow_local
1061 eigvec_idx = row_indices(ii)
1062 eigvec_entry = fm_eigvec%local_data(ii, jj)/sqrt(2.0_dp)
1063 IF (abs(eigvec_entry) > mp2_env%bse%eps_x) THEN
1064 buffer_entries(para_env%mepos)%indx(kk, 1) = (eigvec_idx - 1)/virtual + 1
1065 buffer_entries(para_env%mepos)%indx(kk, 2) = mod(eigvec_idx - 1, virtual) + 1
1066 buffer_entries(para_env%mepos)%msg(kk) = eigvec_entry
1067 kk = kk + 1
1068 END IF
1069 END DO
1070 END DO
1071
1072 DO iproc = 0, para_env%num_pe - 1
1073 CALL para_env%sum(buffer_entries(iproc)%msg)
1074 CALL para_env%sum(buffer_entries(iproc)%indx)
1075 END DO
1076
1077 !Now sum up gathered information
1078 num_entries = sum(num_entries_to_comm)
1079 ALLOCATE (idx_homo(num_entries))
1080 ALLOCATE (idx_virt(num_entries))
1081 ALLOCATE (eigvec_entries(num_entries))
1082
1083 kk = 1
1084 DO iproc = 0, para_env%num_pe - 1
1085 IF (num_entries_to_comm(iproc) /= 0) THEN
1086 DO ii = 1, num_entries_to_comm(iproc)
1087 idx_homo(kk) = buffer_entries(iproc)%indx(ii, 1)
1088 idx_virt(kk) = buffer_entries(iproc)%indx(ii, 2)
1089 eigvec_entries(kk) = buffer_entries(iproc)%msg(ii)
1090 kk = kk + 1
1091 END DO
1092 END IF
1093 END DO
1094
1095 !Deallocate all the used arrays
1096 DO iproc = 0, para_env%num_pe - 1
1097 DEALLOCATE (buffer_entries(iproc)%msg)
1098 DEALLOCATE (buffer_entries(iproc)%indx)
1099 END DO
1100 DEALLOCATE (buffer_entries)
1101 DEALLOCATE (num_entries_to_comm)
1102 NULLIFY (row_indices)
1103 NULLIFY (col_indices)
1104
1105 !Now sort the results according to the involved singleparticle orbitals
1106 ! (homo first, then virtual)
1107 CALL sort_excitations(idx_homo, idx_virt, eigvec_entries)
1108
1109 CALL timestop(handle)
1110
1111 END SUBROUTINE
1112
1113! **************************************************************************************************
1114!> \brief Reads cutoffs for BSE from mp2_env and compares to energies in Eigenval to extract
1115!> reduced homo/virtual and
1116!> \param Eigenval array (1d) with energies, can be e.g. from GW or DFT
1117!> \param homo Total number of occupied orbitals
1118!> \param virtual Total number of unoccupied orbitals
1119!> \param homo_red Total number of occupied orbitals to include after cutoff
1120!> \param virt_red Total number of unoccupied orbitals to include after ctuoff
1121!> \param homo_incl First occupied index to include after cutoff
1122!> \param virt_incl Last unoccupied index to include after cutoff
1123!> \param mp2_env ...
1124! **************************************************************************************************
1125 SUBROUTINE determine_cutoff_indices(Eigenval, &
1126 homo, virtual, &
1127 homo_red, virt_red, &
1128 homo_incl, virt_incl, &
1129 mp2_env)
1130
1131 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
1132 INTEGER, INTENT(IN) :: homo, virtual
1133 INTEGER, INTENT(OUT) :: homo_red, virt_red, homo_incl, virt_incl
1134 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1135
1136 CHARACTER(LEN=*), PARAMETER :: routinen = 'determine_cutoff_indices'
1137
1138 INTEGER :: handle, i_homo, j_virt
1139
1140 CALL timeset(routinen, handle)
1141 ! Determine index in homo and virtual for truncation
1142 ! Uses indices of outermost orbitals within energy range (-mp2_env%bse%bse_cutoff_occ,mp2_env%bse%bse_cutoff_empty)
1143 IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN
1144 IF (-mp2_env%bse%bse_cutoff_occ .LT. eigenval(1) - eigenval(homo) &
1145 .OR. mp2_env%bse%bse_cutoff_occ < 0) THEN
1146 homo_red = homo
1147 homo_incl = 1
1148 ELSE
1149 homo_incl = 1
1150 DO i_homo = 1, homo
1151 IF (eigenval(i_homo) - eigenval(homo) .GT. -mp2_env%bse%bse_cutoff_occ) THEN
1152 homo_incl = i_homo
1153 EXIT
1154 END IF
1155 END DO
1156 homo_red = homo - homo_incl + 1
1157 END IF
1158
1159 IF (mp2_env%bse%bse_cutoff_empty .GT. eigenval(homo + virtual) - eigenval(homo + 1) &
1160 .OR. mp2_env%bse%bse_cutoff_empty < 0) THEN
1161 virt_red = virtual
1162 virt_incl = virtual
1163 ELSE
1164 virt_incl = homo + 1
1165 DO j_virt = 1, virtual
1166 IF (eigenval(homo + j_virt) - eigenval(homo + 1) .GT. mp2_env%bse%bse_cutoff_empty) THEN
1167 virt_incl = j_virt - 1
1168 EXIT
1169 END IF
1170 END DO
1171 virt_red = virt_incl
1172 END IF
1173 ELSE
1174 homo_red = homo
1175 virt_red = virtual
1176 homo_incl = 1
1177 virt_incl = virtual
1178 END IF
1179
1180 CALL timestop(handle)
1181
1182 END SUBROUTINE
1183
1184! **************************************************************************************************
1185!> \brief Determines indices within the given energy cutoffs and truncates Eigenvalues and matrices
1186!> \param fm_mat_S_ia_bse ...
1187!> \param fm_mat_S_ij_bse ...
1188!> \param fm_mat_S_ab_bse ...
1189!> \param fm_mat_S_trunc ...
1190!> \param fm_mat_S_ij_trunc ...
1191!> \param fm_mat_S_ab_trunc ...
1192!> \param Eigenval_scf ...
1193!> \param Eigenval ...
1194!> \param Eigenval_reduced ...
1195!> \param homo ...
1196!> \param virtual ...
1197!> \param dimen_RI ...
1198!> \param unit_nr ...
1199!> \param bse_lev_virt ...
1200!> \param homo_red ...
1201!> \param virt_red ...
1202!> \param mp2_env ...
1203! **************************************************************************************************
1204 SUBROUTINE truncate_bse_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
1205 fm_mat_S_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
1206 Eigenval_scf, Eigenval, Eigenval_reduced, &
1207 homo, virtual, dimen_RI, unit_nr, &
1208 bse_lev_virt, &
1209 homo_red, virt_red, &
1210 mp2_env)
1211
1212 TYPE(cp_fm_type), INTENT(IN) :: fm_mat_s_ia_bse, fm_mat_s_ij_bse, &
1213 fm_mat_s_ab_bse
1214 TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_s_trunc, fm_mat_s_ij_trunc, &
1215 fm_mat_s_ab_trunc
1216 REAL(kind=dp), DIMENSION(:) :: eigenval_scf, eigenval
1217 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_reduced
1218 INTEGER, INTENT(IN) :: homo, virtual, dimen_ri, unit_nr, &
1219 bse_lev_virt
1220 INTEGER, INTENT(OUT) :: homo_red, virt_red
1221 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1222
1223 CHARACTER(LEN=*), PARAMETER :: routinen = 'truncate_BSE_matrices'
1224
1225 INTEGER :: handle, homo_incl, virt_incl
1226 TYPE(cp_blacs_env_type), POINTER :: context
1227 TYPE(cp_fm_struct_type), POINTER :: fm_struct_ab, fm_struct_ia, fm_struct_ij
1228 TYPE(mp_para_env_type), POINTER :: para_env
1229
1230 CALL timeset(routinen, handle)
1231
1232 ! Determine index in homo and virtual for truncation
1233 ! Uses indices of outermost orbitals within energy range (-mp2_env%bse%bse_cutoff_occ,mp2_env%bse%bse_cutoff_empty)
1234
1235 CALL determine_cutoff_indices(eigenval_scf, &
1236 homo, virtual, &
1237 homo_red, virt_red, &
1238 homo_incl, virt_incl, &
1239 mp2_env)
1240
1241 IF (unit_nr > 0) THEN
1242 IF (mp2_env%bse%bse_cutoff_occ > 0) THEN
1243 WRITE (unit_nr, '(T2,A4,T7,A29,T71,F10.3)') 'BSE|', 'Cutoff occupied orbitals [eV]', &
1244 mp2_env%bse%bse_cutoff_occ*evolt
1245 ELSE
1246 WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No cutoff given for occupied orbitals'
1247 END IF
1248 IF (mp2_env%bse%bse_cutoff_empty > 0) THEN
1249 WRITE (unit_nr, '(T2,A4,T7,A26,T71,F10.3)') 'BSE|', 'Cutoff empty orbitals [eV]', &
1250 mp2_env%bse%bse_cutoff_empty*evolt
1251 ELSE
1252 WRITE (unit_nr, '(T2,A4,T7,A34)') 'BSE|', 'No cutoff given for empty orbitals'
1253 END IF
1254 WRITE (unit_nr, '(T2,A4,T7,A20,T71,I10)') 'BSE|', 'First occupied index', homo_incl
1255 WRITE (unit_nr, '(T2,A4,T7,A32,T71,I10)') 'BSE|', 'Last empty index (not MO index!)', virt_incl
1256 WRITE (unit_nr, '(T2,A4,T7,A35,T71,F10.3)') 'BSE|', 'Energy of first occupied index [eV]', eigenval(homo_incl)*evolt
1257 WRITE (unit_nr, '(T2,A4,T7,A31,T71,F10.3)') 'BSE|', 'Energy of last empty index [eV]', eigenval(homo + virt_incl)*evolt
1258 WRITE (unit_nr, '(T2,A4,T7,A54,T71,F10.3)') 'BSE|', 'Energy difference of first occupied index to HOMO [eV]', &
1259 -(eigenval(homo_incl) - eigenval(homo))*evolt
1260 WRITE (unit_nr, '(T2,A4,T7,A50,T71,F10.3)') 'BSE|', 'Energy difference of last empty index to LUMO [eV]', &
1261 (eigenval(homo + virt_incl) - eigenval(homo + 1))*evolt
1262 WRITE (unit_nr, '(T2,A4,T7,A35,T71,I10)') 'BSE|', 'Number of GW-corrected occupied MOs', mp2_env%ri_g0w0%corr_mos_occ
1263 WRITE (unit_nr, '(T2,A4,T7,A32,T71,I10)') 'BSE|', 'Number of GW-corrected empty MOs', mp2_env%ri_g0w0%corr_mos_virt
1264 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1265 END IF
1266 IF (unit_nr > 0) THEN
1267 IF (homo - homo_incl + 1 > mp2_env%ri_g0w0%corr_mos_occ) THEN
1268 cpabort("Number of GW-corrected occupied MOs too small for chosen BSE cutoff")
1269 END IF
1270 IF (virt_incl > mp2_env%ri_g0w0%corr_mos_virt) THEN
1271 cpabort("Number of GW-corrected virtual MOs too small for chosen BSE cutoff")
1272 END IF
1273 END IF
1274 !Truncate full fm_S matrices
1275 !Allocate new truncated matrices of proper size
1276 para_env => fm_mat_s_ia_bse%matrix_struct%para_env
1277 context => fm_mat_s_ia_bse%matrix_struct%context
1278
1279 CALL cp_fm_struct_create(fm_struct_ia, para_env, context, dimen_ri, homo_red*virt_red)
1280 CALL cp_fm_struct_create(fm_struct_ij, para_env, context, dimen_ri, homo_red*homo_red)
1281 CALL cp_fm_struct_create(fm_struct_ab, para_env, context, dimen_ri, virt_red*virt_red)
1282
1283 CALL cp_fm_create(fm_mat_s_trunc, fm_struct_ia, "fm_S_trunc")
1284 CALL cp_fm_create(fm_mat_s_ij_trunc, fm_struct_ij, "fm_S_ij_trunc")
1285 CALL cp_fm_create(fm_mat_s_ab_trunc, fm_struct_ab, "fm_S_ab_trunc")
1286
1287 !Copy parts of original matrices to truncated ones
1288 IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN
1289 !Truncate eigenvals
1290 ALLOCATE (eigenval_reduced(homo_red + virt_red))
1291 ! Include USE_KS_ENERGIES input
1292 IF (mp2_env%bse%use_ks_energies) THEN
1293 eigenval_reduced(:) = eigenval_scf(homo_incl:homo + virt_incl)
1294 ELSE
1295 eigenval_reduced(:) = eigenval(homo_incl:homo + virt_incl)
1296 END IF
1297
1298 CALL truncate_fm(fm_mat_s_trunc, fm_mat_s_ia_bse, virtual, &
1299 homo_red, virt_red, unit_nr, mp2_env, &
1300 nrow_offset=homo_incl)
1301 CALL truncate_fm(fm_mat_s_ij_trunc, fm_mat_s_ij_bse, homo, &
1302 homo_red, homo_red, unit_nr, mp2_env, &
1303 homo_incl, homo_incl)
1304 CALL truncate_fm(fm_mat_s_ab_trunc, fm_mat_s_ab_bse, bse_lev_virt, &
1305 virt_red, virt_red, unit_nr, mp2_env)
1306
1307 ELSE
1308 IF (unit_nr > 0) THEN
1309 WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No truncation of BSE matrices applied'
1310 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1311 END IF
1312 ALLOCATE (eigenval_reduced(homo_red + virt_red))
1313 ! Include USE_KS_ENERGIES input
1314 IF (mp2_env%bse%use_ks_energies) THEN
1315 eigenval_reduced(:) = eigenval_scf(:)
1316 ELSE
1317 eigenval_reduced(:) = eigenval(:)
1318 END IF
1319 CALL cp_fm_to_fm_submat_general(fm_mat_s_ia_bse, fm_mat_s_trunc, dimen_ri, homo_red*virt_red, &
1320 1, 1, 1, 1, context)
1321 CALL cp_fm_to_fm_submat_general(fm_mat_s_ij_bse, fm_mat_s_ij_trunc, dimen_ri, homo_red*homo_red, &
1322 1, 1, 1, 1, context)
1323 CALL cp_fm_to_fm_submat_general(fm_mat_s_ab_bse, fm_mat_s_ab_trunc, dimen_ri, virt_red*virt_red, &
1324 1, 1, 1, 1, context)
1325 END IF
1326
1327 CALL cp_fm_struct_release(fm_struct_ia)
1328 CALL cp_fm_struct_release(fm_struct_ij)
1329 CALL cp_fm_struct_release(fm_struct_ab)
1330
1331 NULLIFY (para_env)
1332 NULLIFY (context)
1333
1334 CALL timestop(handle)
1335
1336 END SUBROUTINE truncate_bse_matrices
1337
1338! **************************************************************************************************
1339!> \brief ...
1340!> \param fm_eigvec ...
1341!> \param fm_eigvec_reshuffled ...
1342!> \param homo ...
1343!> \param virtual ...
1344!> \param n_exc ...
1345!> \param do_transpose ...
1346!> \param unit_nr ...
1347!> \param mp2_env ...
1348! **************************************************************************************************
1349 SUBROUTINE reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, &
1350 unit_nr, mp2_env)
1351
1352 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec
1353 TYPE(cp_fm_type), INTENT(INOUT) :: fm_eigvec_reshuffled
1354 INTEGER, INTENT(IN) :: homo, virtual, n_exc
1355 LOGICAL, INTENT(IN) :: do_transpose
1356 INTEGER, INTENT(IN) :: unit_nr
1357 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1358
1359 CHARACTER(LEN=*), PARAMETER :: routinen = 'reshuffle_eigvec'
1360
1361 INTEGER :: handle, my_m_col, my_n_row
1362 INTEGER, DIMENSION(4) :: reordering
1363 TYPE(cp_fm_struct_type), POINTER :: fm_struct_eigvec_col, &
1364 fm_struct_eigvec_reshuffled
1365 TYPE(cp_fm_type) :: fm_eigvec_col
1366
1367 CALL timeset(routinen, handle)
1368
1369 ! Define reordering:
1370 ! (ia,11) to (a1,i1) for transposition
1371 ! (ia,11) to (i1,a1) for default
1372 IF (do_transpose) THEN
1373 reordering = (/2, 3, 1, 4/)
1374 my_n_row = virtual
1375 my_m_col = homo
1376 ELSE
1377 reordering = (/1, 3, 2, 4/)
1378 my_n_row = homo
1379 my_m_col = virtual
1380 END IF
1381
1382 CALL cp_fm_struct_create(fm_struct_eigvec_col, &
1383 fm_eigvec%matrix_struct%para_env, fm_eigvec%matrix_struct%context, &
1384 homo*virtual, 1)
1385 CALL cp_fm_struct_create(fm_struct_eigvec_reshuffled, &
1386 fm_eigvec%matrix_struct%para_env, fm_eigvec%matrix_struct%context, &
1387 my_n_row, my_m_col)
1388
1389 ! Resort indices
1390 CALL cp_fm_create(fm_eigvec_col, fm_struct_eigvec_col, name="BSE_column_vector")
1391 CALL cp_fm_set_all(fm_eigvec_col, 0.0_dp)
1392 CALL cp_fm_create(fm_eigvec_reshuffled, fm_struct_eigvec_reshuffled, name="BSE_reshuffled_eigenvector")
1393 CALL cp_fm_set_all(fm_eigvec_reshuffled, 0.0_dp)
1394 ! Fill matrix
1395 CALL cp_fm_to_fm_submat(fm_eigvec, fm_eigvec_col, &
1396 homo*virtual, 1, &
1397 1, n_exc, &
1398 1, 1)
1399 ! Reshuffle
1400 CALL fm_general_add_bse(fm_eigvec_reshuffled, fm_eigvec_col, 1.0_dp, &
1401 virtual, 1, &
1402 1, 1, &
1403 unit_nr, reordering, mp2_env)
1404
1405 CALL cp_fm_release(fm_eigvec_col)
1406 CALL cp_fm_struct_release(fm_struct_eigvec_col)
1407 CALL cp_fm_struct_release(fm_struct_eigvec_reshuffled)
1408
1409 CALL timestop(handle)
1410
1411 END SUBROUTINE reshuffle_eigvec
1412
1413! **************************************************************************************************
1414!> \brief Borrowed from the tddfpt module with slight adaptions
1415!> \param qs_env ...
1416!> \param mos ...
1417!> \param istate ...
1418!> \param info_approximation ...
1419!> \param stride ...
1420!> \param append_cube ...
1421!> \param print_section ...
1422! **************************************************************************************************
1423 SUBROUTINE print_bse_nto_cubes(qs_env, mos, istate, info_approximation, &
1424 stride, append_cube, print_section)
1425
1426 TYPE(qs_environment_type), POINTER :: qs_env
1427 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1428 INTEGER, INTENT(IN) :: istate
1429 CHARACTER(LEN=10) :: info_approximation
1430 INTEGER, DIMENSION(:), POINTER :: stride
1431 LOGICAL, INTENT(IN) :: append_cube
1432 TYPE(section_vals_type), POINTER :: print_section
1433
1434 CHARACTER(LEN=*), PARAMETER :: routinen = 'print_bse_nto_cubes'
1435
1436 CHARACTER(LEN=default_path_length) :: filename, info_approx_trunc, &
1437 my_pos_cube, title
1438 INTEGER :: handle, i, iset, nmo, unit_nr_cube
1439 LOGICAL :: mpi_io
1440 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1441 TYPE(cell_type), POINTER :: cell
1442 TYPE(cp_fm_type), POINTER :: mo_coeff
1443 TYPE(cp_logger_type), POINTER :: logger
1444 TYPE(dft_control_type), POINTER :: dft_control
1445 TYPE(particle_list_type), POINTER :: particles
1446 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1447 TYPE(pw_c1d_gs_type) :: wf_g
1448 TYPE(pw_env_type), POINTER :: pw_env
1449 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1450 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1451 TYPE(pw_r3d_rs_type) :: wf_r
1452 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1453 TYPE(qs_subsys_type), POINTER :: subsys
1454
1455 logger => cp_get_default_logger()
1456 CALL timeset(routinen, handle)
1457
1458 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, pw_env=pw_env)
1459 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
1460 CALL auxbas_pw_pool%create_pw(wf_r)
1461 CALL auxbas_pw_pool%create_pw(wf_g)
1462
1463 CALL get_qs_env(qs_env, subsys=subsys)
1464 CALL qs_subsys_get(subsys, particles=particles)
1465
1466 my_pos_cube = "REWIND"
1467 IF (append_cube) THEN
1468 my_pos_cube = "APPEND"
1469 END IF
1470
1471 CALL get_qs_env(qs_env=qs_env, &
1472 atomic_kind_set=atomic_kind_set, &
1473 qs_kind_set=qs_kind_set, &
1474 cell=cell, &
1475 particle_set=particle_set)
1476
1477 DO iset = 1, 2
1478 CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
1479 DO i = 1, nmo
1480 CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1481 cell, dft_control, particle_set, pw_env)
1482 IF (iset == 1) THEN
1483 WRITE (filename, '(A6,I3.3,A5,I2.2,a11)') "_NEXC_", istate, "_NTO_", i, "_Hole_State"
1484 ELSEIF (iset == 2) THEN
1485 WRITE (filename, '(A6,I3.3,A5,I2.2,a15)') "_NEXC_", istate, "_NTO_", i, "_Particle_State"
1486 END IF
1487 info_approx_trunc = trim(adjustl(info_approximation))
1488 info_approx_trunc = info_approx_trunc(2:len_trim(info_approx_trunc) - 1)
1489 filename = trim(info_approx_trunc)//trim(filename)
1490 mpi_io = .true.
1491 unit_nr_cube = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
1492 middle_name=trim(filename), file_position=my_pos_cube, &
1493 log_filename=.false., ignore_should_output=.true., mpi_io=mpi_io)
1494 IF (iset == 1) THEN
1495 WRITE (title, *) "Natural Transition Orbital Hole State", i
1496 ELSEIF (iset == 2) THEN
1497 WRITE (title, *) "Natural Transition Orbital Particle State", i
1498 END IF
1499 CALL cp_pw_to_cube(wf_r, unit_nr_cube, title, particles=particles, stride=stride, mpi_io=mpi_io)
1500 CALL cp_print_key_finished_output(unit_nr_cube, logger, print_section, '', &
1501 ignore_should_output=.true., mpi_io=mpi_io)
1502 END DO
1503 END DO
1504
1505 CALL auxbas_pw_pool%give_back_pw(wf_g)
1506 CALL auxbas_pw_pool%give_back_pw(wf_r)
1507
1508 CALL timestop(handle)
1509 END SUBROUTINE print_bse_nto_cubes
1510
1511! **************************************************************************************************
1512!> \brief Checks BSE input section and adapts them if necessary
1513!> \param homo ...
1514!> \param virtual ...
1515!> \param unit_nr ...
1516!> \param mp2_env ...
1517!> \param qs_env ...
1518! **************************************************************************************************
1519 SUBROUTINE adapt_bse_input_params(homo, virtual, unit_nr, mp2_env, qs_env)
1520
1521 INTEGER, INTENT(IN) :: homo, virtual, unit_nr
1522 TYPE(mp2_type) :: mp2_env
1523 TYPE(qs_environment_type), POINTER :: qs_env
1524
1525 CHARACTER(LEN=*), PARAMETER :: routinen = 'adapt_BSE_input_params'
1526
1527 INTEGER :: handle, i, j, n, ndim_periodic_cell, &
1528 ndim_periodic_poisson, &
1529 num_state_list_exceptions
1530 TYPE(cell_type), POINTER :: cell_ref
1531 TYPE(pw_env_type), POINTER :: pw_env
1532 TYPE(pw_poisson_type), POINTER :: poisson_env
1533
1534 CALL timeset(routinen, handle)
1535 ! Get environment infos for later usage
1536 NULLIFY (pw_env, cell_ref, poisson_env)
1537 CALL get_qs_env(qs_env, pw_env=pw_env, cell_ref=cell_ref)
1538 CALL pw_env_get(pw_env, poisson_env=poisson_env)
1539 ndim_periodic_poisson = count(poisson_env%parameters%periodic == 1)
1540 ndim_periodic_cell = sum(cell_ref%perd(1:3)) ! Borrowed from cell_methods.F/write_cell_low
1541
1542 ! Handle negative NUM_PRINT_EXC
1543 IF (mp2_env%bse%num_print_exc < 0 .OR. &
1544 mp2_env%bse%num_print_exc > homo*virtual) THEN
1545 mp2_env%bse%num_print_exc = homo*virtual
1546 IF (unit_nr > 0) THEN
1547 CALL cp_hint(__location__, &
1548 "Keyword NUM_PRINT_EXC is either negative or too large. "// &
1549 "Printing all computed excitations.")
1550 END IF
1551 END IF
1552
1553 ! Default to NUM_PRINT_EXC if too large or negative,
1554 ! but only if NTOs are called - would be confusing for the user otherwise
1555 ! Prepare and adapt user inputs for NTO analysis
1556 ! Logic: Explicit state list overrides NUM_PRINT_EXC_NTOS
1557 ! If only NUM_PRINT_EXC_NTOS is given, we write the array 1,...,NUM_PRINT_EXC_NTOS to
1558 ! bse_nto_state_list
1559 IF (mp2_env%bse%do_nto_analysis) THEN
1560 IF (mp2_env%bse%explicit_nto_list) THEN
1561 IF (mp2_env%bse%num_print_exc_ntos > 0) THEN
1562 IF (unit_nr > 0) THEN
1563 CALL cp_hint(__location__, &
1564 "Keywords NUM_PRINT_EXC_NTOS and STATE_LIST are both given in input. "// &
1565 "Overriding NUM_PRINT_EXC_NTOS.")
1566 END IF
1567 END IF
1568 ! Check if all states are within the range
1569 ! Count them and initialize new array afterwards
1570 num_state_list_exceptions = 0
1571 DO i = 1, SIZE(mp2_env%bse%bse_nto_state_list)
1572 IF (mp2_env%bse%bse_nto_state_list(i) < 1 .OR. &
1573 mp2_env%bse%bse_nto_state_list(i) > mp2_env%bse%num_print_exc) THEN
1574 num_state_list_exceptions = num_state_list_exceptions + 1
1575 END IF
1576 END DO
1577 IF (num_state_list_exceptions > 0) THEN
1578 IF (unit_nr > 0) THEN
1579 CALL cp_hint(__location__, &
1580 "STATE_LIST contains indices outside the range of included excitation levels. "// &
1581 "Ignoring these states.")
1582 END IF
1583 END IF
1584 n = SIZE(mp2_env%bse%bse_nto_state_list) - num_state_list_exceptions
1585 ALLOCATE (mp2_env%bse%bse_nto_state_list_final(n))
1586 mp2_env%bse%bse_nto_state_list_final(:) = 0
1587 i = 1
1588 DO j = 1, SIZE(mp2_env%bse%bse_nto_state_list)
1589 IF (mp2_env%bse%bse_nto_state_list(j) >= 1 .AND. &
1590 mp2_env%bse%bse_nto_state_list(j) <= mp2_env%bse%num_print_exc) THEN
1591 mp2_env%bse%bse_nto_state_list_final(i) = mp2_env%bse%bse_nto_state_list(j)
1592 i = i + 1
1593 END IF
1594 END DO
1595
1596 mp2_env%bse%num_print_exc_ntos = SIZE(mp2_env%bse%bse_nto_state_list_final)
1597 ELSE
1598 IF (mp2_env%bse%num_print_exc_ntos > mp2_env%bse%num_print_exc .OR. &
1599 mp2_env%bse%num_print_exc_ntos < 0) THEN
1600 mp2_env%bse%num_print_exc_ntos = mp2_env%bse%num_print_exc
1601 END IF
1602 ALLOCATE (mp2_env%bse%bse_nto_state_list_final(mp2_env%bse%num_print_exc_ntos))
1603 DO i = 1, mp2_env%bse%num_print_exc_ntos
1604 mp2_env%bse%bse_nto_state_list_final(i) = i
1605 END DO
1606 END IF
1607 END IF
1608
1609 ! Takes care of triplet states, when oscillator strengths are 0
1610 IF (mp2_env%bse%bse_spin_config /= 0 .AND. &
1611 mp2_env%bse%eps_nto_osc_str > 0) THEN
1612 IF (unit_nr > 0) THEN
1613 CALL cp_warn(__location__, &
1614 "Cannot apply EPS_OSC_STR for Triplet excitations. "// &
1615 "Resetting EPS_OSC_STR to default.")
1616 END IF
1617 mp2_env%bse%eps_nto_osc_str = -1.0_dp
1618 END IF
1619
1620 ! Take care of number for computed exciton descriptors
1621 IF (mp2_env%bse%num_print_exc_descr < 0 .OR. &
1622 mp2_env%bse%num_print_exc_descr > mp2_env%bse%num_print_exc) THEN
1623 IF (unit_nr > 0) THEN
1624 CALL cp_hint(__location__, &
1625 "Keyword NUM_PRINT_EXC_DESCR is either negative or too large. "// &
1626 "Printing exciton descriptors up to NUM_PRINT_EXC.")
1627 END IF
1628 mp2_env%bse%num_print_exc_descr = mp2_env%bse%num_print_exc
1629 END IF
1630
1631 ! Handle screening factor options
1632 IF (mp2_env%BSE%screening_factor > 0.0_dp) THEN
1633 IF (mp2_env%BSE%screening_method /= bse_screening_alpha) THEN
1634 IF (unit_nr > 0) THEN
1635 CALL cp_warn(__location__, &
1636 "Screening factor is only supported for &SCREENING_IN_W ALPHA. "// &
1637 "Resetting SCREENING_IN_W to ALPHA.")
1638 END IF
1639 mp2_env%BSE%screening_method = bse_screening_alpha
1640 END IF
1641 IF (mp2_env%BSE%screening_factor > 1.0_dp) THEN
1642 IF (unit_nr > 0) THEN
1643 CALL cp_warn(__location__, &
1644 "Screening factor is larger than 1.0. ")
1645 END IF
1646 END IF
1647 END IF
1648
1649 IF (mp2_env%BSE%screening_factor < 0.0_dp .AND. &
1650 mp2_env%BSE%screening_method == bse_screening_alpha) THEN
1651 IF (unit_nr > 0) THEN
1652 CALL cp_warn(__location__, &
1653 "Screening factor is negative. Defaulting to 0.25")
1654 END IF
1655 mp2_env%BSE%screening_factor = 0.25_dp
1656 END IF
1657
1658 IF (mp2_env%BSE%screening_factor == 0.0_dp) THEN
1659 ! Use RPA internally in this case
1660 mp2_env%BSE%screening_method = bse_screening_rpa
1661 END IF
1662 IF (mp2_env%BSE%screening_factor == 1.0_dp) THEN
1663 ! Use TDHF internally in this case
1664 mp2_env%BSE%screening_method = bse_screening_tdhf
1665 END IF
1666
1667 ! Add warning for usage of KS energies
1668 IF (mp2_env%bse%use_ks_energies) THEN
1669 IF (unit_nr > 0) THEN
1670 CALL cp_warn(__location__, &
1671 "Using KS energies for BSE calculations. Therefore, no quantities "// &
1672 "of the preceeding GW calculation enter the BSE.")
1673 END IF
1674 END IF
1675
1676 ! Add warning if periodic calculation is invoked
1677 IF (ndim_periodic_poisson /= 0) THEN
1678 IF (unit_nr > 0) THEN
1679 CALL cp_warn(__location__, &
1680 "Poisson solver should be invoked by PERIODIC NONE. "// &
1681 "The applied length gauge might give misleading results for "// &
1682 "oscillator strengths.")
1683 END IF
1684 END IF
1685 IF (ndim_periodic_cell /= 0) THEN
1686 IF (unit_nr > 0) THEN
1687 CALL cp_warn(__location__, &
1688 "CELL in SUBSYS should be invoked with PERIODIC NONE. "// &
1689 "The applied length gauge might give misleading results for "// &
1690 "oscillator strengths.")
1691 END IF
1692 END IF
1693
1694 CALL timestop(handle)
1695 END SUBROUTINE adapt_bse_input_params
1696
1697! **************************************************************************************************
1698
1699! **************************************************************************************************
1700!> \brief ...
1701!> \param fm_multipole_ai_trunc ...
1702!> \param fm_multipole_ij_trunc ...
1703!> \param fm_multipole_ab_trunc ...
1704!> \param qs_env ...
1705!> \param mo_coeff ...
1706!> \param rpoint ...
1707!> \param n_moments ...
1708!> \param homo_red ...
1709!> \param virtual_red ...
1710!> \param context_BSE ...
1711! **************************************************************************************************
1712 SUBROUTINE get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, &
1713 qs_env, mo_coeff, rpoint, n_moments, &
1714 homo_red, virtual_red, context_BSE)
1715
1716 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
1717 INTENT(INOUT) :: fm_multipole_ai_trunc, &
1718 fm_multipole_ij_trunc, &
1719 fm_multipole_ab_trunc
1720 TYPE(qs_environment_type), POINTER :: qs_env
1721 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1722 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: rpoint
1723 INTEGER, INTENT(IN) :: n_moments, homo_red, virtual_red
1724 TYPE(cp_blacs_env_type), POINTER :: context_bse
1725
1726 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_multipoles_mo'
1727
1728 INTEGER :: handle, idir, n_multipole, n_occ, &
1729 n_virt, nao, nmo_mp2
1730 REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
1731 TYPE(cp_fm_struct_type), POINTER :: fm_struct_mp_ab_trunc, fm_struct_mp_ai_trunc, &
1732 fm_struct_mp_ij_trunc, fm_struct_multipoles_ao, fm_struct_nao_nmo, fm_struct_nmo_nmo
1733 TYPE(cp_fm_type) :: fm_work
1734 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_multipole_per_dir
1735 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_multipole, matrix_s
1736 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1737 TYPE(mp_para_env_type), POINTER :: para_env_bse
1738 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1739 POINTER :: sab_orb
1740
1741 CALL timeset(routinen, handle)
1742
1743 !First, we calculate the AO dipoles
1744 NULLIFY (sab_orb, matrix_s)
1745 CALL get_qs_env(qs_env, &
1746 mos=mos, &
1747 matrix_s=matrix_s, &
1748 sab_orb=sab_orb)
1749
1750 ! Use the same blacs environment as for the MO coefficients to ensure correct multiplication dbcsr x fm later on
1751 fm_struct_multipoles_ao => mos(1)%mo_coeff%matrix_struct
1752 ! BSE has different contexts and blacsenvs
1753 para_env_bse => context_bse%para_env
1754 ! Get size of multipole tensor
1755 n_multipole = (6 + 11*n_moments + 6*n_moments**2 + n_moments**3)/6 - 1
1756 NULLIFY (matrix_multipole)
1757 CALL dbcsr_allocate_matrix_set(matrix_multipole, n_multipole)
1758 ALLOCATE (fm_multipole_per_dir(n_multipole))
1759 DO idir = 1, n_multipole
1760 CALL dbcsr_init_p(matrix_multipole(idir)%matrix)
1761 CALL dbcsr_create(matrix_multipole(idir)%matrix, name="ao_multipole", &
1762 template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
1763 CALL cp_dbcsr_alloc_block_from_nbl(matrix_multipole(idir)%matrix, sab_orb)
1764 CALL dbcsr_set(matrix_multipole(idir)%matrix, 0._dp)
1765 END DO
1766
1767 CALL get_reference_point(rpoint, qs_env=qs_env, reference=use_mom_ref_coac, ref_point=ref_point)
1768
1769 CALL build_local_moment_matrix(qs_env, matrix_multipole, n_moments, ref_point=rpoint)
1770
1771 NULLIFY (sab_orb)
1772
1773 ! Now we transform them to MO
1774 ! n_occ is the number of occupied MOs, nao the number of all AOs
1775 ! Writing homo to n_occ instead if nmo,
1776 ! takes care of ADDED_MOS, which would overwrite nmo of qs_env-mos, if invoked
1777 CALL get_mo_set(mo_set=mos(1), homo=n_occ, nao=nao)
1778 ! Takes into account removed nullspace values from SVD
1779 nmo_mp2 = mo_coeff(1)%matrix_struct%ncol_global
1780 n_virt = nmo_mp2 - n_occ
1781
1782 ! At the end, we need four different layouts of matrices in this multiplication, e.g. for a dipole:
1783 ! D_pq = full multipole matrix for occupied and unoccupied
1784 ! Final result:D_pq= C_{mu p} <\mu|\vec{r}|\nu> C_{\nu q} EQ.I
1785 ! \_______/ \___________/ \______/
1786 ! fm_coeff matrix_multipole fm_coeff
1787 ! (EQ.Ia) (EQ.Ib) (EQ.Ia)
1788 ! Intermediate work matrices:
1789 ! fm_work = <\mu|\vec{r}|\nu> C_{\nu q} EQ.II
1790
1791 ! Struct for the full multipole matrix
1792 CALL cp_fm_struct_create(fm_struct_nao_nmo, &
1793 fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, &
1794 nao, nmo_mp2)
1795 CALL cp_fm_struct_create(fm_struct_nmo_nmo, &
1796 fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, &
1797 nmo_mp2, nmo_mp2)
1798
1799 ! At the very end, we copy the multipoles corresponding to truncated BSE indices in i and a
1800 CALL cp_fm_struct_create(fm_struct_mp_ai_trunc, para_env_bse, &
1801 context_bse, virtual_red, homo_red)
1802 CALL cp_fm_struct_create(fm_struct_mp_ij_trunc, para_env_bse, &
1803 context_bse, homo_red, homo_red)
1804 CALL cp_fm_struct_create(fm_struct_mp_ab_trunc, para_env_bse, &
1805 context_bse, virtual_red, virtual_red)
1806 DO idir = 1, n_multipole
1807 CALL cp_fm_create(fm_multipole_ai_trunc(idir), matrix_struct=fm_struct_mp_ai_trunc, &
1808 name="dipoles_mo_ai_trunc")
1809 CALL cp_fm_set_all(fm_multipole_ai_trunc(idir), 0.0_dp)
1810 CALL cp_fm_create(fm_multipole_ij_trunc(idir), matrix_struct=fm_struct_mp_ij_trunc, &
1811 name="dipoles_mo_ij_trunc")
1812 CALL cp_fm_set_all(fm_multipole_ij_trunc(idir), 0.0_dp)
1813 CALL cp_fm_create(fm_multipole_ab_trunc(idir), matrix_struct=fm_struct_mp_ab_trunc, &
1814 name="dipoles_mo_ab_trunc")
1815 CALL cp_fm_set_all(fm_multipole_ab_trunc(idir), 0.0_dp)
1816 END DO
1817
1818 ! Need another temporary matrix to store intermediate result from right multiplication
1819 ! D = C_{mu a} <\mu|\vec{r}|\nu> C_{\nu i}
1820 CALL cp_fm_create(fm_work, matrix_struct=fm_struct_nao_nmo, name="multipole_work")
1821 CALL cp_fm_set_all(fm_work, 0.0_dp)
1822
1823 DO idir = 1, n_multipole
1824 ! Create the full multipole matrix per direction
1825 CALL cp_fm_create(fm_multipole_per_dir(idir), matrix_struct=fm_struct_nmo_nmo, name="multipoles_mo")
1826 CALL cp_fm_set_all(fm_multipole_per_dir(idir), 0.0_dp)
1827 ! Fill final (MO) multipole matrix
1828 CALL cp_dbcsr_sm_fm_multiply(matrix_multipole(idir)%matrix, mo_coeff(1), &
1829 fm_work, ncol=nmo_mp2)
1830 ! Now obtain the multipoles by the final multiplication;
1831 ! We do that inside the loop to obtain multipoles per axis for print
1832 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))
1833
1834 ! Truncate full matrix to the BSE indices
1835 ! D_ai
1836 CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1837 fm_multipole_ai_trunc(idir), &
1838 virtual_red, &
1839 homo_red, &
1840 n_occ + 1, &
1841 n_occ - homo_red + 1, &
1842 1, &
1843 1, &
1844 fm_multipole_per_dir(idir)%matrix_struct%context)
1845 ! D_ij
1846 CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1847 fm_multipole_ij_trunc(idir), &
1848 homo_red, &
1849 homo_red, &
1850 n_occ - homo_red + 1, &
1851 n_occ - homo_red + 1, &
1852 1, &
1853 1, &
1854 fm_multipole_per_dir(idir)%matrix_struct%context)
1855 ! D_ab
1856 CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1857 fm_multipole_ab_trunc(idir), &
1858 virtual_red, &
1859 virtual_red, &
1860 n_occ + 1, &
1861 n_occ + 1, &
1862 1, &
1863 1, &
1864 fm_multipole_per_dir(idir)%matrix_struct%context)
1865 END DO
1866
1867 !Release matrices and structs
1868 NULLIFY (fm_struct_multipoles_ao)
1869 CALL cp_fm_struct_release(fm_struct_mp_ai_trunc)
1870 CALL cp_fm_struct_release(fm_struct_mp_ij_trunc)
1871 CALL cp_fm_struct_release(fm_struct_mp_ab_trunc)
1872 CALL cp_fm_struct_release(fm_struct_nao_nmo)
1873 CALL cp_fm_struct_release(fm_struct_nmo_nmo)
1874 DO idir = 1, n_multipole
1875 CALL cp_fm_release(fm_multipole_per_dir(idir))
1876 END DO
1877 DEALLOCATE (fm_multipole_per_dir)
1878 CALL cp_fm_release(fm_work)
1879 CALL dbcsr_deallocate_matrix_set(matrix_multipole)
1880
1881 CALL timestop(handle)
1882
1883 END SUBROUTINE get_multipoles_mo
1884
1885! **************************************************************************************************
1886!> \brief Computes trace of form Tr{A^T B C} for exciton descriptors
1887!> \param fm_A Full Matrix, typically X or Y, in format homo x virtual
1888!> \param fm_B ...
1889!> \param fm_C ...
1890!> \param alpha ...
1891! **************************************************************************************************
1892 SUBROUTINE trace_exciton_descr(fm_A, fm_B, fm_C, alpha)
1893
1894 TYPE(cp_fm_type), INTENT(IN) :: fm_a, fm_b, fm_c
1895 REAL(kind=dp), INTENT(OUT) :: alpha
1896
1897 CHARACTER(LEN=*), PARAMETER :: routinen = 'trace_exciton_descr'
1898
1899 INTEGER :: handle, ncol_a, ncol_b, ncol_c, nrow_a, &
1900 nrow_b, nrow_c
1901 TYPE(cp_fm_type) :: fm_work_ia
1902
1903 CALL timeset(routinen, handle)
1904
1905 CALL cp_fm_create(fm_work_ia, fm_a%matrix_struct)
1906 CALL cp_fm_get_info(fm_a, nrow_global=nrow_a, ncol_global=ncol_a)
1907 CALL cp_fm_get_info(fm_b, nrow_global=nrow_b, ncol_global=ncol_b)
1908 CALL cp_fm_get_info(fm_c, nrow_global=nrow_c, ncol_global=ncol_c)
1909
1910 ! Check matrix sizes
1911 cpassert(nrow_a == nrow_b .AND. ncol_a == ncol_c .AND. ncol_b == nrow_c)
1912
1913 CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
1914
1915 CALL parallel_gemm("N", "N", nrow_a, ncol_a, nrow_c, 1.0_dp, &
1916 fm_b, fm_c, 0.0_dp, fm_work_ia)
1917
1918 CALL cp_fm_trace(fm_a, fm_work_ia, alpha)
1919
1920 CALL cp_fm_release(fm_work_ia)
1921
1922 CALL timestop(handle)
1923
1924 END SUBROUTINE trace_exciton_descr
1925
1926END MODULE
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:75
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
util::sort_unique
Definition util.F:36
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
bse_util
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_util.F:13
bse_util::filter_eigvec_contrib
subroutine, public filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, i_exc, virtual, num_entries, mp2_env)
Filters eigenvector entries above a given threshold to describe excitations in the singleparticle bas...
Definition bse_util.F:994
bse_util::get_multipoles_mo
subroutine, public get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, qs_env, mo_coeff, rpoint, n_moments, homo_red, virtual_red, context_bse)
...
Definition bse_util.F:1715
bse_util::trace_exciton_descr
subroutine, public trace_exciton_descr(fm_a, fm_b, fm_c, alpha)
Computes trace of form Tr{A^T B C} for exciton descriptors.
Definition bse_util.F:1893
bse_util::determine_cutoff_indices
subroutine, public determine_cutoff_indices(eigenval, homo, virtual, homo_red, virt_red, homo_incl, virt_incl, mp2_env)
Reads cutoffs for BSE from mp2_env and compares to energies in Eigenval to extract reduced homo/virtu...
Definition bse_util.F:1130
bse_util::estimate_bse_resources
subroutine, public estimate_bse_resources(homo_red, virtual_red, unit_nr, bse_abba, para_env, diag_runtime_est)
Roughly estimates the needed runtime and memory during the BSE run.
Definition bse_util.F:936
bse_util::truncate_fm
subroutine, public truncate_fm(fm_out, fm_in, ncol_in, nrow_out, ncol_out, unit_nr, mp2_env, nrow_offset, ncol_offset)
Routine for truncating a full matrix as given by the energy cutoffs in the input file....
Definition bse_util.F:459
bse_util::deallocate_matrices_bse
subroutine, public deallocate_matrices_bse(fm_mat_s_bar_ia_bse, fm_mat_s_bar_ij_bse, fm_mat_s_trunc, fm_mat_s_ij_trunc, fm_mat_s_ab_trunc, fm_mat_q_static_bse_gemm, mp2_env)
...
Definition bse_util.F:752
bse_util::reshuffle_eigvec
subroutine, public reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, unit_nr, mp2_env)
...
Definition bse_util.F:1351
bse_util::mult_b_with_w
subroutine, public mult_b_with_w(fm_mat_s_ij_bse, fm_mat_s_ia_bse, fm_mat_s_bar_ia_bse, fm_mat_s_bar_ij_bse, fm_mat_q_static_bse_gemm, dimen_ri, homo, virtual)
Multiplies B-matrix (RI-3c-Integrals) with W (screening) to obtain \bar{B}.
Definition bse_util.F:113
bse_util::adapt_bse_input_params
subroutine, public adapt_bse_input_params(homo, virtual, unit_nr, mp2_env, qs_env)
Checks BSE input section and adapts them if necessary.
Definition bse_util.F:1520
bse_util::truncate_bse_matrices
subroutine, public truncate_bse_matrices(fm_mat_s_ia_bse, fm_mat_s_ij_bse, fm_mat_s_ab_bse, fm_mat_s_trunc, fm_mat_s_ij_trunc, fm_mat_s_ab_trunc, eigenval_scf, eigenval, eigenval_reduced, homo, virtual, dimen_ri, unit_nr, bse_lev_virt, homo_red, virt_red, mp2_env)
Determines indices within the given energy cutoffs and truncates Eigenvalues and matrices.
Definition bse_util.F:1211
bse_util::comp_eigvec_coeff_bse
subroutine, public comp_eigvec_coeff_bse(fm_work, eig_vals, beta, gamma, do_transpose)
Routine for computing the coefficients of the eigenvectors of the BSE matrix from a multiplication wi...
Definition bse_util.F:787
bse_util::print_bse_nto_cubes
subroutine, public print_bse_nto_cubes(qs_env, mos, istate, info_approximation, stride, append_cube, print_section)
Borrowed from the tddfpt module with slight adaptions.
Definition bse_util.F:1425
bse_util::fm_general_add_bse
subroutine, public fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env)
Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,...
Definition bse_util.F:199
bse_util::sort_excitations
subroutine, public sort_excitations(idx_prim, idx_sec, eigvec_entries)
...
Definition bse_util.F:852
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:205
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:552
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1747
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_invert
subroutine, public cp_fm_cholesky_invert(matrix, n, info_out)
used to replace the cholesky decomposition by the inverse
Definition cp_fm_cholesky.F:122
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:50
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1009
cp_fm_types::cp_fm_to_fm_submat_general
subroutine, public cp_fm_to_fm_submat_general(source, destination, nrows, ncols, s_firstrow, s_firstcol, d_firstrow, d_firstcol, global_context)
General copy of a submatrix of fm matrix to a submatrix of another fm matrix. The two matrices can ha...
Definition cp_fm_types.F:2002
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1472
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:528
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
Definition cp_realspace_grid_cube.F:45
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::use_mom_ref_coac
integer, parameter, public use_mom_ref_coac
Definition input_constants.F:497
input_constants::bse_screening_tdhf
integer, parameter, public bse_screening_tdhf
Definition input_constants.F:1324
input_constants::bse_screening_alpha
integer, parameter, public bse_screening_alpha
Definition input_constants.F:1324
input_constants::bse_screening_rpa
integer, parameter, public bse_screening_rpa
Definition input_constants.F:1324
input_constants::bse_abba
integer, parameter, public bse_abba
Definition input_constants.F:1319
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
moments_utils
Calculates the moment integrals <a|r^m|b>
Definition moments_utils.F:15
moments_utils::get_reference_point
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition moments_utils.F:61
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_wavefunction
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
Definition qs_collocate_density.F:2854
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition qs_moments.F:560
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:134
rpa_communication
Auxiliary routines necessary to redistribute an fm_matrix from a given blacs_env to another.
Definition rpa_communication.F:14
rpa_communication::communicate_buffer
subroutine, public communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array, do_indx, do_msg)
...
Definition rpa_communication.F:434
util
All kind of helpful little routines.
Definition util.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
message_passing::mp_request_type
Definition message_passing.F:623
mp2_types::integ_mat_buffer_type
Definition mp2_types.F:409
mp2_types::mp2_type
Definition mp2_types.F:345
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:122
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56