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