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 Routines for the full diagonalization of GW + Bethe-Salpeter for computing
10 : !> electronic excitations
11 : !> \par History
12 : !> 10.2023 created [Maximilian Graml]
13 : ! **************************************************************************************************
14 : MODULE 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 : get_bse_spin_block_layout,&
29 : get_multipoles_mo,&
30 : occ_of_ia,&
31 : reshuffle_eigvec,&
32 : virt_of_ia
33 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
34 : cp_blacs_env_release,&
35 : cp_blacs_env_type
36 : USE cp_control_types, ONLY: dft_control_type,&
37 : tddfpt2_control_type
38 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
39 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
40 : cp_fm_power
41 : USE cp_fm_redistribution, ONLY: cp_fm_permute
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 :
76 : CONTAINS
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 54 : SUBROUTINE create_A(fm_S_ia, fm_S_bar_ij, fm_S_ab, &
97 54 : fm_A, Eigenval, unit_nr, &
98 54 : 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 54 : INTEGER, ALLOCATABLE, DIMENSION(:) :: eig_offsets, n_ov, offsets
113 : INTEGER, DIMENSION(4) :: reordering
114 54 : 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 54 : CALL timeset(routineN, handle)
125 :
126 54 : nspins = SIZE(homo)
127 270 : ALLOCATE (n_ov(nspins), offsets(nspins), eig_offsets(nspins))
128 54 : 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 54 : eig_offsets(1) = 0
131 68 : DO isp = 2, nspins
132 68 : eig_offsets(isp) = eig_offsets(isp - 1) + homo(isp - 1) + virtual(isp - 1)
133 : END DO
134 :
135 54 : NULLIFY (dft_control, tddfpt_control)
136 54 : dft_control => bse_env%dft_control
137 54 : tddfpt_control => dft_control%tddfpt2_control
138 :
139 54 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
140 3 : 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 108 : SELECT CASE (bse_env%bse_spin_config)
145 : CASE (bse_singlet)
146 54 : alpha = 2.0_dp
147 : CASE (bse_triplet)
148 54 : alpha = 0.0_dp
149 : END SELECT
150 : ! For open-shell (nspins>1): each spin block contributes once; SPIN_CONFIG is ignored.
151 54 : IF (nspins > 1) THEN
152 : CALL cp_warn(__LOCATION__, &
153 14 : "BSE: SPIN_CONFIG ignored for open-shell reference; using alpha=1.")
154 14 : alpha = 1.0_dp
155 : END IF
156 :
157 54 : IF (bse_env%screening_method == bse_screening_alpha) THEN
158 2 : alpha_screening = bse_env%screening_factor
159 : ELSE
160 52 : 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 54 : NULLIFY (blacs_env)
165 54 : 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 54 : para_env=fm_S_ia(1)%matrix_struct%para_env)
175 54 : CALL cp_fm_create(fm_A, fm_struct_A, name="fm_A_iajb")
176 54 : 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 54 : IF (tddfpt_control%do_bse_w_only .AND. nspins == 1) THEN
179 2 : CALL cp_fm_create(fm_A_copy, fm_struct_A, name="fm_A_iajb")
180 2 : 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 54 : IF ((.NOT. tddfpt_control%do_bse) .AND. (.NOT. tddfpt_control%do_bse_w_only)) THEN
186 50 : 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 14 : para_env=fm_S_ia(1)%matrix_struct%para_env)
192 14 : CALL cp_fm_create(fm_S_joint, fm_struct_S_joint, name="fm_S_ia_joint")
193 14 : CALL cp_fm_set_all(fm_S_joint, 0.0_dp)
194 14 : 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 14 : matrix_c=fm_A)
198 14 : CALL cp_fm_release(fm_S_joint)
199 14 : 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 36 : beta=0.0_dp, matrix_c=fm_A)
205 : END IF
206 : END IF
207 :
208 54 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
209 3 : 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 122 : DO isp = 1, nspins
216 122 : 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 66 : para_env=fm_S_ab(isp)%matrix_struct%para_env)
220 66 : CALL cp_fm_create(fm_W, fm_struct_W, name="fm_W_ijab")
221 66 : 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 66 : beta=0.0_dp, matrix_c=fm_W)
227 66 : reordering = [1, 3, 2, 4]
228 : CALL cp_fm_permute(fm_A, fm_W, -1.0_dp, homo(isp), virtual(isp), virtual(isp), virtual(isp), reordering, &
229 66 : row_offset=offsets(isp), col_offset=offsets(isp))
230 66 : IF (nspins == 1 .AND. tddfpt_control%do_bse_w_only) THEN
231 2 : CALL cp_fm_permute(fm_A_copy, fm_W, -1.0_dp, homo(1), virtual(1), virtual(1), virtual(1), reordering)
232 : END IF
233 : ! W and A stash for TDDFPT path (closed-shell only; open-shell deferred)
234 : IF (nspins == 1) THEN
235 38 : IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
236 : tddfpt_control%do_bse_gw_only) THEN
237 4 : NULLIFY (ex_env)
238 4 : ex_env => bse_env%exstate_env
239 4 : IF (.NOT. tddfpt_control%do_bse_gw_only) THEN
240 12 : ALLOCATE (ex_env%bse_w_matrix_MO(1, 1))
241 12 : ALLOCATE (ex_env%bse_a_matrix_MO(1, 1))
242 4 : CALL cp_fm_create(ex_env%bse_w_matrix_MO(1, 1), fm_struct_W)
243 4 : CALL cp_fm_create(ex_env%bse_a_matrix_MO(1, 1), fm_struct_A)
244 4 : CALL cp_fm_to_fm(fm_W, ex_env%bse_w_matrix_MO(1, 1))
245 4 : IF (tddfpt_control%do_bse_w_only) THEN
246 2 : CALL cp_fm_to_fm(fm_A_copy, ex_env%bse_a_matrix_MO(1, 1))
247 : ELSE
248 2 : CALL cp_fm_to_fm(fm_A, ex_env%bse_a_matrix_MO(1, 1))
249 : END IF
250 : END IF
251 : END IF
252 : END IF
253 66 : CALL cp_fm_release(fm_W)
254 132 : CALL cp_fm_struct_release(fm_struct_W)
255 : END IF
256 : END DO
257 :
258 54 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
259 3 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
260 : END IF
261 54 : IF (nspins == 1 .AND. tddfpt_control%do_bse_w_only) CALL cp_fm_release(fm_A_copy)
262 :
263 : ! Get local row/col indices for direct diagonal access
264 : CALL cp_fm_get_info(matrix=fm_A, nrow_local=nrow_local_A, ncol_local=ncol_local_A, &
265 54 : row_indices=row_indices_A, col_indices=col_indices_A)
266 :
267 : !Add (ε_a-ε_i) on the diagonal of each sigma-block; cross-spin blocks have no ε contribution.
268 54 : IF (.NOT. tddfpt_control%do_bse) THEN
269 6744 : DO ii = 1, nrow_local_A
270 6692 : i_row_global = row_indices_A(ii)
271 2133688 : DO jj = 1, ncol_local_A
272 2126944 : j_col_global = col_indices_A(jj)
273 2133636 : IF (i_row_global == j_col_global) THEN
274 : ! Decode spin: isp such that i_row_global in [offsets(isp)+1, offsets(isp)+n_ov(isp)]
275 5170 : isp = nspins
276 5170 : DO k_isp = 1, nspins - 1
277 5170 : IF (i_row_global <= offsets(k_isp) + n_ov(k_isp)) THEN
278 : isp = k_isp
279 : EXIT
280 : END IF
281 : END DO
282 3892 : k_ov = i_row_global - offsets(isp)
283 3892 : i_occ_row = occ_of_ia(k_ov, virtual(isp))
284 3892 : a_virt_row = virt_of_ia(k_ov, virtual(isp))
285 : eigen_diff = Eigenval(eig_offsets(isp) + a_virt_row + homo(isp)) - &
286 3892 : Eigenval(eig_offsets(isp) + i_occ_row)
287 3892 : fm_A%local_data(ii, jj) = fm_A%local_data(ii, jj) + eigen_diff
288 : END IF
289 : END DO
290 : END DO
291 : END IF
292 :
293 : ! GW eigenvalue stash for TDDFPT path (closed-shell only)
294 54 : IF (nspins == 1) THEN
295 40 : IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. &
296 : tddfpt_control%do_bse_gw_only) THEN
297 4 : sizeeigen = SIZE(Eigenval)
298 12 : ALLOCATE (ex_env%gw_eigen(sizeeigen))
299 132 : ex_env%gw_eigen(:) = Eigenval(:)
300 : END IF
301 : END IF
302 :
303 54 : CALL cp_fm_struct_release(fm_struct_A)
304 54 : DEALLOCATE (n_ov, offsets, eig_offsets)
305 :
306 54 : CALL cp_blacs_env_release(blacs_env)
307 :
308 54 : CALL timestop(handle)
309 :
310 162 : END SUBROUTINE create_A
311 :
312 : ! **************************************************************************************************
313 : !> \brief Matrix B constructed from 3c-B-matrices (cf. subroutine screen_slabs)
314 : !> B_ia,jb = α * v_ia,jb - W_ib,aj
315 : !> α is a spin-dependent factor
316 : !> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
317 : !> W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
318 : !> \param fm_S_ia ...
319 : !> \param fm_S_bar_ia ...
320 : !> \param fm_B ...
321 : !> \param homo ...
322 : !> \param virtual ...
323 : !> \param dimen_RI ...
324 : !> \param unit_nr ...
325 : !> \param bse_env the BSE environment (settings)
326 : ! **************************************************************************************************
327 26 : SUBROUTINE create_B(fm_S_ia, fm_S_bar_ia, fm_B, &
328 26 : homo, virtual, dimen_RI, unit_nr, bse_env)
329 :
330 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_S_ia, fm_S_bar_ia
331 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_B
332 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
333 : INTEGER, INTENT(IN) :: dimen_RI, unit_nr
334 : TYPE(bse_env_type), INTENT(IN) :: bse_env
335 :
336 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_B'
337 :
338 : INTEGER :: handle, isp, n_ov_joint, nspins
339 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_ov, offsets
340 : INTEGER, DIMENSION(4) :: reordering
341 : REAL(KIND=dp) :: alpha, alpha_screening
342 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_B, fm_struct_S_joint, &
343 : fm_struct_W
344 : TYPE(cp_fm_type) :: fm_S_joint, fm_W
345 :
346 26 : CALL timeset(routineN, handle)
347 :
348 26 : nspins = SIZE(homo)
349 104 : ALLOCATE (n_ov(nspins), offsets(nspins))
350 26 : CALL get_bse_spin_block_layout(homo, virtual, n_ov, offsets, n_ov_joint)
351 :
352 26 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
353 2 : WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating B'
354 : END IF
355 :
356 : ! Coulomb prefactor: SPIN_CONFIG sector for closed shell; for open shell each spin block
357 : ! contributes once (alpha=1). create_A already emits the SPIN_CONFIG-ignored warning.
358 52 : SELECT CASE (bse_env%bse_spin_config)
359 : CASE (bse_singlet)
360 26 : alpha = 2.0_dp
361 : CASE (bse_triplet)
362 26 : alpha = 0.0_dp
363 : END SELECT
364 26 : IF (nspins > 1) alpha = 1.0_dp
365 :
366 26 : IF (bse_env%screening_method == bse_screening_alpha) THEN
367 2 : alpha_screening = bse_env%screening_factor
368 : ELSE
369 24 : alpha_screening = 1.0_dp
370 : END IF
371 :
372 : ! Joint B over all spin blocks: B_ia,jb = alpha*(ia|bj) - W^sigma_ib,aj (W spin-diagonal)
373 26 : NULLIFY (fm_struct_B)
374 : CALL cp_fm_struct_create(fm_struct_B, context=fm_S_ia(1)%matrix_struct%context, &
375 : nrow_global=n_ov_joint, ncol_global=n_ov_joint, &
376 26 : para_env=fm_S_ia(1)%matrix_struct%para_env)
377 26 : CALL cp_fm_create(fm_B, fm_struct_B, name="fm_B_iajb")
378 26 : CALL cp_fm_set_all(fm_B, 0.0_dp)
379 :
380 : ! Coulomb v_ia,jb = sum_P B^P_ia B^P_jb (= (ia|bj)); cross-spin blocks filled automatically.
381 26 : IF (nspins > 1) THEN
382 4 : NULLIFY (fm_struct_S_joint)
383 : CALL cp_fm_struct_create(fm_struct_S_joint, &
384 : context=fm_S_ia(1)%matrix_struct%context, &
385 : nrow_global=dimen_RI, ncol_global=n_ov_joint, &
386 4 : para_env=fm_S_ia(1)%matrix_struct%para_env)
387 4 : CALL cp_fm_create(fm_S_joint, fm_struct_S_joint, name="fm_S_ia_joint")
388 4 : CALL cp_fm_set_all(fm_S_joint, 0.0_dp)
389 4 : CALL assemble_joint_ov_slab(fm_S_ia, offsets, n_ov, dimen_RI, fm_S_joint)
390 : CALL parallel_gemm(transa="T", transb="N", m=n_ov_joint, n=n_ov_joint, k=dimen_RI, &
391 : alpha=alpha, matrix_a=fm_S_joint, matrix_b=fm_S_joint, beta=0.0_dp, &
392 4 : matrix_c=fm_B)
393 4 : CALL cp_fm_release(fm_S_joint)
394 4 : CALL cp_fm_struct_release(fm_struct_S_joint)
395 : ELSE
396 : CALL parallel_gemm(transa="T", transb="N", m=homo(1)*virtual(1), n=homo(1)*virtual(1), &
397 : k=dimen_RI, alpha=alpha, &
398 : matrix_a=fm_S_ia(1), matrix_b=fm_S_ia(1), &
399 22 : beta=0.0_dp, matrix_c=fm_B)
400 : END IF
401 :
402 26 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
403 2 : WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated B_iajb'
404 : END IF
405 :
406 : ! W^sigma_ib,aj = sum_P barB^P_ib B^P_aj on sigma-diagonal blocks only (offsets place them).
407 : ! reordering [1,4,3,2] maps W_ib,ja -> B_ia,jb. For nspins=1: offsets(1)=0 (original code).
408 26 : IF (bse_env%screening_method /= bse_screening_rpa) THEN
409 52 : DO isp = 1, nspins
410 28 : NULLIFY (fm_struct_W)
411 : CALL cp_fm_struct_create(fm_struct_W, &
412 : context=fm_S_ia(isp)%matrix_struct%context, &
413 : nrow_global=homo(isp)*virtual(isp), &
414 : ncol_global=homo(isp)*virtual(isp), &
415 28 : para_env=fm_S_ia(isp)%matrix_struct%para_env)
416 28 : CALL cp_fm_create(fm_W, fm_struct_W, name="fm_W_ibaj")
417 28 : CALL cp_fm_set_all(fm_W, 0.0_dp)
418 : CALL parallel_gemm(transa="T", transb="N", m=homo(isp)*virtual(isp), &
419 : n=homo(isp)*virtual(isp), k=dimen_RI, alpha=alpha_screening, &
420 : matrix_a=fm_S_bar_ia(isp), matrix_b=fm_S_ia(isp), &
421 28 : beta=0.0_dp, matrix_c=fm_W)
422 28 : reordering = [1, 4, 3, 2]
423 : CALL cp_fm_permute(fm_B, fm_W, -1.0_dp, virtual(isp), virtual(isp), virtual(isp), virtual(isp), &
424 28 : reordering, row_offset=offsets(isp), col_offset=offsets(isp))
425 28 : CALL cp_fm_release(fm_W)
426 80 : CALL cp_fm_struct_release(fm_struct_W)
427 : END DO
428 : END IF
429 :
430 26 : CALL cp_fm_struct_release(fm_struct_B)
431 26 : DEALLOCATE (n_ov, offsets)
432 :
433 26 : CALL timestop(handle)
434 :
435 52 : END SUBROUTINE create_B
436 :
437 : ! **************************************************************************************************
438 : !> \brief Explicit A, and B for ABBA, from the slabs the screening method contracts: the bare ones for
439 : !> TDHF and ALPHA screening, the W-multiplied ones otherwise
440 : !> \param fm_S_ia ...
441 : !> \param fm_S_ij ...
442 : !> \param fm_S_ab ...
443 : !> \param fm_S_bar_ia ...
444 : !> \param fm_S_bar_ij ...
445 : !> \param fm_A ...
446 : !> \param fm_B created only for ABBA
447 : !> \param do_abba ...
448 : !> \param Eigenval ...
449 : !> \param unit_nr ...
450 : !> \param homo ...
451 : !> \param virtual ...
452 : !> \param dimen_RI ...
453 : !> \param bse_env the BSE environment (settings)
454 : ! **************************************************************************************************
455 54 : SUBROUTINE create_A_and_B(fm_S_ia, fm_S_ij, fm_S_ab, fm_S_bar_ia, &
456 108 : fm_S_bar_ij, fm_A, fm_B, do_abba, Eigenval, unit_nr, &
457 54 : homo, virtual, dimen_RI, bse_env)
458 :
459 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_S_ia, fm_S_ij, fm_S_ab, fm_S_bar_ia, &
460 : fm_S_bar_ij
461 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_A, fm_B
462 : LOGICAL, INTENT(IN) :: do_abba
463 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
464 : INTEGER, INTENT(IN) :: unit_nr
465 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
466 : INTEGER, INTENT(IN) :: dimen_RI
467 : TYPE(bse_env_type), INTENT(IN) :: bse_env
468 :
469 54 : IF (bse_env%screening_method == bse_screening_tdhf .OR. &
470 : bse_env%screening_method == bse_screening_alpha) THEN
471 : CALL create_A(fm_S_ia, fm_S_ij, fm_S_ab, fm_A, Eigenval, unit_nr, &
472 4 : homo, virtual, dimen_RI, bse_env)
473 4 : IF (do_abba) THEN
474 : CALL create_B(fm_S_ia, fm_S_ia, fm_B, homo, virtual, dimen_RI, unit_nr, &
475 4 : bse_env)
476 : END IF
477 : ELSE
478 : CALL create_A(fm_S_ia, fm_S_bar_ij, fm_S_ab, fm_A, Eigenval, unit_nr, &
479 50 : homo, virtual, dimen_RI, bse_env)
480 50 : IF (do_abba) THEN
481 : CALL create_B(fm_S_ia, fm_S_bar_ia, fm_B, homo, virtual, dimen_RI, unit_nr, &
482 22 : bse_env)
483 : END IF
484 : END IF
485 :
486 54 : END SUBROUTINE create_A_and_B
487 :
488 : ! **************************************************************************************************
489 : !> \brief Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
490 : !> (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
491 : !> We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
492 : !> of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
493 : !> \param fm_A ...
494 : !> \param fm_B ...
495 : !> \param fm_C ...
496 : !> \param fm_sqrt_A_minus_B ...
497 : !> \param fm_inv_sqrt_A_minus_B ...
498 : !> \param unit_nr ...
499 : !> \param bse_env the BSE environment (settings)
500 : !> \param diag_est ...
501 : ! **************************************************************************************************
502 234 : SUBROUTINE create_hermitian_form_of_ABBA(fm_A, fm_B, fm_C, &
503 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
504 : unit_nr, bse_env, diag_est)
505 :
506 : TYPE(cp_fm_type), INTENT(IN) :: fm_A, fm_B
507 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_C, fm_sqrt_A_minus_B, &
508 : fm_inv_sqrt_A_minus_B
509 : INTEGER, INTENT(IN) :: unit_nr
510 : TYPE(bse_env_type), INTENT(IN) :: bse_env
511 : REAL(KIND=dp), INTENT(IN) :: diag_est
512 :
513 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_hermitian_form_of_ABBA'
514 :
515 : INTEGER :: dim_mat, handle, n_dependent
516 : REAL(KIND=dp), DIMENSION(2) :: eigvals_AB_diff
517 : TYPE(cp_fm_type) :: fm_A_minus_B, fm_A_plus_B, fm_dummy, &
518 : fm_work_product
519 :
520 26 : CALL timeset(routineN, handle)
521 :
522 26 : IF (unit_nr > 0) THEN
523 13 : WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', &
524 26 : ' with size of A. This will take around ', diag_est, " s."
525 : END IF
526 :
527 : ! Create work matrices, which will hold A+B and A-B and their powers
528 : ! C is created afterwards to save memory
529 : ! Final result: C = (A-B)^0.5 (A+B) (A-B)^0.5 EQ.I
530 : ! \_______/ \___/ \______/
531 : ! fm_sqrt_A_minus_B fm_A_plus_B fm_sqrt_A_minus_B
532 : ! (EQ.Ia) (EQ.Ib) (EQ.Ia)
533 : ! Intermediate work matrices:
534 : ! fm_inv_sqrt_A_minus_B: (A-B)^-0.5 EQ.II
535 : ! fm_A_minus_B: (A-B) EQ.III
536 : ! fm_work_product: (A-B)^0.5 (A+B) from (EQ.Ia) and (EQ.Ib) EQ.IV
537 26 : CALL cp_fm_create(fm_A_plus_B, fm_A%matrix_struct)
538 26 : CALL cp_fm_to_fm(fm_A, fm_A_plus_B)
539 26 : CALL cp_fm_create(fm_A_minus_B, fm_A%matrix_struct)
540 26 : CALL cp_fm_to_fm(fm_A, fm_A_minus_B)
541 26 : CALL cp_fm_create(fm_sqrt_A_minus_B, fm_A%matrix_struct)
542 26 : CALL cp_fm_set_all(fm_sqrt_A_minus_B, 0.0_dp)
543 26 : CALL cp_fm_create(fm_inv_sqrt_A_minus_B, fm_A%matrix_struct)
544 26 : CALL cp_fm_set_all(fm_inv_sqrt_A_minus_B, 0.0_dp)
545 :
546 26 : CALL cp_fm_create(fm_work_product, fm_A%matrix_struct)
547 :
548 26 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
549 2 : WRITE (unit_nr, '(T2,A10,T13,A19)') 'BSE|DEBUG|', 'Created work arrays'
550 : END IF
551 :
552 : ! Add/Substract B (cf. EQs. Ib and III)
553 26 : CALL cp_fm_scale_and_add(1.0_dp, fm_A_plus_B, 1.0_dp, fm_B)
554 26 : CALL cp_fm_scale_and_add(1.0_dp, fm_A_minus_B, -1.0_dp, fm_B)
555 :
556 : ! cp_fm_power will overwrite matrix, therefore we create copies
557 26 : CALL cp_fm_to_fm(fm_A_minus_B, fm_inv_sqrt_A_minus_B)
558 :
559 : ! In order to avoid a second diagonalization (cp_fm_power), we create (A-B)^0.5 (EQ.Ia)
560 : ! from (A-B)^-0.5 (EQ.II) by multiplication with (A-B) (EQ.III) afterwards.
561 :
562 : ! Raise A-B to -0.5_dp, no quenching of eigenvectors, hence threshold=0.0_dp
563 26 : CALL cp_fm_create(fm_dummy, fm_A%matrix_struct)
564 : ! Create (A-B)^-0.5 (cf. EQ.II)
565 26 : CALL cp_fm_power(fm_inv_sqrt_A_minus_B, fm_dummy, -0.5_dp, 0.0_dp, n_dependent, eigvals=eigvals_AB_diff)
566 26 : CALL cp_fm_release(fm_dummy)
567 26 : IF (unit_nr > 0) THEN
568 13 : WRITE (unit_nr, '(T2,A4,T7,A,T65,F16.6)') 'BSE|', &
569 26 : 'Smallest eigenvalue of A-B (eV)', eigvals_AB_diff(1)*evolt
570 : END IF
571 : ! Raise an error in case the the matrix A-B is not positive definite (i.e. negative eigenvalues)
572 : ! In this case, the procedure for hermitian form of ABBA is not applicable
573 26 : IF (eigvals_AB_diff(1) < 0) THEN
574 : CALL cp_abort(__LOCATION__, &
575 : "Matrix (A-B) is not positive definite. "// &
576 0 : "Hermitian diagonalization of full ABBA matrix is ill-defined.")
577 : END IF
578 :
579 : ! We keep fm_inv_sqrt_A_minus_B for print of singleparticle transitions of ABBA
580 : ! We further create (A-B)^0.5 for the singleparticle transitions of ABBA
581 : ! Create (A-B)^0.5= (A-B)^-0.5 * (A-B) (EQ.Ia)
582 26 : CALL cp_fm_get_info(fm_A, nrow_global=dim_mat)
583 : 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, &
584 26 : fm_sqrt_A_minus_B)
585 :
586 : ! Compute and store LHS of C, i.e. (A-B)^0.5 (A+B) (EQ.IV)
587 : 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, &
588 26 : fm_work_product)
589 :
590 : ! Release to save memory
591 26 : CALL cp_fm_release(fm_A_plus_B)
592 26 : CALL cp_fm_release(fm_A_minus_B)
593 :
594 : ! Now create full
595 26 : CALL cp_fm_create(fm_C, fm_A%matrix_struct)
596 26 : CALL cp_fm_set_all(fm_C, 0.0_dp)
597 : ! Compute C=(A-B)^0.5 (A+B) (A-B)^0.5 (EQ.I)
598 : 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, &
599 26 : fm_C)
600 26 : CALL cp_fm_release(fm_work_product)
601 :
602 26 : IF (unit_nr > 0 .AND. bse_env%bse_debug_print) THEN
603 2 : WRITE (unit_nr, '(T2,A10,T13,A36)') 'BSE|DEBUG|', 'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5'
604 : END IF
605 :
606 26 : CALL timestop(handle)
607 26 : END SUBROUTINE create_hermitian_form_of_ABBA
608 :
609 : ! **************************************************************************************************
610 : !> \brief Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
611 : !> Here, the eigenvectors Z^n relate to X^n via
612 : !> Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
613 : !> \param fm_C ...
614 : !> \param homo ...
615 : !> \param virtual ...
616 : !> \param homo_irred ...
617 : !> \param fm_sqrt_A_minus_B ...
618 : !> \param fm_inv_sqrt_A_minus_B ...
619 : !> \param unit_nr ...
620 : !> \param diag_est ...
621 : !> \param bse_env the BSE environment (settings, and the state the post-processing reads)
622 : !> \param mo_coeff ...
623 : ! **************************************************************************************************
624 26 : SUBROUTINE diagonalize_C(fm_C, homo, virtual, homo_irred, &
625 : fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
626 26 : unit_nr, diag_est, bse_env, mo_coeff)
627 :
628 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_C
629 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
630 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B
631 : INTEGER, INTENT(IN) :: unit_nr
632 : REAL(KIND=dp), INTENT(IN) :: diag_est
633 : TYPE(bse_env_type), INTENT(IN) :: bse_env
634 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
635 :
636 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_C'
637 :
638 : INTEGER :: diag_info, handle, n_ov_joint, nspins
639 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
640 : TYPE(cp_fm_type) :: fm_eigvec_X, fm_eigvec_Y, fm_eigvec_Z, &
641 : fm_mat_eigvec_transform_diff, &
642 : fm_mat_eigvec_transform_sum
643 :
644 26 : CALL timeset(routineN, handle)
645 :
646 26 : nspins = SIZE(homo)
647 56 : n_ov_joint = SUM(homo*virtual)
648 :
649 26 : IF (unit_nr > 0) THEN
650 13 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', &
651 26 : 'This will take around ', diag_est, ' s.'
652 : END IF
653 :
654 : !We have now the full matrix C=(A-B)^0.5 (A+B) (A-B)^0.5
655 : !Now: Diagonalize it
656 26 : CALL cp_fm_create(fm_eigvec_Z, fm_C%matrix_struct)
657 :
658 78 : ALLOCATE (Exc_ens(n_ov_joint))
659 :
660 26 : CALL choose_eigv_solver(fm_C, fm_eigvec_Z, Exc_ens, diag_info)
661 :
662 26 : IF (diag_info /= 0) THEN
663 : CALL cp_abort(__LOCATION__, &
664 0 : "Diagonalization of C=(A-B)^0.5 (A+B) (A-B)^0.5 failed in BSE")
665 : END IF
666 :
667 : ! C could have negative eigenvalues, since we do not explicitly check A+B
668 : ! for positive definiteness (would make another O(N^6) Diagon. necessary)
669 : ! Instead, we include a check here
670 26 : IF (Exc_ens(1) < 0) THEN
671 0 : IF (unit_nr > 0) THEN
672 : CALL cp_abort(__LOCATION__, &
673 : "Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
674 0 : "(A+B) is not positive definite.")
675 : END IF
676 : END IF
677 2902 : Exc_ens = SQRT(Exc_ens)
678 :
679 : ! Prepare eigenvector for interpretation of singleparticle transitions
680 : ! Compare: F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)
681 : ! We aim for the upper part of the vector (X,Y) for a direct comparison with the TDA result
682 :
683 : ! Following Furche, we basically use Eqs. (A10): First, we multiply
684 : ! the (A-B)^+-0.5 with eigenvectors and then the eigenvalues
685 : ! One has to be careful about the index structure, since the eigenvector matrix is not symmetric anymore!
686 :
687 : ! First, Eq. I from (A10) from Furche: (X+Y)_n = (Ω_n)^-0.5 (A-B)^0.5 T_n
688 26 : CALL cp_fm_create(fm_mat_eigvec_transform_sum, fm_C%matrix_struct)
689 26 : CALL cp_fm_set_all(fm_mat_eigvec_transform_sum, 0.0_dp)
690 : CALL parallel_gemm(transa="N", transb="N", m=n_ov_joint, n=n_ov_joint, k=n_ov_joint, alpha=1.0_dp, &
691 : matrix_a=fm_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
692 26 : matrix_c=fm_mat_eigvec_transform_sum)
693 26 : CALL cp_fm_release(fm_sqrt_A_minus_B)
694 : ! This normalizes the eigenvectors
695 26 : CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_sum, Exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
696 :
697 : ! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n
698 26 : CALL cp_fm_create(fm_mat_eigvec_transform_diff, fm_C%matrix_struct)
699 26 : CALL cp_fm_set_all(fm_mat_eigvec_transform_diff, 0.0_dp)
700 : CALL parallel_gemm(transa="N", transb="N", m=n_ov_joint, n=n_ov_joint, k=n_ov_joint, alpha=1.0_dp, &
701 : matrix_a=fm_inv_sqrt_A_minus_B, matrix_b=fm_eigvec_Z, beta=0.0_dp, &
702 26 : matrix_c=fm_mat_eigvec_transform_diff)
703 26 : CALL cp_fm_release(fm_inv_sqrt_A_minus_B)
704 26 : CALL cp_fm_release(fm_eigvec_Z)
705 :
706 : ! This normalizes the eigenvectors
707 26 : CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_diff, Exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
708 :
709 : ! Now, we add the two equations to obtain X_n
710 : ! Add overwrites the first argument, therefore we copy it beforehand
711 26 : CALL cp_fm_create(fm_eigvec_X, fm_C%matrix_struct)
712 26 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_X)
713 26 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_X, 1.0_dp, fm_mat_eigvec_transform_diff)
714 :
715 : ! Now, we subtract the two equations to obtain Y_n
716 : ! Add overwrites the first argument, therefore we copy it beforehand
717 26 : CALL cp_fm_create(fm_eigvec_Y, fm_C%matrix_struct)
718 26 : CALL cp_fm_to_fm(fm_mat_eigvec_transform_sum, fm_eigvec_Y)
719 26 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_Y, -1.0_dp, fm_mat_eigvec_transform_diff)
720 :
721 : !Cleanup
722 26 : CALL cp_fm_release(fm_mat_eigvec_transform_diff)
723 26 : CALL cp_fm_release(fm_mat_eigvec_transform_sum)
724 :
725 26 : IF (nspins == 1) THEN
726 : CALL postprocess_bse(Exc_ens, fm_eigvec_X, bse_env, mo_coeff, &
727 : homo(1), virtual(1), homo_irred(1), unit_nr, &
728 22 : .FALSE., fm_eigvec_Y)
729 : ELSE
730 : ! Open-shell ABBA: helper forms X+Y internally and prints amplitudes (X and Y).
731 : CALL bse_open_shell_optical(Exc_ens, fm_eigvec_X, homo, virtual, homo_irred, &
732 4 : .FALSE., mo_coeff, bse_env, unit_nr, fm_eigvec_Y)
733 : END IF
734 :
735 26 : DEALLOCATE (Exc_ens)
736 26 : CALL cp_fm_release(fm_eigvec_X)
737 26 : CALL cp_fm_release(fm_eigvec_Y)
738 :
739 26 : CALL timestop(handle)
740 :
741 182 : END SUBROUTINE diagonalize_C
742 :
743 : ! **************************************************************************************************
744 : !> \brief Solving hermitian eigenvalue equation A X^n = Ω^n X^n
745 : !> \param fm_A ...
746 : !> \param homo ...
747 : !> \param virtual ...
748 : !> \param homo_irred ...
749 : !> \param unit_nr ...
750 : !> \param diag_est ...
751 : !> \param bse_env the BSE environment (settings, and the state the post-processing reads)
752 : !> \param mo_coeff ...
753 : ! **************************************************************************************************
754 26 : SUBROUTINE diagonalize_A(fm_A, homo, virtual, homo_irred, &
755 26 : unit_nr, diag_est, bse_env, mo_coeff)
756 :
757 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_A
758 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
759 : INTEGER, INTENT(IN) :: unit_nr
760 : REAL(KIND=dp), INTENT(IN) :: diag_est
761 : TYPE(bse_env_type), INTENT(IN) :: bse_env
762 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
763 :
764 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_A'
765 :
766 : INTEGER :: diag_info, handle, n_ov_joint, nspins
767 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
768 : TYPE(cp_fm_type) :: fm_eigvec
769 :
770 26 : CALL timeset(routineN, handle)
771 :
772 26 : nspins = SIZE(homo)
773 62 : n_ov_joint = SUM(homo*virtual)
774 :
775 26 : IF (unit_nr > 0) THEN
776 13 : WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
777 26 : 'This will take around ', diag_est, ' s.'
778 : END IF
779 :
780 26 : CALL cp_fm_create(fm_eigvec, fm_A%matrix_struct)
781 :
782 78 : ALLOCATE (Exc_ens(n_ov_joint))
783 :
784 26 : CALL choose_eigv_solver(fm_A, fm_eigvec, Exc_ens, diag_info)
785 :
786 26 : IF (diag_info /= 0) THEN
787 : CALL cp_abort(__LOCATION__, &
788 0 : "Diagonalization of A failed in TDA-BSE")
789 : END IF
790 :
791 26 : IF (nspins == 1) THEN
792 : CALL postprocess_bse(Exc_ens, fm_eigvec, bse_env, mo_coeff, &
793 16 : homo(1), virtual(1), homo_irred(1), unit_nr, .TRUE.)
794 : ELSE
795 : CALL bse_open_shell_optical(Exc_ens, fm_eigvec, homo, virtual, homo_irred, &
796 10 : .TRUE., mo_coeff, bse_env, unit_nr)
797 : END IF
798 :
799 26 : CALL cp_fm_release(fm_eigvec)
800 26 : DEALLOCATE (Exc_ens)
801 :
802 26 : CALL timestop(handle)
803 :
804 78 : END SUBROUTINE diagonalize_A
805 :
806 : ! **************************************************************************************************
807 : !> \brief Open-shell (UKS) spin-summed post-processing for the joint spin-block space: joint
808 : !> excitation energies, per-spin transition amplitudes, and oscillator strengths. Mirrors
809 : !> postprocess_bse but spin-summed; exciton descriptors and NTOs are not yet implemented (CPWARN).
810 : !> \param Exc_ens joint excitation energies
811 : !> \param fm_eigvec_X joint X eigenvectors (excitations)
812 : !> \param homo per-spin reduced/active occupied counts
813 : !> \param virtual per-spin reduced/active virtual counts
814 : !> \param homo_irred per-spin full occupied counts (absolute-MO labels; N_e = sum)
815 : !> \param flag_tda .TRUE. -> TDA (coeff=X), .FALSE. -> ABBA (coeff=X+Y)
816 : !> \param mo_coeff per-spin MO coefficients
817 : !> \param bse_env the BSE environment (settings; the AO dipoles, mos)
818 : !> \param unit_nr ...
819 : !> \param fm_eigvec_Y joint Y eigenvectors (deexcitations; ABBA only)
820 : ! **************************************************************************************************
821 14 : SUBROUTINE bse_open_shell_optical(Exc_ens, fm_eigvec_X, homo, virtual, homo_irred, &
822 14 : flag_tda, mo_coeff, bse_env, unit_nr, fm_eigvec_Y)
823 :
824 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
825 : TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
826 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual, homo_irred
827 : LOGICAL, INTENT(IN) :: flag_tda
828 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
829 : TYPE(bse_env_type), INTENT(IN) :: bse_env
830 : INTEGER, INTENT(IN) :: unit_nr
831 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
832 :
833 : CHARACTER(LEN=*), PARAMETER :: routineN = 'bse_open_shell_optical'
834 :
835 : CHARACTER(LEN=10) :: info_approximation, multiplet
836 : INTEGER :: handle, idir, isp, jdir, n, n_ov_joint, &
837 : nspins
838 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_ov_sp, offsets_sp
839 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str_joint
840 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pol_res_joint, trans_mom_joint
841 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_dip_reord, fm_struct_sp, &
842 : fm_struct_tmom
843 : TYPE(cp_fm_type) :: fm_dip_reord_sp, fm_eigvec_sp, &
844 : fm_trans_coeff
845 14 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dip_ab_sp, fm_dip_ai_sp, fm_dip_ij_sp
846 56 : TYPE(cp_fm_type), DIMENSION(3) :: fm_trans_mom_joint
847 :
848 14 : CALL timeset(routineN, handle)
849 :
850 14 : nspins = SIZE(homo)
851 56 : ALLOCATE (n_ov_sp(nspins), offsets_sp(nspins))
852 14 : CALL get_bse_spin_block_layout(homo, virtual, n_ov_sp, offsets_sp, n_ov_joint)
853 :
854 : ! LEN=10 locals auto-pad short literals with spaces (avoids the L-15 short-literal trap);
855 : ! print_excitation_energies prints A6 of these, print_optical_properties prints them as-is.
856 14 : multiplet = "UKS"
857 14 : IF (flag_tda) THEN
858 10 : info_approximation = " -TDA- "
859 : ELSE
860 4 : info_approximation = "-ABBA-"
861 : END IF
862 :
863 14 : IF (unit_nr > 0) THEN
864 7 : WRITE (unit_nr, '(T2,A4,T7,A43)') 'BSE|', 'Joint open-shell BSE excitation energies:'
865 : END IF
866 : CALL print_excitation_energies(Exc_ens, n_ov_joint, 1, flag_tda, multiplet, &
867 14 : info_approximation, bse_env, unit_nr)
868 :
869 : ! Per-spin single-particle transition amplitudes (X via =>, Y via <=).
870 : CALL print_transition_amplitudes(fm_eigvec_X, homo, virtual, homo_irred, &
871 14 : info_approximation, bse_env, unit_nr, fm_eigvec_Y)
872 :
873 : ! Transition coefficient for the spin-summed moment: X (TDA) or X+Y (ABBA).
874 14 : CALL cp_fm_create(fm_trans_coeff, fm_eigvec_X%matrix_struct)
875 14 : CALL cp_fm_to_fm(fm_eigvec_X, fm_trans_coeff)
876 14 : IF (PRESENT(fm_eigvec_Y)) CALL cp_fm_scale_and_add(1.0_dp, fm_trans_coeff, 1.0_dp, fm_eigvec_Y)
877 :
878 : ! Spin-summed transition moments: D^n_dir = sum_σ sum_{ia,σ} D^{dir,σ}_{ai} C_{ia,σ,n}
879 : ! with C = X (TDA) or X+Y (ABBA); explicit spin sum, factor 1.0.
880 140 : ALLOCATE (fm_dip_ai_sp(3), fm_dip_ij_sp(3), fm_dip_ab_sp(3))
881 70 : ALLOCATE (oscill_str_joint(n_ov_joint), trans_mom_joint(3, 1, n_ov_joint))
882 42 : ALLOCATE (pol_res_joint(3, 3, n_ov_joint))
883 14 : trans_mom_joint(:, :, :) = 0.0_dp
884 14 : NULLIFY (fm_struct_dip_reord, fm_struct_sp, fm_struct_tmom)
885 : CALL cp_fm_struct_create(fm_struct_tmom, fm_trans_coeff%matrix_struct%para_env, &
886 14 : fm_trans_coeff%matrix_struct%context, 1, n_ov_joint)
887 56 : DO idir = 1, 3
888 42 : CALL cp_fm_create(fm_trans_mom_joint(idir), fm_struct_tmom)
889 56 : CALL cp_fm_set_all(fm_trans_mom_joint(idir), 0.0_dp)
890 : END DO
891 42 : DO isp = 1, nspins
892 : CALL get_multipoles_mo(fm_dip_ai_sp, fm_dip_ij_sp, fm_dip_ab_sp, &
893 : bse_env%matrix_dipole_ao, bse_env%mos, mo_coeff(isp:isp), &
894 : homo(isp), virtual(isp), fm_trans_coeff%matrix_struct%context, &
895 28 : ispin=isp)
896 28 : NULLIFY (fm_struct_sp, fm_struct_dip_reord)
897 : CALL cp_fm_struct_create(fm_struct_sp, fm_trans_coeff%matrix_struct%para_env, &
898 28 : fm_trans_coeff%matrix_struct%context, n_ov_sp(isp), n_ov_joint)
899 28 : CALL cp_fm_create(fm_eigvec_sp, fm_struct_sp)
900 28 : CALL cp_fm_set_all(fm_eigvec_sp, 0.0_dp)
901 : CALL cp_fm_to_fm_submat(fm_trans_coeff, fm_eigvec_sp, n_ov_sp(isp), n_ov_joint, &
902 28 : offsets_sp(isp) + 1, 1, 1, 1)
903 : CALL cp_fm_struct_create(fm_struct_dip_reord, fm_trans_coeff%matrix_struct%para_env, &
904 28 : fm_trans_coeff%matrix_struct%context, 1, n_ov_sp(isp))
905 112 : DO idir = 1, 3
906 84 : CALL cp_fm_create(fm_dip_reord_sp, fm_struct_dip_reord, name="bse_dip_reord")
907 84 : CALL cp_fm_set_all(fm_dip_reord_sp, 0.0_dp)
908 84 : CALL cp_fm_permute(fm_dip_reord_sp, fm_dip_ai_sp(idir), 1.0_dp, 1, 1, 1, virtual(isp), [2, 4, 3, 1])
909 : CALL parallel_gemm('N', 'N', 1, n_ov_joint, n_ov_sp(isp), 1.0_dp, &
910 84 : fm_dip_reord_sp, fm_eigvec_sp, 1.0_dp, fm_trans_mom_joint(idir))
911 84 : CALL cp_fm_release(fm_dip_reord_sp)
912 84 : CALL cp_fm_release(fm_dip_ai_sp(idir))
913 84 : CALL cp_fm_release(fm_dip_ij_sp(idir))
914 196 : CALL cp_fm_release(fm_dip_ab_sp(idir))
915 : END DO
916 28 : CALL cp_fm_release(fm_eigvec_sp)
917 28 : CALL cp_fm_struct_release(fm_struct_sp)
918 28 : NULLIFY (fm_struct_sp)
919 28 : CALL cp_fm_struct_release(fm_struct_dip_reord)
920 70 : NULLIFY (fm_struct_dip_reord)
921 : END DO
922 56 : DO idir = 1, 3
923 42 : CALL cp_fm_get_submatrix(fm_trans_mom_joint(idir), trans_mom_joint(idir, :, :))
924 56 : CALL cp_fm_release(fm_trans_mom_joint(idir))
925 : END DO
926 14 : CALL cp_fm_struct_release(fm_struct_tmom)
927 5886 : DO n = 1, n_ov_joint
928 23488 : DO idir = 1, 3
929 76336 : DO jdir = 1, 3
930 : pol_res_joint(idir, jdir, n) = 2.0_dp*Exc_ens(n)*trans_mom_joint(idir, 1, n) &
931 70464 : *trans_mom_joint(jdir, 1, n)
932 : END DO
933 : END DO
934 23502 : oscill_str_joint(n) = 2.0_dp/3.0_dp*Exc_ens(n)*SUM(ABS(trans_mom_joint(:, 1, n))**2)
935 : END DO
936 : CALL print_optical_properties(Exc_ens, oscill_str_joint, trans_mom_joint, pol_res_joint, &
937 : n_ov_joint, 1, SUM(homo_irred), flag_tda, info_approximation, &
938 42 : bse_env, unit_nr, open_shell=.TRUE.)
939 : ! Open-shell post-processing is partial: energies, amplitudes, spin-summed oscillator strengths.
940 : CALL cp_warn(__LOCATION__, &
941 : "Open-shell (UKS) BSE: exciton descriptors and NTO analysis are not yet "// &
942 14 : "implemented and have been skipped.")
943 14 : CALL cp_fm_release(fm_trans_coeff)
944 14 : DEALLOCATE (fm_dip_ai_sp, fm_dip_ij_sp, fm_dip_ab_sp)
945 14 : DEALLOCATE (n_ov_sp, offsets_sp, oscill_str_joint, trans_mom_joint, pol_res_joint)
946 :
947 14 : CALL timestop(handle)
948 :
949 56 : END SUBROUTINE bse_open_shell_optical
950 :
951 : ! **************************************************************************************************
952 : !> \brief Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag)
953 : !> \param Exc_ens ...
954 : !> \param fm_eigvec_X ...
955 : !> \param bse_env the BSE environment (settings; the AO multipoles and their reference point, mos, subsys)
956 : !> \param mo_coeff ...
957 : !> \param homo ...
958 : !> \param virtual ...
959 : !> \param homo_irred ...
960 : !> \param unit_nr ...
961 : !> \param flag_TDA ...
962 : !> \param fm_eigvec_Y ...
963 : ! **************************************************************************************************
964 46 : SUBROUTINE postprocess_bse(Exc_ens, fm_eigvec_X, bse_env, mo_coeff, &
965 : homo, virtual, homo_irred, unit_nr, &
966 : flag_TDA, fm_eigvec_Y)
967 :
968 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
969 : TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
970 : TYPE(bse_env_type), INTENT(IN) :: bse_env
971 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
972 : INTEGER :: homo, virtual, homo_irred, unit_nr
973 : LOGICAL, OPTIONAL :: flag_TDA
974 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
975 :
976 : CHARACTER(LEN=*), PARAMETER :: routineN = 'postprocess_bse'
977 :
978 : CHARACTER(LEN=10) :: info_approximation, multiplet
979 : INTEGER :: handle, i_exc, idir, n_moments_di, &
980 : n_moments_quad
981 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str, ref_point_multipole
982 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: polarizability_residues, trans_mom_bse
983 : TYPE(cp_fm_type) :: fm_X_ia, fm_Y_ia
984 46 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_dipole_ab_trunc, fm_dipole_ai_trunc, &
985 46 : fm_dipole_ij_trunc, fm_quadpole_ab_trunc, fm_quadpole_ai_trunc, fm_quadpole_ij_trunc
986 : TYPE(exciton_descr_type), ALLOCATABLE, &
987 46 : DIMENSION(:) :: exc_descr
988 46 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
989 :
990 46 : CALL timeset(routineN, handle)
991 :
992 : !Prepare variables for printing
993 46 : IF (bse_env%bse_spin_config == 0) THEN
994 46 : multiplet = "Singlet"
995 : ELSE
996 0 : multiplet = "Triplet"
997 : END IF
998 46 : IF (.NOT. PRESENT(flag_TDA)) THEN
999 0 : flag_TDA = .FALSE.
1000 : END IF
1001 46 : IF (flag_TDA) THEN
1002 20 : info_approximation = " -TDA- "
1003 : ELSE
1004 26 : info_approximation = "-ABBA-"
1005 : END IF
1006 :
1007 46 : n_moments_di = 3
1008 46 : n_moments_quad = 9
1009 : ! Compute BSE dipoles and oscillator strengths - Keep in memory for later usage
1010 : ! Need dipoles also for spatial expectation values, which are well-defined also for triplets
1011 184 : ALLOCATE (fm_dipole_ij_trunc(n_moments_di))
1012 184 : ALLOCATE (fm_dipole_ab_trunc(n_moments_di))
1013 184 : ALLOCATE (fm_dipole_ai_trunc(n_moments_di))
1014 46 : ALLOCATE (ref_point_multipole(3))
1015 184 : ref_point_multipole(:) = bse_env%multipole_ref_point
1016 : ! Obtain dipoles in MO basis
1017 : CALL get_multipoles_mo(fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc, &
1018 : bse_env%matrix_dipole_ao, bse_env%mos, mo_coeff, &
1019 46 : homo, virtual, fm_eigvec_X%matrix_struct%context)
1020 : ! Compute exciton descriptors from these multipoles
1021 46 : IF (bse_env%num_print_exc_descr > 0) THEN
1022 : ! Obtain quadrupoles in MO basis
1023 40 : ALLOCATE (fm_quadpole_ij_trunc(n_moments_quad))
1024 40 : ALLOCATE (fm_quadpole_ab_trunc(n_moments_quad))
1025 40 : ALLOCATE (fm_quadpole_ai_trunc(n_moments_quad))
1026 : ! the preparer builds the AO quadrupoles whenever NUM_PRINT_EXC_DESCR is nonzero
1027 4 : CPASSERT(ASSOCIATED(bse_env%matrix_quadpole_ao))
1028 : CALL get_multipoles_mo(fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1029 : bse_env%matrix_quadpole_ao, bse_env%mos, mo_coeff, &
1030 4 : homo, virtual, fm_eigvec_X%matrix_struct%context)
1031 : ! Iterate over excitation index outside of routine to make it compatible with tddft module
1032 424 : ALLOCATE (exc_descr(bse_env%num_print_exc_descr))
1033 104 : DO i_exc = 1, bse_env%num_print_exc_descr
1034 100 : CALL reshuffle_eigvec(fm_eigvec_X, fm_X_ia, homo, virtual, i_exc, .FALSE.)
1035 100 : IF (.NOT. flag_TDA) THEN
1036 50 : CALL reshuffle_eigvec(fm_eigvec_Y, fm_Y_ia, homo, virtual, i_exc, .FALSE.)
1037 :
1038 : CALL get_exciton_descriptors(exc_descr, fm_X_ia, &
1039 : fm_quadpole_ij_trunc, fm_quadpole_ab_trunc, &
1040 : fm_quadpole_ai_trunc, &
1041 : i_exc, homo, virtual, &
1042 50 : fm_Y_ia)
1043 : ELSE
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 50 : i_exc, homo, virtual)
1048 : END IF
1049 100 : CALL cp_fm_release(fm_X_ia)
1050 104 : IF (.NOT. flag_TDA) THEN
1051 50 : CALL cp_fm_release(fm_Y_ia)
1052 : END IF
1053 : END DO
1054 : END IF
1055 :
1056 46 : IF (bse_env%bse_spin_config == 0) THEN
1057 : CALL get_oscillator_strengths(fm_eigvec_X, Exc_ens, fm_dipole_ai_trunc, &
1058 : trans_mom_bse, oscill_str, polarizability_residues, &
1059 46 : homo, virtual, fm_eigvec_Y)
1060 : END IF
1061 :
1062 : ! the iterative driver prints the header before its solver runs
1063 46 : IF (bse_env%bse_diag_method /= bse_iterdiag) THEN
1064 38 : CALL print_output_header(homo, virtual, homo_irred, flag_TDA, bse_env, unit_nr)
1065 : END IF
1066 :
1067 : ! Prints excitation energies up to user-specified number
1068 : CALL print_excitation_energies(Exc_ens, homo, virtual, flag_TDA, multiplet, &
1069 46 : info_approximation, bse_env, unit_nr)
1070 :
1071 : ! Print single particle transition amplitudes, i.e. components of eigenvectors X and Y
1072 : CALL print_transition_amplitudes(fm_eigvec_X, [homo], [virtual], [homo_irred], &
1073 184 : info_approximation, bse_env, unit_nr, fm_eigvec_Y)
1074 :
1075 : ! Prints optical properties, if state is a singlet
1076 : CALL print_optical_properties(Exc_ens, oscill_str, trans_mom_bse, polarizability_residues, &
1077 : homo, virtual, homo_irred, flag_TDA, &
1078 46 : info_approximation, bse_env, unit_nr)
1079 : ! Print exciton descriptors if keyword is invoked
1080 46 : IF (bse_env%num_print_exc_descr > 0) THEN
1081 4 : CALL qs_subsys_get(bse_env%subsys, particle_set=particle_set)
1082 : CALL print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, &
1083 : bse_env%num_print_exc_descr, bse_env%bse_debug_print, &
1084 : bse_env%print_directional_exc_descr, &
1085 : bse_env%print_directional_crosscorrelation, &
1086 4 : 'BSE|', particle_set)
1087 : END IF
1088 :
1089 : ! Compute and print excitation wavefunctions
1090 46 : IF (bse_env%do_nto_analysis) THEN
1091 4 : IF (unit_nr > 0) THEN
1092 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1093 : WRITE (unit_nr, '(T2,A4,T7,A47)') &
1094 2 : 'BSE|', "Calculating Natural Transition Orbitals (NTOs)."
1095 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1096 : END IF
1097 : CALL calculate_NTOs(fm_eigvec_X, fm_eigvec_Y, &
1098 : mo_coeff, homo, virtual, &
1099 : info_approximation, &
1100 : oscill_str, &
1101 4 : bse_env, unit_nr)
1102 : END IF
1103 :
1104 184 : DO idir = 1, n_moments_di
1105 138 : CALL cp_fm_release(fm_dipole_ai_trunc(idir))
1106 138 : CALL cp_fm_release(fm_dipole_ij_trunc(idir))
1107 184 : CALL cp_fm_release(fm_dipole_ab_trunc(idir))
1108 : END DO
1109 46 : IF (bse_env%num_print_exc_descr > 0) THEN
1110 40 : DO idir = 1, n_moments_quad
1111 36 : CALL cp_fm_release(fm_quadpole_ai_trunc(idir))
1112 36 : CALL cp_fm_release(fm_quadpole_ij_trunc(idir))
1113 40 : CALL cp_fm_release(fm_quadpole_ab_trunc(idir))
1114 : END DO
1115 4 : DEALLOCATE (fm_quadpole_ai_trunc, fm_quadpole_ij_trunc, fm_quadpole_ab_trunc)
1116 4 : DEALLOCATE (exc_descr)
1117 : END IF
1118 46 : DEALLOCATE (fm_dipole_ai_trunc, fm_dipole_ij_trunc, fm_dipole_ab_trunc)
1119 46 : DEALLOCATE (ref_point_multipole)
1120 46 : IF (bse_env%bse_spin_config == 0) THEN
1121 46 : DEALLOCATE (oscill_str, trans_mom_bse, polarizability_residues)
1122 : END IF
1123 46 : IF (unit_nr > 0) THEN
1124 23 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1125 23 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
1126 : END IF
1127 :
1128 46 : CALL timestop(handle)
1129 :
1130 92 : END SUBROUTINE postprocess_bse
1131 :
1132 : END MODULE bse_full_diag
|
