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