(git:d5a0442)
Loading...
Searching...
No Matches
bse_full_diag.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines for the full diagonalization of GW + Bethe-Salpeter for computing
10!> electronic excitations
11!> \par History
12!> 10.2023 created [Maximilian Graml]
13! **************************************************************************************************
14MODULE bse_full_diag
15
16 USE bse_print, ONLY: print_excitation_energies,&
17 print_exciton_descriptors,&
18 print_optical_properties,&
19 print_output_header,&
20 print_transition_amplitudes
21 USE bse_properties, ONLY: calculate_ntos,&
22 exciton_descr_type,&
23 get_exciton_descriptors,&
24 get_oscillator_strengths
25 USE bse_types, ONLY: bse_env_type
26 USE bse_util, ONLY: assemble_joint_ov_slab,&
27 comp_eigvec_coeff_bse,&
28 fm_general_add_bse,&
29 get_bse_spin_block_layout,&
30 get_multipoles_mo,&
31 occ_of_ia,&
32 reshuffle_eigvec,&
33 virt_of_ia
34 USE cp_blacs_env, ONLY: cp_blacs_env_create,&
35 cp_blacs_env_release,&
36 cp_blacs_env_type
37 USE cp_control_types, ONLY: dft_control_type,&
38 tddfpt2_control_type
39 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
40 USE cp_fm_diag, ONLY: choose_eigv_solver,&
41 cp_fm_power
42 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
43 cp_fm_struct_release,&
44 cp_fm_struct_type
45 USE cp_fm_types, ONLY: cp_fm_create,&
46 cp_fm_get_info,&
47 cp_fm_get_submatrix,&
48 cp_fm_release,&
49 cp_fm_set_all,&
50 cp_fm_to_fm,&
51 cp_fm_to_fm_submat,&
52 cp_fm_type
53 USE exstates_types, ONLY: excited_energy_type
54 USE input_constants, ONLY: bse_iterdiag,&
55 bse_screening_alpha,&
56 bse_screening_rpa,&
57 bse_screening_tdhf,&
58 bse_singlet,&
59 bse_triplet
60 USE kinds, ONLY: dp
61 USE parallel_gemm_api, ONLY: parallel_gemm
62 USE particle_types, ONLY: particle_type
63 USE physcon, ONLY: evolt
64 USE qs_subsys_types, ONLY: qs_subsys_get
65#include "./base/base_uses.f90"
66
67 IMPLICIT NONE
68
69 PRIVATE
70
71 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_full_diag'
72
73 PUBLIC :: create_a, diagonalize_a, create_b, create_a_and_b, create_hermitian_form_of_abba, &
74 diagonalize_c, postprocess_bse
75
76CONTAINS
77
78! **************************************************************************************************
79!> \brief Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine screen_slabs)
80!> A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
81!> ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
82!> α is a spin-dependent factor
83!> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
84!> W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
85!> \param fm_S_ia ...
86!> \param fm_S_bar_ij ...
87!> \param fm_S_ab ...
88!> \param fm_A ...
89!> \param Eigenval ...
90!> \param unit_nr ...
91!> \param homo ...
92!> \param virtual ...
93!> \param dimen_RI ...
94!> \param bse_env the BSE environment (settings; para_env, dft_control, exstate_env)
95! **************************************************************************************************
96 SUBROUTINE create_a(fm_S_ia, fm_S_bar_ij, fm_S_ab, &
97 fm_A, Eigenval, unit_nr, &
98 homo, virtual, dimen_RI, bse_env)
99
100 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_s_ia, fm_s_bar_ij, fm_s_ab
101 TYPE(cp_fm_type), INTENT(INOUT) :: fm_a
102 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
103 INTEGER, INTENT(IN) :: unit_nr
104 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
105 INTEGER, INTENT(IN) :: dimen_ri
106 TYPE(bse_env_type), INTENT(IN) :: bse_env
107
108 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_A'
109
110 INTEGER :: a_virt_row, handle, i_occ_row, i_row_global, ii, isp, j_col_global, jj, k_isp, &
111 k_ov, n_ov_joint, ncol_local_a, nrow_local_a, nspins, sizeeigen
112 INTEGER, ALLOCATABLE, DIMENSION(:) :: eig_offsets, n_ov, offsets
113 INTEGER, DIMENSION(4) :: reordering
114 INTEGER, DIMENSION(:), POINTER :: col_indices_a, row_indices_a
115 REAL(kind=dp) :: alpha, alpha_screening, eigen_diff
116 TYPE(cp_blacs_env_type), POINTER :: blacs_env
117 TYPE(cp_fm_struct_type), POINTER :: fm_struct_a, fm_struct_s_joint, &
118 fm_struct_w
119 TYPE(cp_fm_type) :: fm_a_copy, fm_s_joint, fm_w
120 TYPE(dft_control_type), POINTER :: dft_control
121 TYPE(excited_energy_type), POINTER :: ex_env
122 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
123
124 CALL timeset(routinen, handle)
125
126 nspins = SIZE(homo)
127 ALLOCATE (n_ov(nspins), offsets(nspins), eig_offsets(nspins))
128 CALL get_bse_spin_block_layout(homo, virtual, n_ov, offsets, n_ov_joint)
129 ! Flat Eigenval layout: sigma-block isp at eig_offsets(isp)+1 .. eig_offsets(isp)+homo(isp)+virtual(isp)
130 eig_offsets(1) = 0
131 DO isp = 2, nspins
132 eig_offsets(isp) = eig_offsets(isp - 1) + homo(isp - 1) + virtual(isp - 1)
133 END DO
134
135 NULLIFY (dft_control, tddfpt_control)
136 dft_control => bse_env%dft_control
137 tddfpt_control => dft_control%tddfpt2_control
138
139 IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
140 WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating A'
141 END IF
142
143 !Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
144 SELECT CASE (bse_env%bse_spin_config)
145 CASE (bse_singlet)
146 alpha = 2.0_dp
147 CASE (bse_triplet)
148 alpha = 0.0_dp
149 END SELECT
150 ! For open-shell (nspins>1): each spin block contributes once; SPIN_CONFIG is ignored.
151 IF (nspins > 1) THEN
152 CALL cp_warn(__location__, &
153 "BSE: SPIN_CONFIG ignored for open-shell reference; using alpha=1.")
154 alpha = 1.0_dp
155 END IF
156
157 IF (bse_env%screening_method == bse_screening_alpha) THEN
158 alpha_screening = bse_env%screening_factor
159 ELSE
160 alpha_screening = 1.0_dp
161 END IF
162
163 ! create the blacs env for ij matrices (NOT fm_S_ia%matrix_struct related parallel_gemms!)
164 NULLIFY (blacs_env)
165 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=bse_env%para_env)
166
167 ! We have to use the same blacs_env for A as for the matrices fm_S_ia from RPA
168 ! Logic: A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
169 ! We create v_ia,jb and W_ij,ab, then we communicate entries from local W_ij,ab
170 ! to the full matrix v_ia,jb. By adding these and the energy diffenences: v_ia,jb -> A_ia,jb
171 ! We use the A matrix already from the start instead of v
172 CALL cp_fm_struct_create(fm_struct_a, context=fm_s_ia(1)%matrix_struct%context, &
173 nrow_global=n_ov_joint, ncol_global=n_ov_joint, &
174 para_env=fm_s_ia(1)%matrix_struct%para_env)
175 CALL cp_fm_create(fm_a, fm_struct_a, name="fm_A_iajb")
176 CALL cp_fm_set_all(fm_a, 0.0_dp)
177 ! fm_A_copy only used in the TDDFPT do_bse_w_only path (closed-shell only)
178 IF (tddfpt_control%do_bse_w_only .AND. nspins == 1) THEN
179 CALL cp_fm_create(fm_a_copy, fm_struct_a, name="fm_A_iajb")
180 CALL cp_fm_set_all(fm_a_copy, 0.0_dp)
181 END IF
182
183 ! Create A matrix from GW Energies, v_ia,jb and W_ij,ab
184 ! v_ia,jb = \sum_P B^P_ia B^P_jb (Coulomb)
185 IF ((.NOT. tddfpt_control%do_bse) .AND. (.NOT. tddfpt_control%do_bse_w_only)) THEN
186 IF (nspins > 1) THEN
187 ! Assemble joint ia-slab for a single Coulomb gemm across all spin blocks
188 CALL cp_fm_struct_create(fm_struct_s_joint, &
189 context=fm_s_ia(1)%matrix_struct%context, &
190 nrow_global=dimen_ri, ncol_global=n_ov_joint, &
191 para_env=fm_s_ia(1)%matrix_struct%para_env)
192 CALL cp_fm_create(fm_s_joint, fm_struct_s_joint, name="fm_S_ia_joint")
193 CALL cp_fm_set_all(fm_s_joint, 0.0_dp)
194 CALL assemble_joint_ov_slab(fm_s_ia, offsets, n_ov, dimen_ri, fm_s_joint)
195 CALL parallel_gemm(transa="T", transb="N", m=n_ov_joint, n=n_ov_joint, k=dimen_ri, &
196 alpha=alpha, matrix_a=fm_s_joint, matrix_b=fm_s_joint, beta=0.0_dp, &
197 matrix_c=fm_a)
198 CALL cp_fm_release(fm_s_joint)
199 CALL cp_fm_struct_release(fm_struct_s_joint)
200 ELSE
201 CALL parallel_gemm(transa="T", transb="N", m=homo(1)*virtual(1), n=homo(1)*virtual(1), &
202 k=dimen_ri, alpha=alpha, &
203 matrix_a=fm_s_ia(1), matrix_b=fm_s_ia(1), &
204 beta=0.0_dp, matrix_c=fm_a)
205 END IF
206 END IF
207
208 IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
209 WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated A_iajb'
210 END IF
211
212 ! W term on sigma-diagonal blocks only: W^sigma_ij,ab = sum_P barB^P_ij B^P_ab
213 ! offsets(isp) places each block at the correct position in joint A.
214 ! For nspins=1: offsets(1)=0, equivalent to the original code.
215 DO isp = 1, nspins
216 IF (bse_env%screening_method /= bse_screening_rpa) THEN
217 CALL cp_fm_struct_create(fm_struct_w, context=fm_s_ab(isp)%matrix_struct%context, &
218 nrow_global=homo(isp)**2, ncol_global=virtual(isp)**2, &
219 para_env=fm_s_ab(isp)%matrix_struct%para_env)
220 CALL cp_fm_create(fm_w, fm_struct_w, name="fm_W_ijab")
221 CALL cp_fm_set_all(fm_w, 0.0_dp)
222 !W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab
223 CALL parallel_gemm(transa="T", transb="N", m=homo(isp)**2, n=virtual(isp)**2, &
224 k=dimen_ri, alpha=alpha_screening, &
225 matrix_a=fm_s_bar_ij(isp), matrix_b=fm_s_ab(isp), &
226 beta=0.0_dp, matrix_c=fm_w)
227 reordering = [1, 3, 2, 4]
228 CALL fm_general_add_bse(fm_a, fm_w, -1.0_dp, homo(isp), virtual(isp), &
229 virtual(isp), virtual(isp), unit_nr, reordering, bse_env, &
230 row_offset=offsets(isp), col_offset=offsets(isp))
231 IF (nspins == 1 .AND. tddfpt_control%do_bse_w_only) THEN
232 CALL fm_general_add_bse(fm_a_copy, fm_w, -1.0_dp, homo(1), virtual(1), &
233 virtual(1), virtual(1), unit_nr, reordering, bse_env)
234 END IF
235 ! W and A stash for TDDFPT path (closed-shell only; open-shell deferred)
236 IF (nspins == 1) THEN
237 IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
238 tddfpt_control%do_bse_gw_only) THEN
239 NULLIFY (ex_env)
240 ex_env => bse_env%exstate_env
241 IF (.NOT. tddfpt_control%do_bse_gw_only) THEN
242 ALLOCATE (ex_env%bse_w_matrix_MO(1, 1))
243 ALLOCATE (ex_env%bse_a_matrix_MO(1, 1))
244 CALL cp_fm_create(ex_env%bse_w_matrix_MO(1, 1), fm_struct_w)
245 CALL cp_fm_create(ex_env%bse_a_matrix_MO(1, 1), fm_struct_a)
246 CALL cp_fm_to_fm(fm_w, ex_env%bse_w_matrix_MO(1, 1))
247 IF (tddfpt_control%do_bse_w_only) THEN
248 CALL cp_fm_to_fm(fm_a_copy, ex_env%bse_a_matrix_MO(1, 1))
249 ELSE
250 CALL cp_fm_to_fm(fm_a, ex_env%bse_a_matrix_MO(1, 1))
251 END IF
252 END IF
253 END IF
254 END IF
255 CALL cp_fm_release(fm_w)
256 CALL cp_fm_struct_release(fm_struct_w)
257 END IF
258 END DO
259
260 IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
261 WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
262 END IF
263 IF (nspins == 1 .AND. tddfpt_control%do_bse_w_only) CALL cp_fm_release(fm_a_copy)
264
265 ! Get local row/col indices for direct diagonal access
266 CALL cp_fm_get_info(matrix=fm_a, nrow_local=nrow_local_a, ncol_local=ncol_local_a, &
267 row_indices=row_indices_a, col_indices=col_indices_a)
268
269 !Add (ε_a-ε_i) on the diagonal of each sigma-block; cross-spin blocks have no ε contribution.
270 IF (.NOT. tddfpt_control%do_bse) THEN
271 DO ii = 1, nrow_local_a
272 i_row_global = row_indices_a(ii)
273 DO jj = 1, ncol_local_a
274 j_col_global = col_indices_a(jj)
275 IF (i_row_global == j_col_global) THEN
276 ! Decode spin: isp such that i_row_global in [offsets(isp)+1, offsets(isp)+n_ov(isp)]
277 isp = nspins
278 DO k_isp = 1, nspins - 1
279 IF (i_row_global <= offsets(k_isp) + n_ov(k_isp)) THEN
280 isp = k_isp
281 EXIT
282 END IF
283 END DO
284 k_ov = i_row_global - offsets(isp)
285 i_occ_row = occ_of_ia(k_ov, virtual(isp))
286 a_virt_row = virt_of_ia(k_ov, virtual(isp))
287 eigen_diff = eigenval(eig_offsets(isp) + a_virt_row + homo(isp)) - &
288 eigenval(eig_offsets(isp) + i_occ_row)
289 fm_a%local_data(ii, jj) = fm_a%local_data(ii, jj) + eigen_diff
290 END IF
291 END DO
292 END DO
293 END IF
294
295 ! GW eigenvalue stash for TDDFPT path (closed-shell only)
296 IF (nspins == 1) THEN
297 IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
298 tddfpt_control%do_bse_gw_only) THEN
299 sizeeigen = SIZE(eigenval)
300 ALLOCATE (ex_env%gw_eigen(sizeeigen))
301 ex_env%gw_eigen(:) = eigenval(:)
302 END IF
303 END IF
304
305 CALL cp_fm_struct_release(fm_struct_a)
306 DEALLOCATE (n_ov, offsets, eig_offsets)
307
308 CALL cp_blacs_env_release(blacs_env)
309
310 CALL timestop(handle)
311
312 END SUBROUTINE create_a
313
314! **************************************************************************************************
315!> \brief Matrix B constructed from 3c-B-matrices (cf. subroutine screen_slabs)
316!> B_ia,jb = α * v_ia,jb - W_ib,aj
317!> α is a spin-dependent factor
318!> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
319!> W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
320!> \param fm_S_ia ...
321!> \param fm_S_bar_ia ...
322!> \param fm_B ...
323!> \param homo ...
324!> \param virtual ...
325!> \param dimen_RI ...
326!> \param unit_nr ...
327!> \param bse_env the BSE environment (settings)
328! **************************************************************************************************
329 SUBROUTINE create_b(fm_S_ia, fm_S_bar_ia, fm_B, &
330 homo, virtual, dimen_RI, unit_nr, bse_env)
331
332 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_s_ia, fm_s_bar_ia
333 TYPE(cp_fm_type), INTENT(INOUT) :: fm_b
334 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
335 INTEGER, INTENT(IN) :: dimen_ri, unit_nr
336 TYPE(bse_env_type), INTENT(IN) :: bse_env
337
338 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_B'
339
340 INTEGER :: handle, isp, n_ov_joint, nspins
341 INTEGER, ALLOCATABLE, DIMENSION(:) :: n_ov, offsets
342 INTEGER, DIMENSION(4) :: reordering
343 REAL(kind=dp) :: alpha, alpha_screening
344 TYPE(cp_fm_struct_type), POINTER :: fm_struct_b, fm_struct_s_joint, &
345 fm_struct_w
346 TYPE(cp_fm_type) :: fm_s_joint, fm_w
347
348 CALL timeset(routinen, handle)
349
350 nspins = SIZE(homo)
351 ALLOCATE (n_ov(nspins), offsets(nspins))
352 CALL get_bse_spin_block_layout(homo, virtual, n_ov, offsets, n_ov_joint)
353
354 IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
355 WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating B'
356 END IF
357
358 ! Coulomb prefactor: SPIN_CONFIG sector for closed shell; for open shell each spin block
359 ! contributes once (alpha=1). create_A already emits the SPIN_CONFIG-ignored warning.
360 SELECT CASE (bse_env%bse_spin_config)
361 CASE (bse_singlet)
362 alpha = 2.0_dp
363 CASE (bse_triplet)
364 alpha = 0.0_dp
365 END SELECT
366 IF (nspins > 1) alpha = 1.0_dp
367
368 IF (bse_env%screening_method == bse_screening_alpha) THEN
369 alpha_screening = bse_env%screening_factor
370 ELSE
371 alpha_screening = 1.0_dp
372 END IF
373
374 ! Joint B over all spin blocks: B_ia,jb = alpha*(ia|bj) - W^sigma_ib,aj (W spin-diagonal)
375 NULLIFY (fm_struct_b)
376 CALL cp_fm_struct_create(fm_struct_b, context=fm_s_ia(1)%matrix_struct%context, &
377 nrow_global=n_ov_joint, ncol_global=n_ov_joint, &
378 para_env=fm_s_ia(1)%matrix_struct%para_env)
379 CALL cp_fm_create(fm_b, fm_struct_b, name="fm_B_iajb")
380 CALL cp_fm_set_all(fm_b, 0.0_dp)
381
382 ! Coulomb v_ia,jb = sum_P B^P_ia B^P_jb (= (ia|bj)); cross-spin blocks filled automatically.
383 IF (nspins > 1) THEN
384 NULLIFY (fm_struct_s_joint)
385 CALL cp_fm_struct_create(fm_struct_s_joint, &
386 context=fm_s_ia(1)%matrix_struct%context, &
387 nrow_global=dimen_ri, ncol_global=n_ov_joint, &
388 para_env=fm_s_ia(1)%matrix_struct%para_env)
389 CALL cp_fm_create(fm_s_joint, fm_struct_s_joint, name="fm_S_ia_joint")
390 CALL cp_fm_set_all(fm_s_joint, 0.0_dp)
391 CALL assemble_joint_ov_slab(fm_s_ia, offsets, n_ov, dimen_ri, fm_s_joint)
392 CALL parallel_gemm(transa="T", transb="N", m=n_ov_joint, n=n_ov_joint, k=dimen_ri, &
393 alpha=alpha, matrix_a=fm_s_joint, matrix_b=fm_s_joint, beta=0.0_dp, &
394 matrix_c=fm_b)
395 CALL cp_fm_release(fm_s_joint)
396 CALL cp_fm_struct_release(fm_struct_s_joint)
397 ELSE
398 CALL parallel_gemm(transa="T", transb="N", m=homo(1)*virtual(1), n=homo(1)*virtual(1), &
399 k=dimen_ri, alpha=alpha, &
400 matrix_a=fm_s_ia(1), matrix_b=fm_s_ia(1), &
401 beta=0.0_dp, matrix_c=fm_b)
402 END IF
403
404 IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
405 WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated B_iajb'
406 END IF
407
408 ! W^sigma_ib,aj = sum_P barB^P_ib B^P_aj on sigma-diagonal blocks only (offsets place them).
409 ! reordering [1,4,3,2] maps W_ib,ja -> B_ia,jb. For nspins=1: offsets(1)=0 (original code).
410 IF (bse_env%screening_method /= bse_screening_rpa) THEN
411 DO isp = 1, nspins
412 NULLIFY (fm_struct_w)
413 CALL cp_fm_struct_create(fm_struct_w, &
414 context=fm_s_ia(isp)%matrix_struct%context, &
415 nrow_global=homo(isp)*virtual(isp), &
416 ncol_global=homo(isp)*virtual(isp), &
417 para_env=fm_s_ia(isp)%matrix_struct%para_env)
418 CALL cp_fm_create(fm_w, fm_struct_w, name="fm_W_ibaj")
419 CALL cp_fm_set_all(fm_w, 0.0_dp)
420 CALL parallel_gemm(transa="T", transb="N", m=homo(isp)*virtual(isp), &
421 n=homo(isp)*virtual(isp), k=dimen_ri, alpha=alpha_screening, &
422 matrix_a=fm_s_bar_ia(isp), matrix_b=fm_s_ia(isp), &
423 beta=0.0_dp, matrix_c=fm_w)
424 reordering = [1, 4, 3, 2]
425 CALL fm_general_add_bse(fm_b, fm_w, -1.0_dp, virtual(isp), virtual(isp), &
426 virtual(isp), virtual(isp), unit_nr, reordering, bse_env, &
427 row_offset=offsets(isp), col_offset=offsets(isp))
428 CALL cp_fm_release(fm_w)
429 CALL cp_fm_struct_release(fm_struct_w)
430 END DO
431 END IF
432
433 CALL cp_fm_struct_release(fm_struct_b)
434 DEALLOCATE (n_ov, offsets)
435
436 CALL timestop(handle)
437
438 END SUBROUTINE create_b
439
440! **************************************************************************************************
441!> \brief Explicit A, and B for ABBA, from the slabs the screening method contracts: the bare ones for
442!> TDHF and ALPHA screening, the W-multiplied ones otherwise
443!> \param fm_S_ia ...
444!> \param fm_S_ij ...
445!> \param fm_S_ab ...
446!> \param fm_S_bar_ia ...
447!> \param fm_S_bar_ij ...
448!> \param fm_A ...
449!> \param fm_B created only for ABBA
450!> \param do_abba ...
451!> \param Eigenval ...
452!> \param unit_nr ...
453!> \param homo ...
454!> \param virtual ...
455!> \param dimen_RI ...
456!> \param bse_env the BSE environment (settings)
457! **************************************************************************************************
458 SUBROUTINE create_a_and_b(fm_S_ia, fm_S_ij, fm_S_ab, fm_S_bar_ia, &
459 fm_S_bar_ij, fm_A, fm_B, do_abba, Eigenval, unit_nr, &
460 homo, virtual, dimen_RI, bse_env)
461
462 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_s_ia, fm_s_ij, fm_s_ab, fm_s_bar_ia, &
463 fm_s_bar_ij
464 TYPE(cp_fm_type), INTENT(INOUT) :: fm_a, fm_b
465 LOGICAL, INTENT(IN) :: do_abba
466 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
467 INTEGER, INTENT(IN) :: unit_nr
468 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
469 INTEGER, INTENT(IN) :: dimen_ri
470 TYPE(bse_env_type), INTENT(IN) :: bse_env
471
472 IF (bse_env%screening_method == bse_screening_tdhf .OR. &
473 bse_env%screening_method == bse_screening_alpha) THEN
474 CALL create_a(fm_s_ia, fm_s_ij, fm_s_ab, fm_a, eigenval, unit_nr, &
475 homo, virtual, dimen_ri, bse_env)
476 IF (do_abba) THEN
477 CALL create_b(fm_s_ia, fm_s_ia, fm_b, homo, virtual, dimen_ri, unit_nr, &
478 bse_env)
479 END IF
480 ELSE
481 CALL create_a(fm_s_ia, fm_s_bar_ij, fm_s_ab, fm_a, eigenval, unit_nr, &
482 homo, virtual, dimen_ri, bse_env)
483 IF (do_abba) THEN
484 CALL create_b(fm_s_ia, fm_s_bar_ia, fm_b, homo, virtual, dimen_ri, unit_nr, &
485 bse_env)
486 END IF
487 END IF
488
489 END SUBROUTINE create_a_and_b
490
491 ! **************************************************************************************************
492!> \brief Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
493!> (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
494!> We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
495!> of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
496!> \param fm_A ...
497!> \param fm_B ...
498!> \param fm_C ...
499!> \param fm_sqrt_A_minus_B ...
500!> \param fm_inv_sqrt_A_minus_B ...
501!> \param unit_nr ...
502!> \param bse_env the BSE environment (settings)
503!> \param diag_est ...
504! **************************************************************************************************
505 SUBROUTINE create_hermitian_form_of_abba(fm_A, fm_B, fm_C, &
506 fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
507 unit_nr, bse_env, diag_est)
508
509 TYPE(cp_fm_type), INTENT(IN) :: fm_a, fm_b
510 TYPE(cp_fm_type), INTENT(INOUT) :: fm_c, fm_sqrt_a_minus_b, &
511 fm_inv_sqrt_a_minus_b
512 INTEGER, INTENT(IN) :: unit_nr
513 TYPE(bse_env_type), INTENT(IN) :: bse_env
514 REAL(kind=dp), INTENT(IN) :: diag_est
515
516 CHARACTER(LEN=*), PARAMETER :: routinen = 'create_hermitian_form_of_ABBA'
517
518 INTEGER :: dim_mat, handle, n_dependent
519 REAL(kind=dp), DIMENSION(2) :: eigvals_ab_diff
520 TYPE(cp_fm_type) :: fm_a_minus_b, fm_a_plus_b, fm_dummy, &
521 fm_work_product
522
523 CALL timeset(routinen, handle)
524
525 IF (unit_nr > 0) THEN
526 WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', &
527 ' with size of A. This will take around ', diag_est, " s."
528 END IF
529
530 ! Create work matrices, which will hold A+B and A-B and their powers
531 ! C is created afterwards to save memory
532 ! Final result: C = (A-B)^0.5 (A+B) (A-B)^0.5 EQ.I
533 ! \_______/ \___/ \______/
534 ! fm_sqrt_A_minus_B fm_A_plus_B fm_sqrt_A_minus_B
535 ! (EQ.Ia) (EQ.Ib) (EQ.Ia)
536 ! Intermediate work matrices:
537 ! fm_inv_sqrt_A_minus_B: (A-B)^-0.5 EQ.II
538 ! fm_A_minus_B: (A-B) EQ.III
539 ! fm_work_product: (A-B)^0.5 (A+B) from (EQ.Ia) and (EQ.Ib) EQ.IV
540 CALL cp_fm_create(fm_a_plus_b, fm_a%matrix_struct)
541 CALL cp_fm_to_fm(fm_a, fm_a_plus_b)
542 CALL cp_fm_create(fm_a_minus_b, fm_a%matrix_struct)
543 CALL cp_fm_to_fm(fm_a, fm_a_minus_b)
544 CALL cp_fm_create(fm_sqrt_a_minus_b, fm_a%matrix_struct)
545 CALL cp_fm_set_all(fm_sqrt_a_minus_b, 0.0_dp)
546 CALL cp_fm_create(fm_inv_sqrt_a_minus_b, fm_a%matrix_struct)
547 CALL cp_fm_set_all(fm_inv_sqrt_a_minus_b, 0.0_dp)
548
549 CALL cp_fm_create(fm_work_product, fm_a%matrix_struct)
550
551 IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
552 WRITE (unit_nr, '(T2,A10,T13,A19)') 'BSE|DEBUG|', 'Created work arrays'
553 END IF
554
555 ! Add/Substract B (cf. EQs. Ib and III)
556 CALL cp_fm_scale_and_add(1.0_dp, fm_a_plus_b, 1.0_dp, fm_b)
557 CALL cp_fm_scale_and_add(1.0_dp, fm_a_minus_b, -1.0_dp, fm_b)
558
559 ! cp_fm_power will overwrite matrix, therefore we create copies
560 CALL cp_fm_to_fm(fm_a_minus_b, fm_inv_sqrt_a_minus_b)
561
562 ! In order to avoid a second diagonalization (cp_fm_power), we create (A-B)^0.5 (EQ.Ia)
563 ! from (A-B)^-0.5 (EQ.II) by multiplication with (A-B) (EQ.III) afterwards.
564
565 ! Raise A-B to -0.5_dp, no quenching of eigenvectors, hence threshold=0.0_dp
566 CALL cp_fm_create(fm_dummy, fm_a%matrix_struct)
567 ! Create (A-B)^-0.5 (cf. EQ.II)
568 CALL cp_fm_power(fm_inv_sqrt_a_minus_b, fm_dummy, -0.5_dp, 0.0_dp, n_dependent, eigvals=eigvals_ab_diff)
569 CALL cp_fm_release(fm_dummy)
570 IF (unit_nr > 0) THEN
571 WRITE (unit_nr, '(T2,A4,T7,A,T65,F16.6)') 'BSE|', &
572 'Smallest eigenvalue of A-B (eV)', eigvals_ab_diff(1)*evolt
573 END IF
574 ! Raise an error in case the the matrix A-B is not positive definite (i.e. negative eigenvalues)
575 ! In this case, the procedure for hermitian form of ABBA is not applicable
576 IF (eigvals_ab_diff(1) < 0) THEN
577 CALL cp_abort(__location__, &
578 "Matrix (A-B) is not positive definite. "// &
579 "Hermitian diagonalization of full ABBA matrix is ill-defined.")
580 END IF
581
582 ! We keep fm_inv_sqrt_A_minus_B for print of singleparticle transitions of ABBA
583 ! We further create (A-B)^0.5 for the singleparticle transitions of ABBA
584 ! Create (A-B)^0.5= (A-B)^-0.5 * (A-B) (EQ.Ia)
585 CALL cp_fm_get_info(fm_a, nrow_global=dim_mat)
586 CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_inv_sqrt_a_minus_b, fm_a_minus_b, 0.0_dp, &
587 fm_sqrt_a_minus_b)
588
589 ! Compute and store LHS of C, i.e. (A-B)^0.5 (A+B) (EQ.IV)
590 CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_sqrt_a_minus_b, fm_a_plus_b, 0.0_dp, &
591 fm_work_product)
592
593 ! Release to save memory
594 CALL cp_fm_release(fm_a_plus_b)
595 CALL cp_fm_release(fm_a_minus_b)
596
597 ! Now create full
598 CALL cp_fm_create(fm_c, fm_a%matrix_struct)
599 CALL cp_fm_set_all(fm_c, 0.0_dp)
600 ! Compute C=(A-B)^0.5 (A+B) (A-B)^0.5 (EQ.I)
601 CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_work_product, fm_sqrt_a_minus_b, 0.0_dp, &
602 fm_c)
603 CALL cp_fm_release(fm_work_product)
604
605 IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
606 WRITE (unit_nr, '(T2,A10,T13,A36)') 'BSE|DEBUG|', 'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5'
607 END IF
608
609 CALL timestop(handle)
610 END SUBROUTINE create_hermitian_form_of_abba
611
612! **************************************************************************************************
613!> \brief Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
614!> Here, the eigenvectors Z^n relate to X^n via
615!> Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
616!> \param fm_C ...
617!> \param homo ...
618!> \param virtual ...
619!> \param homo_irred ...
620!> \param fm_sqrt_A_minus_B ...
621!> \param fm_inv_sqrt_A_minus_B ...
622!> \param unit_nr ...
623!> \param diag_est ...
624!> \param bse_env the BSE environment (settings, and the state the post-processing reads)
625!> \param mo_coeff ...
626! **************************************************************************************************
627 SUBROUTINE diagonalize_c(fm_C, homo, virtual, homo_irred, &
628 fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
629 unit_nr, diag_est, bse_env, mo_coeff)
630
631 TYPE(cp_fm_type), INTENT(INOUT) :: fm_c
632 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
633 TYPE(cp_fm_type), INTENT(INOUT) :: fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b
634 INTEGER, INTENT(IN) :: unit_nr
635 REAL(kind=dp), INTENT(IN) :: diag_est
636 TYPE(bse_env_type), INTENT(IN) :: bse_env
637 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
638
639 CHARACTER(LEN=*), PARAMETER :: routinen = 'diagonalize_C'
640
641 INTEGER :: diag_info, handle, n_ov_joint, nspins
642 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
643 TYPE(cp_fm_type) :: fm_eigvec_x, fm_eigvec_y, fm_eigvec_z, &
644 fm_mat_eigvec_transform_diff, &
645 fm_mat_eigvec_transform_sum
646
647 CALL timeset(routinen, handle)
648
649 nspins = SIZE(homo)
650 n_ov_joint = sum(homo*virtual)
651
652 IF (unit_nr > 0) THEN
653 WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', &
654 'This will take around ', diag_est, ' s.'
655 END IF
656
657 !We have now the full matrix C=(A-B)^0.5 (A+B) (A-B)^0.5
658 !Now: Diagonalize it
659 CALL cp_fm_create(fm_eigvec_z, fm_c%matrix_struct)
660
661 ALLOCATE (exc_ens(n_ov_joint))
662
663 CALL choose_eigv_solver(fm_c, fm_eigvec_z, exc_ens, diag_info)
664
665 IF (diag_info /= 0) THEN
666 CALL cp_abort(__location__, &
667 "Diagonalization of C=(A-B)^0.5 (A+B) (A-B)^0.5 failed in BSE")
668 END IF
669
670 ! C could have negative eigenvalues, since we do not explicitly check A+B
671 ! for positive definiteness (would make another O(N^6) Diagon. necessary)
672 ! Instead, we include a check here
673 IF (exc_ens(1) < 0) THEN
674 IF (unit_nr > 0) THEN
675 CALL cp_abort(__location__, &
676 "Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
677 "(A+B) is not positive definite.")
678 END IF
679 END IF
680 exc_ens = sqrt(exc_ens)
681
682 ! Prepare eigenvector for interpretation of singleparticle transitions
683 ! Compare: F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)
684 ! We aim for the upper part of the vector (X,Y) for a direct comparison with the TDA result
685
686 ! Following Furche, we basically use Eqs. (A10): First, we multiply
687 ! the (A-B)^+-0.5 with eigenvectors and then the eigenvalues
688 ! One has to be careful about the index structure, since the eigenvector matrix is not symmetric anymore!
689
690 ! First, Eq. I from (A10) from Furche: (X+Y)_n = (Ω_n)^-0.5 (A-B)^0.5 T_n
691 CALL cp_fm_create(fm_mat_eigvec_transform_sum, fm_c%matrix_struct)
692 CALL cp_fm_set_all(fm_mat_eigvec_transform_sum, 0.0_dp)
693 CALL parallel_gemm(transa="N", transb="N", m=n_ov_joint, n=n_ov_joint, k=n_ov_joint, alpha=1.0_dp, &
694 matrix_a=fm_sqrt_a_minus_b, matrix_b=fm_eigvec_z, beta=0.0_dp, &
695 matrix_c=fm_mat_eigvec_transform_sum)
696 CALL cp_fm_release(fm_sqrt_a_minus_b)
697 ! This normalizes the eigenvectors
698 CALL comp_eigvec_coeff_bse(fm_mat_eigvec_transform_sum, exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.true.)
699
700 ! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n
701 CALL cp_fm_create(fm_mat_eigvec_transform_diff, fm_c%matrix_struct)
702 CALL cp_fm_set_all(fm_mat_eigvec_transform_diff, 0.0_dp)
703 CALL parallel_gemm(transa="N", transb="N", m=n_ov_joint, n=n_ov_joint, k=n_ov_joint, alpha=1.0_dp, &
704 matrix_a=fm_inv_sqrt_a_minus_b, matrix_b=fm_eigvec_z, beta=0.0_dp, &
705 matrix_c=fm_mat_eigvec_transform_diff)
706 CALL cp_fm_release(fm_inv_sqrt_a_minus_b)
707 CALL cp_fm_release(fm_eigvec_z)
708
709 ! This normalizes the eigenvectors
710 CALL comp_eigvec_coeff_bse(fm_mat_eigvec_transform_diff, exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.true.)
711
712 ! Now, we add the two equations to obtain X_n
713 ! Add overwrites the first argument, therefore we copy it beforehand
714 CALL cp_fm_create(fm_eigvec_x, fm_c%matrix_struct)
715 CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_x)
716 CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_x, 1.0_dp, fm_mat_eigvec_transform_diff)
717
718 ! Now, we subtract the two equations to obtain Y_n
719 ! Add overwrites the first argument, therefore we copy it beforehand
720 CALL cp_fm_create(fm_eigvec_y, fm_c%matrix_struct)
721 CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_y)
722 CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_y, -1.0_dp, fm_mat_eigvec_transform_diff)
723
724 !Cleanup
725 CALL cp_fm_release(fm_mat_eigvec_transform_diff)
726 CALL cp_fm_release(fm_mat_eigvec_transform_sum)
727
728 IF (nspins == 1) THEN
729 CALL postprocess_bse(exc_ens, fm_eigvec_x, bse_env, mo_coeff, &
730 homo(1), virtual(1), homo_irred(1), unit_nr, &
731 .false., fm_eigvec_y)
732 ELSE
733 ! Open-shell ABBA: helper forms X+Y internally and prints amplitudes (X and Y).
734 CALL bse_open_shell_optical(exc_ens, fm_eigvec_x, homo, virtual, homo_irred, &
735 .false., mo_coeff, bse_env, unit_nr, fm_eigvec_y)
736 END IF
737
738 DEALLOCATE (exc_ens)
739 CALL cp_fm_release(fm_eigvec_x)
740 CALL cp_fm_release(fm_eigvec_y)
741
742 CALL timestop(handle)
743
744 END SUBROUTINE diagonalize_c
745
746! **************************************************************************************************
747!> \brief Solving hermitian eigenvalue equation A X^n = Ω^n X^n
748!> \param fm_A ...
749!> \param homo ...
750!> \param virtual ...
751!> \param homo_irred ...
752!> \param unit_nr ...
753!> \param diag_est ...
754!> \param bse_env the BSE environment (settings, and the state the post-processing reads)
755!> \param mo_coeff ...
756! **************************************************************************************************
757 SUBROUTINE diagonalize_a(fm_A, homo, virtual, homo_irred, &
758 unit_nr, diag_est, bse_env, mo_coeff)
759
760 TYPE(cp_fm_type), INTENT(INOUT) :: fm_a
761 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
762 INTEGER, INTENT(IN) :: unit_nr
763 REAL(kind=dp), INTENT(IN) :: diag_est
764 TYPE(bse_env_type), INTENT(IN) :: bse_env
765 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
766
767 CHARACTER(LEN=*), PARAMETER :: routinen = 'diagonalize_A'
768
769 INTEGER :: diag_info, handle, n_ov_joint, nspins
770 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
771 TYPE(cp_fm_type) :: fm_eigvec
772
773 CALL timeset(routinen, handle)
774
775 nspins = SIZE(homo)
776 n_ov_joint = sum(homo*virtual)
777
778 IF (unit_nr > 0) THEN
779 WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
780 'This will take around ', diag_est, ' s.'
781 END IF
782
783 CALL cp_fm_create(fm_eigvec, fm_a%matrix_struct)
784
785 ALLOCATE (exc_ens(n_ov_joint))
786
787 CALL choose_eigv_solver(fm_a, fm_eigvec, exc_ens, diag_info)
788
789 IF (diag_info /= 0) THEN
790 CALL cp_abort(__location__, &
791 "Diagonalization of A failed in TDA-BSE")
792 END IF
793
794 IF (nspins == 1) THEN
795 CALL postprocess_bse(exc_ens, fm_eigvec, bse_env, mo_coeff, &
796 homo(1), virtual(1), homo_irred(1), unit_nr, .true.)
797 ELSE
798 CALL bse_open_shell_optical(exc_ens, fm_eigvec, homo, virtual, homo_irred, &
799 .true., mo_coeff, bse_env, unit_nr)
800 END IF
801
802 CALL cp_fm_release(fm_eigvec)
803 DEALLOCATE (exc_ens)
804
805 CALL timestop(handle)
806
807 END SUBROUTINE diagonalize_a
808
809! **************************************************************************************************
810!> \brief Open-shell (UKS) spin-summed post-processing for the joint spin-block space: joint
811!> excitation energies, per-spin transition amplitudes, and oscillator strengths. Mirrors
812!> postprocess_bse but spin-summed; exciton descriptors and NTOs are not yet implemented (CPWARN).
813!> \param Exc_ens joint excitation energies
814!> \param fm_eigvec_X joint X eigenvectors (excitations)
815!> \param homo per-spin reduced/active occupied counts
816!> \param virtual per-spin reduced/active virtual counts
817!> \param homo_irred per-spin full occupied counts (absolute-MO labels; N_e = sum)
818!> \param flag_tda .TRUE. -> TDA (coeff=X), .FALSE. -> ABBA (coeff=X+Y)
819!> \param mo_coeff per-spin MO coefficients
820!> \param bse_env the BSE environment (settings; the AO dipoles, mos)
821!> \param unit_nr ...
822!> \param fm_eigvec_Y joint Y eigenvectors (deexcitations; ABBA only)
823! **************************************************************************************************
824 SUBROUTINE bse_open_shell_optical(Exc_ens, fm_eigvec_X, homo, virtual, homo_irred, &
825 flag_tda, mo_coeff, bse_env, unit_nr, fm_eigvec_Y)
826
827 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
828 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_x
829 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
830 LOGICAL, INTENT(IN) :: flag_tda
831 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
832 TYPE(bse_env_type), INTENT(IN) :: bse_env
833 INTEGER, INTENT(IN) :: unit_nr
834 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_y
835
836 CHARACTER(LEN=*), PARAMETER :: routinen = 'bse_open_shell_optical'
837
838 CHARACTER(LEN=10) :: info_approximation, multiplet
839 INTEGER :: handle, idir, isp, jdir, n, n_ov_joint, &
840 nspins
841 INTEGER, ALLOCATABLE, DIMENSION(:) :: n_ov_sp, offsets_sp
842 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str_joint
843 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pol_res_joint, trans_mom_joint
844 TYPE(cp_fm_struct_type), POINTER :: fm_struct_dip_reord, fm_struct_sp, &
845 fm_struct_tmom
846 TYPE(cp_fm_type) :: fm_dip_reord_sp, fm_eigvec_sp, &
847 fm_trans_coeff
848 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dip_ab_sp, fm_dip_ai_sp, fm_dip_ij_sp
849 TYPE(cp_fm_type), DIMENSION(3) :: fm_trans_mom_joint
850
851 CALL timeset(routinen, handle)
852
853 nspins = SIZE(homo)
854 ALLOCATE (n_ov_sp(nspins), offsets_sp(nspins))
855 CALL get_bse_spin_block_layout(homo, virtual, n_ov_sp, offsets_sp, n_ov_joint)
856
857 ! LEN=10 locals auto-pad short literals with spaces (avoids the L-15 short-literal trap);
858 ! print_excitation_energies prints A6 of these, print_optical_properties prints them as-is.
859 multiplet = "UKS"
860 IF (flag_tda) THEN
861 info_approximation = " -TDA- "
862 ELSE
863 info_approximation = "-ABBA-"
864 END IF
865
866 IF (unit_nr > 0) THEN
867 WRITE (unit_nr, '(T2,A4,T7,A43)') 'BSE|', 'Joint open-shell BSE excitation energies:'
868 END IF
869 CALL print_excitation_energies(exc_ens, n_ov_joint, 1, flag_tda, multiplet, &
870 info_approximation, bse_env, unit_nr)
871
872 ! Per-spin single-particle transition amplitudes (X via =>, Y via <=).
873 CALL print_transition_amplitudes(fm_eigvec_x, homo, virtual, homo_irred, &
874 info_approximation, bse_env, unit_nr, fm_eigvec_y)
875
876 ! Transition coefficient for the spin-summed moment: X (TDA) or X+Y (ABBA).
877 CALL cp_fm_create(fm_trans_coeff, fm_eigvec_x%matrix_struct)
878 CALL cp_fm_to_fm(fm_eigvec_x, fm_trans_coeff)
879 IF (PRESENT(fm_eigvec_y)) CALL cp_fm_scale_and_add(1.0_dp, fm_trans_coeff, 1.0_dp, fm_eigvec_y)
880
881 ! Spin-summed transition moments: D^n_dir = sum_σ sum_{ia,σ} D^{dir,σ}_{ai} C_{ia,σ,n}
882 ! with C = X (TDA) or X+Y (ABBA); explicit spin sum, factor 1.0.
883 ALLOCATE (fm_dip_ai_sp(3), fm_dip_ij_sp(3), fm_dip_ab_sp(3))
884 ALLOCATE (oscill_str_joint(n_ov_joint), trans_mom_joint(3, 1, n_ov_joint))
885 ALLOCATE (pol_res_joint(3, 3, n_ov_joint))
886 trans_mom_joint(:, :, :) = 0.0_dp
887 NULLIFY (fm_struct_dip_reord, fm_struct_sp, fm_struct_tmom)
888 CALL cp_fm_struct_create(fm_struct_tmom, fm_trans_coeff%matrix_struct%para_env, &
889 fm_trans_coeff%matrix_struct%context, 1, n_ov_joint)
890 DO idir = 1, 3
891 CALL cp_fm_create(fm_trans_mom_joint(idir), fm_struct_tmom)
892 CALL cp_fm_set_all(fm_trans_mom_joint(idir), 0.0_dp)
893 END DO
894 DO isp = 1, nspins
895 CALL get_multipoles_mo(fm_dip_ai_sp, fm_dip_ij_sp, fm_dip_ab_sp, &
896 bse_env%matrix_dipole_ao, bse_env%mos, mo_coeff(isp:isp), &
897 homo(isp), virtual(isp), fm_trans_coeff%matrix_struct%context, &
898 ispin=isp)
899 NULLIFY (fm_struct_sp, fm_struct_dip_reord)
900 CALL cp_fm_struct_create(fm_struct_sp, fm_trans_coeff%matrix_struct%para_env, &
901 fm_trans_coeff%matrix_struct%context, n_ov_sp(isp), n_ov_joint)
902 CALL cp_fm_create(fm_eigvec_sp, fm_struct_sp)
903 CALL cp_fm_set_all(fm_eigvec_sp, 0.0_dp)
904 CALL cp_fm_to_fm_submat(fm_trans_coeff, fm_eigvec_sp, n_ov_sp(isp), n_ov_joint, &
905 offsets_sp(isp) + 1, 1, 1, 1)
906 CALL cp_fm_struct_create(fm_struct_dip_reord, fm_trans_coeff%matrix_struct%para_env, &
907 fm_trans_coeff%matrix_struct%context, 1, n_ov_sp(isp))
908 DO idir = 1, 3
909 CALL cp_fm_create(fm_dip_reord_sp, fm_struct_dip_reord, name="bse_dip_reord")
910 CALL cp_fm_set_all(fm_dip_reord_sp, 0.0_dp)
911 CALL fm_general_add_bse(fm_dip_reord_sp, fm_dip_ai_sp(idir), 1.0_dp, &
912 1, 1, 1, virtual(isp), unit_nr, [2, 4, 3, 1], bse_env)
913 CALL parallel_gemm('N', 'N', 1, n_ov_joint, n_ov_sp(isp), 1.0_dp, &
914 fm_dip_reord_sp, fm_eigvec_sp, 1.0_dp, fm_trans_mom_joint(idir))
915 CALL cp_fm_release(fm_dip_reord_sp)
916 CALL cp_fm_release(fm_dip_ai_sp(idir))
917 CALL cp_fm_release(fm_dip_ij_sp(idir))
918 CALL cp_fm_release(fm_dip_ab_sp(idir))
919 END DO
920 CALL cp_fm_release(fm_eigvec_sp)
921 CALL cp_fm_struct_release(fm_struct_sp)
922 NULLIFY (fm_struct_sp)
923 CALL cp_fm_struct_release(fm_struct_dip_reord)
924 NULLIFY (fm_struct_dip_reord)
925 END DO
926 DO idir = 1, 3
927 CALL cp_fm_get_submatrix(fm_trans_mom_joint(idir), trans_mom_joint(idir, :, :))
928 CALL cp_fm_release(fm_trans_mom_joint(idir))
929 END DO
930 CALL cp_fm_struct_release(fm_struct_tmom)
931 DO n = 1, n_ov_joint
932 DO idir = 1, 3
933 DO jdir = 1, 3
934 pol_res_joint(idir, jdir, n) = 2.0_dp*exc_ens(n)*trans_mom_joint(idir, 1, n) &
935 *trans_mom_joint(jdir, 1, n)
936 END DO
937 END DO
938 oscill_str_joint(n) = 2.0_dp/3.0_dp*exc_ens(n)*sum(abs(trans_mom_joint(:, 1, n))**2)
939 END DO
940 CALL print_optical_properties(exc_ens, oscill_str_joint, trans_mom_joint, pol_res_joint, &
941 n_ov_joint, 1, sum(homo_irred), flag_tda, info_approximation, &
942 bse_env, unit_nr, open_shell=.true.)
943 ! Open-shell post-processing is partial: energies, amplitudes, spin-summed oscillator strengths.
944 CALL cp_warn(__location__, &
945 "Open-shell (UKS) BSE: exciton descriptors and NTO analysis are not yet "// &
946 "implemented and have been skipped.")
947 CALL cp_fm_release(fm_trans_coeff)
948 DEALLOCATE (fm_dip_ai_sp, fm_dip_ij_sp, fm_dip_ab_sp)
949 DEALLOCATE (n_ov_sp, offsets_sp, oscill_str_joint, trans_mom_joint, pol_res_joint)
950
951 CALL timestop(handle)
952
953 END SUBROUTINE bse_open_shell_optical
954
955! **************************************************************************************************
956!> \brief Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag)
957!> \param Exc_ens ...
958!> \param fm_eigvec_X ...
959!> \param bse_env the BSE environment (settings; the AO multipoles and their reference point, mos, subsys)
960!> \param mo_coeff ...
961!> \param homo ...
962!> \param virtual ...
963!> \param homo_irred ...
964!> \param unit_nr ...
965!> \param flag_TDA ...
966!> \param fm_eigvec_Y ...
967! **************************************************************************************************
968 SUBROUTINE postprocess_bse(Exc_ens, fm_eigvec_X, bse_env, mo_coeff, &
969 homo, virtual, homo_irred, unit_nr, &
970 flag_TDA, fm_eigvec_Y)
971
972 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: exc_ens
973 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_x
974 TYPE(bse_env_type), INTENT(IN) :: bse_env
975 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
976 INTEGER :: homo, virtual, homo_irred, unit_nr
977 LOGICAL, OPTIONAL :: flag_tda
978 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_y
979
980 CHARACTER(LEN=*), PARAMETER :: routinen = 'postprocess_bse'
981
982 CHARACTER(LEN=10) :: info_approximation, multiplet
983 INTEGER :: handle, i_exc, idir, n_moments_di, &
984 n_moments_quad
985 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str, ref_point_multipole
986 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: polarizability_residues, trans_mom_bse
987 TYPE(cp_fm_type) :: fm_x_ia, fm_y_ia
988 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
989 fm_dipole_ij_trunc, fm_quadpole_ab_trunc, fm_quadpole_ai_trunc, fm_quadpole_ij_trunc
990 TYPE(exciton_descr_type), ALLOCATABLE, &
991 DIMENSION(:) :: exc_descr
992 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
993
994 CALL timeset(routinen, handle)
995
996 !Prepare variables for printing
997 IF (bse_env%bse_spin_config == 0) THEN
998 multiplet = "Singlet"
999 ELSE
1000 multiplet = "Triplet"
1001 END IF
1002 IF (.NOT. PRESENT(flag_tda)) THEN
1003 flag_tda = .false.
1004 END IF
1005 IF (flag_tda) THEN
1006 info_approximation = " -TDA- "
1007 ELSE
1008 info_approximation = "-ABBA-"
1009 END IF
1010
1011 n_moments_di = 3
1012 n_moments_quad = 9
1013 ! Compute BSE dipoles and oscillator strengths - Keep in memory for later usage
1014 ! Need dipoles also for spatial expectation values, which are well-defined also for triplets
1015 ALLOCATE (fm_dipole_ij_trunc(n_moments_di))
1016 ALLOCATE (fm_dipole_ab_trunc(n_moments_di))
1017 ALLOCATE (fm_dipole_ai_trunc(n_moments_di))
1018 ALLOCATE (ref_point_multipole(3))
1019 ref_point_multipole(:) = bse_env%multipole_ref_point
1020 ! Obtain dipoles in MO basis
1021 CALL get_multipoles_mo(fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc, &
1022 bse_env%matrix_dipole_ao, bse_env%mos, mo_coeff, &
1023 homo, virtual, fm_eigvec_x%matrix_struct%context)
1024 ! Compute exciton descriptors from these multipoles
1025 IF (bse_env%num_print_exc_descr > 0) THEN
1026 ! Obtain quadrupoles in MO basis
1027 ALLOCATE (fm_quadpole_ij_trunc(n_moments_quad))
1028 ALLOCATE (fm_quadpole_ab_trunc(n_moments_quad))
1029 ALLOCATE (fm_quadpole_ai_trunc(n_moments_quad))
1030 ! the preparer builds the AO quadrupoles whenever NUM_PRINT_EXC_DESCR is nonzero
1031 cpassert(ASSOCIATED(bse_env%matrix_quadpole_ao))
1032 CALL get_multipoles_mo(fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1033 bse_env%matrix_quadpole_ao, bse_env%mos, mo_coeff, &
1034 homo, virtual, fm_eigvec_x%matrix_struct%context)
1035 ! Iterate over excitation index outside of routine to make it compatible with tddft module
1036 ALLOCATE (exc_descr(bse_env%num_print_exc_descr))
1037 DO i_exc = 1, bse_env%num_print_exc_descr
1038 CALL reshuffle_eigvec(fm_eigvec_x, fm_x_ia, homo, virtual, i_exc, &
1039 .false., unit_nr, bse_env)
1040 IF (.NOT. flag_tda) THEN
1041 CALL reshuffle_eigvec(fm_eigvec_y, fm_y_ia, homo, virtual, i_exc, &
1042 .false., unit_nr, bse_env)
1043
1044 CALL get_exciton_descriptors(exc_descr, fm_x_ia, &
1045 fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1046 fm_quadpole_ai_trunc, &
1047 i_exc, homo, virtual, &
1048 fm_y_ia)
1049 ELSE
1050 CALL get_exciton_descriptors(exc_descr, fm_x_ia, &
1051 fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1052 fm_quadpole_ai_trunc, &
1053 i_exc, homo, virtual)
1054 END IF
1055 CALL cp_fm_release(fm_x_ia)
1056 IF (.NOT. flag_tda) THEN
1057 CALL cp_fm_release(fm_y_ia)
1058 END IF
1059 END DO
1060 END IF
1061
1062 IF (bse_env%bse_spin_config == 0) THEN
1063 CALL get_oscillator_strengths(fm_eigvec_x, exc_ens, fm_dipole_ai_trunc, &
1064 trans_mom_bse, oscill_str, polarizability_residues, &
1065 bse_env, homo, virtual, unit_nr, &
1066 fm_eigvec_y)
1067 END IF
1068
1069 ! the iterative driver prints the header before its solver runs
1070 IF (bse_env%bse_diag_method /= bse_iterdiag) THEN
1071 CALL print_output_header(homo, virtual, homo_irred, flag_tda, bse_env, unit_nr)
1072 END IF
1073
1074 ! Prints excitation energies up to user-specified number
1075 CALL print_excitation_energies(exc_ens, homo, virtual, flag_tda, multiplet, &
1076 info_approximation, bse_env, unit_nr)
1077
1078 ! Print single particle transition amplitudes, i.e. components of eigenvectors X and Y
1079 CALL print_transition_amplitudes(fm_eigvec_x, [homo], [virtual], [homo_irred], &
1080 info_approximation, bse_env, unit_nr, fm_eigvec_y)
1081
1082 ! Prints optical properties, if state is a singlet
1083 CALL print_optical_properties(exc_ens, oscill_str, trans_mom_bse, polarizability_residues, &
1084 homo, virtual, homo_irred, flag_tda, &
1085 info_approximation, bse_env, unit_nr)
1086 ! Print exciton descriptors if keyword is invoked
1087 IF (bse_env%num_print_exc_descr > 0) THEN
1088 CALL qs_subsys_get(bse_env%subsys, particle_set=particle_set)
1089 CALL print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, &
1090 bse_env%num_print_exc_descr, bse_env%bse_debug_print, &
1091 bse_env%print_directional_exc_descr, &
1092 bse_env%print_directional_crosscorrelation, &
1093 'BSE|', particle_set)
1094 END IF
1095
1096 ! Compute and print excitation wavefunctions
1097 IF (bse_env%do_nto_analysis) THEN
1098 IF (unit_nr > 0) THEN
1099 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1100 WRITE (unit_nr, '(T2,A4,T7,A47)') &
1101 'BSE|', "Calculating Natural Transition Orbitals (NTOs)."
1102 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1103 END IF
1104 CALL calculate_ntos(fm_eigvec_x, fm_eigvec_y, &
1105 mo_coeff, homo, virtual, &
1106 info_approximation, &
1107 oscill_str, &
1108 bse_env, unit_nr)
1109 END IF
1110
1111 DO idir = 1, n_moments_di
1112 CALL cp_fm_release(fm_dipole_ai_trunc(idir))
1113 CALL cp_fm_release(fm_dipole_ij_trunc(idir))
1114 CALL cp_fm_release(fm_dipole_ab_trunc(idir))
1115 END DO
1116 IF (bse_env%num_print_exc_descr > 0) THEN
1117 DO idir = 1, n_moments_quad
1118 CALL cp_fm_release(fm_quadpole_ai_trunc(idir))
1119 CALL cp_fm_release(fm_quadpole_ij_trunc(idir))
1120 CALL cp_fm_release(fm_quadpole_ab_trunc(idir))
1121 END DO
1122 DEALLOCATE (fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc)
1123 DEALLOCATE (exc_descr)
1124 END IF
1125 DEALLOCATE (fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc)
1126 DEALLOCATE (ref_point_multipole)
1127 IF (bse_env%bse_spin_config == 0) THEN
1128 DEALLOCATE (oscill_str, trans_mom_bse, polarizability_residues)
1129 END IF
1130 IF (unit_nr > 0) THEN
1131 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1132 WRITE (unit_nr, '(T2,A4)') 'BSE|'
1133 END IF
1134
1135 CALL timestop(handle)
1136
1137 END SUBROUTINE postprocess_bse
1138
1139END MODULE bse_full_diag
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
bse_full_diag
Routines for the full diagonalization of GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_full_diag.F:14
bse_full_diag::postprocess_bse
subroutine, public postprocess_bse(exc_ens, fm_eigvec_x, bse_env, mo_coeff, homo, virtual, homo_irred, unit_nr, flag_tda, fm_eigvec_y)
Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag).
Definition bse_full_diag.F:971
bse_full_diag::create_b
subroutine, public create_b(fm_s_ia, fm_s_bar_ia, fm_b, homo, virtual, dimen_ri, unit_nr, bse_env)
Matrix B constructed from 3c-B-matrices (cf. subroutine screen_slabs) B_ia,jb = α * v_ia,...
Definition bse_full_diag.F:331
bse_full_diag::create_a_and_b
subroutine, public create_a_and_b(fm_s_ia, fm_s_ij, fm_s_ab, fm_s_bar_ia, fm_s_bar_ij, fm_a, fm_b, do_abba, eigenval, unit_nr, homo, virtual, dimen_ri, bse_env)
Explicit A, and B for ABBA, from the slabs the screening method contracts: the bare ones for TDHF and...
Definition bse_full_diag.F:461
bse_full_diag::diagonalize_a
subroutine, public diagonalize_a(fm_a, homo, virtual, homo_irred, unit_nr, diag_est, bse_env, mo_coeff)
Solving hermitian eigenvalue equation A X^n = Ω^n X^n.
Definition bse_full_diag.F:759
bse_full_diag::create_hermitian_form_of_abba
subroutine, public create_hermitian_form_of_abba(fm_a, fm_b, fm_c, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, unit_nr, bse_env, diag_est)
Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem (cf....
Definition bse_full_diag.F:508
bse_full_diag::diagonalize_c
subroutine, public diagonalize_c(fm_c, homo, virtual, homo_irred, fm_sqrt_a_minus_b, fm_inv_sqrt_a_minus_b, unit_nr, diag_est, bse_env, mo_coeff)
Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n . Here, the eigenvectors Z^n relate to X^n via Eq....
Definition bse_full_diag.F:630
bse_full_diag::create_a
subroutine, public create_a(fm_s_ia, fm_s_bar_ij, fm_s_ab, fm_a, eigenval, unit_nr, homo, virtual, dimen_ri, bse_env)
Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine screen_slabs) A_ia,...
Definition bse_full_diag.F:99
bse_print
Routines for printing information in context of the BSE calculation.
Definition bse_print.F:13
bse_print::print_transition_amplitudes
subroutine, public print_transition_amplitudes(fm_eigvec_x, homo, virtual, homo_irred, info_approximation, bse_env, unit_nr, fm_eigvec_y)
...
Definition bse_print.F:363
bse_print::print_excitation_energies
subroutine, public print_excitation_energies(exc_ens, homo, virtual, flag_tda, multiplet, info_approximation, bse_env, unit_nr)
...
Definition bse_print.F:314
bse_print::print_output_header
subroutine, public print_output_header(homo, virtual, homo_irred, flag_tda, bse_env, unit_nr)
Banner, defining equations and legend of one BSE problem, printed before it is solved.
Definition bse_print.F:100
bse_print::print_optical_properties
subroutine, public print_optical_properties(exc_ens, oscill_str, trans_mom_bse, polarizability_residues, homo, virtual, homo_irred, flag_tda, info_approximation, bse_env, unit_nr, open_shell)
...
Definition bse_print.F:458
bse_print::print_exciton_descriptors
subroutine, public print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, num_print_exc_descr, print_checkvalue, print_directional_exc_descr, print_directional_crosscorrelation, prefix_output, particle_set)
Prints exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018).
Definition bse_print.F:669
bse_properties
Routines for computing excitonic properties, e.g. exciton diameter, from the BSE.
Definition bse_properties.F:13
bse_properties::get_oscillator_strengths
subroutine, public get_oscillator_strengths(fm_eigvec_x, exc_ens, fm_dipole_ai_trunc, trans_mom_bse, oscill_str, polarizability_residues, bse_env, homo_red, virtual_red, unit_nr, fm_eigvec_y)
Compute and return BSE dipoles d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n ) and oscil...
Definition bse_properties.F:105
bse_properties::get_exciton_descriptors
subroutine, public get_exciton_descriptors(exc_descr, fm_x_ia, fm_multipole_ij_trunc, fm_multipole_ab_trunc, fm_multipole_ai_trunc, i_exc, homo, virtual, fm_y_ia)
...
Definition bse_properties.F:748
bse_properties::calculate_ntos
subroutine, public calculate_ntos(fm_x, fm_y, mo_coeff, homo, virtual, info_approximation, oscill_str, bse_env, unit_nr)
...
Definition bse_properties.F:419
bse_types
The BSE environment: the settings of the &BSE section and the state a GW path prepares for the solver...
Definition bse_types.F:14
bse_util
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_util.F:13
bse_util::assemble_joint_ov_slab
subroutine, public assemble_joint_ov_slab(fm_s_ia, offsets, n_ov, dimen_ri, fm_s_joint)
Column-concatenate per-spin ia-slabs into the joint dimen_RI x n_ov_joint slab. Sigma-block of spin i...
Definition bse_util.F:1940
bse_util::get_multipoles_mo
subroutine, public get_multipoles_mo(fm_multipole_ai_trunc, fm_multipole_ij_trunc, fm_multipole_ab_trunc, matrix_multipole, mos, mo_coeff, homo_red, virtual_red, context_bse, ispin)
The multipoles in the MO window, D^k_pq = sum_µν C_µp M^k_µν C_νq, cut to the ai, ij and ab blocks of...
Definition bse_util.F:1743
bse_util::occ_of_ia
pure integer function, public occ_of_ia(ia, virt)
Occupied level of the compound transition index ia = (i-1)*virt + a.
Definition bse_util.F:112
bse_util::reshuffle_eigvec
subroutine, public reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, unit_nr, bse_env)
...
Definition bse_util.F:1363
bse_util::get_bse_spin_block_layout
subroutine, public get_bse_spin_block_layout(homo_red, virt_red, n_ov, offsets, n_ov_joint)
Spin-block layout for the open-shell (joint) BSE matrix: per-spin OV-pair counts and the block offset...
Definition bse_util.F:1115
bse_util::virt_of_ia
pure integer function, public virt_of_ia(ia, virt)
Virtual level of the compound transition index ia = (i-1)*virt + a.
Definition bse_util.F:127
bse_util::comp_eigvec_coeff_bse
subroutine, public comp_eigvec_coeff_bse(fm_work, eig_vals, beta, gamma, do_transpose)
Routine for computing the coefficients of the eigenvectors of the BSE matrix from a multiplication wi...
Definition bse_util.F:756
bse_util::fm_general_add_bse
subroutine, public fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, bse_env, row_offset, col_offset)
Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,...
Definition bse_util.F:197
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:171
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:1144
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1044
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1461
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:931
exstates_types
Types for excited states potential energies.
Definition exstates_types.F:14
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::bse_screening_tdhf
integer, parameter, public bse_screening_tdhf
Definition input_constants.F:1401
input_constants::bse_iterdiag
integer, parameter, public bse_iterdiag
Definition input_constants.F:1392
input_constants::bse_singlet
integer, parameter, public bse_singlet
Definition input_constants.F:1388
input_constants::bse_triplet
integer, parameter, public bse_triplet
Definition input_constants.F:1388
input_constants::bse_screening_alpha
integer, parameter, public bse_screening_alpha
Definition input_constants.F:1401
input_constants::bse_screening_rpa
integer, parameter, public bse_screening_rpa
Definition input_constants.F:1401
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:137
bse_properties::exciton_descr_type
Definition bse_properties.F:57
bse_types::bse_env_type
Settings of the &BSE section (read_bse_section, re-read before every BSE run, normalised in place by ...
Definition bse_types.F:43
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:787
cp_control_types::tddfpt2_control_type
Definition cp_control_types.F:663
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
exstates_types::excited_energy_type
Contains information on the excited states energy.
Definition exstates_types.F:53
particle_types::particle_type
Definition particle_types.F:35