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, 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 :
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 41664 : 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 41664 : ia = (i_occ - 1)*virt + a_virt
100 :
101 41664 : 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 18175 : PURE FUNCTION occ_of_ia(ia, virt) RESULT(i_occ)
110 :
111 : INTEGER, INTENT(IN) :: ia, virt
112 : INTEGER :: i_occ
113 :
114 18175 : i_occ = (ia - 1)/virt + 1
115 :
116 18175 : 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 18175 : PURE FUNCTION virt_of_ia(ia, virt) RESULT(a_virt)
125 :
126 : INTEGER, INTENT(IN) :: ia, virt
127 : INTEGER :: a_virt
128 :
129 18175 : a_virt = MOD(ia - 1, virt) + 1
130 :
131 18175 : 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 140 : 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 70 : CALL timeset(routineN, handle)
159 :
160 70 : CALL cp_fm_create(fm_S_bar_ij, fm_S_ij%matrix_struct)
161 70 : 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 70 : matrix_a=fm_eps_inv, matrix_b=fm_S_ij, beta=0.0_dp, matrix_c=fm_S_bar_ij)
166 :
167 70 : IF (do_abba) THEN
168 28 : CALL cp_fm_create(fm_S_bar_ia, fm_S_ia%matrix_struct)
169 28 : 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 28 : 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 70 : CALL timestop(handle)
177 :
178 70 : 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 104 : 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 52 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
203 : LOGICAL :: my_do_transpose
204 : REAL(KIND=dp) :: coeff, my_gamma
205 :
206 52 : CALL timeset(routineN, handle)
207 :
208 52 : IF (PRESENT(gamma)) THEN
209 52 : my_gamma = gamma
210 : ELSE
211 : my_gamma = 2.0_dp
212 : END IF
213 :
214 52 : IF (PRESENT(do_transpose)) THEN
215 52 : 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 52 : col_indices=col_indices)
225 :
226 52 : IF (my_do_transpose) THEN
227 4080 : DO jj = 1, ncol_local
228 4028 : j_col_global = col_indices(jj)
229 1548600 : DO ii = 1, nrow_local
230 1544520 : coeff = (eig_vals(j_col_global)**beta)/my_gamma
231 1548548 : fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff
232 : END DO
233 : END DO
234 : ELSE
235 0 : DO jj = 1, ncol_local
236 0 : DO ii = 1, nrow_local
237 0 : i_row_global = row_indices(ii)
238 0 : coeff = (eig_vals(i_row_global)**beta)/my_gamma
239 0 : fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff
240 : END DO
241 : END DO
242 : END IF
243 :
244 52 : CALL timestop(handle)
245 :
246 52 : 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 1634 : 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 1634 : INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim_work, idx_sec_work, &
266 1634 : idx_spin_work, tmp_index
267 : LOGICAL :: unique_entries
268 1634 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries_work
269 :
270 1634 : CALL timeset(routineN, handle)
271 :
272 1634 : num_entries = SIZE(idx_prim)
273 :
274 4326 : ALLOCATE (tmp_index(num_entries))
275 :
276 1634 : CALL sort(idx_prim, num_entries, tmp_index)
277 :
278 2692 : ALLOCATE (idx_sec_work(num_entries))
279 4326 : ALLOCATE (eigvec_entries_work(num_entries))
280 1830 : IF (PRESENT(idx_spin)) ALLOCATE (idx_spin_work(num_entries))
281 :
282 4664 : DO ii = 1, num_entries
283 3030 : idx_sec_work(ii) = idx_sec(tmp_index(ii))
284 3030 : eigvec_entries_work(ii) = eigvec_entries(tmp_index(ii))
285 4664 : IF (PRESENT(idx_spin)) idx_spin_work(ii) = idx_spin(tmp_index(ii))
286 : END DO
287 :
288 1634 : DEALLOCATE (tmp_index)
289 1634 : DEALLOCATE (idx_sec)
290 1634 : DEALLOCATE (eigvec_entries)
291 :
292 1634 : CALL MOVE_ALLOC(idx_sec_work, idx_sec)
293 1634 : CALL MOVE_ALLOC(eigvec_entries_work, eigvec_entries)
294 1634 : IF (PRESENT(idx_spin)) THEN
295 250 : DEALLOCATE (idx_spin)
296 250 : 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 1634 : CALL sort_unique(idx_prim, unique_entries)
301 1634 : IF (.NOT. unique_entries) THEN
302 928 : ALLOCATE (idx_prim_work(num_entries))
303 2474 : idx_prim_work(:) = idx_prim(:)
304 : ! Find duplicate entries in idx_prim
305 2474 : DO ii = 1, num_entries
306 2010 : IF (idx_prim_work(ii) == 0) CYCLE
307 7180 : num_mults = COUNT(idx_prim_work == idx_prim_work(ii))
308 1178 : IF (num_mults > 1) THEN
309 : !Set all duplicate entries to 0
310 2136 : idx_prim_work(ii:ii + num_mults - 1) = 0
311 : !Start sorting in secondary index
312 1956 : ALLOCATE (idx_sec_work(num_mults))
313 1956 : ALLOCATE (eigvec_entries_work(num_mults))
314 2136 : idx_sec_work(:) = idx_sec(ii:ii + num_mults - 1)
315 2136 : eigvec_entries_work(:) = eigvec_entries(ii:ii + num_mults - 1)
316 652 : IF (PRESENT(idx_spin)) THEN
317 400 : ALLOCATE (idx_spin_work(num_mults))
318 666 : idx_spin_work(:) = idx_spin(ii:ii + num_mults - 1)
319 : END IF
320 1304 : ALLOCATE (tmp_index(num_mults))
321 652 : CALL sort(idx_sec_work, num_mults, tmp_index)
322 :
323 : !Now write newly sorted indices to original arrays
324 2136 : DO kk = ii, ii + num_mults - 1
325 1484 : idx_sec(kk) = idx_sec_work(kk - ii + 1)
326 1484 : eigvec_entries(kk) = eigvec_entries_work(tmp_index(kk - ii + 1))
327 2136 : IF (PRESENT(idx_spin)) idx_spin(kk) = idx_spin_work(tmp_index(kk - ii + 1))
328 : END DO
329 : !Deallocate work arrays
330 652 : DEALLOCATE (tmp_index)
331 652 : DEALLOCATE (idx_sec_work)
332 652 : DEALLOCATE (eigvec_entries_work)
333 652 : IF (PRESENT(idx_spin)) DEALLOCATE (idx_spin_work)
334 : END IF
335 2474 : idx_prim_work(ii) = idx_prim(ii)
336 : END DO
337 464 : DEALLOCATE (idx_prim_work)
338 : END IF
339 :
340 1634 : CALL timestop(handle)
341 :
342 4902 : 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 54 : 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 54 : CALL timeset(routineN, handle)
366 :
367 54 : mem_est = fulldiag_memory_estimate_GB(n_ov_joint, bse_abba)
368 54 : mem_est_per_rank = REAL(mem_est/para_env%num_pe, KIND=dp)
369 :
370 54 : 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 27 : WRITE (unit_nr, '(T2,A4,T7,A40,T68,ES13.3)') 'BSE|', 'Total peak memory estimate from BSE [GB]', &
374 54 : mem_est
375 27 : WRITE (unit_nr, '(T2,A4,T7,A47,T68,F13.3)') 'BSE|', 'Peak memory estimate per MPI rank from BSE [GB]', &
376 54 : mem_est_per_rank
377 27 : 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 54 : 10*32/REAL(para_env%num_pe, KIND=dp)
383 :
384 54 : CALL timestop(handle)
385 :
386 54 : 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 108 : 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 108 : num_BSE_matrices = 2
406 108 : IF (do_abba) num_BSE_matrices = 10
407 108 : 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 108 : mem_GB = 8.0E-9_dp*REAL(full_dim, KIND=dp)
410 :
411 108 : 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 1634 : SUBROUTINE filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, &
429 : i_exc, virtual, num_entries, bse_env, &
430 1634 : 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 1634 : INTEGER, ALLOCATABLE, DIMENSION(:) :: displ, indx_all, indx_local, &
446 1634 : num_entries_to_comm
447 1634 : 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 1634 : CALL timeset(routineN, handle)
453 :
454 1634 : 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 1634 : 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 1634 : num_entries_local = 0
464 104460 : DO jj = 1, ncol_local
465 102826 : IF (col_indices(jj) /= i_exc) CYCLE
466 70140 : DO ii = 1, nrow_local
467 169738 : 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 6606 : ALLOCATE (entries_local(num_entries_local), indx_local(3*num_entries_local))
471 1634 : kk = 0
472 104460 : DO jj = 1, ncol_local
473 102826 : IF (col_indices(jj) /= i_exc) CYCLE
474 70140 : DO ii = 1, nrow_local
475 66912 : eigvec_idx = row_indices(ii)
476 66912 : eigvec_entry = fm_eigvec%local_data(ii, jj)
477 169738 : 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 1515 : isp = 1
481 1515 : r_local = eigvec_idx
482 1515 : v_local = virtual
483 1515 : IF (PRESENT(offsets)) THEN
484 417 : DO ksp = SIZE(offsets), 1, -1
485 417 : IF (eigvec_idx > offsets(ksp)) THEN
486 : isp = ksp
487 : EXIT
488 : END IF
489 : END DO
490 286 : r_local = eigvec_idx - offsets(isp)
491 286 : v_local = virtual_per_spin(isp)
492 : END IF
493 1515 : kk = kk + 1
494 1515 : indx_local(3*kk - 2) = occ_of_ia(r_local, v_local)
495 1515 : indx_local(3*kk - 1) = virt_of_ia(r_local, v_local)
496 1515 : indx_local(3*kk) = isp
497 1515 : 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 6536 : ALLOCATE (num_entries_to_comm(0:para_env%num_pe - 1), displ(0:para_env%num_pe - 1))
505 1634 : CALL para_env%allgather(num_entries_local, num_entries_to_comm)
506 1634 : displ(0) = 0
507 3268 : DO iproc = 1, para_env%num_pe - 1
508 3268 : displ(iproc) = displ(iproc - 1) + num_entries_to_comm(iproc - 1)
509 : END DO
510 4902 : num_entries = SUM(num_entries_to_comm)
511 7018 : ALLOCATE (eigvec_entries(num_entries), indx_all(3*num_entries))
512 1634 : CALL para_env%allgatherv(entries_local, eigvec_entries, num_entries_to_comm, displ)
513 8170 : CALL para_env%allgatherv(indx_local, indx_all, 3*num_entries_to_comm, 3*displ)
514 :
515 5384 : ALLOCATE (idx_homo(num_entries), idx_virt(num_entries))
516 1830 : IF (PRESENT(idx_spin)) ALLOCATE (idx_spin(num_entries))
517 4664 : DO kk = 1, num_entries
518 3030 : idx_homo(kk) = indx_all(3*kk - 2)
519 3030 : idx_virt(kk) = indx_all(3*kk - 1)
520 4664 : IF (PRESENT(idx_spin)) idx_spin(kk) = indx_all(3*kk)
521 : END DO
522 :
523 1634 : DEALLOCATE (entries_local, indx_local, indx_all, num_entries_to_comm, displ)
524 1634 : NULLIFY (row_indices)
525 1634 : 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 1634 : IF (PRESENT(idx_spin)) THEN
530 250 : CALL sort_excitations(idx_homo, idx_virt, eigvec_entries, idx_spin)
531 : ELSE
532 1384 : CALL sort_excitations(idx_homo, idx_virt, eigvec_entries)
533 : END IF
534 :
535 1634 : CALL timestop(handle)
536 :
537 1634 : 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 222 : 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 222 : n_ov_joint = 0
556 504 : DO isp = 1, SIZE(homo_red)
557 282 : offsets(isp) = n_ov_joint
558 282 : n_ov(isp) = homo_red(isp)*virt_red(isp)
559 504 : n_ov_joint = n_ov_joint + n_ov(isp)
560 : END DO
561 :
562 222 : 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 172 : 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 172 : nspins = SIZE(homo)
588 728 : ALLOCATE (spin_windows(nspins))
589 384 : DO ispin = 1, nspins
590 : CALL determine_mo_window(eigenval_scf(:, ispin), n_mo(ispin), homo(ispin), &
591 384 : cutoff_occ, cutoff_empty, spin_windows(ispin))
592 : END DO
593 172 : IF (nspins > 1) THEN
594 40 : CALL combine_mo_windows(spin_windows, window, spins_differ)
595 40 : IF (spins_differ .AND. warn) THEN
596 4 : CPWARN("BSE: spin-resolved active MO cutoff candidates differ; using combined window.")
597 : END IF
598 : ELSE
599 132 : window = spin_windows(1)
600 : END IF
601 172 : DEALLOCATE (spin_windows)
602 :
603 172 : 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 106 : 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 106 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenval_dft
623 106 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
624 :
625 106 : nspins = SIZE(mos)
626 424 : ALLOCATE (homo(nspins), n_mo(nspins))
627 236 : DO ispin = 1, nspins
628 130 : CALL get_mo_set(mos(ispin), nmo=n_mo(ispin), homo=homo(ispin), eigenvalues=eigenvalues)
629 448 : IF (ispin == 1) ALLOCATE (eigenval_dft(n_mo(1), nspins))
630 4326 : eigenval_dft(:, ispin) = eigenvalues(1:n_mo(ispin))
631 : END DO
632 106 : CALL bse_window_from_dft(eigenval_dft, n_mo, homo, cutoff_occ, cutoff_empty, window, warn)
633 106 : DEALLOCATE (homo, n_mo, eigenval_dft)
634 :
635 212 : 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 136 : SUBROUTINE truncate_BSE_matrices(fm_S_ia_full, fm_S_ij_full, fm_S_ab_full, &
657 68 : 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 68 : CALL timeset(routineN, handle)
680 :
681 68 : homo_incl = window%first_mo
682 68 : virt_incl = window%last_mo - homo
683 68 : homo_red = homo - homo_incl + 1
684 68 : virt_red = virt_incl
685 :
686 68 : IF (print_window .AND. unit_nr > 0) THEN
687 22 : IF (bse_env%bse_cutoff_occ > 0) THEN
688 17 : WRITE (unit_nr, '(T2,A4,T7,A29,T71,F10.3)') 'BSE|', 'Cutoff occupied orbitals [eV]', &
689 34 : bse_env%bse_cutoff_occ*evolt
690 : ELSE
691 5 : WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No cutoff given for occupied orbitals'
692 : END IF
693 22 : IF (bse_env%bse_cutoff_empty > 0) THEN
694 18 : WRITE (unit_nr, '(T2,A4,T7,A26,T71,F10.3)') 'BSE|', 'Cutoff empty orbitals [eV]', &
695 36 : bse_env%bse_cutoff_empty*evolt
696 : ELSE
697 4 : WRITE (unit_nr, '(T2,A4,T7,A34)') 'BSE|', 'No cutoff given for empty orbitals'
698 : END IF
699 22 : WRITE (unit_nr, '(T2,A4,T7,A20,T71,I10)') 'BSE|', 'First occupied index', homo_incl
700 22 : WRITE (unit_nr, '(T2,A4,T7,A32,T71,I10)') 'BSE|', 'Last empty index (not MO index!)', virt_incl
701 22 : WRITE (unit_nr, '(T2,A4,T7,A35,T71,F10.3)') 'BSE|', 'Energy of first occupied index [eV]', &
702 44 : Eigenval(homo_incl)*evolt
703 22 : WRITE (unit_nr, '(T2,A4,T7,A31,T71,F10.3)') 'BSE|', 'Energy of last empty index [eV]', &
704 44 : Eigenval(homo + virt_incl)*evolt
705 22 : WRITE (unit_nr, '(T2,A4,T7,A54,T71,F10.3)') 'BSE|', &
706 22 : 'Energy difference of first occupied index to HOMO [eV]', &
707 44 : -(Eigenval(homo_incl) - Eigenval(homo))*evolt
708 22 : WRITE (unit_nr, '(T2,A4,T7,A50,T71,F10.3)') 'BSE|', &
709 22 : 'Energy difference of last empty index to LUMO [eV]', &
710 44 : (Eigenval(homo + virt_incl) - Eigenval(homo + 1))*evolt
711 22 : WRITE (unit_nr, '(T2,A4,T7,A35,T71,I10)') 'BSE|', 'Number of GW-corrected occupied MOs', &
712 44 : bse_env%gw_corr_lev_occ
713 22 : WRITE (unit_nr, '(T2,A4,T7,A32,T71,I10)') 'BSE|', 'Number of GW-corrected empty MOs', &
714 44 : bse_env%gw_corr_lev_virt(ispin)
715 22 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
716 : END IF
717 68 : IF (unit_nr > 0) THEN
718 34 : IF (homo - homo_incl + 1 > bse_env%gw_corr_lev_occ) THEN
719 0 : CPABORT("Number of GW-corrected occupied MOs too small for chosen BSE cutoff")
720 : END IF
721 34 : IF (virt_incl > bse_env%gw_corr_lev_virt(ispin)) THEN
722 0 : 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 68 : para_env => fm_S_ia_full%matrix_struct%para_env
728 68 : context => fm_S_ia_full%matrix_struct%context
729 :
730 68 : CALL cp_fm_struct_create(fm_struct_ia, para_env, context, dimen_RI, homo_red*virt_red)
731 68 : CALL cp_fm_struct_create(fm_struct_ij, para_env, context, dimen_RI, homo_red*homo_red)
732 68 : CALL cp_fm_struct_create(fm_struct_ab, para_env, context, dimen_RI, virt_red*virt_red)
733 :
734 68 : CALL cp_fm_create(bse_env%fm_S_ia(ispin), fm_struct_ia, name="fm_S_trunc", set_zero=.TRUE.)
735 68 : CALL cp_fm_create(bse_env%fm_S_ij(ispin), fm_struct_ij, name="fm_S_ij_trunc", set_zero=.TRUE.)
736 68 : 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 68 : IF (bse_env%bse_cutoff_occ > 0 .OR. bse_env%bse_cutoff_empty > 0) THEN
740 : !Truncate eigenvals
741 208 : ALLOCATE (Eigenval_reduced_scf(homo_red + virt_red), Eigenval_reduced_qp(homo_red + virt_red))
742 1174 : Eigenval_reduced_scf(:) = Eigenval_scf(homo_incl:homo + virt_incl)
743 1174 : 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 468 : 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 468 : 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 260 : last_index=[dimen_RI, 1, virt_red, virt_red])
754 :
755 : ELSE
756 16 : IF (unit_nr > 0) THEN
757 8 : WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No truncation of BSE matrices applied'
758 8 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
759 : END IF
760 64 : ALLOCATE (Eigenval_reduced_scf(homo_red + virt_red), Eigenval_reduced_qp(homo_red + virt_red))
761 696 : Eigenval_reduced_scf(:) = Eigenval_scf(:)
762 696 : 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 16 : 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 16 : 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 16 : 1, 1, 1, 1, context)
769 : END IF
770 :
771 68 : bse_env%homo_red(ispin) = homo_red
772 68 : bse_env%virt_red(ispin) = virt_red
773 :
774 68 : CALL cp_fm_struct_release(fm_struct_ia)
775 68 : CALL cp_fm_struct_release(fm_struct_ij)
776 68 : CALL cp_fm_struct_release(fm_struct_ab)
777 :
778 68 : NULLIFY (para_env)
779 68 : NULLIFY (context)
780 :
781 68 : CALL timestop(handle)
782 :
783 68 : 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 900 : 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 300 : 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 300 : IF (do_transpose) THEN
815 50 : reordering = [2, 3, 1, 4]
816 50 : my_n_row = virtual
817 50 : my_m_col = homo
818 : ELSE
819 250 : reordering = [1, 3, 2, 4]
820 250 : my_n_row = homo
821 250 : 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 300 : 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 300 : my_n_row, my_m_col)
830 :
831 : ! Resort indices
832 300 : CALL cp_fm_create(fm_eigvec_col, fm_struct_eigvec_col, name="BSE_column_vector")
833 300 : CALL cp_fm_set_all(fm_eigvec_col, 0.0_dp)
834 300 : CALL cp_fm_create(fm_eigvec_reshuffled, fm_struct_eigvec_reshuffled, name="BSE_reshuffled_eigenvector")
835 300 : 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 300 : 1, 1)
841 : ! Reshuffle
842 300 : CALL cp_fm_permute(fm_eigvec_reshuffled, fm_eigvec_col, 1.0_dp, virtual, 1, 1, 1, reordering)
843 :
844 300 : CALL cp_fm_release(fm_eigvec_col)
845 300 : CALL cp_fm_struct_release(fm_struct_eigvec_col)
846 300 : CALL cp_fm_struct_release(fm_struct_eigvec_reshuffled)
847 :
848 300 : CALL timestop(handle)
849 :
850 300 : 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 0 : 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 0 : 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 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
886 : TYPE(pw_c1d_gs_type) :: wf_g
887 : TYPE(pw_env_type), POINTER :: pw_env
888 0 : 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 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
892 : TYPE(qs_subsys_type), POINTER :: subsys
893 :
894 0 : logger => cp_get_default_logger()
895 0 : CALL timeset(routineN, handle)
896 :
897 0 : dft_control => bse_env%dft_control
898 : ! PROPERTIES%BSE runs after the GW released the grids; its requirement reader refuses cube files
899 0 : IF (.NOT. ASSOCIATED(bse_env%pw_env)) THEN
900 0 : CPABORT("NTO cube files need the SCF grids: set NTO_ANALYSIS%CUBE_FILES to F")
901 : END IF
902 0 : pw_env => bse_env%pw_env
903 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, pw_pools=pw_pools)
904 0 : CALL auxbas_pw_pool%create_pw(wf_r)
905 0 : CALL auxbas_pw_pool%create_pw(wf_g)
906 :
907 0 : subsys => bse_env%subsys
908 0 : CALL qs_subsys_get(subsys, particles=particles)
909 :
910 0 : my_pos_cube = "REWIND"
911 0 : IF (append_cube) THEN
912 0 : 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 0 : cell=cell, particle_set=particle_set)
917 :
918 0 : DO iset = 1, 2
919 0 : CALL get_mo_set(mo_set=mos(iset), mo_coeff=mo_coeff, nmo=nmo)
920 0 : DO i = 1, nmo
921 : CALL calculate_wavefunction(mo_coeff, i, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
922 0 : cell, dft_control, particle_set, pw_env)
923 0 : IF (iset == 1) THEN
924 0 : WRITE (filename, '(A6,I3.3,A5,I2.2,a11)') "_NEXC_", istate, "_NTO_", i, "_Hole_State"
925 0 : ELSE IF (iset == 2) THEN
926 0 : WRITE (filename, '(A6,I3.3,A5,I2.2,a15)') "_NEXC_", istate, "_NTO_", i, "_Particle_State"
927 : END IF
928 0 : info_approx_trunc = TRIM(ADJUSTL(info_approximation))
929 0 : info_approx_trunc = info_approx_trunc(2:LEN_TRIM(info_approx_trunc) - 1)
930 0 : filename = TRIM(info_approx_trunc)//TRIM(filename)
931 0 : 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 0 : log_filename=.FALSE., ignore_should_output=.TRUE., mpi_io=mpi_io)
935 0 : IF (iset == 1) THEN
936 0 : WRITE (title, *) "Natural Transition Orbital Hole State", i
937 0 : ELSE IF (iset == 2) THEN
938 0 : WRITE (title, *) "Natural Transition Orbital Particle State", i
939 : END IF
940 0 : 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 0 : ignore_should_output=.TRUE., mpi_io=mpi_io)
943 : END DO
944 : END DO
945 :
946 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
947 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
948 :
949 0 : CALL timestop(handle)
950 0 : 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 124 : 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 62 : CALL timeset(routineN, handle)
974 : ! the Poisson periodicity was read from the input (finish_bse_env): the grids may be released
975 62 : NULLIFY (cell_ref)
976 62 : CALL qs_subsys_get(bse_env%subsys, cell_ref=cell_ref)
977 62 : ndim_periodic_poisson = bse_env%n_periodic_poisson
978 248 : 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 62 : n_exc_max = homo*virtual
982 62 : IF (bse_env%bse_diag_method == bse_iterdiag) THEN
983 8 : n_exc_max = MIN(n_exc_max, bse_env%num_exc_en)
984 : END IF
985 :
986 : ! Handle negative NUM_PRINT_EXC
987 62 : IF (bse_env%num_print_exc < 0 .OR. &
988 : bse_env%num_print_exc > n_exc_max) THEN
989 14 : bse_env%num_print_exc = n_exc_max
990 14 : IF (unit_nr > 0) THEN
991 : CALL cp_hint(__LOCATION__, &
992 : "Keyword NUM_PRINT_EXC is either negative or too large. "// &
993 7 : "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 62 : IF (bse_env%do_nto_analysis) THEN
1004 4 : IF (bse_env%explicit_nto_list) THEN
1005 0 : IF (bse_env%num_print_exc_ntos > 0) THEN
1006 0 : 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 0 : "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 0 : num_state_list_exceptions = 0
1015 0 : DO i = 1, SIZE(bse_env%bse_nto_state_list)
1016 0 : IF (bse_env%bse_nto_state_list(i) < 1 .OR. &
1017 0 : bse_env%bse_nto_state_list(i) > bse_env%num_print_exc) THEN
1018 0 : num_state_list_exceptions = num_state_list_exceptions + 1
1019 : END IF
1020 : END DO
1021 0 : IF (num_state_list_exceptions > 0) THEN
1022 0 : IF (unit_nr > 0) THEN
1023 : CALL cp_hint(__LOCATION__, &
1024 : "STATE_LIST contains indices outside the range of included excitation levels. "// &
1025 0 : "Ignoring these states.")
1026 : END IF
1027 : END IF
1028 0 : n = SIZE(bse_env%bse_nto_state_list) - num_state_list_exceptions
1029 0 : ALLOCATE (bse_env%bse_nto_state_list_final(n))
1030 0 : bse_env%bse_nto_state_list_final(:) = 0
1031 : i = 1
1032 0 : DO j = 1, SIZE(bse_env%bse_nto_state_list)
1033 0 : IF (bse_env%bse_nto_state_list(j) >= 1 .AND. &
1034 0 : bse_env%bse_nto_state_list(j) <= bse_env%num_print_exc) THEN
1035 0 : bse_env%bse_nto_state_list_final(i) = bse_env%bse_nto_state_list(j)
1036 0 : i = i + 1
1037 : END IF
1038 : END DO
1039 :
1040 0 : bse_env%num_print_exc_ntos = SIZE(bse_env%bse_nto_state_list_final)
1041 : ELSE
1042 4 : IF (bse_env%num_print_exc_ntos > bse_env%num_print_exc .OR. &
1043 : bse_env%num_print_exc_ntos < 0) THEN
1044 4 : bse_env%num_print_exc_ntos = bse_env%num_print_exc
1045 : END IF
1046 12 : ALLOCATE (bse_env%bse_nto_state_list_final(bse_env%num_print_exc_ntos))
1047 104 : DO i = 1, bse_env%num_print_exc_ntos
1048 104 : 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 62 : IF (bse_env%bse_spin_config /= 0 .AND. &
1055 : bse_env%eps_nto_osc_str > 0) THEN
1056 0 : IF (unit_nr > 0) THEN
1057 : CALL cp_warn(__LOCATION__, &
1058 : "Cannot apply EPS_OSC_STR for Triplet excitations. "// &
1059 0 : "Resetting EPS_OSC_STR to default.")
1060 : END IF
1061 0 : bse_env%eps_nto_osc_str = -1.0_dp
1062 : END IF
1063 :
1064 : ! Take care of number for computed exciton descriptors
1065 62 : IF (bse_env%num_print_exc_descr < 0 .OR. &
1066 : bse_env%num_print_exc_descr > bse_env%num_print_exc) THEN
1067 4 : IF (unit_nr > 0) THEN
1068 : CALL cp_hint(__LOCATION__, &
1069 : "Keyword NUM_PRINT_EXC_DESCR is either negative or too large. "// &
1070 2 : "Printing exciton descriptors up to NUM_PRINT_EXC.")
1071 : END IF
1072 4 : 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 62 : CALL section_vals_val_get(bse_env%bse_section, "PRINT_DIRECTIONAL_CROSSCORRELATION", explicit=explicit)
1077 62 : IF (explicit .AND. bse_env%print_directional_crosscorrelation .AND. &
1078 : .NOT. bse_env%print_directional_exc_descr) THEN
1079 0 : IF (unit_nr > 0) THEN
1080 : CALL cp_warn(__LOCATION__, &
1081 0 : "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 62 : IF (bse_env%use_ks_energies) THEN
1087 4 : IF (unit_nr > 0) THEN
1088 : CALL cp_warn(__LOCATION__, &
1089 : "Using KS energies for BSE calculations. Therefore, no quantities "// &
1090 2 : "of the preceeding GW calculation enter the BSE.")
1091 : END IF
1092 : END IF
1093 :
1094 : ! Add warning if periodic calculation is invoked
1095 62 : IF (ndim_periodic_poisson /= 0) THEN
1096 0 : 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 0 : "oscillator strengths.")
1101 : END IF
1102 : END IF
1103 62 : IF (ndim_periodic_cell /= 0) THEN
1104 0 : 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 0 : "oscillator strengths.")
1109 : END IF
1110 : END IF
1111 :
1112 62 : CALL timestop(handle)
1113 62 : 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 62 : 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 62 : IF (bse_env%screening_factor > 0.0_dp) THEN
1126 2 : IF (bse_env%screening_method /= bse_screening_alpha) THEN
1127 0 : IF (unit_nr > 0) THEN
1128 : CALL cp_warn(__LOCATION__, &
1129 : "Screening factor is only supported for &SCREENING_IN_W ALPHA. "// &
1130 0 : "Resetting SCREENING_IN_W to ALPHA.")
1131 : END IF
1132 0 : bse_env%screening_method = bse_screening_alpha
1133 : END IF
1134 2 : IF (bse_env%screening_factor > 1.0_dp) THEN
1135 0 : IF (unit_nr > 0) THEN
1136 : CALL cp_warn(__LOCATION__, &
1137 0 : "Screening factor is larger than 1.0. ")
1138 : END IF
1139 : END IF
1140 : END IF
1141 :
1142 62 : IF (bse_env%screening_factor < 0.0_dp .AND. &
1143 : bse_env%screening_method == bse_screening_alpha) THEN
1144 0 : IF (unit_nr > 0) THEN
1145 : CALL cp_warn(__LOCATION__, &
1146 0 : "Screening factor is negative. Defaulting to 0.25")
1147 : END IF
1148 0 : bse_env%screening_factor = 0.25_dp
1149 : END IF
1150 :
1151 62 : IF (bse_env%screening_factor == 0.0_dp) THEN
1152 : ! Use RPA internally in this case
1153 0 : bse_env%screening_method = bse_screening_rpa
1154 : END IF
1155 62 : IF (bse_env%screening_factor == 1.0_dp) THEN
1156 : ! Use TDHF internally in this case
1157 0 : bse_env%screening_method = bse_screening_tdhf
1158 : END IF
1159 :
1160 62 : 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 80 : SUBROUTINE get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, &
1179 80 : 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 80 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_multipole_per_dir
1201 : TYPE(mp_para_env_type), POINTER :: para_env_BSE
1202 :
1203 80 : CALL timeset(routineN, handle)
1204 :
1205 80 : my_ispin = 1
1206 80 : 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 80 : fm_struct_multipoles_ao => mos(my_ispin)%mo_coeff%matrix_struct
1210 : ! BSE has different contexts and blacsenvs
1211 80 : para_env_BSE => context_BSE%para_env
1212 80 : n_multipole = SIZE(matrix_multipole)
1213 516 : 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 80 : CALL get_mo_set(mo_set=mos(my_ispin), homo=n_occ, nao=nao)
1220 : ! Takes into account removed nullspace values from SVD
1221 80 : nmo_mp2 = mo_coeff(1)%matrix_struct%ncol_global
1222 80 : 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 80 : 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 80 : 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 80 : context_BSE, virtual_red, homo_red)
1244 : CALL cp_fm_struct_create(fm_struct_mp_ij_trunc, para_env_BSE, &
1245 80 : context_BSE, homo_red, homo_red)
1246 : CALL cp_fm_struct_create(fm_struct_mp_ab_trunc, para_env_BSE, &
1247 80 : context_BSE, virtual_red, virtual_red)
1248 356 : DO idir = 1, n_multipole
1249 : CALL cp_fm_create(fm_multipole_ai_trunc(idir), matrix_struct=fm_struct_mp_ai_trunc, &
1250 276 : name="dipoles_mo_ai_trunc")
1251 276 : 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 276 : name="dipoles_mo_ij_trunc")
1254 276 : 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 276 : name="dipoles_mo_ab_trunc")
1257 356 : 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 80 : CALL cp_fm_create(fm_work, matrix_struct=fm_struct_nao_nmo, name="multipole_work")
1263 80 : CALL cp_fm_set_all(fm_work, 0.0_dp)
1264 :
1265 356 : DO idir = 1, n_multipole
1266 : ! Create the full multipole matrix per direction
1267 276 : CALL cp_fm_create(fm_multipole_per_dir(idir), matrix_struct=fm_struct_nmo_nmo, name="multipoles_mo")
1268 276 : 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 276 : 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 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))
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 276 : 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 276 : 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 356 : fm_multipole_per_dir(idir)%matrix_struct%context)
1307 : END DO
1308 :
1309 : !Release matrices and structs
1310 80 : NULLIFY (fm_struct_multipoles_ao)
1311 80 : CALL cp_fm_struct_release(fm_struct_mp_ai_trunc)
1312 80 : CALL cp_fm_struct_release(fm_struct_mp_ij_trunc)
1313 80 : CALL cp_fm_struct_release(fm_struct_mp_ab_trunc)
1314 80 : CALL cp_fm_struct_release(fm_struct_nao_nmo)
1315 80 : CALL cp_fm_struct_release(fm_struct_nmo_nmo)
1316 356 : DO idir = 1, n_multipole
1317 356 : CALL cp_fm_release(fm_multipole_per_dir(idir))
1318 : END DO
1319 80 : DEALLOCATE (fm_multipole_per_dir)
1320 80 : CALL cp_fm_release(fm_work)
1321 :
1322 80 : CALL timestop(handle)
1323 :
1324 240 : 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 11520 : 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 1920 : CALL timeset(routineN, handle)
1345 :
1346 1920 : CALL cp_fm_create(fm_work_ia, fm_A%matrix_struct)
1347 1920 : CALL cp_fm_get_info(fm_A, nrow_global=nrow_A, ncol_global=ncol_A)
1348 1920 : CALL cp_fm_get_info(fm_B, nrow_global=nrow_B, ncol_global=ncol_B)
1349 1920 : CALL cp_fm_get_info(fm_C, nrow_global=nrow_C, ncol_global=ncol_C)
1350 :
1351 : ! Check matrix sizes
1352 1920 : CPASSERT(nrow_A == nrow_B .AND. ncol_A == ncol_C .AND. ncol_B == nrow_C)
1353 :
1354 1920 : 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 1920 : fm_B, fm_C, 0.0_dp, fm_work_ia)
1358 :
1359 1920 : CALL cp_fm_trace(fm_A, fm_work_ia, alpha)
1360 :
1361 1920 : CALL cp_fm_release(fm_work_ia)
1362 :
1363 1920 : CALL timestop(handle)
1364 :
1365 1920 : 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 18 : 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 18 : CALL timeset(routineN, handle)
1388 54 : DO isp = 1, SIZE(fm_S_ia)
1389 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)
1390 : END DO
1391 18 : CALL timestop(handle)
1392 :
1393 18 : END SUBROUTINE assemble_joint_ov_slab
1394 :
1395 : END MODULE bse_util
|
