(git:3f0d10b)
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 bse_types, ONLY: bse_env_type
16 USE cell_types, ONLY: cell_type
17 USE cp_blacs_env, ONLY: cp_blacs_env_type
18 USE cp_control_types, ONLY: dft_control_type
19 USE cp_dbcsr_api, ONLY: dbcsr_p_type
20 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
21 USE cp_fm_basic_linalg, ONLY: cp_fm_trace
22 USE cp_fm_redistribution, ONLY: cp_fm_permute
23 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
24 cp_fm_struct_release,&
25 cp_fm_struct_type
26 USE cp_fm_types, ONLY: cp_fm_create,&
27 cp_fm_get_info,&
28 cp_fm_release,&
29 cp_fm_set_all,&
30 cp_fm_to_fm_submat,&
31 cp_fm_to_fm_submat_general,&
32 cp_fm_type
33 USE cp_log_handling, ONLY: cp_get_default_logger,&
34 cp_logger_type
35 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
36 cp_print_key_unit_nr
37 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
38 USE input_constants, ONLY: bse_iterdiag,&
39 bse_screening_alpha,&
40 bse_screening_rpa,&
41 bse_screening_tdhf
42 USE input_section_types, ONLY: section_vals_type,&
43 section_vals_val_get
44 USE kinds, ONLY: default_path_length,&
45 dp,&
46 int_8
47 USE message_passing, ONLY: mp_para_env_type
48 USE mo_window, ONLY: combine_mo_windows,&
49 determine_mo_window,&
50 mo_window_type
51 USE parallel_gemm_api, ONLY: parallel_gemm
52 USE particle_list_types, ONLY: particle_list_type
53 USE particle_types, ONLY: particle_type
54 USE physcon, ONLY: evolt
55 USE pw_env_types, ONLY: pw_env_get,&
56 pw_env_type
57 USE pw_pool_types, ONLY: pw_pool_p_type,&
58 pw_pool_type
59 USE pw_types, ONLY: pw_c1d_gs_type,&
60 pw_r3d_rs_type
61 USE qs_collocate_density, ONLY: calculate_wavefunction
62 USE qs_kind_types, ONLY: qs_kind_type
63 USE qs_mo_types, ONLY: get_mo_set,&
64 mo_set_type
65 USE qs_subsys_types, ONLY: qs_subsys_get,&
66 qs_subsys_type
67 USE util, ONLY: sort,&
68 sort_unique
69#include "./base/base_uses.f90"
70
71 IMPLICIT NONE
72
73 PRIVATE
74
75 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_util'
76
77 PUBLIC :: screen_slabs, bse_window_from_dft, bse_window_from_mos, &
78 comp_eigvec_coeff_bse, sort_excitations, &
79 estimate_bse_resources, fulldiag_memory_estimate_gb, filter_eigvec_contrib, truncate_bse_matrices, &
80 adapt_bse_input_params, adapt_bse_screening, get_multipoles_mo, &
81 reshuffle_eigvec, print_bse_nto_cubes, trace_exciton_descr, &
82 get_bse_spin_block_layout, assemble_joint_ov_slab, &
83 ia_of_occ_virt, occ_of_ia, virt_of_ia
84
85CONTAINS
86
87! **************************************************************************************************
88!> \brief Compound transition index ia = (i-1)*virt + a of the pair (i, a): a is the fast index
89!> \param i_occ occupied level
90!> \param a_virt virtual level
91!> \param virt number of virtual levels
92!> \return ia
93! **************************************************************************************************
94 PURE FUNCTION ia_of_occ_virt(i_occ, a_virt, virt) RESULT(ia)
95
96 INTEGER, INTENT(IN) :: i_occ, a_virt, virt
97 INTEGER :: ia
98
99 ia = (i_occ - 1)*virt + a_virt
100
101 END FUNCTION ia_of_occ_virt
102
103! **************************************************************************************************
104!> \brief Occupied level of the compound transition index ia = (i-1)*virt + a
105!> \param ia compound transition index
106!> \param virt number of virtual levels
107!> \return i
108! **************************************************************************************************
109 PURE FUNCTION occ_of_ia(ia, virt) RESULT(i_occ)
110
111 INTEGER, INTENT(IN) :: ia, virt
112 INTEGER :: i_occ
113
114 i_occ = (ia - 1)/virt + 1
115
116 END FUNCTION occ_of_ia
117
118! **************************************************************************************************
119!> \brief Virtual level of the compound transition index ia = (i-1)*virt + a
120!> \param ia compound transition index
121!> \param virt number of virtual levels
122!> \return a
123! **************************************************************************************************
124 PURE FUNCTION virt_of_ia(ia, virt) RESULT(a_virt)
125
126 INTEGER, INTENT(IN) :: ia, virt
127 INTEGER :: a_virt
128
129 a_virt = mod(ia - 1, virt) + 1
130
131 END FUNCTION virt_of_ia
132
133! **************************************************************************************************
134!> \brief The screened slabs \bar{B}^P_ij = sum_Q [1+Q(0)]^-1_PQ B^Q_ij and
135!> \bar{B}^P_ia = sum_Q [1+Q(0)]^-1_PQ B^Q_ia, P, Q ≤ dimen_RI_red
136!> \param fm_eps_inv [1+Q(0)]^-1_PQ
137!> \param fm_S_ij B^P_ij of the window, dimen_RI x homo^2
138!> \param fm_S_ia B^P_ia of the window, dimen_RI x homo*virtual
139!> \param dimen_RI_red reduced RI dimension, the contraction length
140!> \param homo occupied levels in the window
141!> \param virtual virtual levels in the window
142!> \param fm_S_bar_ij \bar{B}^P_ij, created here on fm_S_ij's structure
143!> \param fm_S_bar_ia \bar{B}^P_ia, created here on fm_S_ia's structure for ABBA only
144!> \param do_abba whether the coupling matrix B needs the screened ia slab
145! **************************************************************************************************
146 SUBROUTINE screen_slabs(fm_eps_inv, fm_S_ij, fm_S_ia, dimen_RI_red, homo, virtual, &
147 fm_S_bar_ij, fm_S_bar_ia, do_abba)
148
149 TYPE(cp_fm_type), INTENT(IN) :: fm_eps_inv, fm_s_ij, fm_s_ia
150 INTEGER, INTENT(IN) :: dimen_ri_red, homo, virtual
151 TYPE(cp_fm_type), INTENT(OUT) :: fm_s_bar_ij, fm_s_bar_ia
152 LOGICAL, INTENT(IN) :: do_abba
153
154 CHARACTER(LEN=*), PARAMETER :: routinen = 'screen_slabs'
155
156 INTEGER :: handle
157
158 CALL timeset(routinen, handle)
159
160 CALL cp_fm_create(fm_s_bar_ij, fm_s_ij%matrix_struct)
161 CALL cp_fm_set_all(fm_s_bar_ij, 0.0_dp)
162
163 ! \bar{B}^P_ij = sum_Q [1+Q(0)]^-1_PQ B^Q_ij
164 CALL parallel_gemm(transa="N", transb="N", m=dimen_ri_red, n=homo**2, k=dimen_ri_red, alpha=1.0_dp, &
165 matrix_a=fm_eps_inv, matrix_b=fm_s_ij, beta=0.0_dp, matrix_c=fm_s_bar_ij)
166
167 IF (do_abba) THEN
168 CALL cp_fm_create(fm_s_bar_ia, fm_s_ia%matrix_struct)
169 CALL cp_fm_set_all(fm_s_bar_ia, 0.0_dp)
170 ! \bar{B}^P_ia = sum_Q [1+Q(0)]^-1_PQ B^Q_ia; fm_S_ia comes from RPA and has a different
171 ! blacs_env than fm_S_ij
172 CALL parallel_gemm(transa="N", transb="N", m=dimen_ri_red, n=homo*virtual, k=dimen_ri_red, alpha=1.0_dp, &
173 matrix_a=fm_eps_inv, matrix_b=fm_s_ia, beta=0.0_dp, matrix_c=fm_s_bar_ia)
174 END IF
175
176 CALL timestop(handle)
177
178 END SUBROUTINE screen_slabs
179
180! **************************************************************************************************
181!> \brief Routine for computing the coefficients of the eigenvectors of the BSE matrix from a
182!> multiplication with the eigenvalues
183!> \param fm_work ...
184!> \param eig_vals ...
185!> \param beta ...
186!> \param gamma ...
187!> \param do_transpose ...
188! **************************************************************************************************
189 SUBROUTINE comp_eigvec_coeff_bse(fm_work, eig_vals, beta, gamma, do_transpose)
190
191 TYPE(cp_fm_type), INTENT(INOUT) :: fm_work
192 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
193 INTENT(IN) :: eig_vals
194 REAL(kind=dp), INTENT(IN) :: beta
195 REAL(kind=dp), INTENT(IN), OPTIONAL :: gamma
196 LOGICAL, INTENT(IN), OPTIONAL :: do_transpose
197
198 CHARACTER(LEN=*), PARAMETER :: routinen = 'comp_eigvec_coeff_BSE'
199
200 INTEGER :: handle, i_row_global, ii, j_col_global, &
201 jj, ncol_local, nrow_local
202 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
203 LOGICAL :: my_do_transpose
204 REAL(kind=dp) :: coeff, my_gamma
205
206 CALL timeset(routinen, handle)
207
208 IF (PRESENT(gamma)) THEN
209 my_gamma = gamma
210 ELSE
211 my_gamma = 2.0_dp
212 END IF
213
214 IF (PRESENT(do_transpose)) THEN
215 my_do_transpose = do_transpose
216 ELSE
217 my_do_transpose = .false.
218 END IF
219
220 CALL cp_fm_get_info(matrix=fm_work, &
221 nrow_local=nrow_local, &
222 ncol_local=ncol_local, &
223 row_indices=row_indices, &
224 col_indices=col_indices)
225
226 IF (my_do_transpose) THEN
227 DO jj = 1, ncol_local
228 j_col_global = col_indices(jj)
229 DO ii = 1, nrow_local
230 coeff = (eig_vals(j_col_global)**beta)/my_gamma
231 fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff
232 END DO
233 END DO
234 ELSE
235 DO jj = 1, ncol_local
236 DO ii = 1, nrow_local
237 i_row_global = row_indices(ii)
238 coeff = (eig_vals(i_row_global)**beta)/my_gamma
239 fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff
240 END DO
241 END DO
242 END IF
243
244 CALL timestop(handle)
245
246 END SUBROUTINE comp_eigvec_coeff_bse
247
248! **************************************************************************************************
249!> \brief Sorts excitation entries by ascending primary index, reordering the secondary index,
250!> the eigenvector coefficients and - open shell - the spin index alongside
251!> \param idx_prim Primary index of each entry; sorted in place and used as the sort key
252!> \param idx_sec Secondary index of each entry, reordered to follow idx_prim
253!> \param eigvec_entries Eigenvector coefficients of each entry, reordered to follow idx_prim
254!> \param idx_spin Optional spin index of each entry (open shell), reordered to follow idx_prim
255! **************************************************************************************************
256 SUBROUTINE sort_excitations(idx_prim, idx_sec, eigvec_entries, idx_spin)
257
258 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim, idx_sec
259 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
260 INTEGER, ALLOCATABLE, DIMENSION(:), OPTIONAL :: idx_spin
261
262 CHARACTER(LEN=*), PARAMETER :: routinen = 'sort_excitations'
263
264 INTEGER :: handle, ii, kk, num_entries, num_mults
265 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim_work, idx_sec_work, &
266 idx_spin_work, tmp_index
267 LOGICAL :: unique_entries
268 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries_work
269
270 CALL timeset(routinen, handle)
271
272 num_entries = SIZE(idx_prim)
273
274 ALLOCATE (tmp_index(num_entries))
275
276 CALL sort(idx_prim, num_entries, tmp_index)
277
278 ALLOCATE (idx_sec_work(num_entries))
279 ALLOCATE (eigvec_entries_work(num_entries))
280 IF (PRESENT(idx_spin)) ALLOCATE (idx_spin_work(num_entries))
281
282 DO ii = 1, num_entries
283 idx_sec_work(ii) = idx_sec(tmp_index(ii))
284 eigvec_entries_work(ii) = eigvec_entries(tmp_index(ii))
285 IF (PRESENT(idx_spin)) idx_spin_work(ii) = idx_spin(tmp_index(ii))
286 END DO
287
288 DEALLOCATE (tmp_index)
289 DEALLOCATE (idx_sec)
290 DEALLOCATE (eigvec_entries)
291
292 CALL move_alloc(idx_sec_work, idx_sec)
293 CALL move_alloc(eigvec_entries_work, eigvec_entries)
294 IF (PRESENT(idx_spin)) THEN
295 DEALLOCATE (idx_spin)
296 CALL move_alloc(idx_spin_work, idx_spin)
297 END IF
298
299 !Now check for multiple entries in first idx to check necessity of sorting in second idx
300 CALL sort_unique(idx_prim, unique_entries)
301 IF (.NOT. unique_entries) THEN
302 ALLOCATE (idx_prim_work(num_entries))
303 idx_prim_work(:) = idx_prim(:)
304 ! Find duplicate entries in idx_prim
305 DO ii = 1, num_entries
306 IF (idx_prim_work(ii) == 0) cycle
307 num_mults = count(idx_prim_work == idx_prim_work(ii))
308 IF (num_mults > 1) THEN
309 !Set all duplicate entries to 0
310 idx_prim_work(ii:ii + num_mults - 1) = 0
311 !Start sorting in secondary index
312 ALLOCATE (idx_sec_work(num_mults))
313 ALLOCATE (eigvec_entries_work(num_mults))
314 idx_sec_work(:) = idx_sec(ii:ii + num_mults - 1)
315 eigvec_entries_work(:) = eigvec_entries(ii:ii + num_mults - 1)
316 IF (PRESENT(idx_spin)) THEN
317 ALLOCATE (idx_spin_work(num_mults))
318 idx_spin_work(:) = idx_spin(ii:ii + num_mults - 1)
319 END IF
320 ALLOCATE (tmp_index(num_mults))
321 CALL sort(idx_sec_work, num_mults, tmp_index)
322
323 !Now write newly sorted indices to original arrays
324 DO kk = ii, ii + num_mults - 1
325 idx_sec(kk) = idx_sec_work(kk - ii + 1)
326 eigvec_entries(kk) = eigvec_entries_work(tmp_index(kk - ii + 1))
327 IF (PRESENT(idx_spin)) idx_spin(kk) = idx_spin_work(tmp_index(kk - ii + 1))
328 END DO
329 !Deallocate work arrays
330 DEALLOCATE (tmp_index)
331 DEALLOCATE (idx_sec_work)
332 DEALLOCATE (eigvec_entries_work)
333 IF (PRESENT(idx_spin)) DEALLOCATE (idx_spin_work)
334 END IF
335 idx_prim_work(ii) = idx_prim(ii)
336 END DO
337 DEALLOCATE (idx_prim_work)
338 END IF
339
340 CALL timestop(handle)
341
342 END SUBROUTINE sort_excitations
343
344! **************************************************************************************************
345!> \brief Roughly estimates the needed runtime and memory during the BSE run
346!> \param n_ov_joint ...
347!> \param unit_nr ...
348!> \param bse_abba ...
349!> \param para_env ...
350!> \param diag_runtime_est ...
351! **************************************************************************************************
352 SUBROUTINE estimate_bse_resources(n_ov_joint, unit_nr, bse_abba, &
353 para_env, diag_runtime_est)
354
355 INTEGER, INTENT(IN) :: n_ov_joint, unit_nr
356 LOGICAL :: bse_abba
357 TYPE(mp_para_env_type), POINTER :: para_env
358 REAL(kind=dp) :: diag_runtime_est
359
360 CHARACTER(LEN=*), PARAMETER :: routinen = 'estimate_BSE_resources'
361
362 INTEGER :: handle
363 REAL(kind=dp) :: mem_est, mem_est_per_rank
364
365 CALL timeset(routinen, handle)
366
367 mem_est = fulldiag_memory_estimate_gb(n_ov_joint, bse_abba)
368 mem_est_per_rank = real(mem_est/para_env%num_pe, kind=dp)
369
370 IF (unit_nr > 0) THEN
371 ! WRITE (unit_nr, '(T2,A4,T7,A40,T68,F13.3)') 'BSE|', 'Total peak memory estimate from BSE [GB]', &
372 ! mem_est
373 WRITE (unit_nr, '(T2,A4,T7,A40,T68,ES13.3)') 'BSE|', 'Total peak memory estimate from BSE [GB]', &
374 mem_est
375 WRITE (unit_nr, '(T2,A4,T7,A47,T68,F13.3)') 'BSE|', 'Peak memory estimate per MPI rank from BSE [GB]', &
376 mem_est_per_rank
377 WRITE (unit_nr, '(T2,A4)') 'BSE|'
378 END IF
379 ! Rough estimation of diagonalization runtimes. Baseline was a full BSE Naphthalene
380 ! run with 11000x11000 entries in A/B/C, which took 10s on 32 ranks
381 diag_runtime_est = real(int(n_ov_joint, kind=int_8)/11000_int_8, kind=dp)**3* &
382 10*32/real(para_env%num_pe, kind=dp)
383
384 CALL timestop(handle)
385
386 END SUBROUTINE estimate_bse_resources
387
388! **************************************************************************************************
389!> \brief Peak memory of the full diagonalization in GB, all ranks together: n_ov^2 doubles times
390!> the matrices held at once, 2 for TDA (A and W_ijab), 10 for ABBA (B, C and the six that
391!> build C added)
392!> \param n_ov_joint transition pairs in the active window, summed over spins
393!> \param do_abba full problem, not TDA
394!> \return the estimate in GB
395! **************************************************************************************************
396 PURE FUNCTION fulldiag_memory_estimate_gb(n_ov_joint, do_abba) RESULT(mem_GB)
397
398 INTEGER, INTENT(IN) :: n_ov_joint
399 LOGICAL, INTENT(IN) :: do_abba
400 REAL(kind=dp) :: mem_gb
401
402 INTEGER :: num_bse_matrices
403 INTEGER(KIND=int_8) :: full_dim
404
405 num_bse_matrices = 2
406 IF (do_abba) num_bse_matrices = 10
407 full_dim = int(n_ov_joint, kind=int_8)**2*int(num_bse_matrices, kind=int_8)
408 ! 1 GB = 1E9 B, the unit of the free-memory budget it is compared with
409 mem_gb = 8.0e-9_dp*real(full_dim, kind=dp)
410
411 END FUNCTION fulldiag_memory_estimate_gb
412
413! **************************************************************************************************
414!> \brief Filters eigenvector entries above a given threshold to describe excitations in the
415!> singleparticle basis
416!> \param fm_eigvec ...
417!> \param idx_homo ...
418!> \param idx_virt ...
419!> \param eigvec_entries ...
420!> \param i_exc ...
421!> \param virtual ...
422!> \param num_entries ...
423!> \param bse_env the BSE environment (settings)
424!> \param offsets ...
425!> \param virtual_per_spin ...
426!> \param idx_spin ...
427! **************************************************************************************************
428 SUBROUTINE filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, &
429 i_exc, virtual, num_entries, bse_env, &
430 offsets, virtual_per_spin, idx_spin)
431
432 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec
433 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt
434 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
435 INTEGER :: i_exc, virtual, num_entries
436 TYPE(bse_env_type), INTENT(IN) :: bse_env
437 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: offsets, virtual_per_spin
438 INTEGER, ALLOCATABLE, DIMENSION(:), OPTIONAL :: idx_spin
439
440 CHARACTER(LEN=*), PARAMETER :: routinen = 'filter_eigvec_contrib'
441
442 INTEGER :: eigvec_idx, handle, ii, iproc, isp, jj, &
443 kk, ksp, ncol_local, nrow_local, &
444 num_entries_local, r_local, v_local
445 INTEGER, ALLOCATABLE, DIMENSION(:) :: displ, indx_all, indx_local, &
446 num_entries_to_comm
447 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
448 REAL(kind=dp) :: eigvec_entry
449 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: entries_local
450 TYPE(mp_para_env_type), POINTER :: para_env
451
452 CALL timeset(routinen, handle)
453
454 para_env => fm_eigvec%matrix_struct%para_env
455
456 CALL cp_fm_get_info(matrix=fm_eigvec, &
457 nrow_local=nrow_local, &
458 ncol_local=ncol_local, &
459 row_indices=row_indices, &
460 col_indices=col_indices)
461
462 ! the local entries of column i_exc above EPS_X: count, then collect with their (i, a, sigma)
463 num_entries_local = 0
464 DO jj = 1, ncol_local
465 IF (col_indices(jj) /= i_exc) cycle
466 DO ii = 1, nrow_local
467 IF (abs(fm_eigvec%local_data(ii, jj)) > bse_env%eps_x) num_entries_local = num_entries_local + 1
468 END DO
469 END DO
470 ALLOCATE (entries_local(num_entries_local), indx_local(3*num_entries_local))
471 kk = 0
472 DO jj = 1, ncol_local
473 IF (col_indices(jj) /= i_exc) cycle
474 DO ii = 1, nrow_local
475 eigvec_idx = row_indices(ii)
476 eigvec_entry = fm_eigvec%local_data(ii, jj)
477 IF (abs(eigvec_entry) > bse_env%eps_x) THEN
478 ! Decode spin block from the joint row index (blocks are contiguous; sigma is the
479 ! largest offset strictly below eigvec_idx). offsets absent -> closed shell, sigma=1.
480 isp = 1
481 r_local = eigvec_idx
482 v_local = virtual
483 IF (PRESENT(offsets)) THEN
484 DO ksp = SIZE(offsets), 1, -1
485 IF (eigvec_idx > offsets(ksp)) THEN
486 isp = ksp
487 EXIT
488 END IF
489 END DO
490 r_local = eigvec_idx - offsets(isp)
491 v_local = virtual_per_spin(isp)
492 END IF
493 kk = kk + 1
494 indx_local(3*kk - 2) = occ_of_ia(r_local, v_local)
495 indx_local(3*kk - 1) = virt_of_ia(r_local, v_local)
496 indx_local(3*kk) = isp
497 entries_local(kk) = eigvec_entry
498 END IF
499 END DO
500 END DO
501
502 ! one allgather of the counts, one allgatherv per array: every rank ends with the same
503 ! rank-ordered list
504 ALLOCATE (num_entries_to_comm(0:para_env%num_pe - 1), displ(0:para_env%num_pe - 1))
505 CALL para_env%allgather(num_entries_local, num_entries_to_comm)
506 displ(0) = 0
507 DO iproc = 1, para_env%num_pe - 1
508 displ(iproc) = displ(iproc - 1) + num_entries_to_comm(iproc - 1)
509 END DO
510 num_entries = sum(num_entries_to_comm)
511 ALLOCATE (eigvec_entries(num_entries), indx_all(3*num_entries))
512 CALL para_env%allgatherv(entries_local, eigvec_entries, num_entries_to_comm, displ)
513 CALL para_env%allgatherv(indx_local, indx_all, 3*num_entries_to_comm, 3*displ)
514
515 ALLOCATE (idx_homo(num_entries), idx_virt(num_entries))
516 IF (PRESENT(idx_spin)) ALLOCATE (idx_spin(num_entries))
517 DO kk = 1, num_entries
518 idx_homo(kk) = indx_all(3*kk - 2)
519 idx_virt(kk) = indx_all(3*kk - 1)
520 IF (PRESENT(idx_spin)) idx_spin(kk) = indx_all(3*kk)
521 END DO
522
523 DEALLOCATE (entries_local, indx_local, indx_all, num_entries_to_comm, displ)
524 NULLIFY (row_indices)
525 NULLIFY (col_indices)
526
527 !Now sort the results according to the involved singleparticle orbitals
528 ! (homo first, then virtual). idx_spin is payload, permuted alongside the entries.
529 IF (PRESENT(idx_spin)) THEN
530 CALL sort_excitations(idx_homo, idx_virt, eigvec_entries, idx_spin)
531 ELSE
532 CALL sort_excitations(idx_homo, idx_virt, eigvec_entries)
533 END IF
534
535 CALL timestop(handle)
536
537 END SUBROUTINE filter_eigvec_contrib
538
539! **************************************************************************************************
540!> \brief Spin-block layout for the open-shell (joint) BSE matrix: per-spin OV-pair counts and the
541!> block offsets into the joint matrix of dimension n_ov_joint = sum_sigma homo*virtual.
542!> \param homo_red per-spin (reduced) number of occupied levels
543!> \param virt_red per-spin (reduced) number of virtual levels
544!> \param n_ov per-spin OV-pair count (OUT)
545!> \param offsets per-spin block offset into the joint matrix (OUT)
546!> \param n_ov_joint total joint dimension (OUT)
547! **************************************************************************************************
548 SUBROUTINE get_bse_spin_block_layout(homo_red, virt_red, n_ov, offsets, n_ov_joint)
549 INTEGER, DIMENSION(:), INTENT(IN) :: homo_red, virt_red
550 INTEGER, DIMENSION(:), INTENT(OUT) :: n_ov, offsets
551 INTEGER, INTENT(OUT) :: n_ov_joint
552
553 INTEGER :: isp
554
555 n_ov_joint = 0
556 DO isp = 1, SIZE(homo_red)
557 offsets(isp) = n_ov_joint
558 n_ov(isp) = homo_red(isp)*virt_red(isp)
559 n_ov_joint = n_ov_joint + n_ov(isp)
560 END DO
561
562 END SUBROUTINE get_bse_spin_block_layout
563
564! **************************************************************************************************
565!> \brief The active MO window on the DFT axis: one absolute MO range for every spin; when the spin
566!> windows differ, the smallest range covering all of them
567!> \param eigenval_scf DFT eigenvalues (level, spin)
568!> \param n_mo number of MOs per spin
569!> \param homo occupied MOs per spin
570!> \param cutoff_occ ENERGY_CUTOFF_OCC in Hartree, zero or negative for no cutoff
571!> \param cutoff_empty ENERGY_CUTOFF_EMPTY in Hartree, zero or negative for no cutoff
572!> \param window first and last MO to keep
573!> \param warn warn when the spin windows differ; false for a second call in the same run
574! **************************************************************************************************
575 SUBROUTINE bse_window_from_dft(eigenval_scf, n_mo, homo, cutoff_occ, cutoff_empty, window, warn)
576
577 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval_scf
578 INTEGER, DIMENSION(:), INTENT(IN) :: n_mo, homo
579 REAL(kind=dp), INTENT(IN) :: cutoff_occ, cutoff_empty
580 TYPE(mo_window_type), INTENT(OUT) :: window
581 LOGICAL, INTENT(IN) :: warn
582
583 INTEGER :: ispin, nspins
584 LOGICAL :: spins_differ
585 TYPE(mo_window_type), ALLOCATABLE, DIMENSION(:) :: spin_windows
586
587 nspins = SIZE(homo)
588 ALLOCATE (spin_windows(nspins))
589 DO ispin = 1, nspins
590 CALL determine_mo_window(eigenval_scf(:, ispin), n_mo(ispin), homo(ispin), &
591 cutoff_occ, cutoff_empty, spin_windows(ispin))
592 END DO
593 IF (nspins > 1) THEN
594 CALL combine_mo_windows(spin_windows, window, spins_differ)
595 IF (spins_differ .AND. warn) THEN
596 cpwarn("BSE: spin-resolved active MO cutoff candidates differ; using combined window.")
597 END IF
598 ELSE
599 window = spin_windows(1)
600 END IF
601 DEALLOCATE (spin_windows)
602
603 END SUBROUTINE bse_window_from_dft
604
605! **************************************************************************************************
606!> \brief The active MO window of bse_window_from_dft, from the DFT eigenvalues of the MO sets
607!> \param mos MO sets per spin
608!> \param cutoff_occ ENERGY_CUTOFF_OCC in Hartree, zero or negative for no cutoff
609!> \param cutoff_empty ENERGY_CUTOFF_EMPTY in Hartree, zero or negative for no cutoff
610!> \param window first and last MO to keep, one absolute range for every spin
611!> \param warn warn when the spin windows differ
612! **************************************************************************************************
613 SUBROUTINE bse_window_from_mos(mos, cutoff_occ, cutoff_empty, window, warn)
614
615 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
616 REAL(kind=dp), INTENT(IN) :: cutoff_occ, cutoff_empty
617 TYPE(mo_window_type), INTENT(OUT) :: window
618 LOGICAL, INTENT(IN) :: warn
619
620 INTEGER :: ispin, nspins
621 INTEGER, ALLOCATABLE, DIMENSION(:) :: homo, n_mo
622 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenval_dft
623 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
624
625 nspins = SIZE(mos)
626 ALLOCATE (homo(nspins), n_mo(nspins))
627 DO ispin = 1, nspins
628 CALL get_mo_set(mos(ispin), nmo=n_mo(ispin), homo=homo(ispin), eigenvalues=eigenvalues)
629 IF (ispin == 1) ALLOCATE (eigenval_dft(n_mo(1), nspins))
630 eigenval_dft(:, ispin) = eigenvalues(1:n_mo(ispin))
631 END DO
632 CALL bse_window_from_dft(eigenval_dft, n_mo, homo, cutoff_occ, cutoff_empty, window, warn)
633 DEALLOCATE (homo, n_mo, eigenval_dft)
634
635 END SUBROUTINE bse_window_from_mos
636
637! **************************************************************************************************
638!> \brief Determines indices within the given energy cutoffs and truncates Eigenvalues and matrices
639!> \param fm_S_ia_full ...
640!> \param fm_S_ij_full ...
641!> \param fm_S_ab_full ...
642!> \param Eigenval_scf ...
643!> \param Eigenval ...
644!> \param Eigenval_reduced_scf DFT energies of the window, allocated here
645!> \param Eigenval_reduced_qp quasiparticle energies of the window, allocated here
646!> \param homo ...
647!> \param virtual ...
648!> \param dimen_RI ...
649!> \param unit_nr ...
650!> \param ispin spin of the slabs, selects the GW-corrected virtual count
651!> \param bse_env the BSE environment (settings; the GW-corrected level counts); receives the
652!> window's slabs fm_S_ia, fm_S_ij, fm_S_ab and its level counts homo_red, virt_red of spin ispin
653!> \param window Absolute MO window to retain.
654!> \param print_window Whether to print the cutoff summary.
655! **************************************************************************************************
656 SUBROUTINE truncate_bse_matrices(fm_S_ia_full, fm_S_ij_full, fm_S_ab_full, &
657 Eigenval_scf, Eigenval, Eigenval_reduced_scf, Eigenval_reduced_qp, &
658 homo, virtual, dimen_RI, unit_nr, &
659 ispin, &
660 bse_env, &
661 window, print_window)
662
663 TYPE(cp_fm_type), INTENT(IN) :: fm_s_ia_full, fm_s_ij_full, fm_s_ab_full
664 REAL(kind=dp), DIMENSION(:) :: eigenval_scf, eigenval
665 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenval_reduced_scf, eigenval_reduced_qp
666 INTEGER, INTENT(IN) :: homo, virtual, dimen_ri, unit_nr, ispin
667 TYPE(bse_env_type), INTENT(INOUT) :: bse_env
668 TYPE(mo_window_type), INTENT(IN) :: window
669 LOGICAL, INTENT(IN) :: print_window
670
671 CHARACTER(LEN=*), PARAMETER :: routinen = 'truncate_BSE_matrices'
672
673 INTEGER :: handle, homo_incl, homo_red, virt_incl, &
674 virt_red
675 TYPE(cp_blacs_env_type), POINTER :: context
676 TYPE(cp_fm_struct_type), POINTER :: fm_struct_ab, fm_struct_ia, fm_struct_ij
677 TYPE(mp_para_env_type), POINTER :: para_env
678
679 CALL timeset(routinen, handle)
680
681 homo_incl = window%first_mo
682 virt_incl = window%last_mo - homo
683 homo_red = homo - homo_incl + 1
684 virt_red = virt_incl
685
686 IF (print_window .AND. unit_nr > 0) THEN
687 IF (bse_env%bse_cutoff_occ > 0) THEN
688 WRITE (unit_nr, '(T2,A4,T7,A29,T71,F10.3)') 'BSE|', 'Cutoff occupied orbitals [eV]', &
689 bse_env%bse_cutoff_occ*evolt
690 ELSE
691 WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No cutoff given for occupied orbitals'
692 END IF
693 IF (bse_env%bse_cutoff_empty > 0) THEN
694 WRITE (unit_nr, '(T2,A4,T7,A26,T71,F10.3)') 'BSE|', 'Cutoff empty orbitals [eV]', &
695 bse_env%bse_cutoff_empty*evolt
696 ELSE
697 WRITE (unit_nr, '(T2,A4,T7,A34)') 'BSE|', 'No cutoff given for empty orbitals'
698 END IF
699 WRITE (unit_nr, '(T2,A4,T7,A20,T71,I10)') 'BSE|', 'First occupied index', homo_incl
700 WRITE (unit_nr, '(T2,A4,T7,A32,T71,I10)') 'BSE|', 'Last empty index (not MO index!)', virt_incl
701 WRITE (unit_nr, '(T2,A4,T7,A35,T71,F10.3)') 'BSE|', 'Energy of first occupied index [eV]', &
702 eigenval(homo_incl)*evolt
703 WRITE (unit_nr, '(T2,A4,T7,A31,T71,F10.3)') 'BSE|', 'Energy of last empty index [eV]', &
704 eigenval(homo + virt_incl)*evolt
705 WRITE (unit_nr, '(T2,A4,T7,A54,T71,F10.3)') 'BSE|', &
706 'Energy difference of first occupied index to HOMO [eV]', &
707 -(eigenval(homo_incl) - eigenval(homo))*evolt
708 WRITE (unit_nr, '(T2,A4,T7,A50,T71,F10.3)') 'BSE|', &
709 'Energy difference of last empty index to LUMO [eV]', &
710 (eigenval(homo + virt_incl) - eigenval(homo + 1))*evolt
711 WRITE (unit_nr, '(T2,A4,T7,A35,T71,I10)') 'BSE|', 'Number of GW-corrected occupied MOs', &
712 bse_env%gw_corr_lev_occ
713 WRITE (unit_nr, '(T2,A4,T7,A32,T71,I10)') 'BSE|', 'Number of GW-corrected empty MOs', &
714 bse_env%gw_corr_lev_virt(ispin)
715 WRITE (unit_nr, '(T2,A4)') 'BSE|'
716 END IF
717 IF (unit_nr > 0) THEN
718 IF (homo - homo_incl + 1 > bse_env%gw_corr_lev_occ) THEN
719 cpabort("Number of GW-corrected occupied MOs too small for chosen BSE cutoff")
720 END IF
721 IF (virt_incl > bse_env%gw_corr_lev_virt(ispin)) THEN
722 cpabort("Number of GW-corrected virtual MOs too small for chosen BSE cutoff")
723 END IF
724 END IF
725 !Truncate full fm_S matrices
726 !Allocate new truncated matrices of proper size
727 para_env => fm_s_ia_full%matrix_struct%para_env
728 context => fm_s_ia_full%matrix_struct%context
729
730 CALL cp_fm_struct_create(fm_struct_ia, para_env, context, dimen_ri, homo_red*virt_red)
731 CALL cp_fm_struct_create(fm_struct_ij, para_env, context, dimen_ri, homo_red*homo_red)
732 CALL cp_fm_struct_create(fm_struct_ab, para_env, context, dimen_ri, virt_red*virt_red)
733
734 CALL cp_fm_create(bse_env%fm_S_ia(ispin), fm_struct_ia, name="fm_S_trunc", set_zero=.true.)
735 CALL cp_fm_create(bse_env%fm_S_ij(ispin), fm_struct_ij, name="fm_S_ij_trunc", set_zero=.true.)
736 CALL cp_fm_create(bse_env%fm_S_ab(ispin), fm_struct_ab, name="fm_S_ab_trunc", set_zero=.true.)
737
738 !Copy parts of original matrices to truncated ones
739 IF (bse_env%bse_cutoff_occ > 0 .OR. bse_env%bse_cutoff_empty > 0) THEN
740 !Truncate eigenvals
741 ALLOCATE (eigenval_reduced_scf(homo_red + virt_red), eigenval_reduced_qp(homo_red + virt_red))
742 eigenval_reduced_scf(:) = eigenval_scf(homo_incl:homo + virt_incl)
743 eigenval_reduced_qp(:) = eigenval(homo_incl:homo + virt_incl)
744
745 CALL cp_fm_permute(bse_env%fm_S_ia(ispin), fm_s_ia_full, 1.0_dp, 1, virtual, 1, virt_red, [1, 2, 3, 4], &
746 first_index=[1, 1, homo_incl, 1], &
747 last_index=[dimen_ri, 1, (homo_incl) + homo_red - 1, virt_red])
748 CALL cp_fm_permute(bse_env%fm_S_ij(ispin), fm_s_ij_full, 1.0_dp, 1, homo, 1, homo_red, [1, 2, 3, 4], &
749 first_index=[1, 1, homo_incl, homo_incl], &
750 last_index=[dimen_ri, 1, (homo_incl) + homo_red - 1, (homo_incl) + homo_red - 1])
751 CALL cp_fm_permute(bse_env%fm_S_ab(ispin), fm_s_ab_full, 1.0_dp, 1, bse_env%gw_corr_lev_virt(ispin), 1, &
752 virt_red, [1, 2, 3, 4], first_index=[1, 1, 1, 1], &
753 last_index=[dimen_ri, 1, virt_red, virt_red])
754
755 ELSE
756 IF (unit_nr > 0) THEN
757 WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No truncation of BSE matrices applied'
758 WRITE (unit_nr, '(T2,A4)') 'BSE|'
759 END IF
760 ALLOCATE (eigenval_reduced_scf(homo_red + virt_red), eigenval_reduced_qp(homo_red + virt_red))
761 eigenval_reduced_scf(:) = eigenval_scf(:)
762 eigenval_reduced_qp(:) = eigenval(:)
763 CALL cp_fm_to_fm_submat_general(fm_s_ia_full, bse_env%fm_S_ia(ispin), dimen_ri, homo_red*virt_red, &
764 1, 1, 1, 1, context)
765 CALL cp_fm_to_fm_submat_general(fm_s_ij_full, bse_env%fm_S_ij(ispin), dimen_ri, homo_red*homo_red, &
766 1, 1, 1, 1, context)
767 CALL cp_fm_to_fm_submat_general(fm_s_ab_full, bse_env%fm_S_ab(ispin), dimen_ri, virt_red*virt_red, &
768 1, 1, 1, 1, context)
769 END IF
770
771 bse_env%homo_red(ispin) = homo_red
772 bse_env%virt_red(ispin) = virt_red
773
774 CALL cp_fm_struct_release(fm_struct_ia)
775 CALL cp_fm_struct_release(fm_struct_ij)
776 CALL cp_fm_struct_release(fm_struct_ab)
777
778 NULLIFY (para_env)
779 NULLIFY (context)
780
781 CALL timestop(handle)
782
783 END SUBROUTINE truncate_bse_matrices
784
785! **************************************************************************************************
786!> \brief ...
787!> \param fm_eigvec ...
788!> \param fm_eigvec_reshuffled ...
789!> \param homo ...
790!> \param virtual ...
791!> \param n_exc ...
792!> \param do_transpose ...
793! **************************************************************************************************
794 SUBROUTINE reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose)
795
796 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec
797 TYPE(cp_fm_type), INTENT(INOUT) :: fm_eigvec_reshuffled
798 INTEGER, INTENT(IN) :: homo, virtual, n_exc
799 LOGICAL, INTENT(IN) :: do_transpose
800
801 CHARACTER(LEN=*), PARAMETER :: routinen = 'reshuffle_eigvec'
802
803 INTEGER :: handle, my_m_col, my_n_row
804 INTEGER, DIMENSION(4) :: reordering
805 TYPE(cp_fm_struct_type), POINTER :: fm_struct_eigvec_col, &
806 fm_struct_eigvec_reshuffled
807 TYPE(cp_fm_type) :: fm_eigvec_col
808
809 CALL timeset(routinen, handle)
810
811 ! Define reordering:
812 ! (ia,11) to (a1,i1) for transposition
813 ! (ia,11) to (i1,a1) for default
814 IF (do_transpose) THEN
815 reordering = [2, 3, 1, 4]
816 my_n_row = virtual
817 my_m_col = homo
818 ELSE
819 reordering = [1, 3, 2, 4]
820 my_n_row = homo
821 my_m_col = virtual
822 END IF
823
824 CALL cp_fm_struct_create(fm_struct_eigvec_col, &
825 fm_eigvec%matrix_struct%para_env, fm_eigvec%matrix_struct%context, &
826 homo*virtual, 1)
827 CALL cp_fm_struct_create(fm_struct_eigvec_reshuffled, &
828 fm_eigvec%matrix_struct%para_env, fm_eigvec%matrix_struct%context, &
829 my_n_row, my_m_col)
830
831 ! Resort indices
832 CALL cp_fm_create(fm_eigvec_col, fm_struct_eigvec_col, name="BSE_column_vector")
833 CALL cp_fm_set_all(fm_eigvec_col, 0.0_dp)
834 CALL cp_fm_create(fm_eigvec_reshuffled, fm_struct_eigvec_reshuffled, name="BSE_reshuffled_eigenvector")
835 CALL cp_fm_set_all(fm_eigvec_reshuffled, 0.0_dp)
836 ! Fill matrix
837 CALL cp_fm_to_fm_submat(fm_eigvec, fm_eigvec_col, &
838 homo*virtual, 1, &
839 1, n_exc, &
840 1, 1)
841 ! Reshuffle
842 CALL cp_fm_permute(fm_eigvec_reshuffled, fm_eigvec_col, 1.0_dp, virtual, 1, 1, 1, reordering)
843
844 CALL cp_fm_release(fm_eigvec_col)
845 CALL cp_fm_struct_release(fm_struct_eigvec_col)
846 CALL cp_fm_struct_release(fm_struct_eigvec_reshuffled)
847
848 CALL timestop(handle)
849
850 END SUBROUTINE reshuffle_eigvec
851
852! **************************************************************************************************
853!> \brief Borrowed from the tddfpt module with slight adaptions
854!> \param bse_env the BSE environment (the QS objects the cube writer reads)
855!> \param mos ...
856!> \param istate ...
857!> \param info_approximation ...
858!> \param stride ...
859!> \param append_cube ...
860!> \param print_section ...
861! **************************************************************************************************
862 SUBROUTINE print_bse_nto_cubes(bse_env, mos, istate, info_approximation, &
863 stride, append_cube, print_section)
864
865 TYPE(bse_env_type), INTENT(IN) :: bse_env
866 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
867 INTEGER, INTENT(IN) :: istate
868 CHARACTER(LEN=10) :: info_approximation
869 INTEGER, DIMENSION(:), POINTER :: stride
870 LOGICAL, INTENT(IN) :: append_cube
871 TYPE(section_vals_type), POINTER :: print_section
872
873 CHARACTER(LEN=*), PARAMETER :: routinen = 'print_bse_nto_cubes'
874
875 CHARACTER(LEN=default_path_length) :: filename, info_approx_trunc, &
876 my_pos_cube, title
877 INTEGER :: handle, i, iset, nmo, unit_nr_cube
878 LOGICAL :: mpi_io
879 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
880 TYPE(cell_type), POINTER :: cell
881 TYPE(cp_fm_type), POINTER :: mo_coeff
882 TYPE(cp_logger_type), POINTER :: logger
883 TYPE(dft_control_type), POINTER :: dft_control
884 TYPE(particle_list_type), POINTER :: particles
885 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
886 TYPE(pw_c1d_gs_type) :: wf_g
887 TYPE(pw_env_type), POINTER :: pw_env
888 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
889 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
890 TYPE(pw_r3d_rs_type) :: wf_r
891 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
892 TYPE(qs_subsys_type), POINTER :: subsys
893
894 logger => cp_get_default_logger()
895 CALL timeset(routinen, handle)
896
897 dft_control => bse_env%dft_control
898 ! PROPERTIES%BSE runs after the GW released the grids; its requirement reader refuses cube files
899 IF (.NOT. ASSOCIATED(bse_env%pw_env)) THEN
900 cpabort("NTO cube files need the SCF grids: set NTO_ANALYSIS%CUBE_FILES to F")
901 END IF
902 pw_env => bse_env%pw_env
903 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
904 CALL auxbas_pw_pool%create_pw(wf_r)
905 CALL auxbas_pw_pool%create_pw(wf_g)
906
907 subsys => bse_env%subsys
908 CALL qs_subsys_get(subsys, particles=particles)
909
910 my_pos_cube = "REWIND"
911 IF (append_cube) THEN
912 my_pos_cube = "APPEND"
913 END IF
914
915 CALL qs_subsys_get(subsys, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
916 cell=cell, particle_set=particle_set)
917
918 DO iset = 1, 2
919 CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
920 DO i = 1, nmo
921 CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
922 cell, dft_control, particle_set, pw_env)
923 IF (iset == 1) THEN
924 WRITE (filename, '(A6,I3.3,A5,I2.2,a11)') "_NEXC_", istate, "_NTO_", i, "_Hole_State"
925 ELSE IF (iset == 2) THEN
926 WRITE (filename, '(A6,I3.3,A5,I2.2,a15)') "_NEXC_", istate, "_NTO_", i, "_Particle_State"
927 END IF
928 info_approx_trunc = trim(adjustl(info_approximation))
929 info_approx_trunc = info_approx_trunc(2:len_trim(info_approx_trunc) - 1)
930 filename = trim(info_approx_trunc)//trim(filename)
931 mpi_io = .true.
932 unit_nr_cube = cp_print_key_unit_nr(logger, print_section, '', extension=".cube", &
933 middle_name=trim(filename), file_position=my_pos_cube, &
934 log_filename=.false., ignore_should_output=.true., mpi_io=mpi_io)
935 IF (iset == 1) THEN
936 WRITE (title, *) "Natural Transition Orbital Hole State", i
937 ELSE IF (iset == 2) THEN
938 WRITE (title, *) "Natural Transition Orbital Particle State", i
939 END IF
940 CALL cp_pw_to_cube(wf_r, unit_nr_cube, title, particles=particles, stride=stride, mpi_io=mpi_io)
941 CALL cp_print_key_finished_output(unit_nr_cube, logger, print_section, '', &
942 ignore_should_output=.true., mpi_io=mpi_io)
943 END DO
944 END DO
945
946 CALL auxbas_pw_pool%give_back_pw(wf_g)
947 CALL auxbas_pw_pool%give_back_pw(wf_r)
948
949 CALL timestop(handle)
950 END SUBROUTINE print_bse_nto_cubes
951
952! **************************************************************************************************
953!> \brief Checks BSE input section and adapts them if necessary
954!> \param homo ...
955!> \param virtual ...
956!> \param unit_nr ...
957!> \param bse_env the BSE environment; the settings are normalised in place
958! **************************************************************************************************
959 SUBROUTINE adapt_bse_input_params(homo, virtual, unit_nr, bse_env)
960
961 INTEGER, INTENT(IN) :: homo, virtual, unit_nr
962 TYPE(bse_env_type), INTENT(INOUT) :: bse_env
963
964 CHARACTER(LEN=*), PARAMETER :: routinen = 'adapt_BSE_input_params'
965
966 INTEGER :: handle, i, j, n, n_exc_max, &
967 ndim_periodic_cell, &
968 ndim_periodic_poisson, &
969 num_state_list_exceptions
970 LOGICAL :: explicit
971 TYPE(cell_type), POINTER :: cell_ref
972
973 CALL timeset(routinen, handle)
974 ! the Poisson periodicity was read from the input (finish_bse_env): the grids may be released
975 NULLIFY (cell_ref)
976 CALL qs_subsys_get(bse_env%subsys, cell_ref=cell_ref)
977 ndim_periodic_poisson = bse_env%n_periodic_poisson
978 ndim_periodic_cell = sum(cell_ref%perd(1:3)) ! Borrowed from cell_methods.F/write_cell_low
979
980 ! The iterative solver computes the NUM_EXC_EN lowest excitations only
981 n_exc_max = homo*virtual
982 IF (bse_env%bse_diag_method == bse_iterdiag) THEN
983 n_exc_max = min(n_exc_max, bse_env%num_exc_en)
984 END IF
985
986 ! Handle negative NUM_PRINT_EXC
987 IF (bse_env%num_print_exc < 0 .OR. &
988 bse_env%num_print_exc > n_exc_max) THEN
989 bse_env%num_print_exc = n_exc_max
990 IF (unit_nr > 0) THEN
991 CALL cp_hint(__location__, &
992 "Keyword NUM_PRINT_EXC is either negative or too large. "// &
993 "Printing all computed excitations.")
994 END IF
995 END IF
996
997 ! Default to NUM_PRINT_EXC if too large or negative,
998 ! but only if NTOs are called - would be confusing for the user otherwise
999 ! Prepare and adapt user inputs for NTO analysis
1000 ! Logic: Explicit state list overrides NUM_PRINT_EXC_NTOS
1001 ! If only NUM_PRINT_EXC_NTOS is given, we write the array 1,...,NUM_PRINT_EXC_NTOS to
1002 ! bse_nto_state_list
1003 IF (bse_env%do_nto_analysis) THEN
1004 IF (bse_env%explicit_nto_list) THEN
1005 IF (bse_env%num_print_exc_ntos > 0) THEN
1006 IF (unit_nr > 0) THEN
1007 CALL cp_hint(__location__, &
1008 "Keywords NUM_PRINT_EXC_NTOS and STATE_LIST are both given in input. "// &
1009 "Overriding NUM_PRINT_EXC_NTOS.")
1010 END IF
1011 END IF
1012 ! Check if all states are within the range
1013 ! Count them and initialize new array afterwards
1014 num_state_list_exceptions = 0
1015 DO i = 1, SIZE(bse_env%bse_nto_state_list)
1016 IF (bse_env%bse_nto_state_list(i) < 1 .OR. &
1017 bse_env%bse_nto_state_list(i) > bse_env%num_print_exc) THEN
1018 num_state_list_exceptions = num_state_list_exceptions + 1
1019 END IF
1020 END DO
1021 IF (num_state_list_exceptions > 0) THEN
1022 IF (unit_nr > 0) THEN
1023 CALL cp_hint(__location__, &
1024 "STATE_LIST contains indices outside the range of included excitation levels. "// &
1025 "Ignoring these states.")
1026 END IF
1027 END IF
1028 n = SIZE(bse_env%bse_nto_state_list) - num_state_list_exceptions
1029 ALLOCATE (bse_env%bse_nto_state_list_final(n))
1030 bse_env%bse_nto_state_list_final(:) = 0
1031 i = 1
1032 DO j = 1, SIZE(bse_env%bse_nto_state_list)
1033 IF (bse_env%bse_nto_state_list(j) >= 1 .AND. &
1034 bse_env%bse_nto_state_list(j) <= bse_env%num_print_exc) THEN
1035 bse_env%bse_nto_state_list_final(i) = bse_env%bse_nto_state_list(j)
1036 i = i + 1
1037 END IF
1038 END DO
1039
1040 bse_env%num_print_exc_ntos = SIZE(bse_env%bse_nto_state_list_final)
1041 ELSE
1042 IF (bse_env%num_print_exc_ntos > bse_env%num_print_exc .OR. &
1043 bse_env%num_print_exc_ntos < 0) THEN
1044 bse_env%num_print_exc_ntos = bse_env%num_print_exc
1045 END IF
1046 ALLOCATE (bse_env%bse_nto_state_list_final(bse_env%num_print_exc_ntos))
1047 DO i = 1, bse_env%num_print_exc_ntos
1048 bse_env%bse_nto_state_list_final(i) = i
1049 END DO
1050 END IF
1051 END IF
1052
1053 ! Takes care of triplet states, when oscillator strengths are 0
1054 IF (bse_env%bse_spin_config /= 0 .AND. &
1055 bse_env%eps_nto_osc_str > 0) THEN
1056 IF (unit_nr > 0) THEN
1057 CALL cp_warn(__location__, &
1058 "Cannot apply EPS_OSC_STR for Triplet excitations. "// &
1059 "Resetting EPS_OSC_STR to default.")
1060 END IF
1061 bse_env%eps_nto_osc_str = -1.0_dp
1062 END IF
1063
1064 ! Take care of number for computed exciton descriptors
1065 IF (bse_env%num_print_exc_descr < 0 .OR. &
1066 bse_env%num_print_exc_descr > bse_env%num_print_exc) THEN
1067 IF (unit_nr > 0) THEN
1068 CALL cp_hint(__location__, &
1069 "Keyword NUM_PRINT_EXC_DESCR is either negative or too large. "// &
1070 "Printing exciton descriptors up to NUM_PRINT_EXC.")
1071 END IF
1072 bse_env%num_print_exc_descr = bse_env%num_print_exc
1073 END IF
1074
1075 ! the keyword defaults to true, so only an explicit request is worth a warning
1076 CALL section_vals_val_get(bse_env%bse_section, "PRINT_DIRECTIONAL_CROSSCORRELATION", explicit=explicit)
1077 IF (explicit .AND. bse_env%print_directional_crosscorrelation .AND. &
1078 .NOT. bse_env%print_directional_exc_descr) THEN
1079 IF (unit_nr > 0) THEN
1080 CALL cp_warn(__location__, &
1081 "PRINT_DIRECTIONAL_CROSSCORRELATION has no effect without PRINT_DIRECTIONAL_EXC_DESCR.")
1082 END IF
1083 END IF
1084
1085 ! Add warning for usage of KS energies
1086 IF (bse_env%use_ks_energies) THEN
1087 IF (unit_nr > 0) THEN
1088 CALL cp_warn(__location__, &
1089 "Using KS energies for BSE calculations. Therefore, no quantities "// &
1090 "of the preceeding GW calculation enter the BSE.")
1091 END IF
1092 END IF
1093
1094 ! Add warning if periodic calculation is invoked
1095 IF (ndim_periodic_poisson /= 0) THEN
1096 IF (unit_nr > 0) THEN
1097 CALL cp_warn(__location__, &
1098 "Poisson solver should be invoked by PERIODIC NONE. "// &
1099 "The applied length gauge might give misleading results for "// &
1100 "oscillator strengths.")
1101 END IF
1102 END IF
1103 IF (ndim_periodic_cell /= 0) THEN
1104 IF (unit_nr > 0) THEN
1105 CALL cp_warn(__location__, &
1106 "CELL in SUBSYS should be invoked with PERIODIC NONE. "// &
1107 "The applied length gauge might give misleading results for "// &
1108 "oscillator strengths.")
1109 END IF
1110 END IF
1111
1112 CALL timestop(handle)
1113 END SUBROUTINE adapt_bse_input_params
1114
1115! **************************************************************************************************
1116!> \brief Resolve SCREENING_IN_W and its factor before constructing the dielectric matrix or slabs.
1117!> \param bse_env screening method and factor, normalised in place
1118!> \param unit_nr output unit, positive on the writing rank only
1119! **************************************************************************************************
1120 SUBROUTINE adapt_bse_screening(bse_env, unit_nr)
1121 TYPE(bse_env_type), INTENT(INOUT) :: bse_env
1122 INTEGER, INTENT(IN) :: unit_nr
1123
1124 ! Handle screening factor options
1125 IF (bse_env%screening_factor > 0.0_dp) THEN
1126 IF (bse_env%screening_method /= bse_screening_alpha) THEN
1127 IF (unit_nr > 0) THEN
1128 CALL cp_warn(__location__, &
1129 "Screening factor is only supported for &SCREENING_IN_W ALPHA. "// &
1130 "Resetting SCREENING_IN_W to ALPHA.")
1131 END IF
1132 bse_env%screening_method = bse_screening_alpha
1133 END IF
1134 IF (bse_env%screening_factor > 1.0_dp) THEN
1135 IF (unit_nr > 0) THEN
1136 CALL cp_warn(__location__, &
1137 "Screening factor is larger than 1.0. ")
1138 END IF
1139 END IF
1140 END IF
1141
1142 IF (bse_env%screening_factor < 0.0_dp .AND. &
1143 bse_env%screening_method == bse_screening_alpha) THEN
1144 IF (unit_nr > 0) THEN
1145 CALL cp_warn(__location__, &
1146 "Screening factor is negative. Defaulting to 0.25")
1147 END IF
1148 bse_env%screening_factor = 0.25_dp
1149 END IF
1150
1151 IF (bse_env%screening_factor == 0.0_dp) THEN
1152 ! Use RPA internally in this case
1153 bse_env%screening_method = bse_screening_rpa
1154 END IF
1155 IF (bse_env%screening_factor == 1.0_dp) THEN
1156 ! Use TDHF internally in this case
1157 bse_env%screening_method = bse_screening_tdhf
1158 END IF
1159
1160 END SUBROUTINE adapt_bse_screening
1161
1162! **************************************************************************************************
1163
1164! **************************************************************************************************
1165!> \brief The multipoles in the MO window, D^k_pq = sum_µν C_µp M^k_µν C_νq, cut to the ai, ij and ab
1166!> blocks of the window; M^k_µν are the AO multipoles of get_multipoles_ao
1167!> \param fm_multipole_ai_trunc D^k_ai, one entry per component k
1168!> \param fm_multipole_ij_trunc D^k_ij
1169!> \param fm_multipole_ab_trunc D^k_ab
1170!> \param matrix_multipole M^k_µν, one entry per component k (x, y, z, then xx, xy, xz, yy, yz, zz)
1171!> \param mos MO sets of the reference, for the AO layout and the occupied count
1172!> \param mo_coeff C_µp of the spin
1173!> \param homo_red occupied levels in the window
1174!> \param virtual_red virtual levels in the window
1175!> \param context_BSE ...
1176!> \param ispin spin channel whose mo_set supplies homo/nao (default 1); open-shell beta needs 2
1177! **************************************************************************************************
1178 SUBROUTINE get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, &
1179 matrix_multipole, mos, mo_coeff, &
1180 homo_red, virtual_red, context_BSE, ispin)
1181
1182 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
1183 INTENT(INOUT) :: fm_multipole_ai_trunc, &
1184 fm_multipole_ij_trunc, &
1185 fm_multipole_ab_trunc
1186 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: matrix_multipole
1187 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1188 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1189 INTEGER, INTENT(IN) :: homo_red, virtual_red
1190 TYPE(cp_blacs_env_type), POINTER :: context_bse
1191 INTEGER, INTENT(IN), OPTIONAL :: ispin
1192
1193 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_multipoles_mo'
1194
1195 INTEGER :: handle, idir, my_ispin, n_multipole, &
1196 n_occ, n_virt, nao, nmo_mp2
1197 TYPE(cp_fm_struct_type), POINTER :: fm_struct_mp_ab_trunc, fm_struct_mp_ai_trunc, &
1198 fm_struct_mp_ij_trunc, fm_struct_multipoles_ao, fm_struct_nao_nmo, fm_struct_nmo_nmo
1199 TYPE(cp_fm_type) :: fm_work
1200 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_multipole_per_dir
1201 TYPE(mp_para_env_type), POINTER :: para_env_bse
1202
1203 CALL timeset(routinen, handle)
1204
1205 my_ispin = 1
1206 IF (PRESENT(ispin)) my_ispin = ispin
1207
1208 ! Use the same blacs environment as for the MO coefficients to ensure correct multiplication dbcsr x fm later on
1209 fm_struct_multipoles_ao => mos(my_ispin)%mo_coeff%matrix_struct
1210 ! BSE has different contexts and blacsenvs
1211 para_env_bse => context_bse%para_env
1212 n_multipole = SIZE(matrix_multipole)
1213 ALLOCATE (fm_multipole_per_dir(n_multipole))
1214
1215 ! Now we transform them to MO
1216 ! n_occ is the number of occupied MOs, nao the number of all AOs
1217 ! Writing homo to n_occ instead if nmo,
1218 ! takes care of ADDED_MOS, which would overwrite nmo of qs_env-mos, if invoked
1219 CALL get_mo_set(mo_set=mos(my_ispin), homo=n_occ, nao=nao)
1220 ! Takes into account removed nullspace values from SVD
1221 nmo_mp2 = mo_coeff(1)%matrix_struct%ncol_global
1222 n_virt = nmo_mp2 - n_occ
1223
1224 ! At the end, we need four different layouts of matrices in this multiplication, e.g. for a dipole:
1225 ! D_pq = full multipole matrix for occupied and unoccupied
1226 ! Final result:D_pq= C_{mu p} <\mu|\vec{r}|\nu> C_{\nu q} EQ.I
1227 ! \_______/ \___________/ \______/
1228 ! fm_coeff matrix_multipole fm_coeff
1229 ! (EQ.Ia) (EQ.Ib) (EQ.Ia)
1230 ! Intermediate work matrices:
1231 ! fm_work = <\mu|\vec{r}|\nu> C_{\nu q} EQ.II
1232
1233 ! Struct for the full multipole matrix
1234 CALL cp_fm_struct_create(fm_struct_nao_nmo, &
1235 fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, &
1236 nao, nmo_mp2)
1237 CALL cp_fm_struct_create(fm_struct_nmo_nmo, &
1238 fm_struct_multipoles_ao%para_env, fm_struct_multipoles_ao%context, &
1239 nmo_mp2, nmo_mp2)
1240
1241 ! At the very end, we copy the multipoles corresponding to truncated BSE indices in i and a
1242 CALL cp_fm_struct_create(fm_struct_mp_ai_trunc, para_env_bse, &
1243 context_bse, virtual_red, homo_red)
1244 CALL cp_fm_struct_create(fm_struct_mp_ij_trunc, para_env_bse, &
1245 context_bse, homo_red, homo_red)
1246 CALL cp_fm_struct_create(fm_struct_mp_ab_trunc, para_env_bse, &
1247 context_bse, virtual_red, virtual_red)
1248 DO idir = 1, n_multipole
1249 CALL cp_fm_create(fm_multipole_ai_trunc(idir), matrix_struct=fm_struct_mp_ai_trunc, &
1250 name="dipoles_mo_ai_trunc")
1251 CALL cp_fm_set_all(fm_multipole_ai_trunc(idir), 0.0_dp)
1252 CALL cp_fm_create(fm_multipole_ij_trunc(idir), matrix_struct=fm_struct_mp_ij_trunc, &
1253 name="dipoles_mo_ij_trunc")
1254 CALL cp_fm_set_all(fm_multipole_ij_trunc(idir), 0.0_dp)
1255 CALL cp_fm_create(fm_multipole_ab_trunc(idir), matrix_struct=fm_struct_mp_ab_trunc, &
1256 name="dipoles_mo_ab_trunc")
1257 CALL cp_fm_set_all(fm_multipole_ab_trunc(idir), 0.0_dp)
1258 END DO
1259
1260 ! Need another temporary matrix to store intermediate result from right multiplication
1261 ! D = C_{mu a} <\mu|\vec{r}|\nu> C_{\nu i}
1262 CALL cp_fm_create(fm_work, matrix_struct=fm_struct_nao_nmo, name="multipole_work")
1263 CALL cp_fm_set_all(fm_work, 0.0_dp)
1264
1265 DO idir = 1, n_multipole
1266 ! Create the full multipole matrix per direction
1267 CALL cp_fm_create(fm_multipole_per_dir(idir), matrix_struct=fm_struct_nmo_nmo, name="multipoles_mo")
1268 CALL cp_fm_set_all(fm_multipole_per_dir(idir), 0.0_dp)
1269 ! Fill final (MO) multipole matrix
1270 CALL cp_dbcsr_sm_fm_multiply(matrix_multipole(idir)%matrix, mo_coeff(1), &
1271 fm_work, ncol=nmo_mp2)
1272 ! Now obtain the multipoles by the final multiplication;
1273 ! We do that inside the loop to obtain multipoles per axis for print
1274 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))
1275
1276 ! Truncate full matrix to the BSE indices
1277 ! D_ai
1278 CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1279 fm_multipole_ai_trunc(idir), &
1280 virtual_red, &
1281 homo_red, &
1282 n_occ + 1, &
1283 n_occ - homo_red + 1, &
1284 1, &
1285 1, &
1286 fm_multipole_per_dir(idir)%matrix_struct%context)
1287 ! D_ij
1288 CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1289 fm_multipole_ij_trunc(idir), &
1290 homo_red, &
1291 homo_red, &
1292 n_occ - homo_red + 1, &
1293 n_occ - homo_red + 1, &
1294 1, &
1295 1, &
1296 fm_multipole_per_dir(idir)%matrix_struct%context)
1297 ! D_ab
1298 CALL cp_fm_to_fm_submat_general(fm_multipole_per_dir(idir), &
1299 fm_multipole_ab_trunc(idir), &
1300 virtual_red, &
1301 virtual_red, &
1302 n_occ + 1, &
1303 n_occ + 1, &
1304 1, &
1305 1, &
1306 fm_multipole_per_dir(idir)%matrix_struct%context)
1307 END DO
1308
1309 !Release matrices and structs
1310 NULLIFY (fm_struct_multipoles_ao)
1311 CALL cp_fm_struct_release(fm_struct_mp_ai_trunc)
1312 CALL cp_fm_struct_release(fm_struct_mp_ij_trunc)
1313 CALL cp_fm_struct_release(fm_struct_mp_ab_trunc)
1314 CALL cp_fm_struct_release(fm_struct_nao_nmo)
1315 CALL cp_fm_struct_release(fm_struct_nmo_nmo)
1316 DO idir = 1, n_multipole
1317 CALL cp_fm_release(fm_multipole_per_dir(idir))
1318 END DO
1319 DEALLOCATE (fm_multipole_per_dir)
1320 CALL cp_fm_release(fm_work)
1321
1322 CALL timestop(handle)
1323
1324 END SUBROUTINE get_multipoles_mo
1325
1326! **************************************************************************************************
1327!> \brief Computes trace of form Tr{A^T B C} for exciton descriptors
1328!> \param fm_A Full Matrix, typically X or Y, in format homo x virtual
1329!> \param fm_B ...
1330!> \param fm_C ...
1331!> \param alpha ...
1332! **************************************************************************************************
1333 SUBROUTINE trace_exciton_descr(fm_A, fm_B, fm_C, alpha)
1334
1335 TYPE(cp_fm_type), INTENT(IN) :: fm_a, fm_b, fm_c
1336 REAL(kind=dp), INTENT(OUT) :: alpha
1337
1338 CHARACTER(LEN=*), PARAMETER :: routinen = 'trace_exciton_descr'
1339
1340 INTEGER :: handle, ncol_a, ncol_b, ncol_c, nrow_a, &
1341 nrow_b, nrow_c
1342 TYPE(cp_fm_type) :: fm_work_ia
1343
1344 CALL timeset(routinen, handle)
1345
1346 CALL cp_fm_create(fm_work_ia, fm_a%matrix_struct)
1347 CALL cp_fm_get_info(fm_a, nrow_global=nrow_a, ncol_global=ncol_a)
1348 CALL cp_fm_get_info(fm_b, nrow_global=nrow_b, ncol_global=ncol_b)
1349 CALL cp_fm_get_info(fm_c, nrow_global=nrow_c, ncol_global=ncol_c)
1350
1351 ! Check matrix sizes
1352 cpassert(nrow_a == nrow_b .AND. ncol_a == ncol_c .AND. ncol_b == nrow_c)
1353
1354 CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
1355
1356 CALL parallel_gemm("N", "N", nrow_a, ncol_a, nrow_c, 1.0_dp, &
1357 fm_b, fm_c, 0.0_dp, fm_work_ia)
1358
1359 CALL cp_fm_trace(fm_a, fm_work_ia, alpha)
1360
1361 CALL cp_fm_release(fm_work_ia)
1362
1363 CALL timestop(handle)
1364
1365 END SUBROUTINE trace_exciton_descr
1366
1367! **************************************************************************************************
1368!> \brief Column-concatenate per-spin ia-slabs into the joint dimen_RI x n_ov_joint slab.
1369!> Sigma-block of spin isp occupies columns offsets(isp)+1 .. offsets(isp)+n_ov(isp).
1370!> fm_S_joint must be pre-created and zeroed by the caller.
1371!> \param fm_S_ia per-spin ia-slabs, shape (dimen_RI, n_ov(isp)) per spin
1372!> \param offsets per-spin column offsets into fm_S_joint (0-based)
1373!> \param n_ov per-spin OV-pair counts
1374!> \param dimen_RI RI auxiliary basis dimension (row count)
1375!> \param fm_S_joint pre-created output slab (dimen_RI x n_ov_joint)
1376! **************************************************************************************************
1377 SUBROUTINE assemble_joint_ov_slab(fm_S_ia, offsets, n_ov, dimen_RI, fm_S_joint)
1378 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_s_ia
1379 INTEGER, DIMENSION(:), INTENT(IN) :: offsets, n_ov
1380 INTEGER, INTENT(IN) :: dimen_ri
1381 TYPE(cp_fm_type), INTENT(INOUT) :: fm_s_joint
1382
1383 CHARACTER(LEN=*), PARAMETER :: routinen = 'assemble_joint_ov_slab'
1384
1385 INTEGER :: handle, isp
1386
1387 CALL timeset(routinen, handle)
1388 DO isp = 1, SIZE(fm_s_ia)
1389 CALL cp_fm_to_fm_submat(fm_s_ia(isp), fm_s_joint, dimen_ri, n_ov(isp), 1, 1, 1, offsets(isp) + 1)
1390 END DO
1391 CALL timestop(handle)
1392
1393 END SUBROUTINE assemble_joint_ov_slab
1394
1395END MODULE bse_util
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:79
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
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_types
The BSE environment: the settings of the &BSE section and the state a GW path prepares for the solver...
Definition bse_types.F:14
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:1378
bse_util::adapt_bse_screening
subroutine, public adapt_bse_screening(bse_env, unit_nr)
Resolve SCREENING_IN_W and its factor before constructing the dielectric matrix or slabs.
Definition bse_util.F:1121
bse_util::truncate_bse_matrices
subroutine, public truncate_bse_matrices(fm_s_ia_full, fm_s_ij_full, fm_s_ab_full, eigenval_scf, eigenval, eigenval_reduced_scf, eigenval_reduced_qp, homo, virtual, dimen_ri, unit_nr, ispin, bse_env, window, print_window)
Determines indices within the given energy cutoffs and truncates Eigenvalues and matrices.
Definition bse_util.F:662
bse_util::get_multipoles_mo
subroutine, public get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, matrix_multipole, mos, mo_coeff, homo_red, virtual_red, context_bse, ispin)
The multipoles in the MO window, D^k_pq = sum_µν C_µp M^k_µν C_νq, cut to the ai, ij and ab blocks of...
Definition bse_util.F:1181
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:354
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:1334
bse_util::adapt_bse_input_params
subroutine, public adapt_bse_input_params(homo, virtual, unit_nr, bse_env)
Checks BSE input section and adapts them if necessary.
Definition bse_util.F:960
bse_util::occ_of_ia
pure integer function, public occ_of_ia(ia, virt)
Occupied level of the compound transition index ia = (i-1)*virt + a.
Definition bse_util.F:110
bse_util::fulldiag_memory_estimate_gb
pure real(kind=dp) function, public fulldiag_memory_estimate_gb(n_ov_joint, do_abba)
Peak memory of the full diagonalization in GB, all ranks together: n_ov^2 doubles times the matrices ...
Definition bse_util.F:397
bse_util::filter_eigvec_contrib
subroutine, public filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, i_exc, virtual, num_entries, bse_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:431
bse_util::bse_window_from_dft
subroutine, public bse_window_from_dft(eigenval_scf, n_mo, homo, cutoff_occ, cutoff_empty, window, warn)
The active MO window on the DFT axis: one absolute MO range for every spin; when the spin windows dif...
Definition bse_util.F:576
bse_util::reshuffle_eigvec
subroutine, public reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose)
...
Definition bse_util.F:795
bse_util::ia_of_occ_virt
pure integer function, public ia_of_occ_virt(i_occ, a_virt, virt)
Compound transition index ia = (i-1)*virt + a of the pair (i, a): a is the fast index.
Definition bse_util.F:95
bse_util::screen_slabs
subroutine, public screen_slabs(fm_eps_inv, fm_s_ij, fm_s_ia, dimen_ri_red, homo, virtual, fm_s_bar_ij, fm_s_bar_ia, do_abba)
The screened slabs \bar{B}^P_ij = sum_Q [1+Q(0)]^-1_PQ B^Q_ij and \bar{B}^P_ia = sum_Q [1+Q(0)]^-1_PQ...
Definition bse_util.F:148
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:549
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:257
bse_util::bse_window_from_mos
subroutine, public bse_window_from_mos(mos, cutoff_occ, cutoff_empty, window, warn)
The active MO window of bse_window_from_dft, from the DFT eigenvalues of the MO sets.
Definition bse_util.F:614
bse_util::virt_of_ia
pure integer function, public virt_of_ia(ia, virt)
Virtual level of the compound transition index ia = (i-1)*virt + a.
Definition bse_util.F:125
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:190
bse_util::print_bse_nto_cubes
subroutine, public print_bse_nto_cubes(bse_env, mos, istate, info_approximation, stride, append_cube, print_section)
Borrowed from the tddfpt module with slight adaptions.
Definition bse_util.F:864
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_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:940
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_redistribution
Redistribution of full matrices and distributed arrays.
Definition cp_fm_redistribution.F:11
cp_fm_redistribution::cp_fm_permute
subroutine, public cp_fm_permute(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, nrow_secidx_out, ncol_secidx_out, reordering, row_offset, col_offset, first_index, last_index)
Add a permuted or restricted four-index tensor stored as a full matrix.
Definition cp_fm_redistribution.F:51
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:1046
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:2049
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:1518
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:536
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:167
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_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:1156
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, file_access)
...
Definition cp_output_handling.F:935
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::bse_screening_tdhf
integer, parameter, public bse_screening_tdhf
Definition input_constants.F:1415
input_constants::bse_iterdiag
integer, parameter, public bse_iterdiag
Definition input_constants.F:1406
input_constants::bse_screening_alpha
integer, parameter, public bse_screening_alpha
Definition input_constants.F:1415
input_constants::bse_screening_rpa
integer, parameter, public bse_screening_rpa
Definition input_constants.F:1415
input_constants::bse_abba
integer, parameter, public bse_abba
Definition input_constants.F:1410
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
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
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:27
mo_window
Common selection and union operations for contiguous molecular-orbital windows.
Definition mo_window.F:11
mo_window::combine_mo_windows
subroutine, public combine_mo_windows(spin_windows, combined_window, windows_differ)
Forms the smallest contiguous window covering all supplied spin windows.
Definition mo_window.F:98
mo_window::determine_mo_window
subroutine, public determine_mo_window(eigenvalues, n_mo, n_occ, cutoff_occ, cutoff_empty, window)
Selects one contiguous MO window from an ordered reference spectrum.
Definition mo_window.F:44
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_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:2845
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, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:440
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
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
bse_types::bse_env_type
Settings of the &BSE section (read_bse_section, re-read before every BSE run, normalised in place by ...
Definition bse_types.F:44
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:55
cp_control_types::dft_control_type
Definition cp_control_types.F:786
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:116
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:774
mo_window::mo_window_type
Definition mo_window.F:21
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_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_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:53
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56